1
2# This library is free software; you can redistribute it and/or
3# modify it under the terms of the GNU Library General Public
4# License as published by the Free Software Foundation; either
5# version 2 of the License, or (at your option) any later version.
6#
7# This library is distributed in the hope that it will be useful,
8# but WITHOUT ANY WARRANTY; without even the implied warranty of
9# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10# GNU Library General Public License for more details.
11#
12# You should have received a copy of the GNU Library General
13# Public License along with this library; if not, write to the
14# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
15# MA  02111-1307  USA
16
17
18################################################################################
19# FUNCTION:             DESCRIPTION:
20#  dgh                   Returns density for generalized hyperbolic DF
21#  pgh                   Returns probability for generalized hyperbolic DF
22#  qgh                   Returns quantiles for generalized hyperbolic DF
23#  rgh                   Returns random variates for generalized hyperbolic DF
24################################################################################
25
26
27dgh <-
28function(x, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = -1/2, log = FALSE)
29{
30    # A function implemented by Diethelm Wuertz
31
32    # Description:
33    #   Returns density for the generalized hyperbolic distribution
34
35    # FUNCTION:
36
37    # Parameters:
38    if (length(alpha) == 4) {
39       mu = alpha[4]
40       delta = alpha[3]
41       beta = alpha[2]
42       alpha = alpha[1]
43    }
44
45    # Checks:
46    if (alpha <= 0) stop("alpha must be greater than zero")
47    if (delta <= 0) stop("delta must be greater than zero")
48    if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha")
49
50    # Density:
51    arg = delta*sqrt(alpha^2-beta^2)
52    a = (lambda/2)*log(alpha^2-beta^2) - (
53        log(sqrt(2*pi)) + (lambda-0.5)*log(alpha) + lambda*log(delta) +
54        log(besselK(arg, lambda, expon.scaled = TRUE)) - arg )
55    f = ((lambda-0.5)/2)*log(delta^2+(x - mu)^2)
56
57    # Use exponential scaled form to prevent from overflows:
58    arg = alpha * sqrt(delta^2+(x-mu)^2)
59    k = log(besselK(arg, lambda-0.5, expon.scaled = TRUE)) - arg
60    e = beta*(x-mu)
61
62    # Put all together:
63    ans = a + f + k + e
64    if(!log) ans = exp(ans)
65
66    # Return Value:
67    ans
68}
69
70
71# ------------------------------------------------------------------------------
72
73
74pgh <-
75function(q, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = -1/2)
76{
77    # A function implemented by Diethelm Wuertz
78
79    # Description:
80    #   Returns probability for the generalized hyperbolic distribution
81
82    # FUNCTION:
83
84    # Checks:
85    if (alpha <= 0) stop("alpha must be greater than zero")
86    if (delta <= 0) stop("delta must be greater than zero")
87    if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha")
88
89    # Probability:
90    ans = NULL
91    for (Q in q) {
92        Integral = integrate(dgh, -Inf, Q, stop.on.error = FALSE,
93            alpha = alpha, beta = beta, delta = delta, mu = mu,
94            lambda = lambda)
95        ans = c(ans, as.numeric(unlist(Integral)[1]))
96    }
97
98    # Return Value:
99    ans
100}
101
102
103# ------------------------------------------------------------------------------
104
105
106qgh <-
107function(p, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = -1/2)
108{
109    # A function implemented by Diethelm Wuertz
110
111    # Description:
112    #   Returns quantiles for the generalized hyperbolic distribution
113
114    # FUNCTION:
115
116    # Checks:
117    if (alpha <= 0) stop("alpha must be greater than zero")
118    if (delta <= 0) stop("delta must be greater than zero")
119    if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha")
120
121    # Internal Function:
122    .froot <- function(x, alpha, beta, delta, mu, lambda, p)
123    {
124        pgh(q = x, alpha = alpha, beta = beta, delta = delta,
125            mu = mu, lambda = lambda) - p
126    }
127
128    # Quantiles:
129    result = NULL
130    for (pp in p) {
131        lower = -1
132        upper = +1
133        counter = 0
134        iteration = NA
135        while (is.na(iteration)) {
136            iteration = .unirootNA(f = .froot, interval = c(lower,
137                upper), alpha = alpha, beta = beta, delta = delta,
138                mu = mu, lambda = lambda, p = pp)
139            counter = counter + 1
140            lower = lower - 2^counter
141            upper = upper + 2^counter
142        }
143        result = c(result, iteration)
144    }
145
146    # Return Value:
147    result
148}
149
150
151# ------------------------------------------------------------------------------
152
153
154rgh <-
155function(n, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = -1/2)
156{
157    # A function implemented by Diethelm Wuertz
158
159    # Description:
160    #   Returns random variates for the generalized hyperbolic distribution
161
162    # FUNCTION:
163
164    # Checks:
165    if (alpha <= 0) stop("alpha must be greater than zero")
166    if (delta <= 0) stop("delta must be greater than zero")
167    if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha")
168
169    # Settings:
170    theta = c(lambda, alpha, beta, delta, mu)
171
172    # Random Numbers:
173    ans = .rghyp(n, theta)
174
175    # Attributes:
176    attr(ans, "control") = c(dist = "gh", alpha = alpha, beta = beta,
177    delta = delta, mu = mu, lambda = lambda)
178
179    # Return Value:
180    ans
181}
182
183
184################################################################################
185
186