scripts/Copula_examples.R

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

outdir = "C:/Users/sch/Dropbox/Egne dokumenter/Skole/master/opgave/Figures/Copula/"
library(copula)

## Generate some data
n <- 1000
set.seed(100)
U <- rCopula(n, copula = gumbelCopula(iTau(gumbelCopula(), tau = 0.5)))
X <- qnorm(U)
X. <- qexp(U)

## Plots
pdf(file = paste0(outdir,"mot_norm_exp.pdf"),
    width = 8, height = 6)
layout(rbind(1:2))
plot(X,  xlab = expression(X[1]),     ylab = expression(X[2]))
plot(X., xlab = expression(X[1]*"'"), ylab = expression(X[2]*"'"))
dev.off()

## Same copula for N(0,1) and Exp(1) margins
pdf(file = paste0(outdir,"mot_norm_exp_cop.pdf"),
    width = 8, height = 6)
layout(rbind(1:2))
plot(U,  xlab = expression(U[1]),     ylab = expression(U[2]))
plot(U., xlab = expression(U[1]*"'"), ylab = expression(U[2]*"'"))
dev.off()

pdf(file = paste0(outdir,"mot_norm_exp_ind.pdf"),
    width = 8, height = 6)
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)
dev.off()

#### Sklar_theorem_and_invariance_principle ####
library(mvtnorm)
library(copula)

set.seed(271)
## Sample from a t copula
n <- 1000
d <- 2
rho <- 0.8
P <- matrix(rho, nrow = d, ncol = d)
diag(P) <- 1
nu <- 4
set.seed(271)
X <<- rmvt(n, sigma = P, df = nu)
U <- pt(X, df = nu)
Y <- qexp(U)

## Plot setup
ind <- c(A = 50, B = 322, C = 980)
cols <- c("royalblue3", "maroon3", "darkorange2")

## Scatter plot of a bivariate t distribution
pdf(file = paste0(outdir, "sklar_t_exp.pdf"),
    height = 6, width = 8)
layout(rbind(1:2))
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 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])
layout(1)
dev.off()

U  <- pobs(X)
U. <- pobs(Y)

pdf(file = paste0(outdir, "sklar_t_exp_cop.pdf"),
    height = 6, width = 8)
layout(rbind(1:2))
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])

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])
layout(1)
dev.off()

pdf(file = paste0(outdir,"sklar_t_exp_ind.pdf"),
    width = 8, height = 6)
layout(matrix(1:4, ncol = 2))
plot(U[,1],  ylab = expression(U[1]))
points(which(U[ind["A"],1] == U[,1]), U[ind["A"],1], pch = 21, bg = cols[1], col = cols[1])
points(which(U[ind["B"],1] == U[,1]), U[ind["B"],1], pch = 21, bg = cols[2], col = cols[2])
points(which(U[ind["C"],1] == U[,1]), U[ind["C"],1], pch = 21, bg = cols[3], col = cols[3])
text(which(U[ind["A"],1] == U[,1]), U[ind["A"],1], labels = "A", adj = c(-0.5, 1.75), col = cols[1])
text(which(U[ind["B"],1] == U[,1]), U[ind["B"],1], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(which(U[ind["C"],1] == U[,1]), U[ind["C"],1], labels = "C", adj = c(0.5, -0.75), col = cols[3])

plot(U[,2],  ylab = expression(U[2]))
points(which(U[ind["A"],2] == U[,2]), U[ind["A"],1], pch = 21, bg = cols[1], col = cols[1])
points(which(U[ind["B"],2] == U[,2]), U[ind["B"],1], pch = 21, bg = cols[2], col = cols[2])
points(which(U[ind["C"],2] == U[,2]), U[ind["C"],1], pch = 21, bg = cols[3], col = cols[3])
text(which(U[ind["A"],2] == U[,2]), U[ind["A"],1], labels = "A", adj = c(-0.5, 1.75), col = cols[1])
text(which(U[ind["B"],2] == U[,2]), U[ind["B"],1], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(which(U[ind["C"],2] == U[,2]), U[ind["C"],1], labels = "C", adj = c(0.5, -0.75), col = cols[3])

plot(U.[,1],  ylab = expression(U[1]))
points(which(U.[ind["A"],1] == U.[,1]), U.[ind["A"],1], pch = 21, bg = cols[1], col = cols[1])
points(which(U.[ind["B"],1] == U.[,1]), U.[ind["B"],1], pch = 21, bg = cols[2], col = cols[2])
points(which(U.[ind["C"],1] == U.[,1]), U.[ind["C"],1], pch = 21, bg = cols[3], col = cols[3])
text(which(U.[ind["A"],1] == U.[,1]), U.[ind["A"],1], labels = "A", adj = c(-0.5, 1.75), col = cols[1])
text(which(U.[ind["B"],1] == U.[,1]), U.[ind["B"],1], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(which(U.[ind["C"],1] == U.[,1]), U.[ind["C"],1], labels = "C", adj = c(0.5, -0.75), col = cols[3])

plot(U.[,2],  ylab = expression(U[2]))
points(which(U.[ind["A"],2] == U.[,2]), U.[ind["A"],1], pch = 21, bg = cols[1], col = cols[1])
points(which(U.[ind["B"],2] == U.[,2]), U.[ind["B"],1], pch = 21, bg = cols[2], col = cols[2])
points(which(U.[ind["C"],2] == U.[,2]), U.[ind["C"],1], pch = 21, bg = cols[3], col = cols[3])
text(which(U.[ind["A"],2] == U.[,2]), U.[ind["A"],1], labels = "A", adj = c(-0.5, 1.75), col = cols[1])
text(which(U.[ind["B"],2] == U.[,2]), U.[ind["B"],1], labels = "B", adj = c(0.5, -0.75), col = cols[2])
text(which(U.[ind["C"],2] == U.[,2]), U.[ind["C"],1], labels = "C", adj = c(0.5, -0.75), col = cols[3])
layout(1)
dev.off()

#### Examples_of_copulas ####
library(copula)

n <- 1000
d <- 5
tau <- 0.6
tblack <- function(alpha.f) adjustcolor("black", alpha.f = alpha.f)

set.seed(271)
### Fundamental copulas ##
### Independence copula ##
ic <- indepCopula()

pdf(file = paste0(outdir, "inde_cop.pdf"))
wireframe2(ic, FUN = pCopula)
dev.off()

pdf(file = paste0(outdir, "inde_cop_den.pdf"))
wireframe2(ic, FUN = dCopula, delta = 0.001, zlim = 0:1)
dev.off()

### Frechet--Hoeffding bounds W and M ##
n.grid <- 26
u <- seq(0, 1, length.out = n.grid)
grid <- expand.grid("u[1]" = u, "u[2]" = u)
W <- pmax(grid[,1] + grid[,2] - 1, 0)
M <- pmin(grid[,1], grid[,2])
val.W <- cbind(grid, "W(u[1],u[2])" = W)
val.M <- cbind(grid, "M(u[1],u[2])" = M)

pdf(file = paste0(outdir, "M.pdf"))
wireframe2(val.M, main = "M")
dev.off()

pdf(file = paste0(outdir, "W.pdf"))
wireframe2(val.W, main = "W")
dev.off()

### Implicit copulas ##
### Normal (or: Gauss) copula ##
th <- iTau(normalCopula(), tau = 0.6)
nc <- normalCopula(th)

## Copula (
pdf(file = paste0(outdir, "norm_cop.pdf"))
wireframe2(nc, FUN = pCopula)
dev.off()

## Copula density
pdf(file = paste0(outdir, "norm_cop_den.pdf"))
wireframe2(nc, FUN = dCopula, delta = 0.02)
dev.off()

### t copula ##
nu <- 3
th <- iTau(tCopula(, df = nu), tau = tau)
tc <- tCopula(th, df = nu)

pdf(file = paste0(outdir, "t_cop.pdf"))
wireframe2(tc, FUN = pCopula)
dev.off()

pdf(file = paste0(outdir, "t_cop_den.pdf"))
wireframe2(tc, FUN = dCopula, delta = 0.02)
dev.off()

### Explicit copulas ##

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

pdf(file = paste0(outdir, "cc_cop_den.pdf"))
wireframe2(cc, FUN = dCopula, delta = 0.02)
dev.off()
pdf(file = paste0(outdir, "cc_cop_con.pdf"))
contourplot2(cc, FUN = dCopula)
dev.off()

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

pdf(file = paste0(outdir, "gc_cop_den.pdf"))
wireframe2(gc, FUN = dCopula, delta = 0.02)
dev.off()
pdf(file = paste0(outdir, "gc_cop_con.pdf"))
contourplot2(gc, FUN = dCopula)
dev.off()

#### Exchangeability ####
tau = 0.6
nu <- 3
th <- iTau(tCopula(, df = nu), tau = tau)

tc. <- tCopula(th, dim = 5, df = nu)
U <- rCopula(n, copula = tc.)

pdf(file = paste0(outdir, "exhange_true.pdf"))
pairs2(U, cex = 0.4, col = tblack(0.5))
dev.off()

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

tc.. <- tCopula(P2p(P), dim = d, dispstr = "un", df = 3.5)
U. <- rCopula(n, copula = tc..)
pdf(file = paste0(outdir, "exhange_false.pdf"))
pairs2(U., cex = 0.4, col = tblack(0.5))
dev.off()

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

library(copula)
n <- 1000
set.seed(271)
Z <- rnorm(n)
U <- runif(n)
V <- rep(1, n)
V[U < 1/2] <- -1
X <- cbind(Z, Z*V)
cor(X)[2,1]
Y <- matrix(rnorm(n * 2), ncol = 2)
cor(Y)[2,1]

## Plots
pdf(file = paste0(outdir, "cop_fal_1.pdf"),
    width = 8, height = 6)
layout(matrix(1:2, ncol = 2))
plot(X, xlab = expression(X), ylab = expression(Y))
plot(Y, xlab = expression(phi[1]), ylab = expression(phi[2]))
layout(1)
dev.off()

cor_bound_B <- function(p, method = c("max", "min")) {
  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]))
}

n.grid <- 26
p <- seq(0.01, 0.99, length.out = n.grid)
grid <- expand.grid("p[1]" = p, "p[2]" = p)
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))
pdf(file = paste0(outdir, "cop_fal_2_lower.pdf"))
wireframe2(val.min, main = expression(rho[W]))
dev.off()
pdf(file = paste0(outdir, "cop_fal_2_upper.pdf"))
wireframe2(val.max, main = expression(rho[M]))
dev.off()

#### Spearmans_rho and Kendalls_tau ####
library(copula)
mple_com_ikte_gaus <- function(rho){
  out <- matrix(ncol = 2, nrow = length(rho))
  for (i in 1:length(rho)){
    set.seed(271)
    nc <- normalCopula(rho[i], dim = 3)
    U <- pobs(rCopula(1000, copula = nc))
    fit.N <- fitCopula(normalCopula(, dim = 3, dispstr = "un"), data = U)
    mple_rho <- fit.N@loglik
    P. <- cor(U, method = "kendall")
    ikte_rho <- sum(dCopula(U, copula = normalCopula(P2p(P.), dim = 3, dispstr = "un"), log = TRUE)) # log-likelihood
    out[i,] <- c(mple_rho, ikte_rho)
  }
  colnames(out) <- c("MPLE", "IKTE")
  rownames(out) <- rho
  return(t(out))
}
mple_Compare_ikte_gaus_table <- mple_com_ikte_gaus(c(0.3, 0.4,0.5,0.6, 0.7))
xtable::xtable(mple_Compare_ikte_gaus_table)

mple_com_ikte_t <- function(rho) {
  out <- matrix(ncol = 2, nrow = length(rho))
  for (i in 1:length(rho)){
    set.seed(271)
    tc <- tCopula(rho[i], dim = 3, df = 7)
    U <- pobs(rCopula(1000, copula = tc))
    fit.t <- fitCopula(tCopula(, dim = 3, dispstr = "un"), data = U)
    mple_rho <- fit.t@loglik
    fit.t. <- fitCopula(tCopula(, dim = 3, dispstr = "un"), data = U, method = "itau.mpl")
    ikte <- fit.t.@loglik
    out[i,] <- c(mple_rho, ikte)
  }
  colnames(out) <- c("MPLE", "IKTE")
  rownames(out) <- rho
  return(t(out))
}
mple_Compare_ikte_t_table <- mple_com_ikte_t(c(0.3, 0.4,0.5,0.6, 0.7))
xtable::xtable(mple_Compare_ikte_t_table)

df <- 4
tc <- tCopula(0.8, df = df)
set.seed(271)
U <- lapply(1:3000, function(i) rCopula(100, copula = tc))
X <- lapply(U, function(x) qt(x, df = df))

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]))

ran <- range(cor.pearson, cor.tau)
pdf(file = paste0(outdir, "rho_pear_tau.pdf"),
    height = 6, width = 8)
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(expression(hat(rho)~"by sample correlation coefficient"),
                  expression(hat(rho)~"by inverting Kendall's tau")))
dev.off()

var(cor.pearson)
var(cor.tau)
var(sapply(U, function(x) cor(x)[2,1]))

xtable::xtable(tibble::as_tibble(t(c(rho = var(cor.pearson),
                            tau = var(cor.tau),
                            S = var(sapply(U, function(x) cor(x)[2,1]))))))

var = round(data.frame(t(c(rho = var(cor.pearson),
              tau = var(cor.tau),
              S = var(sapply(U, function(x) cor(x)[2,1]))))), digits = 4)
xtable::xtable(var, digits = 4)