# Copula Constructions for Tail-Dependence Matrices

#### 2017-08-31

require(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
## [1] 0.45532
(lambda. <- lambda1*lambda3) # 2^(-1/4) * 2^(-3/4) = 1/2
## [1] 0.5
## 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)
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)
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)
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()