# Copula generator
library(rmsfuns)
suppressPackageStartupMessages(pacman::p_load("copula", "VineCopula", 
    "magrittr"))
suppressPackageStartupMessages(pacman::p_load("scatterplot3d", 
    "rgl", "ggplot2", "grid"))
set.seed(235)

Intro

In this practical, I will show you how to fit copula models using R.

The practical is not fully completed - but I include it here nonetheless for those interested in seeing how we can apply these techniques.

As this is a completely bonus topic that we will not cover in class, I will not spend too much effort in detailing the steps below.

Probability Distribution function

df_sample <- data.frame(x = c(x = rnorm(10000, mean = 0, sd = 3), 
    y = rnorm(10000, mean = 0, sd = 10), z = rnorm(10000, mean = 0, 
        sd = 2), q = rnorm(10000, mean = 1, sd = 1.5)), g = gl(4, 
    10000))
head(df_sample)
##             x g
## x1  6.4134105 1
## x2  0.8164904 1
## x3 -2.0243966 1
## x4  5.3472214 1
## x5 -0.3749610 1
## x6 -4.3162014 1
ggplot(df_sample, aes(x, colour = g)) + geom_density(adjust = 5)

Cumulative Distribution function

ggplot(df_sample, aes(x, colour = g)) + stat_ecdf()

Copulas

Simulating Copulas

  • Elliptical and Archimedean copulas have nice mathematical proprieties, shape and formulas.
  • They are a good choice for the initial warm up.
  • Assuming you already know the parameters, this is how you would generate a bivariate normal and a t copula
# Generate a bivariate normal copula with rho = 0.7
normal <- tCopula(par = 0.8, dim = 2)
# Generate a bivariate t-copula with rho = 0.8 and df = 2
stc <- tCopula(par = 0.8, dim = 2, df = 1)

# Generate random samples
norm_d <- rCopula(2000, normal)
t_d <- rCopula(2000, stc)

Simulating (cont.)

  • Great simplyfication when trying to model the joint behaviour complex phenomena, mainly because the marginals are easier to model than the joint distribution. Furthermore, if you had to model the joint behaviour of, say more than 20 objects (for instance the life of 20+ capacitors in an electrical component), then you would probably have a hard time by going straight into modelling the joint behaviour while the marginal distributions may be easier to estimate.
# Plot the samples
p1 <- qplot(norm_d[, 1], norm_d[, 2], colour = norm_d[, 1], main = "Norm with rho = 0.8", 
    xlab = "u", ylab = "v")

p1

p2 <- qplot(t_d[, 1], t_d[, 2], colour = t_d[, 1], main = "Student t with rho = 0.8 and df = 1", 
    xlab = "u", ylab = "v")

p2

  • Once the copula has been fitted, you can easily generate random numbers by using rCopula method for the single copula
# Generate random samples
norm_d <- rCopula(2000, normal)
t_d <- rCopula(2000, stc)
  • It also helps to generate the contour plots to see the how your data looks
pacman::p_load(kdecopula)
BiCopKDE(pobs(norm_d)[, 1], pobs(norm_d)[, 2])

BiCopKDE(pobs(t_d)[, 1], pobs(t_d)[, 2])

  • The general form “name” + “Copula()” can be used to build Archimedean copulas as well. Archimedean copulas take only a single parameter \(\theta\)

Copulas with the VineCopula package

# Build a Frank, a Gumbel and a Clayton copula

frank <- BiCop(5, tau = 0.5) %>% BiCopSim(10000, .)
gumbel <- BiCop(4, tau = 0.5) %>% BiCopSim(10000, .)
clayton <- BiCop(3, tau = 0.5) %>% BiCopSim(10000, .)
  • Print information on the Frank copula
print(BiCop(5, tau = 0.5))
## Bivariate copula: Frank (par = 5.82, tau = 0.5)
  • Plot the samples
p1 <- qplot(frank[, 1], frank[, 2], colour = frank[, 1], main = "Frank copula random samples", 
    xlab = "u", ylab = "v")
p2 <- qplot(gumbel[, 1], gumbel[, 2], colour = gumbel[, 1], main = "Gumbel copula random samples", 
    xlab = "u", ylab = "v")
p3 <- qplot(clayton[, 1], clayton[, 2], colour = clayton[, 1], 
    main = "Clayton copula random samples", xlab = "u", ylab = "v")
  • Define grid layout to locate plots and print each graph^(1)
pushViewport(viewport(layout = grid.layout(1, 3)))
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(p3, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))

  • Contours plots
contour(BiCop(5, tau = 0.5))

BiCopKDE(frank[, 1], frank[, 2])

contour(BiCop(4, tau = 0.5))

BiCopKDE(gumbel[, 1], gumbel[, 2])

contour(BiCop(3, tau = 0.5))

BiCopKDE(clayton[, 1], clayton[, 2])

  • Density plots of differing copulas
plot(BiCop(5, tau = 0.5), main = "Frank copula density")

plot(BiCop(4, tau = 0.5), main = "Gumbel copula density")

plot(BiCop(3, tau = 0.5), main = "Clayton copula density")

plot(BiCop(104, par = 3, par2 = 0.5), main = "Tawn type 1 copula ")

plot(BiCop(114, par = 3, par2 = 0.5), main = "Tawn type 1 copula (rot. 180)")

plot(BiCop(204, par = 3, par2 = 0.5), main = " Tawn type 2 copula")

contour(BiCop(family = 104, par = 3, par2 = 0.5))

contour(BiCop(family = 114, par = 3, par2 = 0.5))

contour(BiCop(family = 124, par = -3, par2 = 0.5))

contour(BiCop(family = 134, par = -3, par2 = 0.5))

contour(BiCop(204, par = 3, par2 = 0.5))

  • Changes in dependence measure
t1 <- BiCop(2, tau = 0.2, par2 = 3) %>% BiCopSim(10000, .)

t2 <- BiCop(2, tau = 0.5, par2 = 3) %>% BiCopSim(10000, .)

t3 <- BiCop(2, tau = 0.8, par2 = 3) %>% BiCopSim(10000, .)

p1 <- qplot(t1[, 1], t1[, 2], colour = t1[, 1], main = expression(paste("Student T  where ", 
    tau, "= 0.2")), xlab = "u", ylab = "v")

p2 <- qplot(t2[, 1], t2[, 2], colour = t2[, 1], main = expression(paste("Student T  where ", 
    tau, "= 0.5")), xlab = "u", ylab = "v")

p3 <- qplot(t3[, 1], t3[, 2], colour = t3[, 1], main = expression(paste("Student T  where ", 
    tau, "= 0.8")), xlab = "u", ylab = "v")

pushViewport(viewport(layout = grid.layout(1, 3)))
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(p3, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))

  • Changes in par2 measure (shape and scale parameters)
t1 <- BiCop(2, tau = 0.5, par2 = 3) %>% BiCopSim(10000, .)

t2 <- BiCop(2, tau = 0.5, par2 = 10) %>% BiCopSim(10000, .)

t3 <- BiCop(2, tau = 0.5, par2 = 1e+08) %>% BiCopSim(10000, .)

p1 <- qplot(t1[, 1], t1[, 2], colour = t1[, 1], main = expression(paste("Student T  where ", 
    upsilon, "= 3")), xlab = "u", ylab = "v")

p2 <- qplot(t2[, 1], t2[, 2], colour = t2[, 1], main = expression(paste("Student T  where ", 
    upsilon, "= 10")), xlab = "u", ylab = "v")

p3 <- qplot(t3[, 1], t3[, 2], colour = t3[, 1], main = expression(paste("Student T  where ", 
    upsilon, "= 10e7")), xlab = "u", ylab = "v")

pushViewport(viewport(layout = grid.layout(1, 3)))
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(p3, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))

Fitting Copulas

  • If you know the family of the copula that you are fitting then we can just use the BiCopEst function to get the parameters and test the fit
fit1 <- BiCopEst(clayton[, 1], clayton[, 2], family = 5)
summary(fit1)
## Family
## ------ 
## No:    5
## Name:  Frank
## 
## Parameter(s)
## ------------
## par:  5.67
## 
## Dependence measures
## -------------------
## Kendall's tau:    0.49 (empirical = 0.5, p value < 0.01)
## Upper TD:         0 
## Lower TD:         0 
## 
## Fit statistics
## --------------
## logLik:  3104.25 
## AIC:    -6206.51 
## BIC:    -6199.3
plot(fit1)

  • In terms of testing, there are several goodness of fit tests one can use, Ill just use whites. See goFCopula
# White's information matrix equality (White 1982) as
# introduced by Huang and Prokhorov (2011)
# BiCopGofTest(clayton[, 1], clayton[, 2], fit1)

BiCopKDE(clayton[, 1], clayton[, 2])

contour(fit1, col = 2, add = TRUE, drawlabels = FALSE)

  • If we don’t know the copula type, the VineCopula package provides a function which will determine the appropriate copula and get the parameters of the fit
# Best fit estimation
fit2 <- BiCopSelect(clayton[, 1], clayton[, 2], familyset = c(0:5))

BiCopKDE(clayton[, 1], clayton[, 2])
contour(fit1, col = 2, add = TRUE, drawlabels = FALSE)
contour(fit2, col = 4, add = TRUE, drawlabels = FALSE)

# BiCopGofTest(clayton[, 1], clayton[, 2], fit2)

Now, lets simulate some data from our compula

copsim <- BiCopSim(10000, fit2)

x <- copsim[, 2]
y <- copsim[, 1]
CDF <- BiCopCDF(copsim[, 1], copsim[, 2], BiCop(1, tau = 0.2))

CDF

PDF <- BiCopPDF(copsim[, 1], copsim[, 2], BiCop(1, tau = 0.2))

Lets have a look

# BiCopCompare(u[, 1], u[, 2])
plot(fit2)

# plot3d(x, y, PDF, col = 4)

Vine copula

data(daxreturns)
u <- daxreturns[, 1:4]
pairs.copuladata(u)

vinefit <- RVineStructureSelect(u, familyset = c(1:5))
summary(vinefit)
## tree     edge | family  cop   par   par2 |  tau   utd   ltd 
## ----------------------------------------------------------- 
##    1      2,3 |      2    t  0.61   4.62 | 0.42  0.29  0.29
##           1,2 |      2    t  0.59   4.61 | 0.40  0.28  0.28
##           4,1 |     14   SG  1.61   0.00 | 0.38     -  0.46
##    2    1,3;2 |      2    t  0.18  10.21 | 0.12  0.02  0.02
##         4,2;1 |      2    t  0.27  14.41 | 0.17  0.01  0.01
##    3  4,3;1,2 |      5    F  0.65   0.00 | 0.07     -     -
## ---
## type: D-vine    logLik: 868    AIC: -1716    BIC: -1665.46    
## ---
## 1 <-> ALV.DE,   2 <-> BAS.DE,   3 <-> BAYN.DE,   4 <-> BMW.DE
plot(vinefit, tree = 1, type = 2)

plot(vinefit, tree = 2, type = 2)

plot(vinefit, tree = 3, type = 2)

contour(vinefit)

usim <- RVineSim(1158, vinefit)
pairs(usim, pch = ".", main = "Simulated")

pairs(u, pch = ".", main = "Actuals")