How can I determine the probability of an event in R for a Hyperbolic Distribution?

I would like to know how I can calculate P ≥ -0.0001. I am novice in both Mathematics, Statistics and R.

I am currently modelling some Spread Data in R. I have been able to model the data to a specific distribution (in the case of the below example this is Hyperbolic) and perform a likelihood ratio test to confirm the viability of using the specified distribution.

In my head, once I have obtained the distribution, I am now, more reliably, able to calculate my probabilities as I understand the area underneath the curve more accurately?

I had assumed that I would be able to conduct a similar test to the t-statistic once the distribution was matched however I can't figure out how to do this.

I have seen the pnorm(x) functions and also seen phyper however can't figure out how to simply calculate P ≥ -0.0001.

Below is my code, your help would be hugely appreciated;

`# library (ghyp) library (timeSeries)

# Coverting to Time Series 
E <- timeSeries(A[,"Spread"])

# Fitting
ef <- (density(E))
ghdfit <- fit.ghypuv(E,symmetric = FALSE, control = list(maxit = 1000))
hypfit <- fit.hypuv(E,symmetric = FALSE, control = list(maxit = 1000))
nigfit <- fit.NIGuv(E,symmetric = FALSE, control = list(maxit = 1000))

# Density
ghddens <- dghyp(ef$x, ghdfit)
    hypdens <- dghyp(ef$x, hypfit)
nigdens <- dghyp(ef$x, nigfit)
    nordens <- dnorm(ef$x, mean = mean(E),sd = sd(c(E[,1])))
col.def <- c("black","red","green","orange")
plot(ef, xlab = " Spread ", ylab = expression(f(x)),ylim = c(0,50), main ='CABLE - 3 Day Comparison across 28 Years')
lines(ef$x, ghddens, col = "red")
    lines(ef$x, hypdens, col = "blue")
lines(ef$x, nigdens, col = "green")
    lines(ef$x, nordens, col = "orange")
legend("topleft", legend = c("Empirical","GHD","HYP","NIG","NORM"), col = col.def, lty = 1)



# QQ Plot
        qqghyp(ghdfit, line = TRUE, ghyp.col = "red", plot.legend = FALSE, gaussian = FALSE, main = " ", cex = 0.8)
        qqghyp(hypfit, add = TRUE, ghyp.pch = 2, ghyp.col = "green", gaussian = FALSE, line = FALSE, cex = 0.8)
        qqghyp(nigfit,add = TRUE, ghyp.pch = 3, ghyp.col = "orange", gaussian = FALSE, line = FALSE, cex = 0.8)
        legend("topleft", legend = c("GHD","HYP","NIG"), col = col.def[-c(1,5)], pch = 1:3)

       # Diagnostic
        options(scipen=999)
        AIC <- stepAIC.ghyp(E, dist = c("ghyp","hyp","NIG"), symmetric = FALSE, control = list(maxit=1000))
        LRghdnig <- lik.ratio.test(ghdfit,nigfit)
        LRghdhyp <- lik.ratio.test(ghdfit,hypfit)

LRghdhyp $statistic L 0.07005745

$p.value 1 0.0211198

$df 1 1

$H0 1 FALSE

So, I know the correct distribution and how to fit it. How do I go about determining the probability of > - 0.0001 occurs?

Answer is posted here and involves the pghyp function @fg nu