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

Let's start with something know, the independence case.

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)[1]*"'"), ylab = expression(italic(U)[2]*"'"))
plot(U., xlab = expression(italic(U)[1]*"'"), ylab = expression(italic(U)[2]*"'"))

plot of chunk unnamed-chunk-2

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)[1]), ylab = expression(italic(U)[2]))
plot(U.t., xlab = expression(italic(U)[1]), ylab = expression(italic(U)[2]))

plot of chunk unnamed-chunk-4

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)[1]), ylab = expression(italic(U)[2]))
plot(U.C., xlab = expression(italic(U)[1]), ylab = expression(italic(U)[2]))

plot of chunk unnamed-chunk-6

Marshall–Olkin copula

alpha <- c(0.25, 0.75)
tau <- (alpha[1]*alpha[2]) / (alpha[1]+alpha[2]-alpha[1]*alpha[2])
##' @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]*(1/alpha[2]-1))
    low <- (1-alpha[1])*up
    i1 <- u[,2] <= low
    i3 <- u[,2] >  up
    u2 <- u[,1]^(alpha[1]/alpha[2])
    u2[i1] <- u[i1,1]^alpha[1] * u[i1,2] / (1-alpha[1])
    u2[i3] <- u[i3,2]^(1/(1-alpha[2]))
    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)[1]), ylab = expression(italic(U)[2]))
plot(U.MO., xlab = expression(italic(U)[1]), ylab = expression(italic(U)[2]))

plot of chunk unnamed-chunk-9

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

plot of chunk unnamed-chunk-11

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

plot of chunk unnamed-chunk-12 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)[1]), ylab = expression(italic(U)[2]),
            zlab = expression(italic(U[3])),
            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)[1]), ylab = expression(italic(U)[2]),
            zlab = expression(italic(U[3])),
            par.settings = list(background = list(col = "#ffffff00"),
                                axis.line = list(col = "transparent"),
                                clip = list(panel = "off")))
grid.arrange(p1, p2, ncol = 2)

plot of chunk unnamed-chunk-13

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)

plot of chunk unnamed-chunk-15

Now the quasi-random version.

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

plot of chunk unnamed-chunk-16 Similarly to the 3d t copula case (because of the projections to pairs), not all pairs appear to be 'quasi-random'.

2 Quasi-random numbers for copula models via stochastic representations

For many copula families, it is rarely efficient to sample them via the CDM (or other one-to-one transformations), one typically uses stochastic representations based on simple, easy-to-sample distributions as building blocks. Although, again, not directly visible, quasi-random numbers can also improve the low-discrepancy of the resulting random numbers and thus be used for variance reduction in the context of dependence.

2.1 Colorized scatter plot

To explore this, we sample from a Clayton copula via the CDM (so via a one-to-one transformation) and via the Marshall–Olkin algorithm (so via a stochastic representation in terms of the Gamma frailty distribution and two standard exponentials) based on a three-dimensional Halton sequence.

n <- 1000
family <- "Clayton"
tau <- 0.5
th <- iTau(getAcop(family), tau)
cop <- onacopulaL(family, nacList = list(th, 1:2))
## Generate dependent samples
U <- halton(n, 3)
U_CDM <- cCopula(U[,1:2], copula = cop, inverse = TRUE) # via CDM
U_MO <- copClayton@psi(-log(U[,2:3]) / qgamma(U[,1], 1/th), theta = th) # via Marshall-Olkin (MO)

## Colorization of U[,1:2]
col <- rep("black", n)
col[U[,1] <= 0.5 & U[,2] <= 0.5] <- "maroon3"
col[U[,1] >= 0.5 & U[,2] >= 0.5] <- "royalblue3"

## Colorization of U[,1:3] (= U)
col. <- rep("black", n)
col.[apply(U <= 0.5, 1, all)] <- "maroon3"
col.[apply(U >= 0.5, 1, all)] <- "royalblue3"

Colorized scatter plot (quasi-random numbers and CDM)

par(pty = "s", mfrow = 1:2)
plot(U[,1:2], xlab = expression(italic(U)[1]*"'"), ylab = expression(italic(U)[2]*"'"), col = col)
plot(U_CDM,   xlab = expression(italic(U)[1]), ylab = expression(italic(U)[2]), col = col)