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