R/momIntegrated.R In sjp/GeneralizedHyperbolic: The generalized hyperbolic distribution

Documented in momIntegrated

momIntegrated <- function(densFn = c("ghyp", "hyperb", "gig", "gamma", "invgamma", "vg"),
order, param = NULL, about = 0,
absolute = FALSE) {

if (missing(densFn) | !(is.function(densFn) | is.character(densFn)))
stop("'densFn' must be supplied as a function or name")

## Set default integration limits
low <- -Inf
high <- Inf

if (is.character(densFn)) {
distname <- match.arg(densFn)

if (is.null(densFn))
stop("unsupported distribution")

if (distname == "ghyp" | distname == "generalized hyperbolic") {
if (is.null(param))
param <- c(0, 1, 1, 0, 1)

if (!absolute) {
ddist <- function(x, order, param, about) {
(x - about)^order * dghyp(x, param = param)
}
} else {
ddist <- function(x, order, param, about) {
abs(x - about)^order * dghyp(x, param = param)
}
}
} else if (distname == "hyperb" | distname == "hyperbolic") {
if (is.null(param))
param <- c(0, 1, 1, 0)

if (!absolute) {
ddist <- function(x, order, param, about) {
(x - about)^order * dhyperb(x, param = param)
}
} else {
ddist <- function(x, order, param, about) {
abs(x - about)^order * dhyperb(x, param = param)
}
}
} else  if (distname == "gig" |
distname == "generalized inverse gaussian") {
if (is.null(param))
param <- c(1, 1, 1)

low <- 0

if (!absolute) {
ddist <- function(x, order, param, about) {
(x - about)^order * dgig(x, param = param)
}
} else {
ddist <- function(x, order, param, about) {
abs(x - about)^order * dgig(x, param = param)
}
}
} else if (distname == "gamma") {
if (order <= -param[1])
stop("Order must be greater than shape parameter for gamma")

if (is.null(param))
param <- c(1, 1)

low <- 0

if (!absolute) {
ddist <- function(x, order, param, about) {
dgamma(x, shape = param[1], rate = param[2])
}
} else {
ddist <- function(x, order, param, about) {
dgamma(x, shape = param[1], rate = param[2])
}
}
} else  if (distname == "invgamma" |
distname == "inverse gamma") {
if (is.null(param))
param <- c(-1, 1)

if (param[1] <= order)
stop("Order must be less than shape parameter for inverse gamma")

low <- 0

dinvgamma <- function(x, shape, rate = 1, scale = 1/rate) {
dens <- ifelse(x <= 0, 0,
(scale / x)^shape * exp(-scale / x) / (x * gamma(shape)))
return(dens)
}

if (!absolute) {
ddist <- function(x, order, param, about) {
dinvgamma(x, shape = param[1], rate = param[2])
}
} else {
ddist <- function(x, order, param, about) {
dinvgamma(x, shape = param[1], rate = param[2])
}
}
} else if (distname == "vg" | distname == "variance gamma") {
if (!exists("dvg", mode = "function"))
stop("package Variance Gamma needs to be loaded")

if (is.null(param))
param <- c(0, 1, 0, 1)

if (!absolute) {
ddist <- function(x, order, param, about) {
(x - about)^order * dvg(x, param = param)
}
} else {
ddist <- function(x, order, param, about) {
abs(x - about)^order * dvg(x, param = param)
}
}
}
} else {
if (is.function(densFn))
stop("general density functions code not yet implemented")
}

mom <- integrate(ddist, low, high, param = param,