# Nested Archimedean Lévy Copulas

#### 2022-11-17

require(gsl) # for exponential integral
require(copula)
doPDF <- FALSE

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.

## 1 Auxiliary functions

We start by defining some auxiliary functions used later on.

### Margins

## 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
nu_bar_vargamma <- function(x, th, kap, sig) {
lambda <- (sqrt(th^2+2*sig^2/kap)-th)/sig^2
-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
nu_bar_inv_vargamma <- function(Gamma, th, kap, sig, ...)
{
max.val <- nu_bar_vargamma(.Machine$double.xmin, th=th, kap=kap, sig=sig) res <- numeric(length(Gamma)) large <- Gamma >= max.val res[large] <- 0 # de facto indistinguishable from 0 anyways if(any(!large)) { lambda <- (sqrt(th^2+2*sig^2/kap)-th)/sig^2 nu_bar_vargamma_minus <- function(x, z) -expint_Ei(-lambda*x, give=FALSE)/kap - z res[!large] <- vapply(Gamma[!large], function(Gamma.) uniroot(nu_bar_vargamma_minus, z=Gamma., interval=c(.Machine$double.xmin, 29), ...)$root, NA_real_) } res } ## Transforming Gamma with variance-gamma Levy margins hom_vargamma_Levy <- function(Gamma, th, kap, sig) { U <- runif(nrow(Gamma)) # jump times ord <- order(U) # determine the order of the U's jump_time <- U[ord] # (sorted) jump times jump_size <- apply(Gamma, 2, function(y) nu_bar_inv_vargamma(y, th=th, kap=kap, sig=sig)) # (unsorted) jump sizes (apply inverses of marginal tail integrals) value <- 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 list(jump_time=jump_time, value=value) } ### (Nested) Clayton Lévy copula ## \bar{\psi} for Clayton Levy copulas psi_bar_Clayton <- function(t, theta) t^(-1/theta) ## 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()!) V01_nested_Clayton_Levy <- function(V0, theta0, theta1) copGumbel@V01(V0, theta0=theta0, theta1=theta1) ## 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. Gamma_Clayton_Levy <- function(d, theta, Gamma.star, burn.in) { stopifnot(d >= 1, length(theta) == 1, theta > 0 , Gamma.star > 0, burn.in >= 1) Gamma <- matrix(, nrow=0, ncol=d) count <- 0 W <- 0 repeat { E <- rexp(d+1) W <- W + E[d+1] V <- (W/theta * gamma(1/theta))^theta # generate V = F^{-1}(W) G <- psi_bar_Clayton(E[1:d]/V, theta=theta) # Gamma Gamma <- rbind(Gamma, G) # update Gamma if(count >= burn.in) break # stopping criterion count <- if(any(G <= Gamma.star)) 1 else count + 1 # if there are still Gammas <= Gamma^*, keep generating Gammas } Gamma[Gamma > Gamma.star] <- Inf # => produce \bar{\mu}(.) = 0 (0-height jumps) Gamma } ## Generate Gamma for a 4-dimensional nested Clayton Levy copula Gamma_nested_Clayton_Levy <- function(theta, Gamma.star, burn.in) { stopifnot(d >= 1, length(theta) == 3, theta > 0, min(theta[2:3]) >= theta[1], Gamma.star > 0, burn.in >= 1) d <- 4 # d must be 4 here; obviously, this could be generalized Gamma <- matrix(, nrow=0, ncol=d) count<- 0 W <- 0 repeat { E <- rexp(d+1) W <- W + E[d+1] V0 <- (W/theta[1] * gamma(1/theta[1]))^theta[1] # generate V_0 = F_0^{-1}(W) V01 <- V01_nested_Clayton_Levy(V0, theta0=theta[1], theta1=theta[2]) # generate V_{01} V02 <- V01_nested_Clayton_Levy(V0, theta0=theta[1], theta1=theta[3]) # generate V_{02} G <- c(psi_bar_Clayton(E[1:2]/V01, theta=theta[2]), psi_bar_Clayton(E[3:4]/V02, theta=theta[3])) # Gamma Gamma <- rbind(Gamma, G) # update Gamma if(count >= burn.in) break # stopping criterion count <- if(any(G <= Gamma.star)) 1 else count + 1 # if there are still Gammas <= Gamma^*, keep generating Gammas } Gamma[Gamma > Gamma.star] <- Inf # => produce \bar{\mu}(.) = 0 (0-height jumps) Gamma } ### Plotting ## Plot Gammas plot_Gamma <- function(Gamma, Gamma.star, file=NULL, ...) { stopifnot(is.matrix(Gamma), (d <- ncol(Gamma)) >= 2, is.null(file) || is.character(file)) palette <- colorRampPalette(c("black", "royalblue3", "darkorange2", "maroon3"), space="Lab") cols <- palette(d) ## cols <- adjustcolor(cols, alpha.f=0.1) # no improvement here doPDF <- !is.null(file) 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 plot_Levy <- function(L, file=NULL, ...) { stopifnot(is.matrix(L$value), (d <- ncol(L$value)) >= 2, length(L$jump_time)==nrow(L$value), is.null(file) || is.character(file)) palette <- colorRampPalette(c("black", "royalblue3", "darkorange2", "maroon3"), space="Lab") cols <- palette(d) # d colors x_jump_time <- c(0, L$jump_time, 1) # extended jump times (for nicer plotting)
x_L <- rbind(rep(0, d), L$value, L$value[nrow(L\$value),]) # extended Levy process (for nicer plotting)
doPDF <- !is.null(file)
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()
}

## 2 Sampling paths

Now let’s sample some paths.

## Marginal (variance gamma) parameters
th <- -0.2
kap <- 0.05
sig <- 0.3

## Truncation specification
Gamma.star <- 2000
burn.in <- 500

### 4d positive Clayton Lévy copula

## Gamma
theta <- 4 # theta
d <- 4 # dimension
set.seed(271)
system.time(Gamma <- Gamma_Clayton_Levy(d, theta=theta, Gamma.star=Gamma.star,
burn.in=burn.in))
##        User      System verstrichen
##       0.134       0.000       0.135
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)
L <- hom_vargamma_Levy(Gamma, th=th, kap=kap, sig=sig)
plot_Levy(L, file=if(doPDF) "fig_L_with_positive_Clayton_Levy_copula.pdf" else NULL,
main="Levy process with positive Clayton Levy copula")

### 4d positive nested Clayton Lévy copula

## Gamma
theta <- c(0.7, 4, 2) # theta_0, theta_1, theta_2
set.seed(271)
system.time(Gamma <- Gamma_nested_Clayton_Levy(theta, Gamma.star=Gamma.star,
burn.in=burn.in))
##        User      System verstrichen
##       6.666       0.012       6.749
## 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)
L <- hom_vargamma_Levy(Gamma, th=th, kap=kap, sig=sig)
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")