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.

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)
``` |

`x,q` |
Vector of quantiles. |

`p` |
Vector of probabilities. |

`n` |
Number of random variates to be generated. |

`mu` |
Location parameter |

`delta` |
Scale parameter |

`beta` |
Skewness parameter |

`nu` |
Shape parameter |

`param` |
Specifying the parameters as a vector of the form |

`log,log.p` |
Logical; if |

`method` |
Character. If |

`lower.tail` |
Logical. If |

`tolerance` |
Specified level of tolerance when checking if
parameter |

`subdivisions` |
The maximum number of subdivisions used to integrate the density and determine the accuracy of the distribution function calculation. |

`intTol` |
Value of |

`valueOnly` |
Logical. If |

`nInterpol` |
Number of points used in |

`uniTol` |
Value of |

`...` |
Passes additional arguments to |

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.

`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`

.

David Scott d.scott@auckland.ac.nz, Fiona Grimson

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)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.