The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
require(gsl) # for exponential integral
require(copula)
<- FALSE doPDF
Our goal is to simulate dependent multivariate Lévy processes based on positive (nested) Archimedean Lévy copulas (here: Clayton). As usual, we have to truncate small jumps. We do so by truncating large Gammas (by setting them to \(\infty\) in order for the jump heights to be \(\bar{\nu}^{-1}(\infty) = 0\)). In this sense, we only simulate the largest jumps; see below for more details.
We start by defining some auxiliary functions used later on.
## Tail integral of a variance gamma Levy process
## \bar{\nu}(x) = \int_x^\infty f(z) dz for the Levy density
## f(z) = (c/x)*exp(-lambda*x) for x>0, c=1/kappa and
## lambda=(sqrt(theta^2+2*sigma^2/kappa)-theta)/sigma^2
<- function(x, th, kap, sig) {
nu_bar_vargamma <- (sqrt(th^2+2*sig^2/kap)-th)/sig^2
lambda -expint_Ei(-lambda*x, give=FALSE)/kap
}
## Inverse of the tail integral of a variance gamma Levy process
## \bar{\nu}(x) = \int_x^\infty f(z) dz for the Levy density
## f(z) = (c/x)*exp(-lambda*x) for x>0, c=1/kappa and
## lambda=(sqrt(theta^2+2*sigma^2/kappa)-theta)/sigma^2
<- function(Gamma, th, kap, sig, ...)
nu_bar_inv_vargamma
{<- nu_bar_vargamma(.Machine$double.xmin, th=th, kap=kap, sig=sig)
max.val <- numeric(length(Gamma))
res <- Gamma >= max.val
large <- 0 # de facto indistinguishable from 0 anyways
res[large] if(any(!large)) {
<- (sqrt(th^2+2*sig^2/kap)-th)/sig^2
lambda <- function(x, z)
nu_bar_vargamma_minus -expint_Ei(-lambda*x, give=FALSE)/kap - z
!large] <- vapply(Gamma[!large], function(Gamma.)
res[uniroot(nu_bar_vargamma_minus, z=Gamma.,
interval=c(.Machine$double.xmin, 29), ...)$root, NA_real_)
}
res }
## Transforming Gamma with variance-gamma Levy margins
<- function(Gamma, th, kap, sig)
hom_vargamma_Levy
{<- runif(nrow(Gamma)) # jump times
U <- order(U) # determine the order of the U's
ord <- U[ord] # (sorted) jump times
jump_time <- apply(Gamma, 2, function(y)
jump_size nu_bar_inv_vargamma(y, th=th, kap=kap, sig=sig)) # (unsorted) jump sizes (apply inverses of marginal tail integrals)
<- apply(jump_size, 2, function(x) cumsum(x[ord])) # sort jump sizes according to U's and add them up => (L_t) at jump times
value list(jump_time=jump_time, value=value)
}
## \bar{\psi} for Clayton Levy copulas
<- function(t, theta) t^(-1/theta) psi_bar_Clayton
## V_{01} for nested Clayton Levy copulas
## Note: V_{01,k} | V_{0,k} ~ LS^{-1}[\bar{\psi}_{01}(.; V_{0,k})] with
## \bar{\psi}_{01}(t; V_{0,k}) = \exp(-V_{0,k} t^{\theta_0/\theta_1})
## = copGumbel@V01() (not copClayton@V01()!)
<- function(V0, theta0, theta1)
V01_nested_Clayton_Levy @V01(V0, theta0=theta0, theta1=theta1) copGumbel
## Generate Gamma for a d-dimensional Clayton Levy copula with parameter theta
## Note: - Don't confuse the Clayton parameter theta with the parameter th
## for the marginal tail integral (variance gamma)
## - The advantage of a fixed truncation point Gamma^* is that one can
## correct the bias introduced when cutting off small jumps by adding
## a drift; see Asmussen, Rosinski (2001) for more details
## - The best stopping criterion would be if we are sure that in each
## dimension all generated Gammas which are <= Gamma^* (= Gamma.star)
## form a sample of jump times of a homogeneous Poi(1) process on
## [0, Gamma^*]; this could be tested.
## - We go with a simpler stopping criterion here: Given a burn.in value
## (integer), we stop (only) if in the last burn.in-many generated
## Gammas each had at least one component larger than Gamma^*. So it's
## unlikely that we still get such (uniformly) small Gammas
## (<= Gamma^* in each component); large Gamma => small jump
## => we correctly only truncate (small) jumps.
<- function(d, theta, Gamma.star, burn.in)
Gamma_Clayton_Levy
{stopifnot(d >= 1, length(theta) == 1, theta > 0 , Gamma.star > 0,
>= 1)
burn.in <- matrix(, nrow=0, ncol=d)
Gamma <- 0
count <- 0
W repeat {
<- rexp(d+1)
E <- W + E[d+1]
W <- (W/theta * gamma(1/theta))^theta # generate V = F^{-1}(W)
V <- psi_bar_Clayton(E[1:d]/V, theta=theta) # Gamma
G <- rbind(Gamma, G) # update Gamma
Gamma if(count >= burn.in) break # stopping criterion
<- if(any(G <= Gamma.star)) 1 else count + 1 # if there are still Gammas <= Gamma^*, keep generating Gammas
count
}> Gamma.star] <- Inf # => produce \bar{\mu}(.) = 0 (0-height jumps)
Gamma[Gamma
Gamma }
## Generate Gamma for a 4-dimensional nested Clayton Levy copula
<- function(theta, Gamma.star, burn.in)
Gamma_nested_Clayton_Levy
{stopifnot(d >= 1, length(theta) == 3, theta > 0, min(theta[2:3]) >= theta[1],
> 0, burn.in >= 1)
Gamma.star <- 4 # d must be 4 here; obviously, this could be generalized
d <- matrix(, nrow=0, ncol=d)
Gamma <- 0
count<- 0
W repeat {
<- rexp(d+1)
E <- W + E[d+1]
W <- (W/theta[1] * gamma(1/theta[1]))^theta[1] # generate V_0 = F_0^{-1}(W)
V0 <- V01_nested_Clayton_Levy(V0, theta0=theta[1], theta1=theta[2]) # generate V_{01}
V01 <- V01_nested_Clayton_Levy(V0, theta0=theta[1], theta1=theta[3]) # generate V_{02}
V02 <- c(psi_bar_Clayton(E[1:2]/V01, theta=theta[2]),
G psi_bar_Clayton(E[3:4]/V02, theta=theta[3])) # Gamma
<- rbind(Gamma, G) # update Gamma
Gamma if(count >= burn.in) break # stopping criterion
<- if(any(G <= Gamma.star)) 1 else count + 1 # if there are still Gammas <= Gamma^*, keep generating Gammas
count
}> Gamma.star] <- Inf # => produce \bar{\mu}(.) = 0 (0-height jumps)
Gamma[Gamma
Gamma }
## Plot Gammas
<- function(Gamma, Gamma.star, file=NULL, ...)
plot_Gamma
{stopifnot(is.matrix(Gamma), (d <- ncol(Gamma)) >= 2,
is.null(file) || is.character(file))
<- colorRampPalette(c("black", "royalblue3", "darkorange2",
palette "maroon3"), space="Lab")
<- palette(d)
cols ## cols <- adjustcolor(cols, alpha.f=0.1) # no improvement here
<- !is.null(file)
doPDF if(doPDF) pdf(file=file, width=7, height=7)
plot(Gamma[,1], type="l", ylim=range(Gamma, finite=TRUE), # omit Inf
log="y", xlab="k", ylab="", col=cols[1], ...)
for(j in 2:d) lines(Gamma[,j], col=cols[j])
abline(h=Gamma.star, lty=2, lwd=1.6)
legend("bottomright", bty="n", lty=c(rep(1, d), 2), lwd=c(rep(1,d), 1.6),
col=c(cols, "black"), as.expression( c(lapply(1:d, function(j)
bquote(Gamma[list(k,.(j))])), list(bquote(Gamma*"*")))))
if(doPDF) dev.off()
}
## Plot a multivariate Levy process
<- function(L, file=NULL, ...)
plot_Levy
{stopifnot(is.matrix(L$value), (d <- ncol(L$value)) >= 2,
length(L$jump_time)==nrow(L$value),
is.null(file) || is.character(file))
<- colorRampPalette(c("black", "royalblue3", "darkorange2",
palette "maroon3"), space="Lab")
<- palette(d) # d colors
cols <- c(0, L$jump_time, 1) # extended jump times (for nicer plotting)
x_jump_time <- rbind(rep(0, d), L$value, L$value[nrow(L$value),]) # extended Levy process (for nicer plotting)
x_L <- !is.null(file)
doPDF if(doPDF) pdf(file=file, width=7, height=7)
plot(x_jump_time, x_L[,1], type="s", ylim=range(L),
xlab="t", ylab=expression(bold(L)[t]), col=cols[1], ...)
for(j in 2:d)
lines(x_jump_time, x_L[,j], type="s", col=cols[j])
legend("bottomright", bty="n", lty=rep(1, d), col=cols,
legend=as.expression( lapply(1:d, function(j) bquote(L[list(t,.(j))]))))
if(doPDF) dev.off()
}
Now let’s sample some paths.
## Marginal (variance gamma) parameters
<- -0.2
th <- 0.05
kap <- 0.3
sig
## Truncation specification
<- 2000
Gamma.star <- 500 burn.in
## Gamma
<- 4 # theta
theta <- 4 # dimension
d set.seed(271)
system.time(Gamma <- Gamma_Clayton_Levy(d, theta=theta, Gamma.star=Gamma.star,
burn.in=burn.in))
## User System verstrichen
## 0.131 0.000 0.132
plot_Gamma(Gamma, Gamma.star=Gamma.star,
file=if(doPDF) "fig_Gamma_positive_Clayton_Levy_copula.pdf" else NULL,
main=expression(bold(group("(",Gamma[k],")")~
"for a positive Clayton Levy copula")))
## (L_t)
<- hom_vargamma_Levy(Gamma, th=th, kap=kap, sig=sig)
L plot_Levy(L, file=if(doPDF) "fig_L_with_positive_Clayton_Levy_copula.pdf" else NULL,
main="Levy process with positive Clayton Levy copula")
## Gamma
<- c(0.7, 4, 2) # theta_0, theta_1, theta_2
theta set.seed(271)
system.time(Gamma <- Gamma_nested_Clayton_Levy(theta, Gamma.star=Gamma.star,
burn.in=burn.in))
## User System verstrichen
## 5.293 0.036 5.349
## 15 seconds on 2015-fast platform
plot_Gamma(Gamma, Gamma.star=Gamma.star,
file=if(doPDF) "fig_Gamma_positive_nested_Clayton_Levy_copula.pdf" else NULL,
main=expression(bold(group("(",Gamma[k],")")~
"for a positive nested Clayton Levy copula")))
## (L_t)
<- hom_vargamma_Levy(Gamma, th=th, kap=kap, sig=sig)
L plot_Levy(L, file=if(doPDF) "fig_L_with_positive_nested_Clayton_Levy_copula.pdf" else NULL,
main="Levy process with positive nested Clayton Levy copula")
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.