Examples of Nonstandard Copulas – “Wild Animals”

Marius Hofert and Martin Mächler

2019-04-21

1 Swiss Alps copulas of Hofert, Vrins (2013)

This example implements the Swiss Alps copulas of Hofert, Vrins (2013, “Sibuya copulas”).

Lambda and its inverse

M and its inverse (for \(M_i, i=1,2\)):

S and its inverse (for \(S_i, i=1,2\))

Wrappers for \(p_1\) and \(p_2\) and their inverses:

p and its inverse (for \(p_i(t_k-)\) and \(p_i(t_k), i = 1,2\)):

and the wrappers, which work with p1(), p2(), p1Inv(), p2Inv() as arguments:

Define the copula \(C\)

Compute the singular component (given \(u_1=\) u1, find \(u_2=\) u2 on the singular component) \(u_2 = S2(S1^{-1}(u_1))\):

Generate one bivariate random vector from C:

Draw n vectors of random variates from \(C\)

Wireframe plot to incorporate singular component :

require(lattice)
wf.plot <- function(grid, val.grid, s.comp, val.s.comp, Lambda, S1Inv, S2Inv){
    wireframe(val.grid ~ grid[,1]*grid[,2], xlim=c(0,1), ylim=c(0,1), zlim=c(0,1),
              aspect=1, scales = list(arrows=FALSE, col=1), # remove arrows
              par.settings= list(axis.line = list(col="transparent"), # remove global box
                                 clip = list(panel="off")),
              pts = cbind(s.comp, val.s.comp), # <- add singular component
              panel.3d.wireframe = function(x, y, z, xlim, ylim, zlim, xlim.scaled,
                                            ylim.scaled, zlim.scaled, pts, ...) {
                  panel.3dwire(x=x, y=y, z=z, xlim=xlim, ylim=ylim, zlim=zlim,
                               xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
                               zlim.scaled=zlim.scaled, ...)
                  xx <- xlim.scaled[1]+diff(xlim.scaled)*(pts[,1]-xlim[1])/diff(xlim)
                  yy <- ylim.scaled[1]+diff(ylim.scaled)*(pts[,2]-ylim[1])/diff(ylim)
                  zz <- zlim.scaled[1]+diff(zlim.scaled)*(pts[,3]-zlim[1])/diff(zlim)
                  panel.3dscatter(x=xx, y=yy, z=zz, xlim=xlim, ylim=ylim, zlim=zlim,
                                  xlim.scaled=xlim.scaled, ylim.scaled=ylim.scaled,
                                  zlim.scaled=zlim.scaled, type="l", col=1, ...)
              },
              xlab = expression(italic(u[1])),
              ylab = expression(italic(u[2])),
              zlab = list(expression(italic(C(u[1],u[2]))), rot=90))
}

## Copula plot with singular component
u <- seq(0, 1, length.out=20) # grid points per dimension
grid <- expand.grid(u1=u, u2=u) # grid
val.grid <- apply(grid, 1, C, H=H, Lambda=Lambda, S1Inv=S1Inv, S2Inv=S2Inv) # copula values on grid
s.comp <- cbind(u, sapply(u, s.comp, H=H, S1Inv=S1Inv, S2=S2)) # pairs (u1, u2) on singular component
val.s.comp <- apply(s.comp, 1, C, H=H, Lambda=Lambda, S1Inv=S1Inv, S2Inv=S2Inv) # corresponding z-values
wf.plot(grid=grid, val.grid=val.grid, s.comp=s.comp, val.s.comp=val.s.comp,
        Lambda=Lambda, S1Inv=S1Inv, S2Inv=S2Inv)

2 An example from Wolfgang Trutschnig and Manuela Schreyer

For more details, see Trutschnig, Fernandez Sanchez (2014) “Copulas with continuous, strictly increasing singular conditional distribution functions”

Roughly, one defines an Iterated Function System whose attractor is the word “Copula” and starts the chaos game.

Define the Iterated Function System

Run chaos game B times

\(X :=\) Rotate \(A\) by \(-45^{o} = -\pi/4\) :

Now transform the margings by their marginal ECDF’s so we get uniform margins. Note that, it is equivalent but faster to use rank(*, ties.method="max"):

Now, visually check the margins of U; they are perfectly uniform:

whereas U, the copula sample, indeed is peculiar and contains the word “COPULA” many times if you look closely (well, the “L” is defect …):

3 Sierpinski tetrahedron

This is an implementation of Example 2.3 in https://arxiv.org/pdf/0906.4853.pdf

library(abind) # for merging arrays via abind()
library(lattice) # for cloud()
library(sfsmisc) # for polyn.eval()
## 
## Attaching package: 'sfsmisc'
## The following object is masked _by_ '.GlobalEnv':
## 
##     Sys.memGB

Implement the random number generator:

##' @title Generate samples from the Sierpinski tetrahedron
##' @param n sample size
##' @param N digits in the base-2 expansion
##' @return (n, 3)-matrix
##' @author Marius Hofert
rSierpinskyTetrahedron <- function(n, N)
{
    stopifnot(n >= 1, N >= 1)
    ## Build coefficients in the base-2 expansion
    U12coeff <- array(sample(0:1, size = 2*n*N, replace = TRUE),
                      dim = c(2, n, N), dimnames = list(U12 = c("U1", "U2"),
                                                        sample = 1:n,
                                                        base2digit = 1:N)) # (2, n, N)-array
    U3coeff <- apply(U12coeff, 2:3, function(x) sum(x) %% 2) # (n, N)-matrix
    Ucoeff <- abind(U12coeff, U3 = U3coeff, along = 1)
    ## Convert to U's
    t(apply(Ucoeff, 1:2, function(x)
        polyn.eval(coef = rev(x), x = 2))/2^N) # see sfsmisc::bi2int
}

Draw vectors of random numbers following a “Sierpinski tetrahedron copula”:

set.seed(271)
U <- rSierpinskyTetrahedron(1e4, N = 6)

Use a scatterplot matrix to check all bivariate margins:

pairs(U, gap = 0, cex = 0.25, col = "black",
      labels = as.expression( sapply(1:3, function(j) bquote(U[.(j)])) ))

All pairs “look” independent but, of course, they aren’t:

cloud(U[,3] ~ U[,1] * U[,2], cex = 0.25, col = "black", zoom = 1,
      scales = list(arrows = FALSE, col = "black"), # ticks instead of arrows
      par.settings = list(axis.line = list(col = "transparent"), # to remove box
                          clip = list(panel = "off"),
                          standard.theme(color = FALSE)),
      xlab = expression(U[1]), ylab = expression(U[2]), zlab = expression(U[3]))

Session information

## R version 3.6.0 RC (2019-04-19 r76407)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Fedora 28 (Workstation Edition)
## 
## Matrix products: default
## BLAS:   /data/gannet/ripley/R/R-patched/lib/libRblas.so
## LAPACK: /data/gannet/ripley/R/R-patched/lib/libRlapack.so
## 
## attached base packages:
##  [1] splines   parallel  grid      stats4    tools     stats     graphics 
##  [8] grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] sfsmisc_1.1-3      abind_1.4-5        randtoolbox_1.17.1
##  [4] rngWELL_0.10-5     qrng_0.0-5         gridExtra_2.3     
##  [7] VGAM_1.1-1         rugarch_1.4-1      gsl_2.1-6         
## [10] mev_1.11           lattice_0.20-38    bbmle_1.0.20      
## [13] copula_0.999-19.1 
## 
## loaded via a namespace (and not attached):
##  [1] assertthat_0.2.1            expm_0.999-4               
##  [3] yaml_2.2.0                  GeneralizedHyperbolic_0.8-4
##  [5] numDeriv_2016.8-1           pillar_1.3.1               
##  [7] glue_1.3.1                  DistributionUtils_0.6-0    
##  [9] digest_0.6.18               colorspace_1.4-1           
## [11] sandwich_2.5-1              htmltools_0.3.6            
## [13] Matrix_1.2-17               plyr_1.8.4                 
## [15] pcaPP_1.9-73                pkgconfig_2.0.2            
## [17] ismev_1.42                  purrr_0.3.2                
## [19] mvtnorm_1.0-10              scales_1.0.0               
## [21] rootSolve_1.7               gmm_1.6-2                  
## [23] tibble_2.1.1                nleqslv_3.3.2              
## [25] ADGofTest_0.3               mgcv_1.8-28                
## [27] gmp_0.5-13.5                bayesplot_1.6.0            
## [29] ggplot2_3.1.1               lazyeval_0.2.2             
## [31] magrittr_1.5                crayon_1.3.4               
## [33] mclust_5.4.3                evaluate_0.13              
## [35] ks_1.11.4                   nlme_3.1-139               
## [37] MASS_7.3-51.4               xts_0.11-2                 
## [39] truncnorm_1.0-8             pspline_1.0-18             
## [41] stringr_1.4.0               revdbayes_1.3.3            
## [43] munsell_0.5.0               stabledist_0.7-1           
## [45] compiler_3.6.0              SkewHyperbolic_0.4-0       
## [47] evd_2.3-3                   rlang_0.3.4                
## [49] nloptr_1.2.1                ggridges_0.5.1             
## [51] Rsolnp_1.16                 Runuran_0.27               
## [53] spd_2.0-1                   partitions_1.9-19          
## [55] rmarkdown_1.12              boot_1.3-22                
## [57] gtable_0.3.0                polynom_1.4-0              
## [59] R6_2.4.0                    zoo_1.8-5                  
## [61] knitr_1.22                  dplyr_0.8.0.1              
## [63] KernSmooth_2.23-15          stringi_1.4.3              
## [65] Rcpp_1.0.1                  tidyselect_0.2.5           
## [67] xfun_0.6