R/dghyp.R

In sjp/GeneralizedHyperbolic: The generalized hyperbolic distribution

### Function to calculate the density of the
### generalized hyperbolic distribution
dghyp <- function(x, mu = 0, delta = 1, alpha = 1, beta = 0,
                  lambda = 1, param = c(mu, delta, alpha, beta, lambda)) {

  # Lambda defaults to one if omitted from param vector
  if (length(param) == 4)
    param <- c(param, 1)

  ## check parameters
  parResult <- ghypCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  param <- as.numeric(param)
  mu <- param[1]
  delta <- param[2]
  alpha <- param[3]
  beta <- param[4]
  lambda <- param[5]
  gamma <- sqrt(alpha^2 - beta^2)

  ## Argument of Bessel K function in numerator
  y <- alpha * sqrt(delta^2 + (x - mu)^2)
  bx <- beta * (x - mu)

  ## Deal with underflow in ratio of Bessel K functions
  ## besselK underflows for x > 740
  ## Use exponentially scaled besselK
  if (delta * gamma > 700) {
    # underflow in constant part
    expTerm <- exp(delta * gamma - y + bx)
    besselRatio <- besselK(x = y, nu = lambda - 1 / 2, expon.scaled = TRUE) /
                           besselK(x = delta * gamma, nu = lambda, expon.scaled = TRUE)
    expAndBessel <- expTerm * besselRatio
  } else {
    expAndBessel <- ifelse(y > 700 | bx > 700, # underflow in variable part
                           exp(delta * gamma - y + bx) *
                           besselK(x = y, nu = lambda - 1 / 2,
                                   expon.scaled = TRUE) /
                           besselK(x = delta * gamma, nu = lambda,
                                   expon.scaled = TRUE),
                          exp(bx) * besselK(x = y, nu = lambda - 1 / 2) /
                          besselK(x = delta * gamma, nu = lambda))
  }

  dens <- (y / alpha)^(lambda - 1 / 2) * ((gamma / delta)^lambda) *
          alpha^(1 / 2 - lambda) * expAndBessel / sqrt(2 * pi)

  dens
} ## End of dghyp()


### Cumulative distribution function of the generalized hyperbolic
### New version intended to give guaranteed accuracy
### Uses exponential approximation in tails
### SafeIntegrate() is used over six parts in the middle of the distribution
### Calls ghypBreaks to determine the breaks
###
### DJS 07/01/07
pghyp <- function(q, mu = 0, delta = 1, alpha = 1, beta = 0, lambda = 1,
                  param = c(mu, delta, alpha, beta, lambda), small = 10^(-6),
                  tiny = 10^(-10), deriv = 0.3, subdivisions = 100,
                  accuracy = FALSE, ...) {

  # Lambda defaults to one if omitted from param vector
  if (length(param) == 4)
    param <- c(param, 1)

  ## check parameters
  parResult <- ghypCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  param <- as.numeric(param)
  mu <- param[1]
  delta <- param[2]
  alpha <- param[3]
  beta <- param[4]
  lambda <- param[5]

  bks <- ghypBreaks(param = param, small = small, tiny = tiny, deriv = deriv, ...)
  xTiny <- bks$xTiny
  xSmall <- bks$xSmall
  lowBreak <- bks$lowBreak
  highBreak <- bks$highBreak
  xLarge <- bks$xLarge
  xHuge <- bks$xHuge
  modeDist <- bks$modeDist

  qSort <- sort(q)
  qTiny <- which(qSort < xTiny)
  qSmall <- which(qSort < xSmall)
  qLow <- which(qSort < lowBreak)
  qLessEqMode <- which(qSort <= modeDist)
  qGreatMode <- which(qSort > modeDist)
  qHigh <- which(qSort > highBreak)
  qLarge <- which(qSort > xLarge)
  qHuge <- which(qSort > xHuge)

  ## Break indices into 8 groups: beware of empty groups
  if (length(qLow) > 0) qLessEqMode <- qLessEqMode[qLessEqMode > max(qLow)]
  if (length(qHigh) > 0) qGreatMode <- qGreatMode[qGreatMode < min(qHigh)]
  if (length(qSmall) > 0) qLow <- qLow[qLow > max(qSmall)]
  if (length(qLarge) > 0) qHigh <- qHigh[qHigh < min(qLarge)]
  if (length(qTiny) > 0) qSmall <- qSmall[qSmall > max(qTiny)]
  if (length(qHuge) > 0) qLarge <- qLarge[qLarge < min(qHuge)]
  intFun <- rep(NA, length(q))
  if (length(qTiny) > 0) intFun[qTiny] <- 0
  if (length(qHuge) > 0) intFun[qHuge] <- 1
  intErr <- rep(NA, length(q))
  if (length(qTiny) > 0) intErr[qTiny] <- tiny
  if (length(qHuge) > 0) intErr[qHuge] <- tiny


  ## Use safeIntegrate function between xTiny and xHuge in 6 sections
  dghypInt <- function(q) {
    dghyp(q, param = param)
  }

  ## Calculate integrals and errors to cut points
  resSmall <- safeIntegrate(dghypInt, xTiny, xSmall, subdivisions, ...)
  resLarge <- safeIntegrate(dghypInt, xLarge, xHuge, subdivisions, ...)
  intSmall <- resSmall$value
  intLarge <- resLarge$value
  errSmall <- tiny + resSmall$abs.error
  errLarge <- tiny + resLarge$abs.error
  resLow <- safeIntegrate(dghypInt, xSmall, lowBreak, subdivisions, ...)
  resHigh <- safeIntegrate(dghypInt, highBreak, xLarge, subdivisions, ...)
  intLow <- intSmall + resLow$value
  intHigh <- intLarge + resHigh$value
  errLow <- errSmall + resLow$abs.error
  errHigh <- errLarge + resHigh$abs.error

  for (i in qSmall) {
    intRes <- safeIntegrate(dghypInt, xTiny, qSort[i], subdivisions, ...)
    intFun[i] <- intRes$value
    intErr[i] <- intRes$abs.error + tiny
  }

  for (i in qLarge) {
    intRes <- safeIntegrate(dghypInt, qSort[i], xHuge, subdivisions, ...)
    intFun[i] <-  1- intRes$value
    intErr[i] <- intRes$abs.error + tiny
  }

  for (i in qLow) {
    intRes <- safeIntegrate(dghypInt, xSmall, qSort[i], subdivisions, ...)
    intFun[i] <- intRes$value + intSmall
    intErr[i] <- intRes$abs.error + errSmall
  }

  for (i in qHigh) {
    intRes <- safeIntegrate(dghypInt, qSort[i], xLarge, subdivisions, ...)
    intFun[i] <-  1- intRes$value - intLarge
    intErr[i] <- intRes$abs.error + errLarge
  }

  for (i in qLessEqMode) {
    intRes <- safeIntegrate(dghypInt, lowBreak, qSort[i], subdivisions, ...)
    intFun[i] <- intRes$value + intLow
    intErr[i] <- intRes$abs.error + errLow
  }

  for (i in qGreatMode) {
    intRes <- safeIntegrate(dghypInt, qSort[i], highBreak, subdivisions, ...)
    intFun[i] <-  1- intRes$value - intHigh
    intErr[i] <- intRes$abs.error + errLarge
  }

  if (!accuracy) {
    return(intFun[rank(q)])
  } else {
    return(list(value = intFun[rank(q)], error = intErr[rank(q)]))
  }
} ## End of pghyp()

### Cumulative distribution function of the generalized hyperbolic
### New version intended to give guaranteed accuracy
### qghyp using breaks as for pghyp and splines as in original qghyp
###
### DJS 06/09/06
qghyp <- function(p, mu = 0, delta = 1, alpha = 1, beta = 0, lambda = 1,
                  param = c(mu, delta, alpha, beta, lambda), small = 10^(-6),
                  tiny = 10^(-10), deriv = 0.3, nInterpol = 100,
                  subdivisions = 100, ...) {

  # Lambda defaults to one if omitted from param vector
  if (length(param) == 4)
    param <- c(param, 1)

  ## check parameters
  parResult <- ghypCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  param <- as.numeric(param)
  mu <- param[1]
  delta <- param[2]
  alpha <- param[3]
  beta <- param[4]
  lambda <- param[5]

  bks <- ghypBreaks(param = param, small = small, tiny = tiny, deriv = deriv, ...)
  xTiny <- bks$xTiny
  xSmall <- bks$xSmall
  lowBreak <- bks$lowBreak
  highBreak <- bks$highBreak
  xLarge <- bks$xLarge
  xHuge <- bks$xHuge
  modeDist <- bks$modeDist

  yTiny <- pghyp(xTiny, param = param)
  ySmall <- pghyp(xSmall, param = param)
  yLowBreak <- pghyp(lowBreak, param = param)
  yHighBreak <- pghyp(highBreak, param = param)
  yLarge <- pghyp(xLarge, param = param)
  yHuge <- pghyp(xHuge, param = param)
  yModeDist <- pghyp(modeDist, param = param)

  pSort <- sort(p)
  pSmall <- which(pSort < pghyp(xSmall, param = param))
  pTiny <- which(pSort < pghyp(xTiny, param = param))
  pLarge <- which(pSort > pghyp(xLarge, param = param))
  pHuge <- which(pSort > pghyp(xHuge, param = param))
  pLow <- which(pSort < pghyp(lowBreak, param = param))
  pHigh <- which(pSort > pghyp(highBreak, param = param))
  pLessEqMode <- which(pSort <= pghyp(modeDist, param = param))
  pGreatMode <- which(pSort > pghyp(modeDist, param = param))

  ## Break indices into 8 groups: beware of empty groups
  if (length(pLow) > 0) pLessEqMode <- pLessEqMode[pLessEqMode > max(pLow)]
  if (length(pHigh) > 0) pGreatMode <- pGreatMode[pGreatMode < min(pHigh)]
  if (length(pSmall) > 0) pLow <- pLow[pLow > max(pSmall)]
  if (length(pLarge) > 0) pHigh <- pHigh[pHigh < min(pLarge)]
  if (length(pTiny) > 0) pSmall <- pSmall[pSmall > max(pTiny)]
  if (length(pHuge) > 0) pLarge <- pLarge[pLarge < min(pHuge)]
  qSort <- rep(NA, length(pSort))
  if (length(pTiny) > 0) qSort[pTiny] <- -Inf
  if (length(pHuge) > 0) qSort[pHuge] <- Inf


  if (length(pTiny) > 0) {
    for (i in pTiny) {

      zeroFun <- function(x) {
        pghyp(x, param = param) - pSort[i]
      }

      interval <- c(xTiny - (xSmall - xTiny), xTiny)
      while (zeroFun(interval[1]) * zeroFun(interval[2]) > 0) {
        interval[1] <- interval[1] - (xSmall - xTiny)
      }
      qSort[i] <- uniroot(zeroFun, interval)$root
    }
  }

  if (length(pSmall) > 0) {
    xValues <- seq(xTiny, xSmall, length = nInterpol)
    pghypValues <- pghyp(xValues, param = param, small = small, tiny = tiny, deriv = deriv,
                         subdivisions = subdivisions, accuracy = FALSE)
    pghypSpline <- splinefun(xValues, pghypValues)

    for (i in pSmall) {

      zeroFun <- function(x) {
        pghypSpline(x) - pSort[i]
      }

      if (zeroFun(xTiny) >= 0) {
        qSort[i] <- xTiny
      } else {
        if (zeroFun(xSmall) <= 0) {
          qSort[i] <- xSmall
        } else {
          qSort[i] <- uniroot(zeroFun, interval = c(xTiny, xSmall), ...)$root
        }
      }
    }
  }

  if (length(pLow) > 0) {
    xValues <- seq(xSmall, lowBreak, length = nInterpol)
    pghypValues <- pghyp(xValues, param = param, small = small, tiny = tiny, deriv = deriv,
                         subdivisions = subdivisions, accuracy = FALSE)
    pghypSpline <- splinefun(xValues, pghypValues)

    for (i in pLow) {

      zeroFun <- function(x) {
        pghypSpline(x) - pSort[i]
      }

      if (zeroFun(xSmall) >= 0) {
        qSort[i] <- xSmall
      } else {
        if (zeroFun(lowBreak) <= 0) {
          qSort[i] <- lowBreak
        } else {
          qSort[i] <- uniroot(zeroFun, interval = c(xSmall, lowBreak), ...)$root
        }
      }
    }
  }

  if (length(pLessEqMode) > 0) {
    xValues <- seq(lowBreak, modeDist, length = nInterpol)
    pghypValues <- pghyp(xValues, param = param, small = small, tiny = tiny, deriv = deriv,
                         subdivisions = subdivisions, accuracy = FALSE)
    pghypSpline <- splinefun(xValues, pghypValues)

    for (i in pLessEqMode) {

      zeroFun <- function(x) {
        pghypSpline(x) - pSort[i]
      }

      if (zeroFun(lowBreak) >= 0) {
        qSort[i] <- lowBreak
      } else {
        if (zeroFun(modeDist) <= 0) {
          qSort[i] <- modeDist
        } else {
          qSort[i] <- uniroot(zeroFun, interval = c(lowBreak, modeDist), ...)$root
        }
      }
    }
  }

  if (length(pGreatMode) > 0) {
    xValues <- seq(modeDist, highBreak, length = nInterpol)
    pghypValues <- pghyp(xValues, param = param, small = small, tiny = tiny, deriv = deriv,
                         subdivisions = subdivisions, accuracy = FALSE)
    pghypSpline <- splinefun(xValues, pghypValues)

    for (i in pGreatMode) {

      zeroFun <- function(x) {
        pghypSpline(x) - pSort[i]
      }

      if (zeroFun(modeDist) >= 0) {
        qSort[i] <- modeDist
      } else {
        if (zeroFun(highBreak) <= 0) {
          qSort[i] <- highBreak
        } else {
          qSort[i] <- uniroot(zeroFun, interval = c(modeDist, highBreak), ...)$root
        }
      }
    }
  }

  if (length(pHigh) > 0) {
    xValues <- seq(highBreak, xLarge, length = nInterpol)
    pghypValues <- pghyp(xValues, param = param, small = small, tiny = tiny, deriv = deriv,
                         subdivisions = subdivisions, accuracy = FALSE)
    pghypSpline <- splinefun(xValues, pghypValues)

    for (i in pHigh) {

      zeroFun <- function(x) {
        pghypSpline(x) - pSort[i]
      }

      if (zeroFun(highBreak) >= 0) {
        qSort[i] <- highBreak
      } else {
        if (zeroFun(xLarge) <= 0) {
          qSort[i] <- xLarge
        } else {
          qSort[i] <- uniroot(zeroFun, interval = c(highBreak,xLarge), ...)$root
        }
      }
    }
  }

  if (length(pLarge) > 0) {
    xValues <- seq(xLarge, xHuge, length = nInterpol)
    pghypValues <- pghyp(xValues, param = param, small = small, tiny = tiny, deriv = deriv,
                         subdivisions = subdivisions, accuracy = FALSE)
    pghypSpline <- splinefun(xValues, pghypValues)

    for (i in pLarge) {

      zeroFun <- function(x) {
        pghypSpline(x) - pSort[i]
      }

      if (zeroFun(xLarge) >= 0) {
        qSort[i] <- xLarge
      } else {
        if (zeroFun(xHuge) <= 0) {
          qSort[i] <- xHuge
        } else {
          qSort[i] <- uniroot(zeroFun, interval = c(xLarge,xHuge), ...)$root
        }
      }
    }
  }

  if (length(pHuge) > 0) {
    for (i in pHuge) {

      zeroFun <- function(x) {
        pghyp(x, param = param) - pSort[i]
      }

      interval <- c(xHuge, xHuge + (xHuge - xLarge))

      while (zeroFun(interval[1]) * zeroFun(interval[2]) > 0) {
        interval[1] <- interval[1] + (xHuge - xLarge)
      }
      qSort[i] <- uniroot(zeroFun, interval)$root
    }
  }

  return(qSort[rank(p)])
} # End of qghyp()

### Derivative of the density
ddghyp <- function(x, mu = 0, delta = 1, alpha = 1, beta = 0, lambda = 1,
                   param = c(mu, delta, alpha, beta, lambda)) {

  # Lambda defaults to one if omitted from param vector
  if (length(param) == 4)
    param <- c(param, 1)

  ## check parameters
  parResult <- ghypCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  param <- as.numeric(param)
  mu <- param[1]
  delta <- param[2]
  alpha <- param[3]
  beta <- param[4]
  lambda <- param[5]

  ## Terms for simplification of programming
  t1 <- sqrt(delta^2 + (x - mu)^2)
  t2 <- sqrt(alpha^2 - beta^2)
  t3 <- besselK(x = alpha * t1, nu = lambda - 0.5)
  t4 <- besselK(x = alpha * t1, nu = lambda + 0.5)

  ddghyp <- (t3 * (beta * delta^2 + (2 * lambda - 1) * (x - mu) + beta * (x - mu)^2) -
             t4 * alpha * t1 * (x - mu)) *
             exp(beta * (x - mu)) * t1^(lambda - (5 / 2)) * t2^lambda /
             (sqrt(2 * pi) * alpha^(lambda -  1 / 2) * delta^lambda *
             besselK(x = delta * t2, nu = lambda))

  ddghyp
} ## End of ddghyp()

### Function to set up breaks for pghyp and qghyp
ghypBreaks <- function(mu = 0, delta = 1, alpha = 1, beta = 0, lambda = 1,
                       param = c(mu, delta, alpha, beta, lambda),
                       small = 10^(-6), tiny = 10^(-10), deriv = 0.3, ...) {

  # Lambda defaults to one if omitted from param vector
  if (length(param) == 4)
    param <- c(param, 1)

  ## check parameters
  parResult <- ghypCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  mu <- param[1]
  delta <- param[2]
  alpha <- param[3]
  beta <- param[4]
  lambda <- param[5]

  xTiny <- ghypCalcRange(param = param, tol = tiny, density = TRUE)[1]
  xSmall <- ghypCalcRange(param = param, tol = small, density = TRUE)[1]
  xLarge <- ghypCalcRange(param = param, tol = small, density = TRUE)[2]
  xHuge <- ghypCalcRange(param = param, tol = tiny, density = TRUE)[2]

  modeDist <- ghypMode(param = param)
  ## Determine break points, based on size of derivative
  xDeriv <- seq(xSmall, modeDist, length.out = 101)
  derivVals <- ddghyp(xDeriv, param = param)
  maxDeriv <- max(derivVals)
  minDeriv <- min(derivVals)
  breakSize <- deriv * maxDeriv

  breakFun <- function(x) {
    ddghyp(x, param = param) - breakSize
  }

  if ((maxDeriv < breakSize) | (derivVals[1] > breakSize)) {
    lowBreak <- xSmall
  } else {
    whichMaxDeriv <- which.max(derivVals)
    lowBreak <- uniroot(breakFun, c(xSmall, xDeriv[whichMaxDeriv]))$root
  }

  xDeriv <- seq(modeDist, xLarge, length.out = 101)
  derivVals <- -ddghyp(xDeriv, param = param)
  maxDeriv <- max(derivVals)
  minDeriv <- min(derivVals)
  breakSize <- deriv * maxDeriv

  breakFun <- function(x) {
    -ddghyp(x, param = param) - breakSize
  }

  if ((maxDeriv < breakSize) | (derivVals[101] > breakSize)) {
    highBreak <- xLarge
  } else {
    whichMaxDeriv <- which.max(derivVals)
    highBreak <- uniroot(breakFun, c(xDeriv[whichMaxDeriv], xLarge))$root
  }

  breaks <- c(xTiny, xSmall, lowBreak, highBreak, xLarge, xHuge, modeDist)
  breaks <- list(xTiny = breaks[1], xSmall = breaks[2],
                 lowBreak = breaks[3], highBreak = breaks[4],
                 xLarge = breaks[5], xHuge = breaks[6],
                 modeDist = breaks[7])
  return(breaks)
} ## End of ghypBreaks()


### Function to generate random observations from a
### (generalized) hyperbolic distribution using the
### mixing property of the generalized inverse
### Gaussian distribution.
rghyp <- function(n, mu = 0, delta = 1, alpha = 1, beta = 0, lambda = 1,
                  param = c(mu, delta, alpha, beta, lambda)) {

  # Lambda defaults to one if omitted from param vector
  if (length(param) == 4)
    param <- c(param, 1)

  ## check parameters
  parResult <- ghypCheckPars(param)
  case <- parResult$case
  errMessage <- parResult$errMessage

  if (case == "error")
    stop(errMessage)

  param <- as.numeric(param)
  mu <- param[1]
  delta <- param[2]
  alpha <- param[3]
  beta <- param[4]
  lambda <- param[5]

  chi <- delta^2
  psi <- alpha^2 - beta^2

  if (lambda == 1) {
    X <- rgig1(n, param = c(chi, psi, lambda))
  } else {
    X <- rgig(n, param = c(chi, psi, lambda))
  }

  sigma <- sqrt(X)
  Z <- rnorm(n)
  Y <- mu + beta * sigma^2 + sigma * Z

  Y
} ## End of rghyp()