Quasi-Random Numbers for Copula Models

Marius Hofert

2018-12-20

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”). Note that not all plots are displayed (to keep the tarball small).

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

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

Independence copula

Let’s start with something known, the independence case.

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

Clayton copula

Consider a Clayton copula.

\(t\) copula with three degrees of freedom

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

Marshall–Olkin copula

Now something more fancy, a Marshall–Olkin copula.

3d \(t\) copula with three degrees of freedom

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

First the pseudo-random version.

Now the quasi-random version.

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.

3d R-Vine copula

As another example, consider a three-dimensional R-vine copula. Note that this is not run here to avoid a cyclic dependency (since VineCopula imports copula).

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.

2.2 A variance-reduction example

In this example, we would like to investigate the standard deviation when estimating expected shortfall at \(\alpha=99\%\) confidence level via Monte Carlo simulation based on a Clayton copula with Pareto margins. To this end we consider pseudo-random numbers and quasi-random numbers, as well as two different sampling methods for the Clayton copula (the conditional distribution method and the Marshall–Olkin method (based on a well-known stochastic representation)), hence four different sampling methods.

Here is our setup.

Next, let’s implement a function which can sample the Clayton copula with one of the four approaches.

##' @title Pseudo-/quasi-random number generation for (survival) Clayton copulas
##' @param n Sample size
##' @param d Dimension
##' @param B Number of replications
##' @param theta Clayton parameter
##' @param survival Logicial indicating whether a sample from the survival copula
##'        should be returned
##' @param rng.method Pseudo-/quasi-random number generator
##' @param cop.method Method to construct the pseudo-/quasi-random copula sample
##' @return (n, d, B)-array of pseudo-/quasi-random copula sample
##' @author Marius Hofert
rng_Clayton <- function(n, d, B, theta, survival = FALSE,
                        rng.method = c("runif", "ghalton"),
                        cop.method = c("CDM", "MO"))
{
    ## Sanity checks
    stopifnot(n >= 1, d >= 2, B >= 1, is.logical(survival))
    rng.method <- match.arg(rng.method)
    cop.method <- match.arg(cop.method)

    ## Draw U(0,1) random numbers
    k <- if(cop.method == "CDM") d else d+1
    U. <- switch(rng.method,
    "runif" = {
        array(runif(n*k*B), dim = c(n,k,B)) # (n, k, B)-array
    },
    "ghalton" = {
        replicate(B, expr = ghalton(n, d = k)) # (n, k, B)-array
    },
    stop("Wrong 'rng.method'"))

    ## Convert to pseudo-/quasi-random copula samples
    U <- switch(cop.method, # B-list of (n, d)-matrices
    "CDM" = {
        cop <- onacopulaL("Clayton", nacList = list(theta, 1:d)) # d = k here
        lst <- apply(U., 3, FUN = function(x) list(cCopula(x, copula = cop, inverse = TRUE)))
        lapply(lst, `[[`, 1)
    },
    "MO" = {
        lapply(1:B, function(b) {
            copClayton@psi(-log(U.[,2:k,b]) / qgamma(U.[,1,b], 1/theta), theta = theta)
        })
    },
    stop("Wrong 'cop.method'"))

    ## Return
    if(survival) 1-U else U # B-list of (n, d)-matrices
}

We also need an estimator of expected shortfall; we use the empirical estimator here.

For each of the four methods, each of the chosen sample sizes and the number \(B\) of (bootstrap) replications considered here, generated the samples, aggregate them and compute expected shortfall at the 99% confidence level. To reduce increase comparability, note that samples with smaller sample size are subsets of samples with larger sample size.

Now we can compute the standard deviations, including estimated power curves based on all data stemming from pseudo-random numbers and all data stemming from quasi-random numbers.

And now the results. In a nutshell, in comparison to pseudo-random numbers, quasi-random numbers for copula models can reduce the standard deviations (or variances), and this holds not only for one-to-one transformations such as the conditional distribution method but also for well-known stochastic representations such as Marshall–Olkin’s. Note that the results are more pronounced for larger sample sizes; see also Cambou et al. (2016, “Quasi-random numbers for copula models”).