R Under development (unstable) (2023-12-04 r85659 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2023 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> ## Copyright (C) 2012 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun Yan
> ##
> ## This program is free software; you can redistribute it and/or modify it under
> ## the terms of the GNU General Public License as published by the Free Software
> ## Foundation; either version 3 of the License, or (at your option) any later
> ## version.
> ##
> ## This program is distributed in the hope that it will be useful, but WITHOUT
> ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
> ## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
> ## details.
> ##
> ## You should have received a copy of the GNU General Public License along with
> ## this program; if not, see .
>
>
> require(copula)
Loading required package: copula
> set.seed(1)
> n <- 2000
> N <- 100 # number of variates used for the Kolmogorov-Smirnov test
>
> ## maximal deviation for deciding if sample versions of Kendall's tau are
> ## "close enough" to population versions; NB: depends on 'n'
> eps.tau <- 0.06
>
> (doExtras <- copula:::doExtras())
[1] FALSE
> doPlots <- doExtras
> if(doPlots && !dev.interactive(orNone=TRUE)) pdf("nac-experi.pdf")
>
> options(warn = 2)# warning not allowed!
>
>
> ### 3d check functions #########################################################
>
> ##' correlation check function and run time measuring
> ##'
> ##' @title Check correlation matrix and measure run times of 3d fully nested Archimedean copulas
> ##' @param n number of variates to be drawn
> ##' @param th0 theta0
> ##' @param th1 theta1
> ##' @param cop acopula
> ##' @return a list containing run times for V0 and V01 and Kendall's taus
> ##' @author Marius Hofert, Martin Maechler
> corCheck <- function(n, th0,th1, cop) {
+ mat <- matrix(0,nrow = n,ncol = 3)
+ V0time <- system.time(V0 <- cop@V0(n,th0))
+ V01time <- system.time(V01 <- cop@V01(V0,th0,th1))
+ mat <- cbind(runif(n),
+ exp(-V0*cop@iPsi(cop@psi(rexp(n)/V01,th1),th0)),
+ exp(-V0*cop@iPsi(cop@psi(rexp(n)/V01,th1),th0)))
+ mat[,] <- cop@psi(-log(mat[,])/V0,th0)
+ list(V0time,V01time, name = cop@name, cor = cor(mat, method = "kendall"))
+ }
>
> ##' create output
> prt.tau.diff <- function(c1, c2){
+ stopifnot(is.matrix(c1))
+ delta.c <- 1000 * abs(c1 - c2)
+ cat(sprintf("Max & Mean distance * 1000 to true pairwise Kendall's taus: %7.1f %7.1f\n",
+ max(delta.c), mean(delta.c)))
+ invisible()
+ }
>
> ##' create output
> corCheckout <- function(x, trCorr, famName = x$name) {
+ cat(sprintf("Time [ms] V0 for '%s': %7.1f\n", famName, 1000*x[[1]][1]))
+ cat(sprintf("Time [ms] V01 for '%s': %7.1f\n", famName, 1000*x[[2]][1]))
+ prt.tau.diff(x[["cor"]], trCorr) ; cat("\n")
+ }
>
> ##' Function implementing the chi^2 test
> ##'
> ##' @title The Chi-square test
> ##' @param n [integer] sample size
> ##' @param N [integer] number of replications
> ##' @param cop outer_nacopula to generate from
> ##' @param nInt positive integer: the number of intervals used for each grid
> ##' dimension
> ##' @return an "chiSqChk_cop" object; just a list(...) and a print method
> ##' @author Marius Hofert, Martin Maechler
> chiSq_check_cop <- function(n,N,cop,nInt, verbose = interactive()){
+ copName <- deparse(substitute(cop)) # copula name
+ d <- dim(cop) # copula dimension
+ stopifnot(is.numeric(d), d >= 1, is.numeric(nInt), nInt >= 1)
+ pts <- (1:nInt)/nInt # (upper) division points [lower = upper - h
+ # = {0, ..., 1-h}; h=1/nInt]
+ mygrid <- do.call("expand.grid", rep.int(list(pts), d)) # build grid
+ m <- nInt^d # == grid length == nrow(mygrid)
+ v.cube <- nInt^(0:(d-1))
+
+ ## build a function that returns the number of the cube in which each row
+ ## of U falls
+ cube <- function(U, pts) {
+ di <- dim(U)
+ intervals <- array(cut(U, breaks = pts, include.lowest = TRUE,
+ labels = FALSE), # find "interval number" for
+ # each component of U; these numbers
+ # are in {NA,1,2,3,...,nInt-1}
+ dim = di)
+ intervals[is.na(intervals)] <- 0 # NAs correspond to smallest interval
+ as.vector(intervals %*% v.cube) + 1
+ }
+
+ ## determine the expected number of observations in each cube
+ prob_up <- function(u) {
+ ## probability mass in cube with upper corner 'u'
+ mesh <- 1/nInt
+ l <- u - mesh
+ prob(cop, l, u)
+ }
+ masscube <- apply(mygrid, 1, prob_up)
+ E_nobs <- n * masscube # expected number of observations in each cube
+
+ ## now simulate data, count observations in each cube, and compute test
+ ## statistic
+ k <- 0
+ CPU <- system.time({
+ T <- replicate(N, {
+ if(verbose) cat(sprintf("%2d%1s",{k <<- k+1}, if(k %% 20)""
+ else "\n"))
+ U <- rnacopula(n,cop) # generate data
+ cubenumbers <- cube(U,pts) # for each row vector of U, find the
+ # number of the cube in which the
+ # vector falls
+ nobs <- tabulate(cubenumbers, nbins = m) # number of observations
+ # in each cube
+ sum((nobs - E_nobs)^2 / E_nobs) # chi^2 test statistic
+ }); if(verbose) cat("done\n")
+ })[1]
+
+ structure(class = "chiSqChk_cop",
+ list( ## compute the result of the Kolmogorov-Smirnov test based
+ ## on the N realizations of the chi^2 test statistics:
+ ks = ks.test(T, "pchisq", df = m-1),
+ T = T, CPU = CPU,
+ n=n, N=N, copName = copName, m = m,
+ ## percentage of cubes that fulfill the rule of thumb:
+ percentrot = (sum(E_nobs >= 5)/m)*100
+ ))
+ }## end{chiSq_check_cop}
>
> ##' a print method for this class:
> print.chiSqChk_cop <- function(x, ...) {
+ stopifnot(is.list(x), all(c("ks","T","CPU","n","N") %in% names(x)),
+ is.numeric(pv <- x$ks[[2]]))
+ cat(sprintf("%s (3d)NAcopula (n=%d):\n %s (N=%d): %s\n ",
+ x$copName, x$n,
+ "P-value of the chi-square test", x$N, format.pval(pv)),
+ sprintf("Percentage fulfilling chi^2 rule of thumb: %4.1f\n",
+ x$percentrot),
+ sprintf("Time (user) needed = c(N,n; cop) = %8.1f [ms]\n",
+ 1000 * x$CPU), sep="")
+ if(pv < 0.05) {
+ if(pv < 0.01)
+ cat("\n*************** P-value < 0.01 <<<<<<<<<<<<<<<<<<<<<<<<\n",
+ "\n*************** ============== <<<<<<<<<<<<<<<<<<<<<<<<\n\n")
+ else cat("\n*** > > > P-value < 0.05 <<<<<<<<<<<<<<<<<<<<<<<<\n\n")
+ stopifnot(pv > 0.001)
+ }
+ invisible(x)
+ }
>
> ##' compute the probability to fall in a cube with
> ##' lower point l and upper point u for d=3
> probin3dcube <- function(cop,l,u) {
+ pCopula(u, cop)+
+ - pCopula(c(l[1],u[2],u[3]), cop)+
+ - pCopula(c(u[1],l[2],u[3]), cop)+
+ - pCopula(c(u[1],u[2],l[3]), cop)+
+ + pCopula(c(l[1],l[2],u[3]), cop)+
+ + pCopula(c(l[1],u[2],l[3]), cop)+
+ + pCopula(c(u[1],l[2],l[3]), cop)+
+ - pCopula(l, cop)
+ }
>
>
> ### 3d examples ################################################################
>
> ### AMH ########################################################################
>
> theta0 <- 0.7135 # tau_{12}=tau_{13}=0.2, tau_{23}=0.3
> theta1 <- 0.9430
>
> ## check 1
> corCheckAMH <- corCheck(n,theta0,theta1,copAMH)
> trCorr <- rbind(c(1,0.2,0.2),
+ c(0.2,1,0.3),
+ c(0.2,0.3,1))
> corCheckout(corCheckAMH,trCorr)
Time [ms] V0 for 'AMH': 0.0
Time [ms] V01 for 'AMH': 0.0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 6.5 2.2
> stopifnot(max(abs(corCheckAMH[["cor"]]-trCorr)) < eps.tau)
>
> ## check 2
> AMH3d <-
+ new("outer_nacopula", copula = setTheta(copAMH, theta0),
+ comp = 1L,
+ childCops = list(new("nacopula",
+ copula = setTheta(copAMH, theta1),
+ comp = 2:3)) # no childCops
+ )
>
> ## constructor forms of the above:
> rr <- onacopula("A", C(0.7135, 1, list(C(0.943, 2:3, NULL))))
> r0 <- onacopula("A", C(0.7135, 1, C(0.943, 2:3, NULL)))
> r1 <- onacopula("A", C(0.7135, 1, C(0.943, 2:3, )))
> r2 <- onacopula("AMH", C(0.7135, 1, C(0.943, 2:3 )))
> stopifnot(identical(AMH3d, rr), identical(rr, r0),
+ identical(r0, r1), identical(r1, r2))
>
> ## check
> (chkAMH <- chiSq_check_cop(n,N,AMH3d,5))
AMH3d (3d)NAcopula (n=2000):
P-value of the chi-square test (N=100): 0.61377
Percentage fulfilling chi^2 rule of thumb: 92.0
Time (user) needed = c(N,n; cop) = 170.0 [ms]
>
> ## check probability
> l <- c(.1, .05, .3)
> u <- c(.4, .7, .6)
> stopifnot(all.equal(print( prob(AMH3d,l,u)),
+ probin3dcube(AMH3d,l,u), tolerance=1e-14))
[1] 0.06694093
>
> ### Clayton ####################################################################
>
> theta0 <- 0.5 # tau_{12}=tau_{13}=0.2, tau_{23}=0.5
> theta1 <- 2
>
> ## check 1
> corCheckClayton <- corCheck(n,theta0,theta1,copClayton)
> trCorr <- rbind(c(1,0.2,0.2),
+ c(0.2,1,0.5),
+ c(0.2,0.5,1))
> corCheckout(corCheckClayton,trCorr)
Time [ms] V0 for 'Clayton': 0.0
Time [ms] V01 for 'Clayton': 0.0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 25.5 9.6
> stopifnot(max(abs(corCheckClayton[["cor"]]-trCorr)) < eps.tau)
>
> ## check 2
> Clayton3d <- onacopula("Clayton", C(theta0, 1, C(theta1, 2:3)))
> (chkClayton <- chiSq_check_cop(512,100,Clayton3d,5))
Clayton3d (3d)NAcopula (n=512):
P-value of the chi-square test (N=100): 0.35079
Percentage fulfilling chi^2 rule of thumb: 33.6
Time (user) needed = c(N,n; cop) = 230.0 [ms]
>
> ## check probability
> stopifnot(all.equal(print( prob(Clayton3d,l,u)),
+ probin3dcube(Clayton3d,l,u), tolerance=1e-14))
[1] 0.0709828
>
> ### Frank ######################################################################
>
> theta0 <- 1.8609 # tau_{12}=tau_{13}=0.2, tau_{23}=0.5
> theta1 <- 5.7363
>
> ## check 1
> corCheckFrank <- corCheck(n,theta0,theta1,copFrank)
> corCheckout(corCheckFrank,trCorr)
Time [ms] V0 for 'Frank': 0.0
Time [ms] V01 for 'Frank': 10.0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 3.5 1.5
> stopifnot(max(abs(corCheckFrank[["cor"]]-trCorr)) < eps.tau)
>
> ## check 2
> Frank3d <- onacopula("F", C(theta0, 1, C(theta1, 2:3)))
> (chkFrank <- chiSq_check_cop(n,N,Frank3d,5))
Frank3d (3d)NAcopula (n=2000):
P-value of the chi-square test (N=100): 0.25145
Percentage fulfilling chi^2 rule of thumb: 76.0
Time (user) needed = c(N,n; cop) = 720.0 [ms]
>
> ## check probability
> stopifnot(all.equal(print( prob(Frank3d,l,u)),
+ probin3dcube(Frank3d,l,u), tolerance=1e-14))
[1] 0.07750099
>
> ### Gumbel #####################################################################
>
> theta0 <- 1.25
> theta1 <- 2 #--> tau_{12}=tau_{13}=0.2, tau_{23}=0.5
> trCorr <- rbind(c(1,0.2,0.2),
+ c(0.2,1,0.5),
+ c(0.2,0.5,1))
> ## check 1
> corCheckGumbel <- corCheck(n,theta0,theta1,copGumbel)
> corCheckout(corCheckGumbel,trCorr)
Time [ms] V0 for 'Gumbel': 0.0
Time [ms] V01 for 'Gumbel': 0.0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 27.5 15.0
> stopifnot(max(abs(corCheckGumbel[["cor"]]-trCorr)) < eps.tau)
>
> ## check 2
> Gumbel3d <- onacopula("Gumbel", C(theta0, 1, C(theta1, 2:3)))
> (chkGumbel <- chiSq_check_cop(n,N,Gumbel3d,5))
Gumbel3d (3d)NAcopula (n=2000):
P-value of the chi-square test (N=100): 0.69828
Percentage fulfilling chi^2 rule of thumb: 80.8
Time (user) needed = c(N,n; cop) = 460.0 [ms]
>
> ## check probability
> stopifnot(all.equal(print( prob(Gumbel3d,l,u)),
+ probin3dcube(Gumbel3d,l,u), tolerance=1e-14))
[1] 0.08035166
>
> ### Joe ########################################################################
>
> theta0 <- 1.4438#tau_{12}=tau_{13}=0.2, tau_{23}=0.5
> theta1 <- 2.8562
>
> ## check 1
> corCheckJoe <- corCheck(n,theta0,theta1,copJoe)
> corCheckout(corCheckJoe,trCorr)
Time [ms] V0 for 'Joe': 0.0
Time [ms] V01 for 'Joe': 20.0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 41.2 15.3
> stopifnot(max(abs(corCheckJoe[["cor"]]-trCorr)) < eps.tau)
>
> ## check 2
> Joe3d <- onacopula("J", C(theta0, 1, C(theta1, 2:3)))
> (chkJoe <- chiSq_check_cop(n,N,Joe3d,5))
Joe3d (3d)NAcopula (n=2000):
P-value of the chi-square test (N=100): 0.63188
Percentage fulfilling chi^2 rule of thumb: 74.4
Time (user) needed = c(N,n; cop) = 700.0 [ms]
>
> ## check probability
> stopifnot(all.equal(print( prob(Joe3d,l,u)),
+ probin3dcube(Joe3d,l,u), tolerance=1e-14))
[1] 0.08367711
>
> ### Examples that check pnacopula() and rnacopula() ############################
>
> ## generate output for the examples
> prt.stats <- function(c1,c2, rt) {
+ cat("Time [ms] for generating", n,
+ "vectors of variates: ", round(1000*rt[1],1), "\n")
+ prt.tau.diff(c1, c2) ; cat("\n")
+ }
>
> ### 3d Ali-Mikhail-Haq copula example ##########################################
>
> c3 <- onacopula("A", C(0.7135, 1, list(C(0.943, 2:3))))
>
> ## basic check
> d <- dim(c3)
> stopifnot(d == 3,
+ allComp(c3) == 1:3,
+ allComp(c3@childCops[[1]]) == 2:3)
>
> ## test pCopula(., ) {was pnacopula()}
> u <- c(.3, .4, .5)
> ## with function:
> v <- pCopula(u, c3)
> ## by hand
> psi <- function(t,theta) { (1-theta)/(exp(t)-theta) }
> iPsi <- function(t,theta) { log((1-theta*(1-t))/t) }
> th0 <- 0.7135
> th1 <- 0.9430
> level1 <- psi(iPsi(u[2],th1) + iPsi(u[3],th1), th1)
> level0 <- psi(iPsi(u[1],th0) + iPsi(level1, th0), th0)
> stopifnot(all.equal(v, level0, tolerance = 1e-14))
>
> ## test rnacopula()
> rt <- system.time(rC3 <- rnacopula(n,c3))
> C3 <- cor(rC3,method = "kendall")
> trCorr <- rbind(c(1,0.2,0.2),
+ c(0.2,1,0.3),
+ c(0.2,0.3,1)) # tau_{12}=tau_{13}=0.2, tau_{23}=0.3
> stopifnot(is.numeric(rC3), is.matrix(rC3),
+ dim(rC3) == c(n, 3),max(abs(C3-trCorr)) < eps.tau)
> prt.stats(C3,trCorr,rt)
Time [ms] for generating 2000 vectors of variates: 0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 13.0 6.5
> if(doPlots) {
+ stopifnot(require("KernSmooth"))## for smoothScatter()
+ pairs2(rC3, panel = function(...) { par(new = TRUE); smoothScatter(...) })
+ }
>
> ### 2d Clayton copula example ##################################################
>
> c2 <- onacopula("Clayton", C(0.5, c(1,2))) # no childCops
> ## or simply c2 <- onacopula("Clayton", C(0.5, 1:2))
>
> ## basic check
> d <- dim(c2)
> stopifnot(d == 2,
+ allComp(c2) == 1:2)
>
> ## test pCopula()
> v <- pCopula(c(.3, .4), c2)
> stopifnot(all.equal(v,
+ local( { u1 <- .3; u2 <- .4
+ (u1^(-1/2)+u2^(-1/2)-1)^(-2) }),
+ tolerance = 1e-14))
>
> ## test rnacopula()
> racopula <- copula:::racopula
> set.seed(17) ; rt <- system.time(rC2 <- rnacopula(n,c2))
> set.seed(17) ; rt. <- system.time(rc2 <- racopula (n, c2@copula, d=2))
> stopifnot(identical(rC2, rc2))
>
> C2 <- cor(rC2,method = "kendall")
> trCorr <- rbind(c(1,0.2),
+ c(0.2,1)) # tau_{12}=0.2
> stopifnot(is.numeric(rC2), is.matrix(rC2),
+ dim(rC2) == c(n, 2), max(abs(C2-trCorr)) < eps.tau)
> prt.stats(C2,trCorr,rt)
Time [ms] for generating 2000 vectors of variates: 0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 2.1 1.1
> if(doPlots)
+ smoothScatter(rC2)
>
> ### 3d Clayton copula example ##################################################
>
> c3 <- onacopula("C", C(0.5, 1, C(2., c(2,3))))
>
> ## basic check
> d <- dim(c3)
> stopifnot(d == 3,
+ allComp(c3) == 1:3,
+ allComp(c3@childCops[[1]]) == 2:3)
>
> ## test pCopula()
> v <- pCopula(c(.3, .4, .5), c3)
> stopifnot(all.equal(v,
+ local( { u1 <- .3; u2 <- .4; u3 <- .5
+ 1/((1/u2^2 +1/u3^2 -1)^(1/4) -1 +1/sqrt(u1))^2 }),
+ tolerance = 1e-14))
>
> ## test rnacopula()
> rt <- system.time(rC3 <- rnacopula(n,c3))
> C3 <- cor(rC3,method = "kendall")
> trCorr <- matrix(c(1,0.2,0.2,0.2,1,0.5,0.2,0.5,1),nrow = 3,byrow = TRUE)
> # tau_{12}=tau_{13}=0.2, tau_{23}=0.5
> stopifnot(is.numeric(rC3), is.matrix(rC3),
+ dim(rC3) == c(n, 3),max(abs(C3-trCorr)) < eps.tau)
> prt.stats(C3,trCorr,rt)
Time [ms] for generating 2000 vectors of variates: 20
Max & Mean distance * 1000 to true pairwise Kendall's taus: 2.2 1.3
>
> if(doPlots)
+ pairs2(rC3, panel = function(...) { par(new = TRUE); smoothScatter(...) })
>
> ### 9d Clayton copula example ##################################################
>
> c9 <- onacopula("Clayton", C(0.5, c(3,6,1),
+ C(2., c(9,2,7,5),
+ C(3., c(8,4)))))
> c9Lis <- list(0.5, c(3,6,1),
+ list(list(2., c(9,2,7,5),
+ list(list(3., c(8,4))))))
> ## consistency onacopula() <-> onacopulaL() :
> stopifnot(identical(c9, onacopulaL("Clayton", c9Lis)))
>
>
> ## basic check
> d <- dim(c9)
> stopifnot(d == 9,
+ allComp(c9) == c(3,6,1,9,2,7,5,8,4),
+ allComp(c9@childCops[[1]]) == c(9,2,7,5,8,4),
+ allComp(c9@childCops[[1]]@childCops[[1]]) == c(8,4))
>
> ## test pCopula()
> u <- seq(0.1,0.9,by = 0.1)
> v <- pCopula(u, c9)
> ## by hand
> psi <- function(t,theta) { (1+t)^(-1/theta) }
> iPsi <- function(t,theta) { t^(-theta) - 1 }
> th0 <- 0.5
> th1 <- 2
> th2 <- 3
> level2 <- psi(iPsi(u[8],th2) + iPsi(u[4],th2), th2)
> level1 <- psi(iPsi(u[9],th1)+
+ iPsi(u[2],th1)+
+ iPsi(u[7],th1)+
+ iPsi(u[5],th1) +
+ iPsi(level2, th1), th1)
> level0 <- psi(iPsi(u[3],th0)+
+ iPsi(u[6],th0)+
+ iPsi(u[1],th0)+
+ iPsi(level1, th0), th0)
> stopifnot(all.equal(v, level0, tolerance = 1e-14))
>
> ## test rnacopula()
> rt <- system.time(rC9 <- rnacopula(n,c9))
> C9 <- cor(rC9,method = "kendall")
>
> ## Theoretical values:
> ## (11,12,13,14,15,16,17,18,19)=(1,0.2,0.2,0.2,0.2,0.2,0.2,0.2,0.2)
> ## (22,23,24,25,26,27,28,29)=(1,0.2,0.5,0.5,0.2,0.5,0.5,0.5)
> ## (33,34,35,36,37,38,39)=(1,0.2,0.2,0.2,0.2,0.2,0.2)
> ## (44,45,46,47,48,49)=(1,0.5,0.2,0.5,0.6,0.5)
> ## (55,56,57,58,59)=(1,0.2,0.5,0.5,0.5)
> ## (66,67,68,69)=(1,0.2,0.2,0.2)
> ## (77,78,79)=(1,0.5,0.5)
> ## (88,89)=(1,0.5)
>
> C9.true <- rbind(c(1. ,rep(0.2,8)),
+ c(0.2,1. ,0.2,0.5,0.5,0.2, rep(0.5,3)),
+ c(0.2,0.2,1. , rep(0.2,6)),
+ c(0.2,0.5,0.2,1. ,0.5,0.2,0.5,0.6,0.5),
+ c(0.2,0.5,0.2,0.5,1. ,0.2, rep(0.5,3)),
+ c(rep(0.2,5), 1. , rep(0.2,3)),
+ c(0.2,0.5,0.2,0.5,0.5,0.2,1. ,0.5,0.5),
+ c(0.2,0.5,0.2,0.6,0.5,0.2,0.5,1. ,0.5),
+ c(0.2,0.5,0.2,0.5,0.5,0.2,0.5,0.5,1. ))
> stopifnot(dim(rC9) == c(n, 9),
+ max(abs(C9-C9.true)) < eps.tau)
> prt.stats(C9,C9.true,rt)
Time [ms] for generating 2000 vectors of variates: 0
Max & Mean distance * 1000 to true pairwise Kendall's taus: 33.3 11.2
> if(doPlots && dev.interactive(orNone=TRUE)) # "large"
+ pairs2(rC9, gap = .1, pch = 20, cex = 0.2, col = rgb(.2,.1,.7, alpha = .5),
+ main = paste0(n," vectors of a ", d,"-dimensional nested Clayton copula"))
>
> ### 25d Clayton ==============================================
> c25 <- onacopula("Clayton", C(0.5, 17,
+ list(C(2, 20:18),
+ C(2.5, c(25,23, 8:12)),
+ C(2.25,c(24,21), C(4, 3:5)),
+ C(1.5, c(22,15:16), C(1.7, 1:2)),
+ C(3, c(6:7, 14:13)))))
> stopifnot(identical(sort(allComp(c25)), 1:25))
> c25
Nested Archimedean copula ("outer_nacopula" of dim. 25), with slot
'comp' = (17) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (0.5)
and 5 child copulas
1) Nested Archimedean copula ("nacopula"), with slot
'comp' = (20, 19, 18) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (2)
and *no* child copulas
2) Nested Archimedean copula ("nacopula"), with slot
'comp' = (25, 23, 8, 9, 10, 11, 12) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (2.5)
and *no* child copulas
3) Nested Archimedean copula ("nacopula"), with slot
'comp' = (24, 21) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (2.25)
and 1 child copula
Nested Archimedean copula ("nacopula"), with slot
'comp' = (3, 4, 5) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (4)
and *no* child copulas
4) Nested Archimedean copula ("nacopula"), with slot
'comp' = (22, 15, 16) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (1.5)
and 1 child copula
Nested Archimedean copula ("nacopula"), with slot
'comp' = (1, 2) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (1.7)
and *no* child copulas
5) Nested Archimedean copula ("nacopula"), with slot
'comp' = (6, 7, 14, 13) and root
'copula' = Archimedean copula ("acopula"), family "Clayton", theta= (3)
and *no* child copulas
> stopifnot(
+ all.equal(pCopula(rep(.01, 25), c25), 0.0001734511294041, tol = 9e-9)# 1.84e-13
+ ,
+ all.equal(pCopula(rep(.99, 25), c25), 0.79506048556858, tol = 9e-9)# 3.77e-15
+ )
>
> ### 125d Clayton copula example ################################################
>
> c125 <- onacopula("Clayton", C(0.5, , # no direct components
+ list(C(2, 1:10),
+ C(3, 11:40),
+ C(2, 41:60),
+ C(2, 61:85),
+ C(3, 86:105),
+ C(2,106:125))))
> c125Lis <- list(0.5, integer(0), # <- could use NULL
+ list(list(2, 1:10),
+ list(3, 11:40),
+ list(2, 41:60),
+ list(2, 61:85),
+ list(3, 86:105),
+ list(2,106:125)))
> ## consistency onacopula() <-> onacopulaL() :
> stopifnot(identical(c125, onacopulaL("C", c125Lis)))
>
>
>
> ## basic check
> d <- dim(c125)
> stopifnot(d == 125,
+ allComp(c125) == 1:125,
+ allComp(c125@childCops[[1]]) == 1:10,
+ allComp(c125@childCops[[2]]) == 11:40,
+ allComp(c125@childCops[[3]]) == 41:60,
+ allComp(c125@childCops[[4]]) == 61:85,
+ allComp(c125@childCops[[5]]) == 86:105,
+ allComp(c125@childCops[[6]]) == 106:125
+ )
>
> ## test rnacopula()
> rt <- system.time(rC125 <- rnacopula(n,c125))
> stopifnot(is.numeric(rC125), is.matrix(rC125), dim(rC125) == c(n, 125))
> cat("Time elapsed for generating ",n," vectors of variates:\n",sep = "")
Time elapsed for generating 2000 vectors of variates:
> rt
user system elapsed
0.08 0.00 0.08
> summary(p125 <- pCopula(rC125, c125))
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.100e-07 2.423e-03 5.434e-03 7.021e-03 9.965e-03 3.967e-02
> stopifnot(is.finite(p125),
+ all.equal(range(p125), c(4.13323986e-07, 0.03966520023))# tol=4.9e-11
+ ,
+ all.equal(quantile(p125, (1:3)/4, names=FALSE),
+ c(0.002422534035, 0.005433890506, 0.009965210908)) # tol=1.5e-11
+ )
>
>
> theta <- c(2,8)
> copG4 <- onacopulaL("Gumbel",
+ list(theta[1], NULL,
+ list(list(theta[2], c(1,2)),
+ list(theta[2], c(3,4)))))
> set.seed(11)
> uG4 <- rCopula(1000, copG4)
> u4. <- rbind(c(1, 0.5, 1, 0.5),
+ c(0.8, 0.4, 0.8, 0.5),
+ c(0.9, 0.5, 0.9, 0.5))
> (pu4 <- pCopula(u4., copula = copG4))
[1] 0.3752142 0.3169733 0.3752142
>
> stopifnot(
+ all.equal(
+ pu4, c(0.375214227246, 0.316973265328, 0.375214214144), tol = 9e-9)# 8e-13
+ ,
+ all.equal(print(
+ prob(copG4, l = rep(0.9, 4), u = rep(1, 4))
+ ), 0.056602230621, tol = 1e-9)
+ ,
+ ## and the probability of an empty corner is very small
+ all.equal(print(
+ prob(copG4,
+ l = c(0, 0.7, 0, 0.7),
+ u = c(0.4, 1 , 0.4,1 ))
+ , digits = 15), 8.832425e-10, tol = 1e-5)
+ )
[1] 0.05660223
[1] 8.83242534666095e-10
>
> ## Less "even" example:
> theta <- c(2, 4, 8, 25)
> copG8 <- onacopulaL("Gumbel",
+ list(theta[1], 5:6,
+ list(list(theta[2], c(1,2,7)),
+ list(theta[3], 3, list(list(theta[4], c(8,4)))))))
> uu <- rCopula(500, copG8)
> if(doPlots)
+ splom2(uu, col = adjustcolor("darkseagreen4", 0.5), cex = 0.25)
> summary(puu <- pCopula(uu, copG8))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000054 0.0305491 0.1466765 0.2474256 0.4116874 0.9984074
> stopifnot(is.finite(puu),
+ all.equal(min(puu), 5.353253447e-06, tol= 9e-9) # 2.33e-12
+ ,
+ all.equal(max(puu), 0.99840742713, tol= 9e-9) # 2.34e-12
+ )
>
>
> cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons''
Time elapsed: 9.21 0.2 9.42 NA NA
>
> proc.time()
user system elapsed
9.21 0.20 9.42