demo/tail_compatibility.R

In copula: Multivariate Dependence with Copulas

## Copyright (C) 2012 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan
##
## 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/>.


### This demo accompanies the paper Embrechts, Hofert, Wang ("Bernoulli and
### Tail-Dependence Compatibility", 2016)

library(copula)
doPDF <- FALSE


### 1 Plot Liebscher copula ####################################################


## The example we consider is
## \[
## 	\mathbf{U}=(\mathrm{max}\{U_{11}^2, U_{21}^2\},\ \mathrm{max}\{U_{12}^2, U_{22}^2\})\sim C\quad\text{for}\quad C(u_1,u_2) = C_1(\sqrt{u_1},\sqrt{u_2}) C_2(\sqrt{u_1},\sqrt{u_2}).
## \]

## General example setup
n <- 2000
set.seed(271)

## Define copula/df of U1
## lambda_l = 2^{-1/th}; th = 2^{-(1-tau)/(2*tau)}
family <- "Clayton"
th <- 4
copU1 <- onacopulaL(family, nacList=list(th, 1:2))
U1 <- rCopula(n, copula=copU1)
lambda1 <- copClayton@lambdaL(th)

## Define copula/df of U2
## lambda_l = 2*t_{n+1}(-sqrt(nu+1)*sqrt(1-th)/sqrt(1+th)); th = sin((pi/2) * tau)
family <- "t" # copula family
nu <- 3 # degrees of freedom
th <- 0.8 # parameter
copU2 <- ellipCopula(family, param=th, dim=2, df=nu) # define copula object
U2 <- rCopula(n, copula=copU2)
lambda2 <- 2*pt(-sqrt(nu+1)*sqrt(1-th)/sqrt(1+th), df=nu+1) # = 2t_4(-2/3)

## Sample a survival MO copula (U3)
alpha <- c(2^(-3/4), 0.8) # ~= 0.6, 0.8
U <- matrix(runif(3*n), ncol=3) # U'_1, U'_2, U'_12
U3 <- 1 - cbind(pmax(U[,1]^(1/(1-alpha[1])), U[,3]^(1/alpha[1])),
                pmax(U[,2]^(1/(1-alpha[2])), U[,3]^(1/alpha[2])))
lambda3 <- min(alpha) # lambda_l > 0, lambda_u = 0

## Define U
U <- cbind(pmax(U1[,1], U2[,1]), pmax(U1[,2], U2[,2]))^2 # Liebscher based on C, t3
U. <- cbind(pmax(U1[,1], U3[,1]), pmax(U1[,2], U3[,2]))^2 # Liebscher based on C, survival MO
## => off-diagonal entry of Lambda(_l) is:
(lambda <- lambda1*lambda2) # 2^(-1/4) * 2 * pt(-2/3, df=4) ~= 0.4553
(lambda. <- lambda1*lambda3) # 2^(-1/4) * 2^(-3/4) = 1/2

## Plots
if(doPDF) pdf(file=(file <- "U_Liebscher_n=2000_U1=C_th=4_U2=t3_th=0.5.pdf"), width=6, height=6)
par(pty="s")
plot(U, xlab=expression(italic(U[1])), ylab=expression(italic(U[2])))
if(doPDF) dev.off()
if(doPDF) pdf(file=(file <- "U_Liebscher_n=2000_U1=C_th=4_U2=sMO_a1=0.6_a2=0.8.pdf"), width=6, height=6)
par(pty="s")
plot(U., xlab=expression(italic(U[1])), ylab=expression(italic(U[2])))
if(doPDF) dev.off()


### 2 Sample $\mathbf{Y} = X\mathbf{U} + \mathbf{Z}\circ\mathbf{V}$ ############

## Note that the margins are standard uniform by construction.

##' @title Sampling Y = XU + Z circ V
##' @param n sample size
##' @param copU copula of U (m-dimensional) or an (n, m) matrix of samples
##' @param copV copula of V (d-dimensional) or an (n, d) matrix of samples
##'             lower tail-dependence matrix must be the identity!
##' @param X.method method for generating the sub-stochastic (d,m)-matrix X
##'        (no row has more than one 1)
##' @param ... additional arguments passed to rmultinom()
##' @return (n, d)-matrix containing samples from Y
##' @author Marius Hofert
##' @note does not allow for m=1 (due to copU being at least bivariate)
rY <- function(n, copU, copV, X.method=c("random", "multinomial"), ...)
{
    stopifnot(n >= 1)

    ## 1) sample/deal with U
    U <- if(length(dim(copU)) == 1) {
        rCopula(n, copula=copU)
    } else {
        stopifnot(is.matrix(copU), dim(copU) == c(n,m))
        copU
    }
    m <- ncol(U)
    stopifnot(m >= 2)

    ## 2) sample/deal with V
    if(length(dim(copV)) == 1) { # copV is a copula object
        d <- dim(copV)
        V <- rCopula(n, copula=copV)
    } else { # copV is a matrix of samples
        stopifnot(is.matrix(copV), dim(copV) == c(n, d))
        d <- ncol(copV)
        V <- copV # must have the identity matrix as lower tail-dependence matrix
    }
    stopifnot(d >= 2)

    ## 3) deal with argument 'X.method'
    X.method <- match.arg(X.method)
    if(X.method == "multinomial") {
        stopifnot(hasArg(prob))
        prob <- list(...)$prob
        stopifnot(length(prob) == m)
    }

    ## 4) fill Y
    Y <- matrix(, nrow=n, ncol=d)
    for(i in 1:n) {
        ## 4.1) sample the sub-stochastic (d,m)-matrix X (no row has more than one 1)
        switch(X.method,
               "random" = {
                   X <- matrix(0, nrow=d, ncol=m)
                   num.1 <- sample(0:d, size=1) # sample number of rows with one 1
                   if(num.1 > 0) {
                       ind.row <- sample(1:d, size=num.1) # indices of rows with one 1
                       ind.col <- sample(1:m, size=length(ind.row)) # randomly pick corresponding columns
                       for(k in seq_along(ind.row)) X[ind.row[k], ind.col[k]] <- 1 # set 1s in X
                   }
               },
               "multinomial" = {
                   X <- t(rmultinom(d, size=1, prob=prob)) # (d,m)-matrix; generates precisely one 1 in each row (no 0 matrix)
               },
               stop("wrong X.method"))
        ## 4.2) build Z
        Z <- 1-rowSums(X)
        stopifnot(Z >= 0, Z <= 1) # defensive programming
        ## 4.3) multiply together
        Y[i,] <- X %*% U[i,] + Z * V[i,] # = X %*% t(U[i,, drop=FALSE]) + Z * V[i,]
    }
    Y
}


### 2.1 $\mathbf{U}\sim t_3$, $\mathbf{V}\sim\Pi$ or $\mathbf{V}\sim G$ ########

d <- 2
set.seed(271)

## Setup for U
m <- 2
family.U <- "t" # copula family
nu <- 3 # degrees of freedom
tau.U <- 0.75 # Kendall's tau
th.U <- iTau(ellipCopula(family.U, df=nu), tau.U) # corresponding parameter
copU <- ellipCopula(family.U, param=th.U, dim=m, df=nu) # define copula object

## Setup for V
V <- matrix(runif(n*d), ncol=2) # V ~ Pi
family.V2 <- "Gumbel"
tau.V2 <- 0.75
th.V2 <- iTau(archmCopula(family.V2), tau.V2)
copV2 <- archmCopula(family.V2, param=th.V2, dim=d) # V ~ G

## Sample Y
Y1 <- rY(n, copU=copU, copV=V, X.method="random")
Y2 <- rY(n, copU=copU, copV=V, X.method="multinomial", prob=c(0.5, 0.5))
Y3 <- rY(n, copU=copU, copV=copV2, X.method="random")
Y4 <- rY(n, copU=copU, copV=copV2, X.method="multinomial", prob=c(0.5, 0.5))

## Plot (left: X ~ random; right: X ~ multinomial; top: V ~ Pi; bottom: V ~ G;
##       this plot is not contained in the paper)
opar <- par(no.readonly=TRUE)
par(pty="s", mar=c(2.6, 2, 1, 1) + 0.1)
lay <- matrix(1:4, ncol=2, byrow=TRUE) # layout matrix
layout(lay, widths=c(1, 1), heights=c(1, 1)) # layout
plot(Y1, xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
plot(Y2, xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
plot(Y3, xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
plot(Y4, xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
par(opar)
## Note:
## - the shape of Y is mainly determined by U, X
## - for X ~ multinomial, the choice of copula for V is barely visible


### 2.2 $\mathbf{U}\sim$ survival MO, $\mathbf{V}\sim\Pi$ or $\mathbf{V}\sim G$

## Sample U
alpha <- c(0.25, 0.75) # parameters
U. <- matrix(runif(3*n), ncol=3) # U'_1, U'_2, U'_12
U <- 1 - cbind(pmax(U.[,1]^(1/(1-alpha[1])), U.[,3]^(1/alpha[1])),
               pmax(U.[,2]^(1/(1-alpha[2])), U.[,3]^(1/alpha[2])))

## Setup for V (V ~ Pi as before)
family.V3 <- "normal"
tau.V3 <- 0.9
th.V3 <- iTau(ellipCopula(family.V3), tau.V3)
copV3 <- ellipCopula(family.V3, param=th.V3, dim=d) # V ~ Ga

## Sample Y
Y1. <- rY(n, copU=U, copV=V, X.method="random")
Y2. <- rY(n, copU=U, copV=V, X.method="multinomial", prob=c(0.5, 0.5))
Y3. <- rY(n, copU=U, copV=copV3, X.method="random")
Y4. <- rY(n, copU=U, copV=copV3, X.method="multinomial", prob=c(0.5, 0.5))

## Plot (left: X ~ random; right: X ~ multinomial; top: V ~ Pi; bottom: V ~ G;
##       this plot is not contained in the paper)
opar <- par(no.readonly=TRUE)
par(pty="s", mar=c(2.6, 2, 1, 1) + 0.1)
lay <- matrix(1:4, ncol=2, byrow=TRUE) # layout matrix
layout(lay, widths=c(1, 1), heights=c(1, 1)) # layout
plot(Y1., xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
plot(Y2., xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
plot(Y3., xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
plot(Y4., xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)
par(opar)
## Note:
## - the shape of Y is mainly determined by U, X
## - for X ~ multinomial, the choice of copula for V is barely visible

## Plot for paper (V ~ Pi; left: X ~ random; right: X ~ multinomial; top: U ~ t3; bottom: U ~ sMO)
file <- "Y_n=2000_U=t3_tau=0.75_or_sMO_a1=0.25_a2=0.75_X=random_or_multinomial_V=Pi.pdf"
if(doPDF) pdf(file=file, width=8, height=8)
opar <- par(no.readonly=TRUE)
par(pty="s", mar=c(2.6, 2, 1, 1) + 0.1)
lay <- matrix(1:4, ncol=2, byrow=TRUE) # layout matrix
layout(lay, widths=c(1, 1), heights=c(1, 1)) # layout
plot(Y1, xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)  # U=t3,  X=random, V=Pi
plot(Y2, xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5)  # U=t3,  X=multinomial, V=Pi
plot(Y1., xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5) # U=sMO, X=random, V=Pi
plot(Y2., xlab=expression(italic(Y[1])), ylab=expression(italic(Y[2])), cex=0.5) # U=sMO, X=multinomial, V=Pi
par(opar)
if(doPDF) dev.off()


### 3 Example of Federico Degen ################################################

## Y for the example of Federico Degen (note: m=1)
rY.FD <- function(n, copV, alpha)
{
    stopifnot(n >= 1)

    ## 1) sample U
    U <- runif(n)

    ## 2) sample/deal with V
    if(length(dim(copV)) == 1) { # copV is a copula object
        d <- dim(copV)
        V <- rCopula(n, copula=copV)
    } else { # copV is a matrix of samples
        stopifnot(is.matrix(copV), dim(copV) == c(n, d))
        d <- ncol(copV)
        V <- copV # must have the identity matrix as lower tail-dependence matrix
    }
    stopifnot(d >= 2, 0 <= alpha, alpha <= 1/(d-1))

    ## 3) sample X
    ##    interpret each *row* of X as 1 *column* vector of the
    ##    (mathematical) X (in {0,1}^{d,m=1})
    ##    => sanity check: the row sum has to be in {0,1,2}
    X <- matrix(0, nrow=n, ncol=d)
    X[,d] <- 1 # last element is 1; the d-1 before are 1 each with probability alpha (but at most one can be 1)
    W <- sample(1:d, size=n, replace=TRUE, prob=c(rep(alpha, d-1), 1-(d-1)*alpha))
    ## note: the following is *not* correct as sample() scales prob
    ## W <- if(alpha>0) sample(1:(d-1), size=n, replace=TRUE, prob=rep(alpha, d-1)) else rep(0, n)
    X[cbind(1:n, W)] <- 1 # if W[i] == d, then we don't change X[,d] => fine
    stopifnot(rowSums(X) <= 2) # sanity check

    ## 4) build Z
    Z <- 1-X # true for m=1
    stopifnot(Z >= 0, Z <= 1) # defensive programming

    ## 5) build Y
    X*U + Z*V
}

## General example setup
d <- 4
alpha <- 1/(d-1) # maximal allowed alpha
set.seed(271)

## Setup for V
V <- matrix(runif(n*d), ncol=d) # V ~ Pi
family <- "normal"
tau <- 0.8
th <- iTau(ellipCopula(family), tau)
copV <- ellipCopula(family, param=th, dim=d) # V ~ Ga

## Sample Y
Y  <- rY.FD(n, copV=V, alpha=alpha)
Y. <- rY.FD(n, copV=copV, alpha=alpha)

## Plots
if(doPDF) pdf(file=(file <- "Y_n=2000_FD_max_alpha_V=Pi.pdf"), width=6, height=6)
par(pty="s")
pairs(Y, labels=as.expression( sapply(1:d, function(j) bquote(italic(Y[.(j)]))) ), gap=0, cex=0.25)
if(doPDF) dev.off()

if(doPDF) pdf(file=(file <- "Y_n=2000_FD_max_alpha_V=Ga_tau=0.8.pdf"), width=6, height=6)
par(pty="s")
pairs(Y., labels=as.expression( sapply(1:d, function(j) bquote(italic(Y[.(j)]))) ), gap=0, cex=0.25)
if(doPDF) dev.off()