Notebook: rvinecopulib basics

library('rvinecopulib')

set.seed(867-5309)

Bivariate copulas

Constructing and displaying

Note that the name parameter for the t-copula is given in the documentation as “student”, but the function actually takes the string “t”.

tcop <- bicop_dist('t', parameters=c(0.25, 4))
print(tcop)

Bivariate copula ('bicop_dist'): family = t, rotation = 0, parameters = 0.25, 4, var_types = c,c

plot(tcop)

contour(tcop)

Contour plots are generated with standard normal margins by default, but you can choose from standard uniform (“unif”), standard exponential (“exp”), or flipped exponential (“flexp”)

gcop <- bicop_dist('gumbel', rotation=90, parameters=2.5)
contour(gcop)

contour(gcop, margins='unif')

Density, distribution, conditional distribution, and RNG functions

Density and distribution functions are computed on (standard uniform) pseudo-observations.

u1 <- matrix(c(0.8, 0.2, 0.9, 0.2, 0.8, 0.9), ncol=2)
u1

     [,1] [,2]
[1,]  0.8  0.2
[2,]  0.2  0.8
[3,]  0.9  0.9

Density:

dbicop(u1, gcop)

[1] 2.220596 3.087756 0.017942

Probability (distribution function):

pbicop(u1, gcop)

[1] 0.08040761 0.05505112 0.80004124

The “h-functions” are conditional distributions. The cond_var argument determines which is the conditioning variable. So, \[h_1(u,v) = P(V\leq v | U = u),\] and \[h_2(u,v) = P(U \leq u | V = v).\]

u2 <- matrix(c(0.2, 0.4, 0.6, 0.8, 0.8, 0.8), ncol=2)
u2

     [,1] [,2]
[1,]  0.2  0.8
[2,]  0.4  0.8
[3,]  0.6  0.8

Unconditional (joint):

pbicop(u2, gcop)

[1] 0.05505112 0.21472948 0.40423109

Conditional:

hbicop(u2, gcop, cond_var = 1)

[1] 0.6143537 0.9083536 0.9724446

The RNG just takes a number of points to generate. The margins of the outputs are standard uniform.

ur <- rbicop(1000, gcop)
plot(ur)

You can also specify quasi-random instead of pseudorandom points.

uq <- rbicop(1000, gcop, qrng=TRUE)
plot(uq)

Vine copulas

Vine structure

Vines are represented in rvinecopulib using a special structure. The easiest way to create one of those structures is to convert it from the matrix representation. The vine matrices used here are flipped left to right from the version in Czado’s book, but they are otherwise compatible.

Here is a vine constructed from Example 5.3 from Czado (sort of. \(T_1\) is the same, but I constructed the subsequent trees a little differently):

m5.3 <- matrix(c(5,1,2,3,4,6, 1,2,3,4,5,0, 3,1,2,4,0,0, 1,2,3,rep(0,3), 1,2,rep(0,4), 1,rep(0,5)), ncol=6)
m5.3

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    5    1    3    1    1    1
[2,]    1    2    1    2    2    0
[3,]    2    3    2    3    0    0
[4,]    3    4    4    0    0    0
[5,]    4    5    0    0    0    0
[6,]    6    0    0    0    0    0

v5.3 <- as_rvine_structure(m5.3)
v5.3

6-dimensional R-vine structure ('rvine_structure')
5 1 3 1 1 1
1 2 1 2 2  
2 3 2 3    
3 4 4      
4 5        
6          

plot(v5.3, tree=1:5)

There are convenience functions for constructing D-Vines and C-Vines. (For the first parameter you could also just put the number of dimensions, and the labels would be chosen in a default order)

vd5 <- dvine_structure(seq(5,1))
vd5

5-dimensional R-vine structure ('rvine_structure')
4 3 2 1 1
3 2 1 2  
2 1 3    
1 4      
5        

plot(vd5, tree='ALL')

vc5 <- cvine_structure(5)
plot(vc5, tree='ALL')

Specifying a “truncation level” allows you to suppress higher-order interactions.

vd5_2 <- dvine_structure(seq(5,1), 2)
print(vd5_2)

5-dimensional R-vine structure ('rvine_structure'), 2-truncated
4 3 2 1 1
3 2 1 2  
    3    
  4      
5        

dim(vd5_2)

      dim trunc_lvl 
        5         2 

plot(vd5_2, tree='ALL')

Vine copula distributions

You can build a vine copula distribution with a vine copula structure and a bunch of bivariate copulas. The bivariate copulas are organized into a list of lists, with each element of the top-level list corresponding to a tree (organized as \(T_1\) through \(T_n\)), and each element of the second-level lists corresponding to one of the edges in that tree. Within each tree the edges are taken left-to-right from the columns of the structure matrix.

ci <- bicop_dist()
cnp <- bicop_dist('gaussian', parameters=0.7)
cnn <- bicop_dist('gaussian', parameters=-0.7)
pclist <- list(
  list(cnp, ci, ci, cnn),
  list(ci,ci,ci),
  list(ci,ci),
  list(ci)
)
vcop <- vinecop_dist(pclist, vd5)
contour(vcop)

Density, distribution, and RNG functions

Like the bivariate copulas, there are density and distribution functions for vine copulas. (There is, however, no equivalent of the h-functions.)

u4 <- matrix(c(
  rep(0.5, 5),           # center of the distribution
  rep(0.75, 5),          
  c(0.25, rep(0.75, 4)), # Higher density than the previous b/c 1,2 are anticorrelated
  c(rep(0.75,4), 0.25)   # Low density b/c 1,2 and 4,5 have the wrong association.
  ),
  ncol=5,
  byrow=TRUE)
u4

     [,1] [,2] [,3] [,4] [,5]
[1,] 0.50 0.50 0.50 0.50 0.50
[2,] 0.75 0.75 0.75 0.75 0.75
[3,] 0.25 0.75 0.75 0.75 0.75
[4,] 0.75 0.75 0.75 0.75 0.25

dvinecop(u4, vcop, cores=2)  # Supports multicore evaluation!

[1] 1.9607843 0.8180376 2.8519360 0.2346425

pvinecop(u4, vcop, cores=2)

[1] 0.0234 0.2459 0.0491 0.0924

The RNG function allows for a QRNG option. With Gaussian bivariate copulas across the board, The D-Vine with truncation level 2 should be equivalent to an AR(2) model, albeit without the ability to easily extend the series to arbitrary length. On the other hand, we can introduce tail dependence and other such using the copula model.

vcop2 <- vinecop_dist(
  list(rep(list(tcop), 4),
       rep(list(cnn), 3)),   # after controlling for the lag-1 var, the lag-2 var is anticorrelated
  vd5_2
)
contour(vcop2)

z <- qnorm(rvinecop(500, vcop2, qrng=TRUE, cores=2))

GGally::ggpairs(as.data.frame(z), progress=FALSE)   # plot on std. normal margins