### R code for slides
### "Dependence modeling with copulas"
library(copula)
library(mnormt)               
library(lattice)               
library(mvtnorm)

## Motivating example
# Generate some data 

n <- 1000 # sample size
set.seed(271) # for reproducibility
U <- rCopula(n, copula = gumbelCopula(iTau(gumbelCopula(), tau = 0.5)))
X <- qnorm(U) # quantile transformation
X. <- qexp(U) # quantile transformation
plot(X,  xlab = expression(X[1]), ylab = expression(X[2]))
plot(X., xlab = expression(Y[1]), ylab = expression(Y[2]))

#build and plot pseudo-observations (scaled ranks)
U  <- pobs(X)
U. <- pobs(X.)
plot(U,  xlab = expression(U[1]), 
     ylab = expression(U[2]))
plot(U., xlab = expression(U[1]*"'"), 
     ylab = expression(U[2]*"'"))

layout(matrix(1:4, ncol=2))
par(mar=c(3,4,2,2))
plot(U[,1],  ylab = expression(U[1]), xlab="")
plot(U[,2],  ylab = expression(U[2]), xlab="")
plot(U.[,1], ylab = expression(U[1]*"'"), xlab="")
plot(U.[,2], ylab = expression(U[2]*"'"), xlab="")

###################################################
## Basic copulas
## Independence Copula

d<-2
ic<-indepCopula(dim=d)
par(mar=c(1,1,1,1))
wireframe2(ic, FUN=pCopula,
           col.4=adjustcolor("black", alpha.f=0.25), col=3)
par(mar=c(3,3,3,3))
contourplot2(ic, FUN=pCopula, col=3)

# c-volume 
n<-1000
U<-rCopula(n, copula=ic)
par(mar=c(4,4,1,1))
plot(U, xlab=quote(U[1]), ylab=quote(U[2]))
segments(1/4, 1/2, 1/4, 1, lwd = 3, col=2)
segments(1/3, 1/2, 1/3, 1, lwd = 3, col=2)
segments(1/4, 1/2, 1/3, 1/2, lwd = 3, col=2)
segments(1/4, 1, 1/3, 1, lwd = 3, col=2)

## Fréchet–Hoeffding Bounds

U<-runif(100)
plot(cbind(U, 1-U), 
     xlab=quote(U[1]), ylab=quote(U[2]))
plot(cbind(U, U), 
     xlab=quote(U[1]), ylab=quote(U[2]))
u<-seq(0,1, length.out=40)
u12<-expand.grid("u[1]"=u, "u[2]"=u)
W<-pmax(u12[,1]+u12[,2]-1, 0)
M<-pmin(u12[,1], u12[,2])
val.W<-cbind(u12, "W(u[1], u[2])"= W)
val.M<-cbind(u12, "M(u[1], u[2])"= M)
wireframe2(val.W)
wireframe2(val.M)


## Marshall-Olkin copula
C<-function(u, alpha)
  pmin(u[,1]*u[,2]^(1-alpha[2]), 
       u[,1]^(1-alpha[1])*u[,2])
alpha<-c(0.2, 0.8)
u<-seq(0,1, length.out=40)
u12<-expand.grid("u[1]"=u, "u[2]"=u)
val<-cbind(u12, "C(u[1], u[2])"=C(u12, alpha=alpha))
set.seed(34)
V<-matrix(runif(1000*3), ncol=3)
U<-cbind(pmax(V[,1]^(1/(1-alpha[1])), V[,3]^(1/alpha[1])), 
         pmax(V[,2]^(1/(1-alpha[2])), V[,3]^(1/alpha[2])))
wireframe2(val)
plot(U, xlab=quote(U[1]), ylab=quote(U[2]))

## Invariance property: Example
#######################################

# From bivariate normal to normal copula
set.seed(332)
d<-2
rho<-0.7
P<-matrix(rho, nrow=d, ncol=d)
diag(P)<-1
u<-runif(d)
x<-qnorm(u)
nc<-normalCopula(rho)

set.seed(332)
par(mfrow=c(1,2))
X<-rmvnorm(1000, mean=c(0,0), sigma = P)
plot(X, xlim=c(-4, 4), ylim=c(-4, 4), 
     xlab=quote(X[1]), ylab=quote(X[2]))
U<-pnorm(X)
#same as
set.seed(332)
U.<-rCopula(1000, normalCopula(rho,dim=2))
stopifnot(all.equal(U, U.))
plot(U, xlim=c(0, 1), ylim=c(0, 1), 
     xlab=quote(U[1]), ylab=quote(U[2]))

#From normal copula to meta-Gaussian sample with exp. margins
par(mfrow=c(1,2))
plot(U, xlim=c(0, 1), ylim=c(0, 1), 
     xlab=quote(U[1]), ylab=quote(U[2]))
Y<-qexp(U, 2)
plot(Y, xlim=c(0, 4), ylim=c(0, 4), 
     xlab=quote(Y[1]), ylab=quote(Y[2]))


## Gaussian copula

r<-0.5
mu <- c(0,0)
sigma <- matrix(c(1,r,r,1),2,2)                                
x<-seq(-5,5,0.2)
y<-seq(-5,5,0.2)
f<-function(x,y){dmnorm(cbind(x,y), mu, sigma)}
z<-outer(x,y,f)
persp(x,y,z, theta=30, axes=T, ticktype="detailed",
      cex.axis=0.8)

nc <- normalCopula(r, dim=2)
## Copula (wire frame and level curves)
wireframe2(nc, FUN = pCopula, draw.4.pCoplines=FALSE)
# contourplot2(nc, FUN = pCopula)
## Copula density (wire frame and level curves)
wireframe2(nc, FUN = dCopula, delta = 0.002)

## scatter plots of biv. simulated data form the normal copula 
set.seed(999)
n<-1000
nc
U2<-rCopula(n, copula=setTheta(nc, value=0.1))
plot(U2, xlab=quote(U[1]), ylab=quote(U[2]))
U<-rCopula(n, copula = nc)
plot(U, xlab=quote(U[1]), ylab=quote(U[2]))
U3<-rCopula(n, copula=setTheta(nc, value=0.7))
plot(U3, xlab=quote(U[1]), ylab=quote(U[2]))

#Student t copulas

nu <- 4 
tc <- tCopula(iTau(tCopula(df = nu), tau = 0.6), df = nu)
tc
set.seed(271)
U <- rCopula(1000, copula = tc) # sample from the t copula
wireframe2(tc, FUN = dCopula, delta = 0.025) # density
contourplot2(tc, FUN = pCopula) # copula
contourplot2(tc, FUN = dCopula, n.grid = 42, cuts = 27) # density
plot(U, xlab = quote(U[1]), ylab = quote(U[2])) # scatter plot

### Comparison between copulas

library(gridExtra) # for grid.arrange()
n <- 1000 # sample size
set.seed(27)
# Define the copula parameters chosen such 
# that Pearson's correlation coefficient,
# when computed with N(0,1) margins, is 0.7 
# tc <- frankCopula(iRho(frankCopula(), rho = 0.7))
th.f<- 5.821 # Frank 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
nf <- frankCopula(th.f) # Frank copula
gc <- gumbelCopula(th.g) # Gumbel copula
cc <- claytonCopula(th.c) # Clayton copula
tc <- tCopula(th.t) # t_4 copula

## Generate copula data
U.nf <- rCopula(n, copula = nf)
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.nf <- qnorm(U.nf)
X.gc <- qnorm(U.gc)
X.cc <- qnorm(U.cc)
X.tc <- qnorm(U.tc)

par(mfrow=c(2,2))
plot(U.nf, xlab = expression(U[1]), ylab = expression(U[2]), # Frank copula
     cex = 0.4, main = "Frank 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
par(mfrow=c(2,2))
m <- max(abs(X.nf), abs(X.gc), abs(X.cc), abs(X.tc))
lim <- c(-m, m)
plot(X.nf, xlim = lim, ylim = lim, # meta-Frank
     xlab = expression(X[1]), ylab = expression(X[2]),
     cex = 0.4, main = "Meta-Frank sample")
plot(X.gc, xlim = lim, ylim = lim,
     xlab = expression(X[1]), ylab = expression(X[2]), # meta-Gumbel
     cex = 0.4, main = "Meta-Gumbel sample")
plot(X.cc, xlim = lim, ylim = lim, # meta-Clayton
     xlab = expression(X[1]), ylab = expression(X[2]),
     cex = 0.4, main = "Meta-Clayton sample")
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")))
par(opar) # restore graphical parameters

## Density plots ############################################################
## Define the corresponding (meta-C) densities 
## 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)) )
dMetaCop<- function(x, copula)
  exp(dCopula(pnorm(x), copula = copula, log = TRUE) + rowSums(dnorm(x, log = TRUE)))

par(mfrow=c(2,2))
## Wire frame plots
n.grid <- 25 # number of grid points
lim <- c(-2.6, 2.6)
s <- seq(lim[1], lim[2], length.out = n.grid)
grid <- as.matrix(expand.grid("x[1]" = s, "x[2]" = s))

val.f <- cbind(grid, "f(x[1],x[2])" = dMetaCop(grid, copula = nf))
val.g <- cbind(grid, "f(x[1],x[2])" = dMetaCop(grid, copula = gc))
val.c <- cbind(grid, "f(x[1],x[2])" = dMetaCop(grid, copula = cc))
val.t <- cbind(grid, "f(x[1],x[2])" = dMetaCop(grid, copula = tc))

zlim <- c(0, max(val.f[,3], val.g[,3], val.c[,3], val.t[,3]))

## Density plots
w.f <- wireframe2(val.f, xlim = lim, ylim = lim, zlim = zlim, 
            main = list(label="Meta-Frank density - N(0,1) margins", 
                   cex=0.8), scales=list(arrows = FALSE, cex = 0.6)) 
w.g <- wireframe2(val.g, xlim = lim, ylim = lim, zlim = zlim, 
                  main = list(label="Meta-Gumbel density - N(0,1) margins", 
                              cex=0.8), 
                  scales=list(arrows = FALSE, cex = 0.6))
w.c <- wireframe2(val.c, xlim = lim, ylim = lim, zlim = zlim, 
          main = list(label="Meta-Clayton density - N(0,1) margins", cex=0.8),
          scales=list(arrows = FALSE, cex = 0.6))
w.t <- wireframe2(val.t, xlim = lim, ylim = lim, zlim = zlim, 
        main = list(label= expression(bold("Meta-"*italic(t)[4]~"density - N(0,1) margins")), 
                    cex=0.8),scales=list(arrows = FALSE, cex = 0.6))
grid.arrange(w.f, w.g, w.c, w.t, ncol = 2)

## Contour plots
lay <- matrix(1:4, ncol = 2, byrow = TRUE) # layout matrix
layout(lay) # layout
c.f <- contourplot2(val.f, xlim = lim, ylim = lim,
        main = list(label = "Meta-Frank contours - N(0,1) margins", cex = 0.8))
c.g <- contourplot2(val.g, xlim = lim, ylim = lim,
        main = list(label = "Meta-Gumbel contours - N(0,1) margins", cex =0.8))
c.c <- contourplot2(val.c, xlim = lim, ylim = lim,
        main = list(label = "Meta-Clayton contours - N(0,1) margins", cex = 0.8))
c.t <- contourplot2(val.t, xlim = lim, ylim = lim,
      main = list(label = expression(bold("Meta-"*italic(t)[4]~"contours - N(0,1) margins")),
      cex = 0.8))
grid.arrange(c.f, c.g, c.c, c.t, ncol = 2)



## 3-dim t copula object
nu <- 5 # degrees of freedom
th <- iTau(tCopula(df = nu), tau = 0.5) # correlation parameter
tc. <- tCopula(th, dim = 3, df = nu)
U <- rCopula(n, copula = tc.)
cloud2(U, xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))
pairs2(U, cex = 0.4) 

### Sampling from survival copulas
## survival copula 
cop<-claytonCopula(2)
set.seed(332)
U<-rCopula(1000, copula=cop)
V<-1-U # sample from the survival Clayton copula
par(mfrow=c(1,2))
plot(U, xlab = quote(U[1]),  ylab = quote(U[2]))
plot(V, xlab = quote(V[1]),  ylab = quote(V[2]))


### Symmetry and Exchangeability
par(pty = "s")
plot(rCopula(1000, copula=tCopula(0.7, df=3.5)),
     xlab = quote(U[1]),  ylab = quote(U[2]))
contourplot2(tCopula(0.7, df=3.5),
             FUN=dCopula, n.grid=64, lwd=1/2)
plot(rCopula(1000, copula=gumbelCopula(2)),
     xlab = quote(U[1]),  ylab = quote(U[2]))
contourplot2(gumbelCopula(2),
             FUN=dCopula, n.grid=64, lwd=1/4, 
             pretty=FALSE, cuts=42, 
             col.regions=gray(seq(0.5,1, length.out=128)))

#########################################################
### (II) Dependence concepts and measures
### Two models with N(0; 1) margins and zero correlation.
n<-1000
set.seed(314)
Z<-rnorm(n)
U<-runif(n)
V<-rep(1,n)
V[U<1/2]<--1
Y<-cbind(Z, Z*V)
X<-matrix(rnorm(n*2), ncol=2)

par(mfrow=c(1,2),mar=c(4,4,2,2))
plot(X, xlab=quote(X[1]), ylab=quote(X[2])) # indep copula
plot(Y, xlab=quote(Y[1]), ylab=quote(Y[2])) # mixture

#### VaR for aggregated risks (corr. fallacies)
layout(1:1)
q<-seq(0.9, 0.99, by=0.001)
plot(q, sqrt(2)*qnorm(q), type="l", 
     ylim=c(1.5, 4), lty=2, 
     ylab="Value-at-Risk")
curve(2*qnorm(2*x-1), add=TRUE)
legend("topleft", legend=c(expression(X[1]+X[2]), 
                           expression(Y[1]+Y[2])), 
       lty=c(2,1), bty="n")


## Define the Function to compute the correlation bounds for LN(0, sigma_.^2) margins

cor_bound <- 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(grid, method = "min"))
val.max <- cbind(grid, "bar(Cor)(sigma[1],sigma[2])" =
                   cor_bound(grid))
par(mar=c(4,4,3,3))
plot(s,cor_bound(cbind(rep(1, n.grid), s), method="min"), 
     ylim=c(-1,1),xlim=c(0,5), type="l", 
     xlab=quote(sigma), col="coral", lwd=2, ylab="correlation")
points(s, cor_bound(cbind(rep(1, n.grid), s), method="max"),
       type="l", col="blue", lwd=2)
abline(h=0)
legend("topright", legend=c("minimum corr.", "maximum corr."), 
       lwd=c(2,2), bty="n", col=c("coral", "blue"), cex=0.7)


### rank correlation: Gaussian copula
rho<-seq(-1, 1, by= 0.01)
rho.s<-(6/pi)*asin(rho/2)
tau<-(2/pi)*asin(rho)
plot(rho, rho.s, 
     type="l", col=2, xlab=quote(rho), 
     lwd=2, ylab="rank correlation")
abline(0,1)
lines(rho, tau, col=4, lwd=2, lty=2)
legend("topleft", bty="n", col=c(2,4),
       lwd=c(2,2), lty=c(1,2), 
       legend=c(expression(rho[s]),expression(tau)))

## Tail dependence

par(mfrow=c(2,2), mar=c(4,4,3,3))
tau<-0.6
th.n<-iTau(normalCopula(), tau=tau)
th.t<-iTau(tCopula(df=3), tau=tau)
th.c<-iTau(claytonCopula(), tau=tau)
th.g<-iTau(gumbelCopula(), tau=tau)
n<-10000
marg<-list(list(mean=0, sd=1),list(mean=0, sd=1))
X.n<-rMvdc(n, mvdc=mvdc(normalCopula(th.n), c("norm", "norm"), marg))
X.t<-rMvdc(n, mvdc=mvdc(tCopula(th.t), c("norm", "norm"), marg))
X.c<-rMvdc(n, mvdc=mvdc(claytonCopula(th.c), c("norm", "norm"), marg))
X.g<-rMvdc(n, mvdc=mvdc(gumbelCopula(th.g), c("norm", "norm"), marg))
a<-0.005
q<-qnorm(c(a, 1-a))
lim<-c(floor(range(q, X.n, X.g, X.c, X.t)[1]), floor(range(q, X.n, X.g, X.c, X.t)[2]))
plot(X.n, xlim=lim, ylim=lim, cex=0.4, 
     xlab=quote(X[1]), ylab=quote(X[2]), 
     main="Normal copula") 
abline(h=q, v=q, lty=2, col="gray")
plot(X.t, xlim=lim, ylim=lim, cex=0.4, 
     xlab=quote(X[1]), ylab=quote(X[2]), 
     main="t copula (3 df)") 
abline(h=q, v=q, lty=2, col="gray")
plot(X.c, xlim=lim, ylim=lim, cex=0.4, 
     xlab=quote(X[1]), ylab=quote(X[2]), 
     main="Clayton copula") 
abline(h=q, v=q, lty=2, col="gray")
plot(X.g, xlim=lim, ylim=lim, cex=0.4, 
     xlab=quote(X[1]), ylab=quote(X[2]), 
     main="Gumbel-Hougaard copula") 
abline(h=q, v=q, lty=2, col="gray")

## (Tail Dependence of t Copulas)
pdf("tcopTD.pdf", width = 10, height = 4)
par(mfrow=c(1,2), mar=c(4,4,3,3))
rho<-seq(-1, 1, by=0.01)
nu<-c(3, 5, 7, Inf)
n.nu<-length(nu)
lam.rho<-sapply(nu, function(nu.)
  sapply(rho, function(rho.) lambda(tCopula(rho., df=nu.))[["lower"]]))
expr.rho<-as.expression(lapply(1:n.nu, function(j)
  bquote(nu==.(if(nu[j]==Inf) quote (infinity) else nu[j]))))
matplot(rho, lam.rho, type="l", lty=1, lwd=2, col=1:n.nu,
        xlab=quote(rho), ylab=quote(lambda))
legend("topleft", legend=expr.rho, bty="n", lwd=2, col=1:n.nu)

rho.<-c(-1, -0.5, 0, 0.5, 1)
nu.<-c(seq(3, 12, by=0.2), Inf)
n.rho<-length(rho.)
lam.nu<-sapply(rho., function(rh)
  sapply(nu., function(nu) lambda(tCopula(rh, df=nu))[["lower"]]))
expr<-as.expression(lapply(1:n.rho, function(j)
  bquote(rho==.(rho.[j]))))
matplot(nu., lam.nu, type="l", lty=1, lwd=2, col=1:n.rho,
        xlab=quote(nu), ylab=quote(lambda))
legend("topright", legend=expr, 
       bty="n", lwd=2, col=1:n.rho)
dev.off()