https://github.com/cran/nacopula
Tip revision: f64e14b64fb7985b724358c00e7a2493729e6047 authored by Martin Maechler on 21 September 2011, 00:00:00 UTC
version 0.7-9
version 0.7-9
Tip revision: f64e14b
GIG.R
## Copyright (C) 2010 Marius Hofert and Martin Maechler
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
source(system.file("Rsource", "utils.R", package="nacopula"))
##--> tryCatch.W.E(), canGet()
stopifnot(require(Runuran), require(optimx))
## ==== GIG generator and related functions ====================================
## generator (parameters theta_1 in IR, theta_2 in (0,Inf))
## note: theta_1!=0 and theta_2=0 reduces to copClayton@psi(t, 1/theta_1)
psi.GIG <- function(t, theta)
(1+t)^(-theta[1]/2)*besselK(theta[2]*sqrt(1+t),nu=theta[1])/
besselK(theta[2],nu=theta[1])
## generator inverse
## note: the initial interval is only valid for theta_1 > -1/2
## it can be large, but it's finite; it is based on an inequality about
## modified Bessel functions of the third kind (besselK) given in
## Paris (1984) ("An inequality for the Bessel function J_v(vx)",
## SIAM J. Math. Anal. 15, 203--205)
psiInv.GIG <- function(t, theta, upper=NULL, ...){
if(is.null(upper)){
if(theta[1] > -0.5) upper <- function(x) (1-log(x)/theta[2])^2-1 else
stop("initial interval for psiInv.GIG fails")
}
res <- numeric(length(t))
is0 <- t==0
is1 <- t==1
n01 <- !(is0 | is1)
res[is0] <- Inf
res[is1] <- 0
t. <- t[n01]
up <- upper(t.)
res[n01] <- unlist(lapply(seq_along(t.), function(i){
uniroot(function(t..) psi.GIG(t.., theta)-t.[i],
interval=c(0, up[i]), ...)$root
}))
if(is.matrix(t)) matrix(res, ncol=ncol(t)) else res
}
## generator derivatives
psiDabs.GIG <- function(t, theta, degree=1, n.MC=0, log=FALSE){
res <- numeric(length(t))
iInf <- is.infinite(t)
res[iInf] <- -Inf # log(0)
if(any(!iInf)){
t. <- t[!iInf]
if(n.MC>0){
V <- V0.GIG(n.MC, theta)
res[!iInf] <- nacopula:::lsum(-V %*% t(t.) + degree*log(V) - log(n.MC))
}else{
res[!iInf] <- degree*log(theta[2]/2)-((theta[1]+degree)/2)*log1p(t.) +
log(besselK(theta[2]*sqrt(1+t.), nu=theta[1]+degree, expon.scaled=TRUE)) -
log(besselK(theta[2], nu=theta[1], expon.scaled=TRUE)) -
(sqrt(1+t.)-1)*theta[2]
}
}
r <- if(log) res else exp(res)
if(is.matrix(t)) matrix(r, ncol=ncol(t)) else r
}
## derivative of the generator inverse
psiInvD1abs.GIG <- function(t, theta, log=FALSE){
res <- -psiDabs.GIG(psiInv.GIG(t, theta=theta), theta=theta, log=TRUE)
if(log) res else exp(res)
}
## density of the GIG copula
dacopula.GIG <- function(u, theta, n.MC=0, log=FALSE){
if(!is.matrix(u)) u <- rbind(u)
if((d <- ncol(u)) < 2) stop("u should be at least bivariate") # check that d >= 2
## f() := NaN outside and on the boundary of the unit hypercube
res <- rep.int(NaN, n <- nrow(u))
n01 <- apply(u,1,function(x) all(0 < x, x < 1)) # indices for which density has to be evaluated
if(!any(n01)) return(res)
## auxiliary results
u. <- u[n01,, drop=FALSE]
psiI <- psiInv.GIG(u., theta=theta) # this is costly if d is large
psiI.sum <- rowSums(psiI)
## main part
if(n.MC > 0){ # Monte Carlo
res[n01] <- psiDabs.GIG(psiI.sum, theta=theta, degree=d, n.MC=n.MC, log=TRUE) -
rowSums(psiDabs.GIG(psiI, theta=theta, log=TRUE))
}else{ # explicit
## auxiliary function for density evaluation
h <- function(t, k, theta, log=FALSE){
s1pt <- sqrt(1+t)
B <- besselK(theta[2]*s1pt, nu=theta[1]+k)
if(log) log(B)-(theta[1]+k)*log(s1pt) else B/s1pt^(theta[1]+k)
}
## result
res[n01] <- h(psiI.sum, k=d, theta=theta, log=TRUE) + (d-1)*
log(besselK(theta[2], nu=theta[1])) - rowSums(h(psiI, k=1, theta=theta,
log=TRUE))
}
if(log) res else exp(res)
}
## V0 random number generator
V0.GIG <- function(n, theta){
dens <- udgig(theta[1], 1, theta[2]^2) # the three args lambda, psi, chi for the GIG density as on page 497 in McNeil, Frey, Embrechts (2005)
gen <- pinvd.new(dens) # works fast, via approximation of the quantile function by piecewise polynomials
ur(gen, n)/2
}
## density of V0
dV0.GIG <- function(x, theta, log=FALSE){
dens <- udgig(theta[1], 1, theta[2]^2)
if(log) log(2*ud(dens, 2*x)) else 2*ud(dens, 2*x)
}
## Kendall's tau for fixed theta_1 as a function in theta_2
## quiet==TRUE means that neither warnings nor errors are returned, just NAs
tau.GIG. <- function(theta, upper=c(200,200), quiet=TRUE, ...){
if(theta[1]<0) stop("tau.GIG.: only implemented for nu >= 0")
stopifnot(theta[1]<=upper[1], 0<theta[2], theta[2]<=upper[2])
## determine minimal theta such that integration still works (before
## asymptotic is used)
theta.min <- if(theta[1]<50){
1e-4
}else{
if(theta[1]<100){
0.1
}else{
l10 <- log(10)
exp(l10/50*theta[1]-3*l10)
}
}
## either use heuristic or integration
if(theta[2]<theta.min){
1/(1+2*theta[1])
}else{
tau.integrand <- function(t, theta){
st <- sqrt(1+t)
t*(theta[2]*besselK(theta[2]*st, nu=theta[1]+1)/(st^(theta[1]+1)*
besselK(theta[2], nu=theta[1])))^2
}
int <- tryCatch.W.E(integrate(tau.integrand, lower=0, upper=Inf, theta=theta,
subdivisions=1000, ...)$value)
warn.err <- inherits(int$value, what="simpleError") | inherits(int$value, what="simpleWarning")
if(warn.err && quiet) NA else 1-int$value
}
}
## wrapper for tau.GIG. [vectorized]
tau.GIG <- function(theta, upper=c(200,200), quiet=TRUE, ...){
if(!is.matrix(theta)) theta <- matrix(theta, ncol=2)
apply(theta, 1, tau.GIG., upper=upper, quiet=quiet, ...)
}
## inverse of Kendall's tau for all parameters fixed except the one given as NA
## note: initial interval has to be specified [e.g., c(1e-30,200)] since there
## are no adequate bounds known
tauInv.GIG <- function(tau, theta, upper=c(200,200), quiet=TRUE, iargs=list(), ...){
stopifnot(length(i <- which(is.na(theta))) == 1)
tau.min <- tau.GIG(upper, upper=upper, quiet=quiet)
if(tau<=tau.min) stop(paste("tauInv.GIG: tau must be greater than tau.min=",tau.min,sep=""))
if(i==1){ # wanted: nu; given: theta
interval <- c(0, upper[1])
}else{ # wanted: theta; given: nu
if(theta[1]<0) stop("tauInv.GIG: only implemented for nu >= 0")
tau.max <- 1/(1+2*theta[1])
## make sure we are in the range of attainable Kendall's tau
if(tau>=tau.max) stop(paste("tauInv.GIG: the supremum of attainable taus is",
round(tau.max,8))) # this assumes tau to be falling in theta (numerical experiments show this behavior)
interval <- c(1e-100, upper[2])
}
replace(theta, i, uniroot(function(th) do.call(tau.GIG, c(list(replace(theta,
i, th),
upper=upper,
quiet=quiet),
iargs))-tau,
interval=interval, ...)$root)
}
## generate vectors of random variates from a GIG copula
rnacopula.GIG <- function(n, d, theta)
psi.GIG(matrix(rexp(n*d), nrow=n, ncol=d)/V0.GIG(n, theta), theta)