Hierarchical Archimax Copulas

Marius Hofert and Avinash Prasad

2023-12-07

library(copula) # for hierarchical frailties, generators etc.
library(mev) # for rmev()

Setup and auxiliary functions

The (frozen) parameters considered are as follows.

## Parameters
d. <- c(2, 3) # copula sector dimensions
d <- sum(d.) # copula dimension
stopifnot(d >= 3) # for 3d plots
n <- 777 # sample size  (was 1000; save space for *.html)

We also implement an auxiliary function to generate the hierarchical frailties involved.

##' @title Generated two hierarchical frailties (V0, V01, V02)
##' @param n sample size
##' @param family copula family
##' @param tau Kendall's tau
##' @return 3-column matrix containing the hierarchical frailties
##' @author Marius Hofert and Avinash Prasad
rV012 <- function(n, family, tau)
{
    stopifnot(n >= 1, is.character(family), length(tau) == 3, tau > 0)
    cop <- getAcop(family) # corresponding AC
    th <- iTau(cop, tau = tau) # copula parameters
    V0 <- cop@V0(n, theta = th[1]) # sample frailties V_0
    V01 <- cop@V01(V0, theta0 = th[1], theta1 = th[2]) # sample frailties V_01
    V02 <- cop@V01(V0, theta0 = th[1], theta1 = th[3]) # sample frailties V_02
    cbind(V0  = V0, V01 = V01, V02 = V02)
}

And a short helper function for plotting (possibly to PDF).

## Plot with possible export to PDF
mypairs <- function(x, pch = ".", file = NULL, width = 6, height = 6, ...)
{
    opar <- par(pty = "s")
    on.exit(par(opar))
    doPDF <- !is.null(file)
    if(doPDF) {
        pdf(file = file, width = width, height = height)
        stopifnot(require(crop))
    }
    pairs2(x, pch = pch, ...)
    if(doPDF) dev.off.crop(file)
}

1 ACs vs AXCs vs NACs vs (different) HAXCs

In this example, we compare Archimedean copulas (ACs), Archimax copulas (AXCs), nested Archimedean copulas (NACs) and hierarchical Archimax copulas (HAXCs) with hierarchical frailties only, with hierarchical frailties and hierarchical extreme-value copula (HEVC) but both of the same hierarchical structure, and with hierarchical frailties and HEVC but of different hierarchical structure.

1.1 AC (Clayton copula)

To begin with, we sample from a standard Clayton copula.

## Sample frailties (to recycle the frailties, we apply the Marshall--Olkin (MO)
## algorithm manually here)
tau <- 0.4 # Kendall's tau
family <- "Clayton" # frailty family
cop <- getAcop(family) # corresponding AC
th <- iTau(cop, tau = tau) # copula parameter
set.seed(271) # for reproducibility
V <- cop@V0(n, theta = th) # sample frailty
E <- matrix(rexp(n * d), ncol = d) # sample Exp(1)
U.AC <- cop@psi(E/V, theta = th) # MO

## Plot
mypairs(U.AC)

1.2 AXC (Clayton frailties and Gumbel EVC)

Now we sample from an Archimax copula with gamma frailties (as underlying the Clayton family; we recycle the frailties from 1.1) and Gumbel EVC (parameters chosen such that Kendall’s tau equals 0.5).

## Generate Gumbel EVC sample with Exp(1) margins
tau.EVC <- 0.5 # Kendall's tau
family.EVC <- "Gumbel" # EVC family
th.EVC <- iTau(getAcop(family.EVC), tau = tau.EVC) # EVC parameter
cop.EVC <- onacopulaL(family.EVC, list(th.EVC, 1:d)) # EVC
set.seed(271) # for reproducibility
U.EVC <- rCopula(n, copula = cop.EVC) # sample EVC
E.EVC <- -log(U.EVC) # map the EVC to Exp(1) margins

## Combine with the (same) Clayton frailties (as before)
U.AXC <- cop@psi(E.EVC/V, theta = th)

## Plot
mypairs(U.AXC)

1.3 NAC (nested Clayton)

Sampling from a NAC based on Clayton’s family with parameters such that Kendall’s tau equals 0.2 between the two sectors, 0.4 within the first sector and 0.6 within the second sector, can be done as follows.

## Generate hierarchical frailties
tau.N <- c(0.2, 0.4, 0.6) # Kendall's taus
family.N <- "Clayton" # frailty family
cop.N <- getAcop(family.N) # corresponding AC
th.N <- iTau(cop.N, tau = tau.N) # copula parameters
set.seed(271) # for reproducibility
V.N <- rV012(n, family = family.N, tau = tau.N) # generate corresponding frailties
V0  <- V.N[,"V0"]
V01 <- V.N[,"V01"]
V02 <- V.N[,"V02"]

## Combine with independent Exp(1)
U.NAC <- cbind(cop.N@psi(E[,1:d.[1]]/V01,     theta = th.N[2]),
               cop.N@psi(E[,(d.[1]+1):d]/V02, theta = th.N[3]))

## Plot
mypairs(U.NAC)

1.4 HAXC (hierarchical Clayton frailties and Gumbel EVC)

Sampling from a NAXC based on hierarchical Clayton frailties (recycled from 1.3) and Gumbel EVC (realizations recycled from 1.2). Note that hierarchies only take place at the levels of the frailties in this case.

## Generate samples
U.HAXC <- cbind(cop.N@psi(E.EVC[,1:d.[1]]/V01,     theta = th.N[2]),
                cop.N@psi(E.EVC[,(d.[1]+1):d]/V02, theta = th.N[3]))

## Plot
mypairs(U.HAXC)

1.5 HAXC (hierarchical Clayton frailties and hierarchical Gumbel EVC, same hierarchical structure)

After this warm-up, we can consider sampling from a HAXC based on hierarchical Clayton frailties (recycled from 1.3) and a hierarchical Gumbel EVC (with sector sizes 2 and 3 and parameters such that Kendall’s tau equals 0.2 between the two sectors, 0.5 within the first sector and 0.7 within the second sector). Note that there are two types of hierarchies involved, at the level of the (hierarchical) frailties and at the level of the (hierarchical) EVC. Furthermore, their hierarchical structure is the same.

## Sampling from a hierarchical Gumbel copula
tau.HEVC <- c(0.2, 0.5, 0.7) # Kendall's taus
family.HEVC <- "Gumbel" # EVC family
cop.HEVC <- getAcop(family.HEVC) # corresponding base EVC
th.HEVC <- iTau(cop.HEVC, tau = tau.HEVC) # copula parameters
cop.HEVC <- onacopulaL(family.HEVC, list(th.HEVC[1], NULL,
                                     list(list(th.HEVC[2], 1:d.[1]),
                                          list(th.HEVC[3], (d.[1]+1):d)))) # hierarchical EVC
set.seed(271) # for reproducibility
U.HEVC <- rCopula(n, copula = cop.HEVC) # sample the HEVC
E.HEVC <- -log(U.HEVC) # map the HEVC samples to Exp(1) margins

## Combine the hierarchical Gumbel EVC with Exp(1) margins with the hierarchical frailties
U.HAXC.HEVC.same <- cbind(cop.N@psi(E.HEVC[,1:d.[1]]/V01,     theta = th.N[2]),
                          cop.N@psi(E.HEVC[,(d.[1]+1):d]/V02, theta = th.N[3]))

## Plot
mypairs(U.HAXC.HEVC.same)