# Quasi-Random Numbers for Copula Models

#### 2017-08-31

In this vignette, we present some findings around quasi-random numbers for copula models; see also Cambou et al. (2016, ``Quasi-random numbers for copula models’’).

``````library(lattice)
library(copula)
library(VGAM)
library(VineCopula)
library(gridExtra)
library(qrng)
library(randtoolbox)``````

## 1 Quasi-random numbers for copula models via conditional distribution method

### Independence copula

``````n <- 1000 # sample size
set.seed(271) # set the seed (for reproducibility)
U <- matrix(runif(n*2), ncol = 2) # pseudo-random numbers
U. <- halton(n, dim = 2) # quasi-random numbers
par(pty = "s", mfrow = 1:2)
plot(U,  xlab = expression(italic(U)*"'"), ylab = expression(italic(U)*"'"))
plot(U., xlab = expression(italic(U)*"'"), ylab = expression(italic(U)*"'"))`````` Let’s check if the more equally spaced points (less gaps, less clusters) are preserved in the copula world when determined with one-to-one transformations (such as the conditional distribution method (CDM); this can be obtained via `cCopula(, inverse=TRUE)`).

### \(t\) copula with three degrees of freedom

Consider a \(t\) copula with three degrees of freedom.

``````family <- "t"
nu <- 3 # degrees of freedom
tau <- 0.5 # Kendall's tau (determines the copula parameter rho)
th <- iTau(ellipCopula(family, df = nu), tau)
cop <- ellipCopula(family, param = th, df = nu)``````
``````U.t  <- cCopula(U,  copula = cop, inverse = TRUE) # via PRNG
U.t. <- cCopula(U., copula = cop, inverse = TRUE) # via QRNG
par(pty = "s", mfrow = 1:2)
plot(U.t,  xlab = expression(italic(U)), ylab = expression(italic(U)))
plot(U.t., xlab = expression(italic(U)), ylab = expression(italic(U)))`````` ### Clayton copula

``````family <- "Clayton"
tau <- 0.5
th <- iTau(getAcop(family), tau)
cop <- onacopulaL(family, nacList = list(th, 1:2))``````
``````U.C  <- cCopula(U,  copula = cop, inverse = TRUE) # via PRNG
U.C. <- cCopula(U., copula = cop, inverse = TRUE) # via QRNG
par(pty = "s", mfrow = 1:2)
plot(U.C,  xlab = expression(italic(U)), ylab = expression(italic(U)))
plot(U.C., xlab = expression(italic(U)), ylab = expression(italic(U)))`````` ### Marshall–Olkin copula

``````alpha <- c(0.25, 0.75)
tau <- (alpha*alpha) / (alpha+alpha-alpha*alpha)``````
``````##' @title Inverse of the bivariate conditional Marshall--Olkin copula
##' @param u (n,2) matrix of U[0,1] random numbers to be transformed to
##'        (u[,1], C^-(u[,2]|u[,1]))
##' @param alpha bivariate parameter vector
##' @return (u[,1], C^-(u[,2]|u[,1])) for C being a MO copula
##' @author Marius Hofert
inv_cond_cop_MO <- function(u, alpha)
{
stopifnot(is.matrix(u), 0 <= alpha, alpha <= 1)
up <- u[,1]^(alpha*(1/alpha-1))
low <- (1-alpha)*up
i1 <- u[,2] <= low
i3 <- u[,2] >  up
u2 <- u[,1]^(alpha/alpha)
u2[i1] <- u[i1,1]^alpha * u[i1,2] / (1-alpha)
u2[i3] <- u[i3,2]^(1/(1-alpha))
cbind(u[,1], u2)
}``````
``````U.MO  <- inv_cond_cop_MO(U,  alpha = alpha)
U.MO. <- inv_cond_cop_MO(U., alpha = alpha)
par(pty = "s", mfrow = 1:2)
plot(U.MO,  xlab = expression(italic(U)), ylab = expression(italic(U)))
plot(U.MO., xlab = expression(italic(U)), ylab = expression(italic(U)))`````` ### 3d \(t\) copula with three degrees of freedom

Let’s consider a three-dimensional \(t\) copula with three degrees of freedom.

``````family <- "t"
nu <- 3 # degrees of freedom
tau <- 0.5 # Kendall's tau (determines the copula parameter rho)
th <- iTau(ellipCopula(family, df = nu), tau)
cop <- ellipCopula(family, param = th, dim = 3, df = nu)``````

First the pseudo-random version.

``````U.3d <- matrix(runif(n*3), ncol = 3)
U.t.3d <- cCopula(U.3d, copula = cop, inverse = TRUE)
par(pty = "s")
pairs(U.t.3d, gap = 0,
labels = as.expression(sapply(1:3, function(j) bquote(italic(U[.(j)])))))`````` Now the quasi-random version.

``````U.3d. <- halton(n, dim = 3)
U.t.3d. <- cCopula(U.3d., copula = cop, inverse = TRUE)
par(pty = "s")
pairs(U.t.3d., gap = 0,
labels = as.expression(sapply(1:3, function(j) bquote(italic(U[.(j)])))))`````` Note that projections (here: to pairs) can appear not to be `quasi-random’ (or appear not to possess a lower discrepancy), but see Section 2.2 below! Visualization in more than two dimensions seems difficult; we have just seen bivariate projections and ‘quasi-randomness’ is also not easily visible from a 3d cloud plot.

``````p1 <- cloud(U.t.3d[,3]~U.t.3d[,1]+U.t.3d[,2], scales = list(col = 1, arrows = FALSE),
col = 1, xlab = expression(italic(U)), ylab = expression(italic(U)),
zlab = expression(italic(U)),
par.settings = list(background = list(col = "#ffffff00"),
axis.line = list(col = "transparent"),
clip = list(panel = "off")))
p2 <- cloud(U.t.3d.[,3]~U.t.3d.[,1]+U.t.3d.[,2], scales = list(col = 1, arrows = FALSE),
col = 1, xlab = expression(italic(U)), ylab = expression(italic(U)),
zlab = expression(italic(U)),
par.settings = list(background = list(col = "#ffffff00"),
axis.line = list(col = "transparent"),
clip = list(panel = "off")))
grid.arrange(p1, p2, ncol = 2)`````` ### 3d R-Vine copula

As another example, consider a three-dimensional R-vine copula.

``````M <- matrix(c(3, 1, 2,
0, 2, 1,
0, 0, 1), ncol = 3) # R-vine tree structure matrix
family <- matrix(c(0, 3, 3, # C, C
0, 0, 3, #    C
0, 0, 0), ncol = 3) # R-vine pair-copula family matrix (0 = Pi)
param <- matrix(c(0, 1, 1,
0, 0, 1,
0, 0, 0), ncol = 3) # R-vine pair-copula parameter matrix
param. <- matrix(0, nrow = 3, ncol = 3) # 2nd R-vine pair-copula parameter matrix
RVM <- RVineMatrix(Matrix = M, family = family, par = param, par2 = param.) # RVineMatrix object``````

First the pseudo-random version.

``````U <- RVineSim(n, RVM) # PRNG
par(pty = "s")
pairs(U, labels = as.expression( sapply(1:3, function(j) bquote(italic(U[.(j)]))) ),
gap = 0, cex = 0.3)``````