scripts/Q7_Copulas_and_Dependence.R

In 3schwartz/SpecialeScrAndFun: Script and functions used for master thesis.

#### A_motivating_example ####
library(copula)

## Generate some data (just skip this step, let's pretend we don't see it)
n <- 1000 # sample size
set.seed(271) # for reproducibility
U <- rCopula(n, copula = gumbelCopula(iTau(gumbelCopula(), tau = 0.5))) # geek zone
X <- qnorm(U) # quantile transformation
X. <- qexp(U) # quantile transformation

## Plots
layout(rbind(1:2))
plot(X,  xlab = expression(X[1]),     ylab = expression(X[2]))
plot(X., xlab = expression(X[1]*"'"), ylab = expression(X[2]*"'"))

## Q: For which data is the dependence between the two variables "larger"?

## This is difficult to answer, as both axis show different scales/ranges etc.
## We can transform the data to have the same marginals by using the probability
## transformation. Instead of ("cheating" by) using the above chosen marginal
## distribution functions (dfs), we make the realistic assumption that we don't
## know them, so we would estimate them by their empirical dfs (edfs) and then
## apply these marginal empirical dfs to the respective data columns. The function
## pobs() of the R package indirectly does that for you; see (*) for the 'why'.
## By the probability transformation, we then get the the marginal dfs of both
## data sets X and X. are (approximately) U(0,1), so allow for a comparison.
U  <- pobs(X)
U. <- pobs(X.)
## (*) Applying the edf \hat{F}_{n,j} of the jth margin to the data X_{ij},
## i in {1,..,n}, we obtain U_{ij} = \hat{F}_{n,j}(X_ij)
## = (1/n)\sum_{k=1}^n \I_{\{X_{kj} \le X_{ij}\}} = R_{ij}/n, where R_{ij}
## denotes the rank of X_{ij} among X_{1j},..,X_{nj}. One then typically uses
## the slightly scaled version R_{ij}/(n+1) to avoid that that U_{ij} can be 1
## (so on the boundary of [0,1]^d). The function pobs() computes exactly that.

## Now plot the pseudo-observations
plot(U,  xlab = expression(U[1]),     ylab = expression(U[2]))
plot(U., xlab = expression(U[1]*"'"), ylab = expression(U[2]*"'"))

## A: The dependence between the components (= columns) of X and X. is the same!
##    Therefore, in order to study dependence (independently of the margins),
##    we need to study the distributions of random vectors with U(0,1) marginals,
##    hence copulas.

## Visually check whether the margins are indeed (approximately) U(0,1)
layout(matrix(1:4, ncol = 2))
plot(U[,1],  ylab = expression(U[1]))
plot(U[,2],  ylab = expression(U[2]))
plot(U.[,1], ylab = expression(U[1]*"'"))
plot(U.[,2], ylab = expression(U[2]*"'"))
layout(1)
## Hence U and U. above are "pseudo"-observations (hence the name "pobs") of the
## underlying copulas. They give us insight in the actual dependence structure
## (copula) underlying our data sets X and X.

#### First_principles ####
## Generate and transform multivariate data (with the probability and quantile
## transformations)
library(mvtnorm) # for sampling from N(0, P)
library(copula) # for pairs2()
library(qrmtools) # for qPar()

set.seed(271) # for reproducibility

### 1 Normal copula sample ##

## Sample from a bivariate normal distribution
P <- matrix(0.7, nrow = 2, ncol = 2) # correlation matrix; play with the entry!
diag(P) <- 1 # set diagonals to 1
n <- 1000 # sample size
set.seed(271) # for reproducibility
X <- rmvnorm(n, sigma = P) # sample from N(0, P)
plot(X, xlab = expression(X[1]), ylab = expression(X[2])) # scatter plot

## Marginally apply probability transforms (with the corresponding dfs)
U <- pnorm(X) # probability transformation
plot(U, xlab = expression(U[1]), ylab = expression(U[2])) # sample from the Gauss copula

## (Visually) check that the margins of our generated data are indeed U(0,1)
plot(U[,1], ylab = expression(U[1])) # scatter plot
hist(U[,1], xlab = expression(U[1]), probability = TRUE, main = "") # histogram
plot(U[,2], ylab = expression(U[2])) # scatter plot
hist(U[,2], xlab = expression(U[2]), probability = TRUE, main = "") # histogram

## Map the Gauss copula sample to one with t_3.5 margins
nu <- 3.5 # degrees of freedom
Y <- qt(U, df = nu) # quantile transformation
plot(Y, xlab = expression(Y[1]), ylab = expression(Y[2]), col = "royalblue3") # scatter plot
points(X)
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
       legend = c(expression(N(0,P)),
                  substitute("With"~italic(t)[nu.]~"margins", list(nu. = nu))))
## => Sample is stretched out by heavier tailed marginal distributions


### 2 t copula sample ##

## Build a block correlation matrix
d <- 5
rho <- c(0.3, 0.6, 0.9)
P <- matrix(rho[1], nrow = d, ncol = d)
P[1:2, 1:2] <- rho[2]
P[3:5, 3:5] <- rho[3]
diag(P) <- 1
P

## Sample from a multivariate t distribution
X <- rmvt(n, sigma = P, df = nu) # sample from t_{P, nu}
pairs2(X, labels.null.lab = "X", cex = 0.4, col = adjustcolor("black", alpha.f = 0.5))

## Marginally apply probability transforms (here: empirically estimated dfs)
U <- pobs(X) # probability transformation with the empirically estimated margins
pairs2(U, cex = 0.4, col = adjustcolor("black", alpha.f = 0.5))

## Map the t copula sample to one with N(0,1) and Par(3) margins
Y <- cbind(qPar(U[,1], shape = 3), qnorm(U[,2:5])) # quantile transformation
pairs2(Y, cex = 0.4, col = adjustcolor("black", alpha.f = 0.5),
       labels.null.lab = "Y")
## ... and the components have the same dependence structure as those of X above!

plot(X[,1:2], xlab = expression(Y[1]), ylab = expression(Y[2]), col = "royalblue3") # scatter plot
points(Y[,1:2])
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
       legend = c(expression(Par(3)),
                  substitute("With"~italic(t)[nu.]~"margins", list(nu. = nu))))

plot(X[,4:5], xlab = expression(Y[4]), ylab = expression(Y[5]), col = "royalblue3") # scatter plot
points(Y[,4:5])
legend("bottomright", bty = "n", pch = c(1,1), col = c("black", "royalblue3"),
       legend = c(expression(N(0,1)),
                  substitute("With"~italic(t)[nu.]~"margins", list(nu. = nu))))

#### Sklar_theorem_and_invariance_principle ####
## Visually explaining Sklar's Theorem and the invariance principle

library(mvtnorm)

set.seed(271)

## Sample from a t copula
n <- 1000 # sample size
d <- 2 # dimension
rho <- 0.7 # off-diagonal entry in the correlation matrix P
P <- matrix(rho, nrow = d, ncol = d) # build the correlation matrix P
diag(P) <- 1
nu <- 3.5 # degrees of freedom
set.seed(271)
X <- rmvt(n, sigma = P, df = nu) # n multiv. t observations
U <- pt(X, df = nu) # n ind. realizations from the corresponding t copula
Y <- qexp(U) # transform U (t copula) to Exp(1) margins

## Plot setup
ind <- c(A = 119, B = 516, C = 53) # use 'plot(X); identify(X)' to identify these points
cols <- c("royalblue3", "maroon3", "darkorange2")
par(pty = "s")

## Scatter plot of a bivariate t distribution
plot(X, xlab = quote(X[1]), ylab = quote(X[2]))
points(X[ind["A"],1], X[ind["A"],2], pch = 21, bg = cols[1], col = cols[1])
points(X[ind["B"],1], X[ind["B"],2], pch = 21, bg = cols[2], col = cols[2])
points(X[ind["C"],1], X[ind["C"],2], pch = 21, bg = cols[3], col = cols[3])
text(X[ind["A"],1], X[ind["A"],2], labels = "A", adj = c(0.5, -0.75), col = cols[1])
text(X[ind["B"],1], X[ind["B"],2], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(X[ind["C"],1], X[ind["C"],2], labels = "C", adj = c(0.5, -0.75), col = cols[3])

## Scatter plot of the corresponding t copula (1st part of Sklar's Theorem)
plot(U, xlab = quote(U[1]), ylab = quote(U[2]))
points(U[ind["A"],1], U[ind["A"],2], pch = 21, bg = cols[1], col = cols[1])
points(U[ind["B"],1], U[ind["B"],2], pch = 21, bg = cols[2], col = cols[2])
points(U[ind["C"],1], U[ind["C"],2], pch = 21, bg = cols[3], col = cols[3])
text(U[ind["A"],1], U[ind["A"],2], labels = "A", adj = c(0.5, -0.75), col = cols[1])
text(U[ind["B"],1], U[ind["B"],2], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(U[ind["C"],1], U[ind["C"],2], labels = "C", adj = c(0.5, -0.75), col = cols[3])

## Scatter plot of the meta-t distribution with Exp(1) margins (2nd part of Sklar's Theorem)
plot(Y, xlab = quote(Y[1]), ylab = quote(Y[2]))
points(Y[ind["A"],1], Y[ind["A"],2], pch = 21, bg = cols[1], col = cols[1])
points(Y[ind["B"],1], Y[ind["B"],2], pch = 21, bg = cols[2], col = cols[2])
points(Y[ind["C"],1], Y[ind["C"],2], pch = 21, bg = cols[3], col = cols[3])
text(Y[ind["A"],1], Y[ind["A"],2], labels = "A", adj = c(0.5, -0.75), col = cols[1])
text(Y[ind["B"],1], Y[ind["B"],2], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(Y[ind["C"],1], Y[ind["C"],2], labels = "C", adj = c(0.5, -0.75), col = cols[3])
## We see that the *relative* locations of all points remains the same.
## We thus only change the marginal distributions, but not the dependence
## (by applying the probability and quantile transformations). This also
## visually confirms the invariance principle.

#### Examples_of_copulas ####
## Selected copulas and copula families

library(copula)

n <- 1000 # sample size
d <- 5 # max. considered dimension
tau <- 0.5 # Kendall's tau
tblack <- function(alpha.f) adjustcolor("black", alpha.f = alpha.f) # color (transparent black)

set.seed(271) # for reproducibility


### 1 Fundamental copulas ##

### 1.1 Independence copula ##

## Define the independence copula object
ic <- indepCopula()

## Copula (wire frame and level curves)
wireframe2(ic, FUN = pCopula)
contourplot2(ic, FUN = pCopula)

## Copula density (wire frame)
wireframe2(ic, FUN = dCopula, delta = 0.001, zlim = 0:1)
contourplot2(ic, FUN = dCopula, delta = 0.001, zlim = 0:1)

## Scatter plots
U <- matrix(runif(n*d), ncol = 5)
plot(U[,1:2], xlab = expression(U[1]), ylab = expression(U[2])) # d = 2
cloud2(U[,1:3], # d = 3
       xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(U) # d = 5


### 1.2 Frechet--Hoeffding bounds W and M ##

## Copulas (wire frame and level curves)
n.grid <- 26 # number of grid points
u <- seq(0, 1, length.out = n.grid) # subdivison points in each dimension
grid <- expand.grid("u[1]" = u, "u[2]" = u) # build a grid
W <- pmax(grid[,1] + grid[,2] - 1, 0) # values of W on grid
M <- pmin(grid[,1], grid[,2]) # values of M on grid
val.W <- cbind(grid, "W(u[1],u[2])" = W) # append grid
val.M <- cbind(grid, "M(u[1],u[2])" = M) # append grid
wireframe2(val.W) # wire frame plot of W

contourplot2(val.W, xlim = 0:1, ylim = 0:1) # level curves of W
wireframe2(val.M) # wire frame plot of M
contourplot2(val.M, xlim = 0:1, ylim = 0:1) # level curves of M

## Scatter plots
U <- runif(n)
plot(cbind(U, 1-U), xlab = expression(U[1]), ylab = expression(U[2]), col = tblack(0.5)) # sample of W for d = 2
plot(cbind(U, U), xlab = expression(U[1]), ylab = expression(U[2]), col = tblack(0.5)) # sample of M for d = 2
cloud2(do.call(cbind, rep(list(U), 3)), col = tblack(0.05), # sample of M for d = 3
       xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(do.call(cbind, rep(list(U), d)), col = tblack(0.05)) # sample of M for d = 5


### 2 Implicit copulas ##

### 2.1 Normal (or: Gauss) copula ##

## Define the normal copula object
th <- iTau(normalCopula(), tau = 0.5)
nc <- normalCopula(th)

## Copula (wire frame and level curves)
wireframe2(nc, FUN = pCopula)
contourplot2(nc, FUN = pCopula)

## Copula density (wire frame and level curves)
wireframe2(nc, FUN = dCopula, delta = 0.02)
contourplot2(nc, FUN = dCopula)

## Scatter plots
nc. <- normalCopula(th, dim = 5)
U <- rCopula(n, copula = nc.)
plot(U[,1:2], xlab = expression(U[1]), ylab = expression(U[2])) # d = 2
cloud2(U[,1:3], # d = 3
       xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(U, cex = 0.4, col = tblack(0.5)) # d = 5


### 2.2 t copula ##

## Define the t copula object
nu <- 3 # degrees of freedom
th <- iTau(tCopula(, df = nu), tau = tau) # correlation parameter
tc <- tCopula(th, df = nu)

## Copula (wire frame and level curves)
## Note: pCopula() is only available for integer degrees of freedom
wireframe2(tc, FUN = pCopula)
contourplot2(tc, FUN = pCopula)

## Copula density (wire frame and level curves)
wireframe2(tc, FUN = dCopula, delta = 0.02)
contourplot2(tc, FUN = dCopula)

## Scatter plots
tc. <- tCopula(th, dim = 5, df = nu)
U <- rCopula(n, copula = tc.)
plot(U[,1:2], xlab = expression(U[1]), ylab = expression(U[2])) # d = 2
cloud2(U[,1:3], # d = 3
       xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(U, cex = 0.4, col = tblack(0.5)) # d = 5

## A nonexchangeable example
## Build a block correlation matrix
rho <- c(0.3, 0.6, 0.9)
P <- matrix(rho[1], nrow = d, ncol = d)
P[1:2, 1:2] <- rho[2]
P[3:5, 3:5] <- rho[3]
diag(P) <- 1
## Define, sample and plot a t copula
tc.. <- tCopula(P2p(P), dim = d, dispstr = "un", df = 3.5)
U. <- rCopula(n, copula = tc..)
pairs2(U., cex = 0.4, col = tblack(0.5)) # d = 5

## A model more flexible than the multivariate t (more degrees of freedom)
X <- sapply(1:ncol(U.), function(j) qt(U.[,j], df = j))
pairs2(X, labels.null.lab = "X", cex = 0.4, col = tblack(0.5))


### 3 Explicit copulas ##

### 3.1 Clayton copula ##

## Define the Clayton copula object
th <- iTau(claytonCopula(), tau = tau)
cc <- claytonCopula(th)

## Copula (wire frame and level curves)
wireframe2(cc, FUN = pCopula)
contourplot2(cc, FUN = pCopula)

## Copula density (wire frame and level curves)
wireframe2(cc, FUN = dCopula, delta = 0.02)
contourplot2(cc, FUN = dCopula)

## Scatter plots
cc. <- claytonCopula(th, dim = 5)
U <- rCopula(n, copula = cc.)
plot(U[,1:2], xlab = expression(U[1]), ylab = expression(U[2])) # d = 2
cloud2(U[,1:3], # d = 3
       xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(U, cex = 0.4, col = tblack(0.5)) # d = 5


### 3.2 Gumbel copula ##

## Define the Gumbel copula object
th <- iTau(gumbelCopula(), tau = tau)
gc <- gumbelCopula(th)

## Copula (wire frame and level curves)
wireframe2(gc, FUN = pCopula)
contourplot2(gc, FUN = pCopula)

## Copula density (wire frame and level curves)
wireframe2(gc, FUN = dCopula, delta = 0.02)
contourplot2(gc, FUN = dCopula)

## Scatter plots
gc. <- gumbelCopula(th, dim = 5)
U <- rCopula(n, copula = gc.)
plot(U[,1:2], xlab = expression(U[1]), ylab = expression(U[2])) # d = 2
cloud2(U[,1:3], # d = 3
       xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(U, cex = 0.4, col = tblack(0.5)) # d = 5

## Scatter plot of a survival Gumbel copula
pairs2(1-U, cex = 0.4, col = tblack(0.5))

## Only the first component flipped
pairs2(U, cex = 0.4, col = tblack(0.5)) # original
pairs2(cbind(1-U[,1], U[,2:5]), cex = 0.4, col = tblack(0.5)) # first flipped


### 3.3 Outer power Clayton copula ##

## Note: Outer power Clayton copulas can have both lower and upper tail dependence

## Define the outer power Clayton copula
tauC <- 0.3 # Kendall's tau for the underlying Clayton copula
thC <- copClayton@iTau(tauC) # choose Clayton's generator s.t. Kendall's tau is tauC
opC <- opower(copClayton, thC) # define an outer power Clayton copula (its parameter theta is not specified yet)
th <- opC@iTau(tau) # define the outer power Clayton copula s.t. has Kendall's tau is tau
opcc <- onacopulaL(opC, list(th, 1:2)) # define the outer power Clayton copula

## Scatter plot
U <- rCopula(n, copula = opcc)
plot(U[,1:2], xlab = expression(U[1]), ylab = expression(U[2])) # d = 2


### 3.4 Marshall--Olkin copulas ##

## Note: Marshall--Olkin copulas have a singular component, i.e., a set of
##       Lebesgue measure 0 with positive probability mass assigned.

## Define the MO copula
C <- function(u, alpha) {
  if(!is.matrix(u)) u <- rbind(u)
  stopifnot(0 <= alpha, alpha <= 1, length(alpha) == 2)
  pmin(u[,1]*u[,2]^(1-alpha[2]), u[,1]^(1-alpha[1])*u[,2])
}

## Sampling
rMO <- function(n, alpha) {
  stopifnot(n >= 1, 0 <= alpha, alpha <= 1, length(alpha) == 2)
  U. <- matrix(runif(n*3), ncol = 3)
  U <- 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])))
}

## Define the singular component
S.C <- function(u, alpha) {
  stopifnot(0 <= u, u <= 1,
            0 <= alpha, alpha <= 1, length(alpha) == 2)
  tau <- alpha[1] * alpha[2] / (alpha[1] + alpha[2] - alpha[1] * alpha[2])
  tau * pmin(u[,1]^alpha[1], u[,2]^alpha[2])^(1/tau)
}

## Define the absolutely continuous component
A.C <- function(u, alpha) C(u, alpha) - S.C(u, alpha)


## Scatter plot
alpha <- c(0.2, 0.8)
U <- rMO(n, alpha = alpha)
plot(U, xlab = expression(U[1]), ylab = expression(U[2]))
## Interpretation: Given u_1 = 0.6, for example, u_2 lies with a non-zero
## probability p on the curve and with the remaining probability 1-p
## anywhere else (*not*: uniform) along the line u_1 = 0.6.

## Check the margins
plot(U[,1], ylab = expression(U[1]))
plot(U[,2], ylab = expression(U[2]))

## Evaluate the copula (and its singular and absolutely continuous components)
grid <- expand.grid("u[1]" = u, "u[2]" = u) # build a grid
val.C <- cbind(grid, "C(u[1],u[2])" = C(grid, alpha = alpha)) # append C
val.S <- cbind(grid, "S[C](u[1],u[2])" = S.C(grid, alpha = alpha)) # append S.C
val.A <- cbind(grid, "A[C](u[1],u[2])" = A.C(grid, alpha = alpha)) # append S.C

## Copula (wire frame and level curves)
wireframe2(val.C) # wire frame plot
contourplot2(val.C, xlim = 0:1, ylim = 0:1) # level curves
## The copula has a kink, that is, a so-called singular component.
## A singular component is a set of Lebesgue measure 0 where C puts
## mass at (here: a curve). This is better visible from the scatter plots below.

## Singular and absolutely continuous component (wire frame)
wireframe2(val.S) # singular component
wireframe2(val.A) # absolutely continuous component

#### Meta_copula_models ####
## Illustrating meta-C models

library(copula)
library(gridExtra) # for grid.arrange()

n <- 1000 # sample size
set.seed(271)


### 1 Generate the data ##

## Define the copula parameters
## Note: They are chosen such that Pearson's correlation coefficient,
##       when computed with N(0,1) margins, is 0.7 (see below)
th.n <- 0.7 # Gauss copula parameter
th.g <- 2 # Gumbel copula parameter
th.c <- 2.2 # Clayton copula parameter
th.t <- 0.71 # t_4 copula parameter

## Define the copulas
nc <- normalCopula(th.n) # Gauss copula
gc <- gumbelCopula(th.g) # Gumbel copula
cc <- claytonCopula(th.c) # Clayton copula
tc <- tCopula(th.t) # t_4 copula

## Generate copula data
U.nc <- rCopula(n, copula = nc)
U.gc <- rCopula(n, copula = gc)
U.cc <- rCopula(n, copula = cc)
U.tc <- rCopula(n, copula = tc)

## Map to N(0,1) margins (meta-copula data)
X.nc <- qnorm(U.nc)
X.gc <- qnorm(U.gc)
X.cc <- qnorm(U.cc)
X.tc <- qnorm(U.tc)

## Correlations
cors <- vapply(list(X.nc, X.gc, X.cc, X.tc), function(x) cor(x)[1,2], NA_real_)
stopifnot(all.equal(cors, rep(0.7, 4), tol = 0.015))

## Define the corresponding (meta-C) densities (via Sklar's Theorem)
## Note: Density f(x_1, x_2) = c(F_1(x_1), F_2(x_2)) * f_1(x_1) * f_2(x_2)
##                           = exp( log(c(F_1(x_1), F_2(x_2))) + log(f_1(x_1)) + log(f_2(x_2)) )
dMetaCopulaN01 <- function(x, copula)
  exp(dCopula(pnorm(x), copula = copula, log = TRUE) + rowSums(dnorm(x, log = TRUE)))
## Alternatively, we could work with dMvdc() here


### 2 Scatter plots ##

## Copula samples (see Figure 7.3)
opar <- par(pty = "s", mar = c(5.1, 4.1, 4.1, 2.1) - 1)
lay <- matrix(1:4, ncol = 2, byrow = TRUE) # layout matrix
layout(lay) # layout
plot(U.nc, xlab = expression(U[1]), ylab = expression(U[2]), # Gauss copula
     cex = 0.4, main = "Gauss copula sample")
plot(U.gc, xlab = expression(U[1]), ylab = expression(U[2]), # Gumbel copula
     cex = 0.4, main = "Gumbel copula sample")
plot(U.cc, xlab = expression(U[1]), ylab = expression(U[2]), # Clayton copula
     cex = 0.4, main = "Clayton copula sample")
plot(U.tc, xlab = expression(U[1]), ylab = expression(U[2]), # t_4 copula
     cex = 0.4, main = expression(bold(italic(t)[4]~"copula sample")))

## Meta-copula samples with N(0,1) margins (see Figure 7.4)
m <- max(abs(X.nc), abs(X.gc), abs(X.cc), abs(X.tc))
lim <- c(-m, m)
plot(X.nc, xlim = lim, ylim = lim, # meta-Gauss
     xlab = expression(X[1]), ylab = expression(X[2]),
     cex = 0.4, main = "Meta-Gauss sample")
mtext("N(0,1) margins", side = 4, line = 0.25, adj = 0, cex = 0.7)
plot(X.gc, xlim = lim, ylim = lim,
     xlab = expression(X[1]), ylab = expression(X[2]), # meta-Gumbel
     cex = 0.4, main = "Meta-Gumbel sample")
mtext("N(0,1) margins", side = 4, line = 0.25, adj = 0, cex = 0.7)
plot(X.cc, xlim = lim, ylim = lim, # meta-Clayton
     xlab = expression(X[1]), ylab = expression(X[2]),
     cex = 0.4, main = "Meta-Clayton sample")
mtext("N(0,1) margins", side = 4, line = 0.25, adj = 0, cex = 0.7)
plot(X.tc, xlim = lim, ylim = lim, # meta-t_4
     xlab = expression(X[1]), ylab = expression(X[2]),
     cex = 0.4, main = expression(bold("Meta-"*italic(t)[4]~"sample")))
mtext("N(0,1) margins", side = 4, line = 0.25, adj = 0, cex = 0.7)
par(opar) # restore graphical parameters


### 3 Density plots ##

## Wire frame plots
n.grid <- 26 # number of grid points
lim <- c(-2.5, 2.5)
s <- seq(lim[1], lim[2], length.out = n.grid)
grid <- as.matrix(expand.grid("x[1]" = s, "x[2]" = s))
val.n <- cbind(grid, "f(x[1],x[2])" = dMetaCopulaN01(grid, copula = nc))
val.g <- cbind(grid, "f(x[1],x[2])" = dMetaCopulaN01(grid, copula = gc))
val.c <- cbind(grid, "f(x[1],x[2])" = dMetaCopulaN01(grid, copula = cc))
val.t <- cbind(grid, "f(x[1],x[2])" = dMetaCopulaN01(grid, copula = tc))
zlim <- c(0, max(val.n[,3], val.g[,3], val.c[,3], val.t[,3]))

## Density plots
zm <- 1.1 # zoom in a bit more
scs <- list(arrows = FALSE, col = "black", cex = 0.6) # scale down ticks
tcx <- 0.95 # scale back titles
w.n <- wireframe2(val.n, xlim = lim, ylim = lim, zlim = zlim, zoom = zm, scales = scs,
                  main = list(label = "Meta-Gauss density with N(0,1) margins", cex = tcx))
w.g <- wireframe2(val.g, xlim = lim, ylim = lim, zlim = zlim, zoom = zm, scales = scs,
                  main = list(label = "Meta-Gumbel density with N(0,1) margins", cex = tcx))
w.c <- wireframe2(val.c, xlim = lim, ylim = lim, zlim = zlim, zoom = zm, scales = scs,
                  main = list(label = "Meta-Clayton density with N(0,1) margins", cex = tcx))
w.t <- wireframe2(val.t, xlim = lim, ylim = lim, zlim = zlim, zoom = zm, scales = scs,
                  main = list(label = expression(bold("Meta-"*italic(t)[4]~"density with N(0,1) margins")),
                              cex = tcx))
grid.arrange(w.n, w.g, w.c, w.t, ncol = 2)

## Contour plots
tcx <- 0.9
c.n <- contourplot2(val.n, xlim = lim, ylim = lim,
                    main = list(label = "Meta-Gauss contours with N(0,1) margins", cex = tcx))
c.g <- contourplot2(val.g, xlim = lim, ylim = lim,
                    main = list(label = "Meta-Gumbel contours with N(0,1) margins", cex = tcx))
c.c <- contourplot2(val.c, xlim = lim, ylim = lim,
                    main = list(label = "Meta-Clayton contours with N(0,1) margins", cex = tcx))
c.t <- contourplot2(val.t, xlim = lim, ylim = lim,
                    main = list(label = expression(bold("Meta-"*italic(t)[4]~"contours with N(0,1) margins")),
                                cex = tcx))
grid.arrange(c.n, c.g, c.c, c.t, ncol = 2)

#### Correlation_pitfalls ####
## Selected correlation pitfalls

library(copula)

### Fallacy 1 (Ex. 1): F_1, F_2 and correlation rho uniquely determine F ##

## Simple example with same margins, rho = 0, but (obviously) different models
n <- 1000
set.seed(271)
Z <- rnorm(n)
U <- runif(n)
V <- rep(1, n)
V[U < 1/2] <- -1 # => V in {-1,1}, each with probability 1/2
X <- cbind(Z, Z*V) # convex combination of W and M (each with prob. 1/2 => rho = 0)
cor(X)[2,1]
Y <- matrix(rnorm(n * 2), ncol = 2) # independent N(0,1) (=> rho = 0)
cor(Y)[2,1]

## Plots
plot(X, xlab = expression(X[1]), ylab = expression(X[2]))
plot(Y, xlab = expression(Y[1]), ylab = expression(Y[2]))


### Fallacy 1 (Ex. 2): F_1, F_2 and correlation rho uniquely determine F ##

## We construct a parametric copula family here which *always* produces
## correlation 0. To this end, let C(u1,u2) = u1 * u2 + f1(u1) * f2(u2),
## for f1, f2 continuously differentiable s.t. f1(0) = f2(0) = f1(1) = f2(1) = 0
## and s.t. 1 + f1'(u1) * f2'(u2) >= 0 for all u1, u2 in [0,1].
## By computing C's density, C can be seen to be a copula. Let (U_1,U_2) ~ C.
## As the margins are U(0,1), Var(U_1) = Var(U_2) = 1/12 and thus Pearson's
## correlation coefficient rho equals rho = 12 Cov(U_1, U_2). By Hoeffding's
## identity, rho = 12 \int_0^1 \int_0^1 (C(u1,u2) - u1 * u2) du1 du2. Plugging
## in our form of C, we obtain rho = 12 \int_0^1 f1(u1) du1 * \int_0^1 f2(u2) du2.
## So if one of the f's is point symmetric about 1/2, rho equals 0. For example,
## take f1(x) = 2x(x-1/2)(x-1) and f2(x) = theta * x * (1-x). Then, for any
## theta in [-1,1], f1 and f2 fulfill all conditions and will generate a
## proper copula with Pearson's correlation coefficient equal to 0.
## Let's do that!

## Define auxiliary functions
f1 <- function(x) 2*x*(x-1/2)*(x-1) # f1
f1. <- function(x) 6*x^2-6*x+1 # derivative of f1
f2 <- function(x, th) th * x * (1-x) # f2
f2. <- function(x, th) th * (-2 * x + 1) # derivative of f2

## Define the conditional copula C(u2 | u_1)
cCop <- function(u, th) {
  if(!is.matrix(u)) u <- rbind(u)
  stopifnot(-1 <= th, th <= 1)
  u[,2] + f1.(u[,1]) * f2(u[,2], th = th)
}

## Density
dCop <- function(u, th) {
  if(!is.matrix(u)) u <- rbind(u)
  stopifnot(-1 <= th, th <= 1)
  1 + f1.(u[,1]) * f2.(u[,2], th = th)
}

## Conditional distribution method for this copula
## Note; We numerically invert C(u2 | u_1) to get C^{-1}(u2 | u_1)
CDM <- function(u, th) {
  if(!is.matrix(u)) u <- rbind(u)
  cbind(u[,1], apply(u, 1, function(u.)
    uniroot(function(u2) cCop(c(u.[1], u2), th = th) - u.[2],
            interval = 0:1, maxiter = 1000)$root))
}

## Sample
th <- 1
n <- 10000
set.seed(1)
U <- CDM(matrix(runif(2*n), ncol = 2), th = th)
plot(U, xlab = expression(U[1]), ylab = expression(U[2]))

## 'Show' that rho is (close to) 0
stopifnot(all.equal(cor(U)[2,1], 0.0016, tol = 1e-2))

## Visualize the copula density
n.grid <- 26 # number of grid points
u <- seq(0, 1, length.out = n.grid) # subdivison points in each dimension
grid <- expand.grid("u[1]" = u, "u[2]" = u) # build a grid
val.c <- cbind(grid, "c(u[1],u[2])" = dCop(grid, th = th)) # evaluate density on the grid
wireframe2(val.c) # wire frame plot of the density
contourplot2(val.c, xlim = 0:1, ylim = 0:1) # level curves


### Fallacy 2: Given F_1, F_2, any rho in [-1,1] is attainable ##

## Function to compute the correlation bounds for LN(0, sigma_.^2) margins
## Note: The derivation can be done with the moment-generating function of
##       the standard normal distribution
cor_bound_LN <- function(s, method = c("max", "min")) {
  ## s = (sigma_1, sigma_2)
  if(!is.matrix(s)) s <- rbind(s)
  method <- match.arg(method)
  if(method == "min") s[,2] <- -s[,2]
  (exp((s[,1]+s[,2])^2/2)-exp((s[,1]^2+s[,2]^2)/2)) /
    sqrt(expm1(s[,1]^2)*exp(s[,1]^2)*expm1(s[,2]^2)*exp(s[,2]^2))
}

## Evaluate and plot correlation bounds on a grid
n.grid <- 26 # number of grid points in each dimension
s <- seq(0.01, 5, length.out = n.grid) # subdivision points in each dimension
grid <- expand.grid("sigma[1]" = s, "sigma[2]" = s) # build a grid
val.min <- cbind(grid, "underline(Cor)(sigma[1],sigma[2])" =
                   cor_bound_LN(grid, method = "min"))
val.max <- cbind(grid, "bar(Cor)(sigma[1],sigma[2])" =
                   cor_bound_LN(grid))
wireframe2(val.min)
wireframe2(val.max)


## Another example of this type: X_1, X_2 being (e.g., default) indicators

## Function to compute the correlation bounds for B(1, p_.) margins
## Note: Using F_j (B(1,p_j) df), one obtains F_j^{-1}(u) = I_{(1-p_j,1]}(u) for
##       all u in (0, 1] and thus X_j = I_{(1-p_j,1]}(U_j) for (U_1, U_2) ~ C.
##       Hence, E(X_1 X_2) = E(I_{(1-p_1,1] x (1-p_2,1]}(U_1, U_2))
##       = P(U_1 > 1-p_1, U_2 > 1-p_2) = P(1-U_1 <= p_1, 1-U_2 <= p_2)
##       = bar(C)(p_1, p_2) with bounds W(p_1, p_2) and M(p_1, p_2).
##       Together with E(X_j) = P(U_j > 1-p_j) = p_j and Var(X_j) = p_j*(1-p_j),
##       one thus obtains the correlation bounds.
cor_bound_B <- function(p, method = c("max", "min")) {
  ## p = (p_1, p_2)
  if(!is.matrix(p)) p <- rbind(p)
  method <- match.arg(method)
  C.bounds <- if(method == "min") pmax(p[,1]+p[,2]-1, 0) else pmin(p[,1], p[,2])
  (C.bounds - p[,1]*p[,2]) / sqrt(p[,1]*(1-p[,1]) * p[,2]*(1-p[,2]))
}

## Evaluate and plot correlation bounds on a grid
n.grid <- 26 # number of grid points in each dimension
p <- seq(0.01, 0.99, length.out = n.grid) # subdivision points in each dimension
grid <- expand.grid("p[1]" = p, "p[2]" = p) # build a grid
val.min <- cbind(grid, "underline(Cor)(p[1],p[2])" =
                   cor_bound_B(grid, method = "min"))
val.max <- cbind(grid, "bar(Cor)(p[1],p[2])" =
                   cor_bound_B(grid))
wireframe2(val.min)
wireframe2(val.max)

#### Copula_tail_probabilities ####
## Computing and plotting tail probabilities

library(copula)

set.seed(271) # due to the Monte Carlo (MC) error for evaluating Ga/t copulas


### Example 7.40 (Financial example with five stocks) ##

d <- 5
cop <- ellipCopula("normal", param=0.5, dim=d)
pCopula(rep(0.01, d), copula=cop) # ~= 7.44e-05 (for our seed)


### Table 7.2 (Joint tail probabilities P(U_1 > u, U_2 > u) as a plot ##

n <- 128
u <- seq(0.95, to=0.9999, length.out=n) # quantiles u where to evaluate P(U_1 > u, U_2 > u)
rho <- c(0.75, 0.5) # correlation parameter rho
nu <- c(3, 4, 8, Inf) # degrees of freedom
len <- length(rho)*length(nu)
col <- c("maroon3", "darkorange2", "royalblue3", "black") # colors for different nu
tail.prob <- matrix(, nrow=length(u), ncol=1+len, # tail probabilities
                    dimnames=list(rownames=seq_along(u),
                                  colnames=c("u", "rho.075.nu.3", "rho.075.nu.4",
                                             "rho.075.nu.8", "rho.075.nu.Inf",
                                             "rho.05.nu.3", "rho.05.nu.4",
                                             "rho.05.nu.8", "rho.05.nu.Inf")))
tail.prob[,1] <- u
expr <- vector("expression", len) # vector of expressions
ltys <- numeric(len) # line types
cols <- numeric(len) # colors
for(i in seq_along(rho)) { # rho
  for(j in seq_along(nu)) { # degrees of freedom
    k <- length(nu)*(i-1)+j
    ## Create the copula
    cop <- ellipCopula("t", param=rho[i], df=nu[j])
    ## Evaluate P(U_1 > u, U_2 > u) = P(U_1 <= 1-u, U_2 <= 1-u)
    tail.prob[,k+1] <- pCopula(cbind(1-u, 1-u), copula=cop)
    ## Create plot information
    expr[k] <- as.expression(
      substitute(group("(",list(rho, nu),")")==group("(", list(r., n.), ")"),
                 list(r.=rho[i], n.=if(nu[j]==Inf) bquote(infinity) else nu[j])))
    ltys[k] <- length(rho)-i+1
    cols[k] <- col[j]
  }
}

## Standardize w.r.t. Gauss case
tail.prob.fact <- tail.prob # tail probabilities as factors in comparison to Gauss case
tail.prob.fact[,2:5] <- tail.prob.fact[,2:5]/tail.prob.fact[,5]
tail.prob.fact[,6:9] <- tail.prob.fact[,6:9]/tail.prob.fact[,9]

## Plot tail probabilities
matplot(tail.prob[,1], tail.prob[,-1], type="l", lty=ltys, col=cols,
        xlab="Quantile u", ylab=expression("Tail probability"~P(U[1]>u,U[2]>u)))
legend("topright", inset=0.01, bty="n", lty=ltys, col=cols, legend=expr)

## Plot standardized tail probabilities
matplot(tail.prob.fact[,1], tail.prob.fact[,-1], log="y", type="l", lty=ltys, col=cols,
        lwd=c(1,1,1,1.6,1,1,1,1), xlab="Quantile u",
        ylab=expression("Tail probability"~P(U[1]>u,U[2]>u)~"standardized by Gauss case"))
legend("topleft", inset=0.01, bty="n", lty=ltys, col=cols, legend=expr)


### Table 7.3 (Joint tail probabilities P(U_1 > u, ..., U_d > u) as a plot ##

## Note: Due to radial symmetry, we can compute probabilities in the lower tail
d <- 2:20
u <- 0.99 # quantile u where to evaluate P(U_1 > u, ..., U_d > u)
rho <- c(0.75, 0.5) # correlation parameter rho
nu <- c(3, 4, 8, Inf) # degrees of freedom
len <- length(rho)*length(nu)
col <- c("maroon3", "darkorange2", "royalblue3", "black") # colors for different nu
tail.prob <- matrix(, nrow=length(d), ncol=1+len, # tail probabilities
                    dimnames=list(rownames=seq_along(d),
                                  colnames=c("d", "rho.075.nu.3", "rho.075.nu.4",
                                             "rho.075.nu.8", "rho.075.nu.Inf",
                                             "rho.05.nu.3", "rho.05.nu.4",
                                             "rho.05.nu.8", "rho.05.nu.Inf")))
tail.prob[,1] <- d
expr <- vector("expression", len) # vector of expressions
ltys <- numeric(len) # line types
cols <- numeric(len) # colors
for(i in seq_along(rho)) { # rho
  for(j in seq_along(nu)) { # degrees of freedom
    k <- length(nu)*(i-1)+j
    tail.prob <- cbind(tail.prob, rep(NA, length(d)))
    for(l in seq_along(d)) { # dimension
      ## Create the copula
      cop <- ellipCopula("t", param=rho[i], dim=d[l], df=nu[j])
      ## Evaluate P(U_1 > u, ..., U_d > u) = P(U_1 <= 1-u, ..., U_d <= 1-u)
      tail.prob[l,k+1] <- pCopula(rep(1-u, d[l]), copula=cop)
    }
    ## Create plot information
    expr[k] <- as.expression(
      substitute(group("(", list(rho, nu), ")")==group("(", list(r., n.), ")"),
                 list(r.=rho[i], n.=if(nu[j]==Inf) bquote(infinity) else nu[j])))
    ltys[k] <- length(rho)-i+1
    cols[k] <- col[j]
  }
}

## Standardize w.r.t. Gauss case
tail.prob.fact <- tail.prob # tail probabilities as factors in comparison to Gauss case
tail.prob.fact[,2:5] <- tail.prob.fact[,2:5]/tail.prob.fact[,5]
tail.prob.fact[,6:9] <- tail.prob.fact[,6:9]/tail.prob.fact[,9]

## Plot tail probabilities
matplot(tail.prob[,1], tail.prob[,-1], type="l", lty=ltys, col=cols,
        xlab="Dimension d", ylab=expression("Tail probability"~P(U[1]>u,...,U[d]>u)))
legend("topright", inset=0.01, bty="n", lty=ltys, col=cols, legend=expr)

## Plot standardized tail probabilities
matplot(tail.prob.fact[,1], tail.prob.fact[,-1], log="y", type="l", lty=ltys, col=cols,
        lwd=c(1,1,1,1.6,1,1,1,1), xlab="Dimension d",
        ylab=expression("Tail probability"~P(U[1]>u,..,U[d]>u)~"standardized by Gauss case"))
legend("topleft", inset=0.01, bty="n", lty=ltys, col=cols, legend=expr)

## Results:
##
## 'Table 7.2' case:
## - Obviously, P(U_1 > u, U_2 > u) is decreasing in u.
## - The larger rho or the smaller nu, the larger P(U_1 > u, U_2 > u).
## - Divided by the tail probability for the Gauss copula, P(U_1 > u, U_2 > u) is...
##   + ... increasing in u (the further we are in the tail, the more pronounced the difference);
##   + ... larger the smaller nu (the more we deviate from the Gauss copula);
##   + ... larger the smaller rho.
##
## 'Table 7.3' case:
## - Obviously, P(U_1 > u, ..., U_d > u) is decreasing in d.
## - The larger rho or the smaller nu, the larger P(U_1 > u, ..., U_d > u).
## - Divided by the tail probability for the Gauss copula, P(U_1 > u, ..., U_d > u) is...
##   + ... increasing in d (the larger the dimension, the more pronounced the difference);
##   + ... larger the smaller nu (the more we deviate from the Gauss copula);
##   + ... larger the smaller rho.

#### Copula_estimation ( estimation) ####
## Basic copula estimation (2d example, 4 copulas)

library(xts)
library(copula)
library(qrmdata)
library(qrmtools)


### 1 Working with the data ##

## Load the time series
data("SP500")
data("FTSE")

## Build negative log-returns
X.SP500 <- -returns(SP500)
X.FTSE  <- -returns(FTSE)

## Merge
X <- merge(X.SP500, X.FTSE, all = FALSE)
time <- c("2003-01-01", "2012-12-31")
X <- X[paste0(time, collapese = "/"),]

## Basic plot
plot.zoo(X)

## Observe that there are zero values caused by market closures
apply(X == 0, 2, sum)
rexcl <- (X[,1] == 0) | (X[,2] == 0)
X. <- X[!rexcl,]

## Aggregating by week
dim(X.w <- apply.weekly(X., FUN = colSums))
plot.zoo(X.w, type = "h")
plot(as.numeric(X.w[,1]), as.numeric(X.w[,2]))
# plot(qexp(U[,1]), qexp(U[,2]))

## Compute pseudo-observations
U <- as.matrix(pobs(X.w))
plot(U, xlab = expression(U[1]), ylab = expression(U[2]))

### 2 Fitting copulas ##

## Compare various bivariate copulas
fit.N <- fitCopula(normalCopula(),  data = U)
fit.t <- fitCopula(tCopula(),       data = U) # df of freedom are estimated, too
fit.C <- fitCopula(claytonCopula(), data = U)
fit.G <- fitCopula(gumbelCopula(),  data = U)

## Comparing the likelihoods
sort(c(N = fit.N@loglik, t = fit.t@loglik, C = fit.C@loglik, G = fit.G@loglik),
     decreasing = TRUE)
iTau(fit.N, )
fit.N@estimate
fit.t@estimate
fit.C@estimate
fit.G@estimate
## Maximum pseudo-likelihood estimation of P and the d.o.f.
fit.N <- fitCopula(normalCopula(, dispstr = "un"),  data = U)
fit.t <- fitCopula(tCopula(, dispstr = "un"),       data = U) # df of freedom are estimated, too
fit.C <- fitCopula(claytonCopula(), data = U)
fit.G <- fitCopula(gumbelCopula(),  data = U)
sort(c(N = fit.N@loglik, t = fit.t@loglik, C = fit.C@loglik, G = fit.G@loglik),
     decreasing = TRUE)
fit.N@estimate
fit.t@estimate
fit.C@estimate
fit.G@estimate
fit.N@var.est
fit.t@var.est
fit.C@var.est
fit.G@var.est



#
#### Copula_estimation_and_goodness-of-fit (estimation) ####
## Copula estimation and goodness-of-fit (5d example, 2 copulas)

library(rugarch)
library(xts)
library(ADGofTest)
library(qqtest)
library(copula)
library(qrmdata)
library(qrmtools)


### 1 Working with the data ##

## Select the data we work with
data("SP500_const") # load the constituents data of the S&P 500
stocks <- c("INTC", "QCOM", "GOOGL", "AAPL", "MSFT") # Intel, Qualcomm, Google, Apple, Microsoft
time <- c("2007-01-03", "2009-12-31") # time period
S <- SP500_const[paste0(time, collapse = "/"), stocks] # data

## Check for missing data
stopifnot(all(!is.na(S))) # na.fill(, fill = "extend") is often helpful


### 2 Fitting marginal ARMA(1,1)-GARCH(1,1) models with standardized t residuals

## Build negative log-returns
X <- -returns(S) # -log-returns

## Basic plot
plot.zoo(X, main = "Negative log-returns")

## Fit marginal time series
uspec <- rep(list(ugarchspec(distribution.model = "std")), ncol(X))
fit.ARMA.GARCH <- fit_ARMA_GARCH(X, ugarchspec.list = uspec)
stopifnot(sapply(fit.ARMA.GARCH$error, is.null)) # NULL = no error
if(FALSE)
  fit.ARMA.GARCH$warning
## => Warning comes from finding initial values and can be ignored here
fits <- fit.ARMA.GARCH$fit # fitted models
Z <- as.matrix(do.call(merge, lapply(fits, residuals, standardize = TRUE))) # grab out standardized residuals
colnames(Z) <- colnames(S)
(nu.mar <- vapply(fits, function(x) x@fit$coef[["shape"]], NA_real_)) # vector of estimated df
n <- nrow(X) # sample size
d <- ncol(X) # dimension


### 3 Fitting copulas ##

## Compute pseudo-observations from the standardized t residuals
U <- pobs(Z)
pairs2(U, cex = 0.4, col = adjustcolor("black", alpha.f = 0.5))

## Fitting a Gumbel copula
fit.gc <- fitCopula(gumbelCopula(dim = d),
                    data = U, method = "mpl")
fit.gc@estimate # estimated copula parameter
gc <- fit.gc@copula # fitted copula

## Compute matrices of pairwise Kendall's tau and upper tail-dependence coefficients
p2P(tau(gc), d = d)
p2P(lambda(gc)["upper"], d = d)

## Fitting a t copula
fit.tc <- fitCopula(tCopula(dim = d, dispstr = "un"),
                    data = U, method = "itau.mpl")
(nu <- tail(fit.tc@estimate, n = 1)) # estimated degrees of freedom nu
(P <- p2P(head(fit.tc@estimate, n = -1))) # estimated correlation matrix
tc <- fit.tc@copula # fitted copula

## Compute matrices of pairwise Kendall's tau and upper tail-dependence coefficients
p2P(tau(tc))
p2P(lambda(tc)[(choose(d,2)+1):(d*(d-1))])


### 4 Goodness-of-fit ##

## We use the parametric bootstrap here, based on the maximum pseudo-likelihood
## estimator.
set.seed(271) # for reproducibility
N <- 100 # this is to save run time; it should be larger!

## Check the Gumbel copula
gof.gc <- gofCopula(gc, x = U, N = N) # warning also because the copula doesn't fit well here
stopifnot(gof.gc$p.value < 0.05) # => rejection

## Check the t copula
## Note: - This can currently only be done for fixed and integer degrees of
##         freedom as there is no algorithm to evaluate the multivariate t df for
##         non-integer degrees of freedom.
##       - ... yet it's still quite slow. We thus check the model graphically
##         after mapping the variates to a U(0,1)^d distribution with the
##         Rosenblatt transform.
U.Rsnbl <- cCopula(U, copula = tc)
pairs2(U.Rsnbl, cex = 0.4, col = adjustcolor("black", alpha.f = 0.5)) # looks ok

## Map it to a K_d distribution ("kay") and do a Q-Q plot
U.Rsnbl.K <- sqrt(rowMeans(qnorm(U.Rsnbl)^2)) # map to a kay distribution
pK <- function(q, d) pchisq(d*q*q, df = d) # df of a K_d distribution
AD <- ad.test(U.Rsnbl.K, distr.fun = pK, d = d) # compute an AD test
stopifnot(AD$p.value >= 0.05)
## Note: The AD test here does not take into account parameter estimation
##       (that would require a parametric bootstrap, for example)

## A (sophisticated) Q-Q plot
qqtest(U.Rsnbl.K, dist = "kay", df = d, nreps = 1000, pch = 1,
       col = adjustcolor("black", alpha.f = 0.5), main = "",
       xlab = substitute(K[dof]~"quantiles", list(dof = d))) # => looks ok


### 5 Simulating paths from the full model ##

## Simulate paths
B <- 200
m <- ceiling(n/10) # length of the simulates paths
X.lst <- lapply(1:B, function(b) {
  ## 1) Simulate from the fitted copula
  U. <- rCopula(m, copula = tc)
  ## 2) Quantile-transform to standardized t distributions (for ugarchsim())
  Z. <- sapply(1:d, function(j) sqrt((nu.mar[j]-2)/nu.mar[j]) * qt(U.[,j], df = nu.mar[j]))
  ## 3) Use these multivariate dependent t innovations to sample from the
  ##    time series
  sim <- lapply(1:d, function(j)
    ugarchsim(fits[[j]], n.sim = m, m.sim = 1,
              custom.dist = list(name = "sample",
                                 distfit = Z.[,j, drop = FALSE])))
  ## 4) Extract simulated series
  sapply(sim, function(x) fitted(x)) # simulated multivariate series X_t (= x@simulation$seriesSim)
})
## => List of length B containing (n x d)-matrices


### 6 Predict the aggregated loss and VaR_0.99 ##

## Note: - This is merely a demo of what can be done with the simulated data.
##       - See also the vignette "ARMA_GARCH_VaR" in qrmtools

## Predict VaR of the aggregated loss nonparametrically
Xs <- rowSums(X) # aggregated loss; n-vector
Xs. <- sapply(X.lst, rowSums) # simulated aggregated losses; (m, B)-matrix
Xs.mean <- rowMeans(Xs.) # predicted aggregated loss; m-vector
Xs.CI <- apply(Xs., 1, function(x) quantile(x, probs = c(0.025, 0.975))) # CIs; (2, m)-matrix
alpha <- 0.99 # confidence level
VaR <- apply(Xs., 1, function(x) quantile(x, probs = alpha)) # VaR_alpha; m-vector

## Plot
tm <- index(SP500_const)
start <- match(time[1], as.character(tm))
past <- tm[start:(start+n-1)]
future <- tm[(start+n):(start+n+m-1)]
plot(past, Xs, type = "l", xlim = range(c(past, future)), xlab = "", ylab = "") # actual (past) losses
polygon(c(future, rev(future)), c(Xs.CI[1,], rev(Xs.CI[2,])),
        border = NA, col = "grey80") # CI region
lines(future, Xs.mean, col = "royalblue3") # predicted aggregated loss
lines(future, Xs.CI[1,], col = "grey50") # lower CI
lines(future, Xs.CI[2,], col = "grey50") # upper CI
lines(future, VaR, col = "maroon3") # VaR_alpha
legend("bottomright", bty = "n", lty = rep(1, 4),
       col = c("black", "royalblue3", "grey50", "maroon3"),
       legend = c("(Aggregated) loss", "(Simulated) predicted loss",
                  "95% CIs", as.expression(substitute("Simulated"~VaR[a], list(a = alpha)))))

#### Quasi-random_numbers_for_copulas ####
## Pseudo-random numbers vs quasi-random numbers

library(qrng)
library(copula)

n <- 1000 # sample size


### 1 Examples ##

### 1.1 Independence copula ##

## Sample
set.seed(271) # for reproducibility
U <- matrix(runif(n*2), ncol = 2) # pseudo-random numbers
U. <- ghalton(n, d = 2) # quasi-random numbers (faster alternative: sobol())

## Plot
layout(rbind(1:2))
plot(U,  xlab = expression(U[1]), ylab = expression(U[2]))
plot(U., xlab = expression(U[1]), ylab = expression(U[2]))


### 1.2 t copula ##

## Define the copula
nu <- 3 # degrees of freedom
th <- iTau(tCopula(df = nu), tau = 0.5) # correlation parameter
cop.t3 <- tCopula(param = th, df = nu) # t copula

## Sample
U.t <-  cCopula(U,  copula = cop.t3, inverse = TRUE)
U.t. <- cCopula(U., copula = cop.t3, inverse = TRUE)

## Plot
layout(rbind(1:2))
plot(U.t,  xlab = expression(U[1]), ylab = expression(U[2]))
plot(U.t., xlab = expression(U[1]), ylab = expression(U[2]))
layout(1) # reset layout

## 3d plot
U.3d. <- ghalton(n, d = 3)
cop <- tCopula(param = th, dim = 3, df = nu)
U.t.3d. <- cCopula(U.3d., copula = cop, inverse = TRUE)
cloud2(U.t.3d., xlab = expression(U[1]), ylab = expression(U[2]),
       zlab = expression(U[3])) # not much visible
pairs2(U.t.3d., labels.null.lab = "U")
## Structure doesn't *look* very convincing, but it is a low-discrepancy sequence


### 1.3 Clayton copula ##

## Define the copula
th <- iTau(claytonCopula(), tau = 0.5)
cop.C <- claytonCopula(th)

## Sample
U.C <-  cCopula(U,  copula = cop.C, inverse = TRUE)
U.C. <- cCopula(U., copula = cop.C, inverse = TRUE)

## Plot
layout(rbind(1:2))
plot(U.C,  xlab = expression(U[1]), ylab = expression(U[2]))
plot(U.C., xlab = expression(U[1]), ylab = expression(U[2]))
layout(1) # reset layout


### 2 Why considering quasi-random numbers? ##

## Computing P(U_1 > u_1,..., U_d > u_d) by MC based on PRNG, QRNG
tail_prob <- function(n, copula, u) # sample size, copula, lower-left endpoint
{
  d <- length(u)
  stopifnot(n >= 1, inherits(copula, "Copula"), 0 < u, u < 1,
            d == dim(copula))
  umat <- rep(u, each = n)

  ## Pseudo-random numbers
  clock <- proc.time() # start watch
  U <- rCopula(n, copula = copula) # unfair comparison but *still* slower than Sobol' + cCopula()
  prob.PRNG <- mean(rowSums(U > umat) == d) # = sum(<rows with all entries TRUE>) / n
  usr.time.PRNG <- 1000 * (proc.time() - clock)[[1]] # time in ms

  ## (Randomized) quasi-random numbers
  clock <- proc.time() # start watch
  U. <- cCopula(sobol(n, d = d, randomize = TRUE), # needed for unbiasedness
                copula = copula, inverse = TRUE)
  prob.QRNG <- mean(rowSums(U. > umat) == d) # = sum(<rows with all entries TRUE>) / n
  usr.time.QRNG <- 1000 * (proc.time() - clock)[[1]] # time in ms

  ## Return
  list(PRNG = c(prob = prob.PRNG, rt = usr.time.PRNG),
       QRNG = c(prob = prob.QRNG, rt = usr.time.QRNG))
}

## Simulate the probabilities of falling in (u_1, 1] x ... x (u_d, 1]
N <- 200 # number of replications
n <- 2000 # sample size
d <- 2 # dimension
u <- rep(0.99, d) # lower-left endpoint of the considered cube
nu <- 3 # degrees of freedom
tau <- 0.5 # Kendall's tau
rho <- iTau(tCopula(df = nu), tau = 0.5) # correlation parameter
cop <- tCopula(param = rho, dim = d, df = nu) # t copula

## Run
set.seed(271) # for reproducibility
system.time(res <- lapply(1:N, function(N.) tail_prob(n, copula = cop, u = u)))

## Grab out the results
prob.PRNG <- sapply(res, function(x) x[["PRNG"]][["prob"]])
prob.QRNG <- sapply(res, function(x) x[["QRNG"]][["prob"]])
rt.PRNG   <- sapply(res, function(x) x[["PRNG"]][["rt"]])
rt.QRNG   <- sapply(res, function(x) x[["QRNG"]][["rt"]])

## Boxplot of computed exceedance probabilities
boxplot(list(PRNG = prob.PRNG, QRNG = prob.QRNG),
        main = substitute("Simulated"~
                            P(bold(U) > bold(u))~~"for a"~t[nu.]~"copula",
                          list(nu. = nu)))
mtext(sprintf("N = %d replications with n = %d and d = %d", N, n, d),
      side = 4, line = 1, adj = 0, las = 0)
## QRNGs provide a smaller variance

## Variance reduction factor and % improvement
(vrf <- var(prob.PRNG)/var(prob.QRNG)) # estimated variance reduction factor

## Boxplot of measured run times (in milliseconds)
boxplot(list(PRNG = rt.PRNG, QRNG = rt.QRNG), outline = FALSE,
        main = substitute("Run times for a" ~ t[nu.]~"copula",
                          list(nu. = nu)))
mtext(sprintf("N = %d replications with n = %d and d = %d", N, n, d),
      side = 4, line = 1, adj = 0, las = 0)
## Clearly (the way we sampled here!), QRNG is slower here in this comparison

## Run-time reduction factor
(rrf <- mean(rt.PRNG)/mean(rt.QRNG)) # estimated run-time factor PRNG w.r.t. QRNG
## < 1 (PRNG faster) here because rCopula() is *way* faster than cCopula()

## Effort comparison (variance per unit of time)
vrf * rrf # = (Var('PRNG') * time('PRNG')) / (Var('QRNG') * time('QRNG'))
## Still better to use QRNG! (if > 1, QRNG is preferable)

## Remark:
## - QRNGs can estimate tail probabilities with smaller variance than PRNGs
##   (=> need less random variates to obtain the same precision => good for
##    memory/storage-intensive methods).
## - QRNGs can be faster than PRNGs (but it depends: Sobol': yes; generalized
##   Halton: no); this may depend on the dimension, too.
## - If more efficient sampling methods for copula-QRNGs are used (so similar
##   to rCopula() instead of cCopula()), run time for QRNGs is significantly
##   reduced; see Cambou, Hofert, Lemieux ("Quasi-random numbers for copula models").
## - Even if slower, it can still be advantageous to use QRNGs instead of PRNGs
##   (because of memory/storage limitations).

## Just comparing a PRNG and QRNGs
n. <- 2e7 # 20 Mio
set.seed(271)
system.time(runif(n.)) # PRNG
system.time(ghalton(n., d = 1)) # QRNG
system.time(sobol(n.,   d = 1, randomize = TRUE)) # Faster QRNG
## Run time also depends on the *type* of QRNG

if(TRUE) {
  ## Profiling: See where run time is spent
  Rprof(profiling <- tempfile(), line.profiling = TRUE) # enable profiling
  res <- lapply(1:N, function(N.) tail_prob(n, copula = cop, u = u))
  Rprof(NULL) # disable profiling
  (profile <- summaryRprof(profiling, lines = "both")) # get a summary
  ## => by.total => most of the run time is spent in cCopula()
}

#### Computing_volumes ####
## Computing C-volumes

library(copula)

## Volume we consider
a <- c(1/4, 1/2) # lower left end point
b <- c(1/3, 1) # upper right end point
stopifnot(0 <= a, a <= 1, 0 <= b, b <= 1, a <= b)


### 1 Independence copula ##

## Define the independence copula
ic <- indepCopula()

## Compute the Pi_{(a,b]} volume
vol <- prob(ic, l = a, u = b)

## Let's manually compute this volume. Under independence, the probability of falling
## in the two-dimensional interval (a, b] is (b[1] - a[1]) * (b[2] - a[2])
p <- (b[1] - a[1]) * (b[2] - a[2])

## Check whether they are (numerically) the same
stopifnot(all.equal(vol, p))

## d = 3
vol. <- prob(indepCopula(dim = 3), l = c(a, 1/8), u = c(b, 3/8))
p. <- p * (3/8 - 1/8)
stopifnot(all.equal(vol., p.))


### 2 Frechet--Hoeffding lower bound W ##

## Define W
W <- function(u1, u2) max(u1 + u2 - 1, 0)

## Compute the W_{(a,b]} volume
vol <- W(b[1], b[2]) - W(b[1], a[2]) - W(a[1], b[2]) + W(a[1], a[2])

## Let's manually compute this volume. Under W, the probability of falling
## in the two-dimensional interval (a, b] is given by the length of the line
## segment cut out of the secondary diagonal by (a, b]. Given how a and b
## are located here, this is:
p <- sqrt(2*(b[1] - a[1])^2) / sqrt(2); sqrt((b[1] - a[1])^2)

## Check whether they are (numerically) the same
stopifnot(all.equal(vol, p))

## Compute the volume of a rectangle below or above the secondary diagonal
## Rectangle below the secondary diagonal
a. <- a
b. <- c(1/3, 1/2)
stopifnot(0 <= a., a. <= 1, 0 <= b., b. <= 1, a.[1]+a.[2] <= 1, b. >= a.)
W(b.[1], b.[2]) - W(b.[1], a.[2]) - W(a.[1], b.[2]) + W(a.[1], a.[2])
## Rectangle above the secondary diagonal
a. <- c(2/3, 1/2)
b. <- c(3/4, 3/4)
stopifnot(0 <= a., a. <= 1, 0 <= b., b. <= 1, a.[1]+a.[2] >= 1, b. >= a.)
W(b.[1], b.[2]) - W(b.[1], a.[2]) - W(a.[1], b.[2]) + W(a.[1], a.[2])


### 3 t copula ##

## Define the t copula
nu <- 3 # note: pCopula() is currently only available for integer degrees of freedom
tc <- tCopula(iTau(tCopula(df = nu), tau = 0.5), df = nu)

## Compute the C_{(a,b]} volume
vol <- prob(tc, l = a, u = b)

## Let's manually compute this volume and check by Monte-Carlo
p <- pCopula(b, copula = tc) - pCopula(c(b[1], a[2]), copula = tc) -
  pCopula(c(a[1], b[2]), copula = tc) + pCopula(a, copula = tc)

## Check whether they are (numerically) the same
stopifnot(all.equal(vol, p))

## Let's also compute it by MC by approximating the volume by the relative
## frequency of samples falling in the two-dimensional interval (a, b]
N <- 1e6 # (large) sample size
set.seed(271) # set a seed (for reproducibility)
U <- rCopula(N, copula = tc) # sample the t copula
volMC <- mean(a[1] < U[,1] & U[,1] <= b[1] & a[2] < U[,2] & U[,2] <= b[2])

## Check whether we get (numerically) the same
stopifnot(all.equal(volMC, p, tol = 1e-2))

#### Spearmans_rho_Kendalls_tau ####
## Various experiments with Spearman's rho and Kendall's tau

library(copula)


### 1 How close do we get to the maximum log-likelihood for estimating P of a
###   Gaussian copula if Spearman's rho is used to estimate P? ##

## Generate pseudo-observations
set.seed(271)
nc <- normalCopula(0.6, dim = 3) # define a Gauss copula
U <- pobs(rCopula(1000, copula = nc)) # sample pseudo-observations from it

## Maximum pseudo-likelihood estimation of P
fit.N <- fitCopula(normalCopula(, dim = 3, dispstr = "un"), data = U)
fit.N@loglik # log-likelihood

## Estimation of P via Spearman's rho
P. <- cor(U, method = "spearman")
sum(dCopula(U, copula = normalCopula(P2p(P.), dim = 3, dispstr = "un"), log = TRUE)) # log-likelihood


### 2 How close do we get to the maximum log-likelihood for estimating P of a
###   t copula if Kendall's tau is used to estimate P? ##

## Generate pseudo-observations
set.seed(271)
tc <- tCopula(0.6, dim = 3, df = 10) # define a t_10 copula
U <- pobs(rCopula(1000, copula = tc)) # sample pseudo-observations from it

## Maximum pseudo-likelihood estimation of P and the d.o.f.
fit.t <- fitCopula(tCopula(, dim = 3, dispstr = "un"), data = U)
fit.t@loglik # log-likelihood

## Estimation of P via Kendall's tau (and of the d.o.f. via the likelihood)
fit.t. <- fitCopula(tCopula(, dim = 3, dispstr = "un"), data = U, method = "itau.mpl")
fit.t.@loglik # log-likelihood


### 3 How well does a Gauss copula fit t copula data? ##

fit.N <- fitCopula(normalCopula(, dim = 3, dispstr = "un"), data = U)
fit.N@loglik


### 4 Shall we estimate correlations with the sample correlation or by inverting
###   Kendall's tau? ##

## Generate data
df <- 3
tc <- tCopula(0.5, df = df) # define a t_3 copula
set.seed(271)
U <- lapply(1:3000, function(i) rCopula(90, copula = tc)) # (90, 2, 3000)-array
X <- lapply(U, function(x) qt(x, df = df)) # map to multivariate t

## Compute the correlation coefficient and Kendall's tau
cor.pearson <- sapply(X, function(x) cor(x)[2,1])
cor.tau <- iTau(tc, tau = sapply(X, function(x) cor(x, method = "kendall")[2,1]))

## Plot
ran <- range(cor.pearson, cor.tau)
plot(cor.pearson, ylim = ran, xlab = "Sample",
     ylab = expression("Estimates of"~rho))
points(cor.tau, col = "royalblue3")
legend("bottomleft", pch = c(1,1), bty = "n", col = c("black", "royalblue"),
       legend = c("Via sample correlation coefficient", "Via inverting Kendall's tau"))

## Variance
var(cor.pearson)
var(cor.tau) # smaller
var(sapply(U, function(x) cor(x)[2,1])) # equally small (= Spearman's rho; based on U)


### 5 Approximation of Spearman’s rho for t copulas by Pearson's rho ##

## Note: If it exists, Pearson's rho equals the correlation parameter of the
##       t copula

## Random angle for computing Spearman's rho for t copulas
Theta <- function(nu, n.MC = 1e5)
{
  W1 <- nu/rchisq(n.MC, df=nu)
  W2 <- nu/rchisq(n.MC, df=nu)
  W3 <- nu/rchisq(n.MC, df=nu)
  2 * asin(sqrt(W1 * W2 / ((W1 + W3) * (W2 + W3))))
}

## Spearman's rho for t copulas
## (see McNeil, Frey, Embrechts (2015, Prop. 7.44 and Eq. (7.42)))
rho_t <- function(rho, nu, n.MC = 1e5)
{
  stopifnot(-1 <= rho, rho <= 1, nu > 0, n.MC > 0)
  (6 / pi) * colMeans(asin(outer(sin(Theta(nu, n.MC = n.MC) / 2), rho, FUN = "*"))) # length(rho)
}

## Generate data
set.seed(271)
nu <- 10^seq(-1, 1, length.out = 3)
rho <- seq(-1, 1, by = 0.01)
rho.t <- sapply(nu, function(nu.) rho_t(rho, nu = nu.)) # (length(rho), length(nu)-matrix

## Plot of Spearman's rho vs Pearson's rho for t copulas
cols <- c("black", "maroon3", "darkorange2", "royalblue3")
plot(rho, rho.t[,1], type = "l", col = cols[2],
     xlab = expression("Pearson's"~~rho), ylab = expression("Spearman's"~~rho[S]))
lines(rho, rho.t[,2], col = cols[3])
lines(rho, rho.t[,3], col = cols[4])
abline(0, 1)
legend("bottomright", bty = "n", lty = rep(1, 4), col = cols,
       legend = c(expression(rho[S] == rho),
                  substitute(nu==nu., list(nu.=nu[1])),
                  substitute(nu==nu., list(nu.=nu[2])),
                  substitute(nu==nu., list(nu.=nu[3]))))