Examples of Nonstandard Copulas – “Wild Animals”

M1 <- function(t) 0.01*t M1Inv <- function(t) t/0.01 M2 <- function(t) { f <- function(t.) { ct <- ceiling(t.) if(ct %% 2) t. - (ct-1)/2 else ct/2 } unlist(lapply(t, f)) } M2Inv <- function(t) { f <- function(t.) { if(t. == 0) return(0) t. + ceiling(t) - 1 } unlist(lapply(t, f)) }

S1 <- function(t, H) exp(-(M1(t) + Lambda(t)*(1-exp(-H)))) S1Inv <- function(u, H, upper=1e6) unlist(lapply(u, function(u.) uniroot(function(x) S1(x,H=H)-u., interval=c(0, upper))$root)) S2 <- function(t, H) exp(-(M2(t) + Lambda(t)*(1-exp(-H)))) S2Inv <- function(u, H, upper=1e6) unlist(lapply(u, function(u.) uniroot(function(x) S2(x,H=H)-u., interval=c(0, upper))$root))

p1 <- function(t, k, H) exp(-M1(t)-H*k) p1Inv <- function(u, k, H, upper=1e6) unlist(lapply(u, function(u.) uniroot(function(x) p1(x,k=k,H=H)-u., interval=c(0, upper))$root)) p2 <- function(t, k, H) exp(-M2(t)-H*k) p2Inv <- function(u, k, H, upper=1e6) unlist(lapply(u, function(u.) uniroot(function(x) p2(x,k=k,H=H)-u., interval=c(0, upper))$root))

p <- function(t, k, H, I, p1, p2){ if((lI <- length(I)) == 0){ stop("error in p") }else if(lI==1){ if(I==1) p1(t, k=k, H=H) else p2(t, k=k, H=H) }else{ # lI == 2 c(p1(t, k=k, H=H), p2(t, k=k, H=H)) } } pInv <- function(u, k, H, I, p1Inv, p2Inv){ if((lI <- length(I)) == 0){ stop("error in pInv") }else if(lI==1){ if(I==1) p1Inv(u, k=k, H=H) else p2Inv(u, k=k, H=H) }else{ # lI == 2 c(p1Inv(u[1], k=k, H=H), p2Inv(u[2], k=k, H=H)) } }

rC1 <- function(H, LambdaInv, S1, S2, p1, p2, p1Inv, p2Inv) { d <- 2 # dim = 2 ## (1) U <- runif(d) # for determining the default times of all components ## (2) -- t_{h,0} := initial value for the occurrence of the homogeneous ## Poisson process with unit intensity t_h <- 0 t <- 0 # t_0; initial value for the occurrence of the jump process k <- 1 # indices for the sets I I_ <- list(1:d) # I_k; indices for which tau has to be determined tau <- rep(Inf,d) # in the beginning, set all default times to Inf (= "no default") ## (3) repeat{ ## (4) ## k-th occurrence of a homogeneous Poisson process with unit intensity: t_h[k] <- rexp(1) + if(k == 1) 0 else t_h[k-1] ## k-th occ. of non-homogeneous Poisson proc. with integrated rate function Lambda: t[k] <- LambdaInv(t_h[k]) ## (5) pvec <- p(t[k], k=k, H=H, I=c(1,2), p1=p1, p2=p2) # c(p_1(t_k), p_2(t_k)) I. <- (1:d)[U >= pvec] # determine all i in I I <- intersect(I., I_[[k]]) # determine I ## (6)--(10) if(length(I) > 0){ default.at.t <- U[I] <= p(t[k], k=k-1, H=H, I=I, p1=p1, p2=p2) tau[I[default.at.t]] <- t[k] Ic <- I[!default.at.t] # I complement if(length(Ic) > 0) tau[Ic] <- pInv(U[Ic], k=k-1, H=H, I=Ic, p1Inv=p1Inv, p2Inv=p2Inv) } ## (11) -- define I_{k+1} := I_k \ I I_[[k+1]] <- setdiff(I_[[k]], I) ## (12) if(length(I_[[k+1]]) == 0) break else k <- k+1 } ## (14) c(S1(tau[1], H=H), S2(tau[2], H=H)) }

rC <- function(n, H, LambdaInv, S1, S2, p1, p2, p1Inv, p2Inv){ mat <- t(sapply(rep(H, n), rC1, LambdaInv=LambdaInv, S1=S1, S2=S2, p1=p1, p2=p2, p1Inv=p1Inv, p2Inv=p2Inv)) row.names(mat) <- NULL mat } ## Generate copula data n <- 2e4 H <- 10 set.seed(271) U <- rC(n, H=H, LambdaInv=LambdaInv, S1=S1, S2=S2, p1=p1, p2=p2, p1Inv=p1Inv, p2Inv=p2Inv) ## Check margins of U par(pty="s") hist(U[,1], probability=TRUE, main="Histogram of the first component", xlab=expression(italic(U[1])))

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)

IFS <- local({ ## Using `local`, so `n` is part of IFS n <- 23 list(function(x) c(3*x[1]/n, x[2]/4), function(x) c(-(x[2]-1)/n, x[1]/2+1/4), function(x) c(3*x[1]/n, x[2]/4+3/4), function(x) c((3*x[1]+4)/n, x[2]/4), function(x) c(-(x[2]-5)/n, x[1]/2+1/4), function(x) c((3*x[1]+4)/n, x[2]/4+3/4), function(x) c(-(x[2]-7)/n, x[1]/2+1/4), function(x) c(-x[2]/n+9/n, 3*x[1]/4), function(x) c((3*x[1]+8)/n, x[2]/4+3/4), function(x) c(x[1]/n+10/n, x[2]/8+1/2+1/8), function(x) c(2*x[1]/n+9/n, x[2]/4+1/4+1/8), function(x) c(-x[2]/n+13/n, (3*x[1]+1)/4), function(x) c((3*x[1]+12)/n, x[2]/4), function(x) c((3*x[1]+12)/n, x[2]/4), function(x) c(-x[2]/n+15/n, (3*x[1]+1)/4), function(x) c((3*x[1]+16)/n, x[2]/4), function(x) c(-(x[2]-21)/n, 3*x[1]/4), function(x) c((3*x[1]+20)/n, x[2]/4+3/4), function(x) c((x[1]+21)/n, x[2]/4+1/4+1/8), function(x) c(-(x[2]-23)/n, 3*x[1]/4)) })

B <- 20 # replications n.steps <- 20000 # number of steps AA <- vector("list", length=B) set.seed(271) for(i in 1:B) { ind <- sample(length(IFS), size=n.steps, replace=TRUE) # (randomly) 'bootstrap' functions res <- matrix(0, nrow=n.steps+1, ncol=2) # result matrix (for each i) pt <- c(0, 0) # initial point for(r in seq_len(n.steps)) { res[r+1,] <- IFS[[ind[r]]](pt) # evaluate randomly chosen functions at pt pt <- res[r+1,] # redefine point } AA[[i]] <- res # keep this matrix } A <- do.call(rbind, AA) # rbind (n.steps+1, 2)-matrices n <- nrow(A) stopifnot(ncol(A) == 2, n == B*(n.steps+1)) # sanity check

par(pty="s") sfsmisc::mult.fig(mfcol = c(1,2), main = "Margins are uniform") hist(U[,1], probability=TRUE, main="Histogram of U[,1]", xlab=quote(italic(U[1]))) hist(U[,2], probability=TRUE, main="Histogram of U[,2]", xlab=quote(italic(U[2])))

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.5.2 RC (2018-12-12 r75846)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Fedora 28 (Twenty Eight)
## 
## Matrix products: default
## BLAS: /sfs/u/maechler/R/D/r-patched/F28-64-inst/lib/libRblas.so
## LAPACK: /sfs/u/maechler/R/D/r-patched/F28-64-inst/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.0-6         rugarch_1.4-0      gsl_1.9-10.3      
## [10] mev_1.11           lattice_0.20-38    bbmle_1.0.20      
## [13] copula_0.999-19   
## 
## loaded via a namespace (and not attached):
##  [1] assertthat_0.2.0            expm_0.999-3               
##  [3] yaml_2.2.0                  GeneralizedHyperbolic_0.8-4
##  [5] numDeriv_2016.8-1           pillar_1.3.1               
##  [7] glue_1.3.0                  DistributionUtils_0.6-0    
##  [9] digest_0.6.18               colorspace_1.3-2           
## [11] sandwich_2.5-0              htmltools_0.3.6            
## [13] Matrix_1.2-15               plyr_1.8.4                 
## [15] pcaPP_1.9-73                pkgconfig_2.0.2            
## [17] ismev_1.42                  purrr_0.2.5                
## [19] mvtnorm_1.0-8               scales_1.0.0               
## [21] rootSolve_1.7               gmm_1.6-2                  
## [23] tibble_1.4.2                nleqslv_3.3.2              
## [25] ADGofTest_0.3               gmp_0.5-13.2               
## [27] mgcv_1.8-26                 bayesplot_1.6.0            
## [29] ggplot2_3.1.0               lazyeval_0.2.1             
## [31] magrittr_1.5                crayon_1.3.4               
## [33] mclust_5.4.2                evaluate_0.12              
## [35] ks_1.11.3                   nlme_3.1-137               
## [37] MASS_7.3-51.1               xts_0.11-2                 
## [39] truncnorm_1.0-8             pspline_1.0-18             
## [41] stringr_1.3.1               revdbayes_1.3.2            
## [43] munsell_0.5.0               stabledist_0.7-1           
## [45] bindrcpp_0.2.2              compiler_3.5.2             
## [47] SkewHyperbolic_0.4-0        evd_2.3-3                  
## [49] rlang_0.3.0.1               nloptr_1.2.1               
## [51] ggridges_0.5.1              Rsolnp_1.16                
## [53] Runuran_0.24                spd_2.0-1                  
## [55] partitions_1.9-19           rmarkdown_1.11             
## [57] boot_1.3-20                 gtable_0.2.0               
## [59] polynom_1.3-9               R6_2.3.0                   
## [61] zoo_1.8-4                   knitr_1.21                 
## [63] dplyr_0.7.7                 bindr_0.1.1                
## [65] KernSmooth_2.23-15          stringi_1.2.4              
## [67] Rcpp_1.0.0                  tidyselect_0.2.5           
## [69] xfun_0.4

Examples of Nonstandard Copulas – “Wild Animals”

Marius Hofert and Martin Mächler

2018-12-20

1 Swiss Alps copulas of Hofert, Vrins (2013)

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:

Define the copula \(C\)

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

2 An example from Wolfgang Trutschnig and Manuela Schreyer

Define the Iterated Function System

Run chaos game B times

3 Sierpinski tetrahedron

Session information