# Function to generate random observations from a # hyperbolic distribution using the # mixing property of the generalized inverse # Gaussian distribution and Dagpunar's algorithm # for the generalized inverse Gaussian rhypJD <- function(n,Theta){ hyperbPi <- Theta[1] zeta <- Theta[2] delta <- Theta[3] mu <- Theta[4] if(zeta<=0) stop("zeta must be positive") if(delta<=0) stop("delta must be positive") lowerupper <- hyperbCalcRange(c(hyperbPi,zeta,1,0)) lower <- lowerupper[1] upper <- lowerupper[2] alpha <- zeta * sqrt(1 + hyperbPi^2)/delta hyperbBeta <- zeta * hyperbPi/delta chi <- delta^2 psi <- alpha^2 - hyperbBeta^2 gigAlpha <- sqrt(psi/chi) gigBeta <- sqrt(psi*chi) m <- abs(gigBeta)/gigBeta g <- function(y){ 0.5*gigBeta*y^3 - y^2*(0.5*gigBeta*m+2) + y*(-0.5*gigBeta) + 0.5*gigBeta*m } upper <- m while(g(upper)<=0) upper <- 2*upper yM <- uniroot(g,interval=c(0,m))\$root yP <- uniroot(g,interval=c(m,upper))\$root a <- (yP-m)*exp(-0.25*gigBeta*(yP+1/yP-m-1/m)) b <- (yM-m)*exp(-0.25*gigBeta*(yM+1/yM-m-1/m)) c <- -0.25*gigBeta*(m+1/m) X <- numeric(n) for(i in 1:n){ need.value <- TRUE while(need.value==TRUE){ R1 <- runif(1) R2 <- runif(1) Y <- m + a*R2/R1 + b*(1-R2)/R1 if(Y>0){ if(-log(R1)>=0.25*gigBeta*(Y+1/Y)+c){ need.value <- FALSE } } } X[i] <- Y/gigAlpha } sigma <- sqrt(X) Z <- rnorm(n) output <- mu + hyperbBeta*sigma^2 + sigma*Z output }