# Function to generate random observations from a # hyperbolic distribution using the # ratio of uniforms method rhypRoU <- 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) hyperbBeta <- zeta * hyperbPi quasiDens <- function(x,alpha,hyperbBeta){ quasiDens <- exp(-alpha*sqrt(1+x^2)+hyperbBeta*x) } vfun <- function(x,alpha,hyperbBeta){ vfun <- -alpha*x/sqrt(1+x^2)+hyperbBeta } vplusloc <- uniroot(function(x)vfun(x,alpha,hyperbBeta), lower=lower,upper=upper)\$root vplus <- sqrt(quasiDens(vplusloc,alpha,hyperbBeta)) ufun <- function(x,alpha,hyperbBeta){ ufun <- x*exp(-(1/2)*alpha*sqrt(1+x^2)+(1/2)*hyperbBeta*x) } dufun <- function(x,alpha,hyperbBeta){ dufun <- (-(1/2)*exp(-alpha*sqrt(1+x^2)/2+hyperbBeta*x/2)* (-2*sqrt(1+x^2)+alpha*x*x-x*hyperbBeta*sqrt(1+x^2))/ sqrt(1+x^2)) } uminusloc <- uniroot(function(x)dufun(x,alpha,hyperbBeta), lower=lower,upper=0)\$root uplusloc <- uniroot(function(x)dufun(x,alpha,hyperbBeta), lower=0,upper=upper)\$root uminus <- uminusloc*sqrt(quasiDens(uminusloc,alpha,hyperbBeta)) uplus <- uplusloc*sqrt(quasiDens(uplusloc,alpha,hyperbBeta)) # Make sure the maximum and minimum are the right way around if (uminus > uplus){ temp <- uminus uminus <- uplus uplus <- temp } output <- numeric(n) for(i in 1:n){ need.value <- TRUE while(need.value==TRUE){ U <- uminus+(uplus-uminus)*runif(1) V <- vplus*runif(1) X <- U/V fX <- quasiDens(X,alpha,hyperbBeta) if(V^2