# SkewHyperbolicDistribution: Skewed Hyperbolic Student t-Distribution In SkewHyperbolic: The Skew Hyperbolic Student t-Distribution

## Description

Density function, distribution function, quantiles and random number generation for the skew hyperbolic Student t-distribution, with parameters beta (skewness), delta (scale), mu (location) and nu (shape). Also a function for the derivative of the density function.

## Usage

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18``` ```dskewhyp(x, mu = 0, delta = 1, beta = 1, nu = 1, param = c(mu,delta,beta,nu), log = FALSE, tolerance = .Machine\$double.eps^0.5) pskewhyp(q, mu = 0, delta = 1, beta = 1, nu = 1, param = c(mu, delta, beta, nu), log.p = FALSE, lower.tail = TRUE, subdivisions = 100, intTol = .Machine\$double.eps^0.25, valueOnly = TRUE, ...) qskewhyp(p, mu = 0, delta = 1, beta = 1, nu = 1, param = c(mu,delta, beta, nu), lower.tail = TRUE, log.p = FALSE, method = c("spline","integrate"), nInterpol = 501, uniTol = .Machine\$double.eps^0.25, subdivisions = 100, intTol = uniTol, ...) rskewhyp(n, mu = 0, delta = 1, beta = 1, nu = 1, param = c(mu,delta,beta,nu), log = FALSE) ddskewhyp(x, mu = 0, delta = 1, beta = 1, nu = 1, param = c(mu,delta,beta,nu),log = FALSE, tolerance = .Machine\$double.eps^0.5) ```

## Arguments

 `x,q` Vector of quantiles. `p` Vector of probabilities. `n` Number of random variates to be generated. `mu` Location parameter mu, default is 0. `delta` Scale parameter delta, default is 1. `beta` Skewness parameter beta, default is 1. `nu` Shape parameter nu, default is 1. `param` Specifying the parameters as a vector of the form `c(mu,delta,beta,nu)`. `log,log.p` Logical; if `log = TRUE`, probabilities are given as log(p). `method` Character. If `"spline"` quantiles are found from a spline approximation to the distribution function. If `"integrate"`, the distribution function used is always obtained by integration. `lower.tail` Logical. If `lower.tail = TRUE`, the cumulative density is taken from the lower tail. `tolerance` Specified level of tolerance when checking if parameter `beta` is equal to 0. `subdivisions` The maximum number of subdivisions used to integrate the density and determine the accuracy of the distribution function calculation. `intTol` Value of `rel.tol` and hence `abs.tol` in calls to `integrate`. See `integrate`. `valueOnly` Logical. If `valueOnly = TRUE` calls to `pskewhyp` only return the value obtained for the integral. If `valueOnly = FALSE` an estimate of the accuracy of the numerical integration is also returned. `nInterpol` Number of points used in `qskewhyp` for cubic spline interpolation of the distribution function. `uniTol` Value of `tol` in calls to `uniroot`. See `uniroot`. `...` Passes additional arguments to `integrate` in `pskewhyp` and `qskewhyp`, and to `uniroot` in `qskewhyp`.

## Details

Users may either specify the values of the parameters individually or as a vector. If both forms are specified, then the values specified by the vector `param` will overwrite the other ones. In addition the parameter values are examined by calling the function `skewhypCheckPars` to see if they are valid.

The density function is

f(x) = (2^((1-nu)/2) delta^nu abs(beta)^((nu+1)/2) K_((nu+1)/2)(sqrt(beta^2 (delta^2+(x-mu)^2)) ) exp(beta (x-mu)))/ (gamma(nu/2) sqrt(pi) (sqrt(delta^2+(x-mu)^2))^((nu+1)/2))

when beta != 0, and

f(x)=(gamma((nu+1)/2)/(sqrt(pi) delta gamma(nu/2)))(1+((x-mu)^2)/(delta^2))^(-(nu+1)/2)

when beta = 0, where K_nu(.) is the modified Bessel function of the third kind with order nu, and gamma(.) is the gamma function.

`pskewhyp` uses the function `integrate` to numerically integrate the density function. The integration is from `-Inf` to `x` if `x` is to the left of the mode, and from `x` to `Inf` if `x` is to the right of the mode. The probability calculated this way is subtracted from 1 if required. Integration in this manner appears to make calculation of the quantile function more stable in extreme cases.

Calculation of quantiles using `qhyperb` permits the use of two different methods. Both methods use `uniroot` to find the value of x for which a given q is equal F(x) where F denotes the cumulative distribution function. The difference is in how the numerical approximation to F is obtained. The obvious and more accurate method is to calculate the value of F(x) whenever it is required using a call to `phyperb`. This is what is done if the method is specified as `"integrate"`. It is clear that the time required for this approach is roughly linear in the number of quantiles being calculated. A Q-Q plot of a large data set will clearly take some time. The alternative (and default) method is that for the major part of the distribution a spline approximation to F(x) is calculated and quantiles found using `uniroot` with this approximation. For extreme values (for which the tail probability is less than 10^(-7)), the integration method is still used even when the method specifed is `"spline"`.

If accurate probabilities or quantiles are required, tolerances (`intTol` and `uniTol`) should be set to small values, say 10^(-10) or 10^(-12) with ```method = "integrate"```. Generally then accuracy might be expected to be at least 10^(-9). If the default values of the functions are used, accuracy can only be expected to be around 10^(-4). Note that on 32-bit systems `.Machine\$double.eps^0.25 = 0.0001220703` is a typical value.

Note that when small values of nu are used, and the density is skewed, there are often some extreme values generated by `rskewhyp`. These look like outliers, but are caused by the heaviness of the skewed tail, see Examples.

The extreme skewness of the distribution when beta is large in absolute value and nu is small make this distribution very challenging numerically.

## Value

`dskewhyp` gives the density function, `pskewhyp` gives the distribution function, `qskewhyp` gives the quantile function and `rskewhyp` generates random variates.

An estimate of the accuracy of the approximation to the distribution function can be found by setting `valueOnly = FALSE` in the call to `pskewyhp` which returns a list with components `value` and `error`.

`ddskewhyp` gives the derivative of `dskewhyp`.

## Author(s)

David Scott [email protected], Fiona Grimson

## References

Aas, K. and Haff, I. H. (2006). The Generalised Hyperbolic Skew Student's t-distribution, Journal of Financial Econometrics, 4, 275–309.

`safeIntegrate`, `integrate` for its shortfalls, `skewhypCheckPars`, `logHist`. Also `skewhypMean` for information on moments and mode, and `skewhypFit` for fitting to data.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43``` ```param <- c(0,1,40,10) par(mfrow = c(1,2)) range <- skewhypCalcRange(param = param, tol = 10^(-2)) ### curves of density and distribution curve(dskewhyp(x, param = param), range[1], range[2], n = 1000) title("Density of the \n Skew Hyperbolic Distribution") curve(pskewhyp(x, param = param), range[1], range[2], n = 500) title("Distribution Function of the \n Skew Hyperbolic Distribution") ### curves of density and log density par(mfrow = c(1,2)) data <- rskewhyp(1000, param = param) curve(dskewhyp(x, param = param), range(data)[1], range(data)[2], n = 1000, col = 2) hist(data, freq = FALSE, add = TRUE) title("Density and Histogram of the\n Skew Hyperbolic Distribution") logHist(data, main = "Log-Density and Log-Histogram of\n the Skew Hyperbolic Distribution") curve(dskewhyp(x, param = param, log = TRUE), range(data)[1], range(data)[2], n = 500, add = TRUE, col = 2) ##plots of density and derivative par(mfrow = c(2,1)) curve(dskewhyp(x, param = param), range[1], range[2], n = 1000) title("Density of the Skew\n Hyperbolic Distribution") curve(ddskewhyp(x, param = param), range[1], range[2], n = 1000) title("Derivative of the Density\n of the Skew Hyperbolic Distribution") ##example of density and random numbers for beta large and nu small par(mfrow = c(1,2)) param1 <- c(0,1,10,1) data1 <- rskewhyp(1000, param = param1) curve(dskewhyp(x, param = param1), range(data1)[1], range(data1)[2], n = 1000, col = 2) hist(data1, freq = FALSE, add = TRUE) title("Density and Histogram\n when nu is small") logHist(data1, main = "Log-Density and Log-Histogram\n when nu is small") curve(dskewhyp(x, param = param1, log = TRUE), range(data1)[1], range(data1)[2], n = 500, add = TRUE, col = 2) ```