Nested Archimedean Lévy Copulas

Marius Hofert

2019-04-21

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

(Nested) Clayton Lévy copula

## 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
}

Plotting

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

##    user  system elapsed 
##   0.441   0.024   0.469

4d positive nested Clayton Lévy copula

##    user  system elapsed 
##  16.408   0.070  16.592