R Under development (unstable) (2026-02-16 r89426 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 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 . > > > ### (Nested) Archimedean Copulas ----------------------------------------------- > > require(copula) Loading required package: copula > > source(system.file("Rsource", "utils.R", package="copula", mustWork=TRUE)) Loading required package: tools > ## tryCatch.W.E(), showProc.time() etc > (doExtras <- copula:::doExtras()) [1] FALSE > > if(!dev.interactive(orNone=TRUE)) pdf("copula-play.pdf") > > ### testing psi > > myCop <- setTheta(copAMH, value = 0.5) # is maybe more natural > > ## Care: copula *does* define psi() already! > setGeneric("psi.", function(cop) standardGeneric("psi.")) [1] "psi." > setMethod(psi., "acopula", + function(cop) { function(t) cop@psi(t, theta = cop@theta) }) > psi.(myCop) # is a function function (t) cop@psi(t, theta = cop@theta) > psi.(myCop)(0:4) [1] 1.00000000 0.22539967 0.07257888 0.02552904 0.00924246 > curve(psi.(myCop)(x), 0, 4) > ##' but this can also be done directly [ => same curve "on top" :] > curve(myCop@psi(x, theta = myCop@theta), 0, 4, col = 2, add = TRUE) > showProc.time() Time (user system elapsed): 0.1 0 0.1 > > ### testing Kendall's tau > > p.Tau <- function(cop, n = 201, xlim = pmin(paraI, 50), ...) { + stopifnot(is(cop, "acopula")) + paraI <- cop@paraInterval + theta <- seq(xlim[1], xlim[2], length.out = n) + tit <- substitute(tau[NAME](theta), list(NAME = cop@name)) + plot(theta, cop@tau(theta), type = "l", main = tit, ...) + abline(h = c(0,1), lty = 3, col = "gray20") + } > > p.Tau(copAMH) > p.Tau(copClayton) > p.Tau(copFrank, xlim = c(0, 80), ylim= 0:1) # fast via debye_1() > p.Tau(copGumbel) > p.Tau(copJoe, ylim = 0:1, yaxs="i") > showProc.time() Time (user system elapsed): 0.04 0 0.05 > > ### test function ############################################################## > > ##' @title stopifnot() plus output > ##' @param expr > ##' @param prefix > ##' @param true > ##' @return > ##' @author Martin Maechler > checkifnot <- function(expr, prefix = "check if ", true = "[Ok]") + { + c0 <- function(...) cat(..., sep = "") + ## match.call(): not "calling" expr too early: + c0(prefix, deparse(match.call()[[2]])[1],": ") + stopifnot(expr) + c0(true,"\n") + } > > > ##' @title Perform a set of checks on a Archimedean copula object (with theta set) > ##' @param cop acopula > ##' @param theta1 parameter theta1 > ##' @param thetavec vector of parameters > ##' @param x values where psi is evaluated; x' := unique(round(x)) is used > ##' for dV0(), dV01() as some of these (AMH,Frank,Joe) are defined only on the natural numbers > ##' @param nRnd number of generated V0's and V01's > ##' @param nu2e n := 2^nu2e = #{u01} used by default > ##' @param u01 values where psiinv is evaluated > ##' @param lambdaLvec vector of lower tail-dependence coefficients > ##' @param lambdaUvec vector of upper tail-dependence coefficients > ##' @return list of measurements > ##' @author Marius Hofert, Martin Maechler > tstCop <- function(cop, theta1 = cop@theta, thetavec = cop@theta, + x = seq(1, 10, length.out = if(doExtras) 33 else 5), + nRnd = if(doExtras) 100 else 32, + nu2e = if(doExtras) 8 else 5, + u01 = local({n <- 2^nu2e; seq_len(n-1L)/n }),# exact binary fractions + K.MC.do = doExtras, + iTau.do = doExtras || cop@name != "Joe", + lambdaLvec = NA_real_, lambdaUvec = NA_real_) + { + stopifnot(is(cop, "acopula")) + cat0 <- function(...) cat(..., "\n", sep = "") + theta0 <- cop@theta + CT <- list() + + ### (1) cop name + + cat0(sprintf("(1) copula family: %10s, theta0 = %g", + cop@name, theta0)) + + ### (2) generator + + ### (2.1) psi and iPsi + + cat("\n(2) values of psi at x:\n") + CT[["psi"]] <- system.time( p.i <- cop@psi(x,theta = theta0)) + print(p.i) + checkifnot(identical(numeric(0), cop@iPsi(numeric(0), theta = theta0))) + checkifnot(cop@iPsi(0, theta = theta0) == Inf) + cat0("\nvalues of iPsi at u01:") + CT[["psiI"]] <- system.time( pi.t <- cop@iPsi(u01, theta = theta0)) + print(pi.t) + CT[["psiI"]] <- CT[["psiI"]] + + system.time(pi.pi <- cop@iPsi(p.i,theta = theta0)) + CT[["psi" ]] <- CT[["psi" ]] + + system.time(p.pit <- cop@psi(pi.t, theta = theta0)) + cat0("check if iPsi(psi(x))==x: ", all.equal(pi.pi, x)) + cat0("check if psi(iPsi(u01))==u01: ", all.equal(p.pit, u01)) + + ### (2.2) absdPsi + + ## absdPsi with degree = 10 + cat0("\nvalues of absdPsi with degree=10 at x:") + CT[["absdPsi"]] <- system.time( + p.D <- cop@absdPsi(x,theta = theta0, degree = 10)) + print(p.D) + cat0("check if all values are nonnegative") + stopifnot(is.vector(p.D), all(p.D >= 0)) + cat("check absdPsi(Inf,theta,degree=10) = 0 and the class of absdPsi(0,theta,degree=10): ") + at.0 <- cop@absdPsi(0, theta = theta0, degree = 10) + stopifnot(cop@absdPsi(Inf, theta = theta0, degree = 10) == 0, + is.numeric(at.0), !is.nan(at.0)) + cat0("[Ok]") + ## absdPsi with degree = 10 and MC + cat("\nvalues of absdPsi with degree=10 and MC at x:\n") + CT[["absdPsi"]] <- system.time( p.D <- cop@absdPsi(x,theta = theta0, degree = 10, n.MC = 1000)) + print(p.D) + cat0("check if all values are nonnegative") + stopifnot(all(p.D >= 0)) + cat("check absdPsi(Inf,theta,degree=10,n.MC=1000) = 0 and the class of absdPsi(0,theta,degree=10,n.MC=1000): ") + at.0 <- cop@absdPsi(0, theta = theta0, degree = 10, n.MC = 1000) + stopifnot(cop@absdPsi(Inf, theta = theta0, degree = 10, n.MC = 1000)==0, + is.numeric(at.0), !is.nan(at.0)) + cat0("[Ok]") + + ### (2.3) absdiPsi + + cat0("\nvalues of absdiPsi at u01:") + CT[["absdiPsi."]] <- system.time( absdiPsi. <- cop@absdiPsi(u01, theta = theta0)) + print(absdiPsi.) + stopifnot(all(absdiPsi. >= 0, is.numeric(absdiPsi.), !is.nan(absdiPsi.))) + cat("check the class of absdiPsi(0,theta): ") + at.0 <- cop@absdiPsi(0, theta = theta0) + stopifnot(is.numeric(at.0),!is.nan(at.0)) + cat0("[Ok]") + + ### (3) parameter interval + + cat("\n(3) parameter interval:\n") + print(cop@paraInterval) + cat0("theta1=",theta1) + cat0("nesting condition for theta0 and theta1 fulfilled: ", + cop@nestConstr(theta0,theta1)) + + ### (4) V0, dV0, V01, dV01 + + ## V0 + CT[["V0"]] <- system.time(V0 <- cop@V0(nRnd,theta0)) + cat0("\n(4) ",nRnd," generated V0's:") + print(summary(V0)) + ## dV0 + x. <- unique(round(x)) + cat("\nvalues of dV0 at x':\n") + CT[["dV0"]] <- system.time(dV0.i <- cop@dV0(x.,theta0)) + print(dV0.i) + ## V01 + CT[["V01"]] <- system.time(V01 <- cop@V01(V0,theta0,theta1)) + cat0("\n",nRnd," generated V01's:") + print(summary(V01)) + nt <- length(thetavec) + ## dV01 + cat("\nvalues of dV01 at x':\n") + CT[["dV01"]] <- system.time( dV01.i <- cop@dV01(x.,V0=1,theta0=theta0, theta1=theta1)) + print(dV01.i) + + ### (5) cCopula {was "cacopula"} + cat("\n(5) values of cCopula(cbind(v,rev(v)), copula = cop) for v=u01:\n") + cop. <- onacopulaL(cop@name, list(theta0, 1:2)) + CT[["cCopula."]] <- system.time( + cac <- cCopula(cbind(u01, rev(u01)), copula = cop., indices = 2)) + stopifnot(identical(dim(cac), c(length(u01),1L)), 0 <= cac, cac <= 1) + print(c(cac)) + + ### (6) dCopula (log = TRUE) {was dnacopula()} + + u <- matrix(runif(400),ncol=20) + ocop.2d <- onacopulaL(cop@name,list(theta0,1:2)) + ocop.20d <- onacopulaL(cop@name,list(theta0,1:20)) + + ## d = 2 + cat("\n(6) check dCopula(*, log = TRUE) for u being a random (20x2)-matrix:\n") + CT[["dCopula."]] <- system.time(lD <- dCopula(u[,1:2], ocop.2d, log = TRUE)) + print(lD); stopifnot(is.numeric(lD), is.finite(lD)); cat0("[Ok]") + cat("check at (0,0.5) and (1,0.5):\n") + stopifnot(dCopula(cbind(0:1,0.5), ocop.2d, log = FALSE) == 0, + dCopula(cbind(0:1,0.5), ocop.2d, log = TRUE ) == -Inf) + cat0("[Ok]") + + ## d = 20, n.MC = 0 + cat("\n check dCopula(*, log = TRUE) for u being a random (20x20)-matrix:\n") + CT[["dCopula.20"]] <- system.time(lD. <- dCopula(u, ocop.20d, log = TRUE)) + print(lD.); stopifnot(is.numeric(lD.), is.finite(lD.)); cat0("[Ok]") + + ## d = 20, n.MC > 0 + cat("\n check dCopula(*, log = TRUE) and MC for u being a random (20x20)-matrix:\n") + CT[["dCopula.MC"]] <- system.time(lD.. <- dCopula(u, ocop.20d, n.MC = 1000, log = TRUE)) + print(lD..); stopifnot(is.numeric(lD..), is.finite(lD..)); cat0("[Ok]") + + ## d = 20, check if n.MC > 0 is close to n.MC = 0 + stopifnot(all.equal(lD., lD.., tolerance=0.5)) + + ### (7) K + + check.K.u01 <- function(K){ + d.K <- diff(K) + if(any(neg <- d.K < 0)){ # happens for AMH, Clayton, and Frank (near 1) + if(any(Neg <- abs(d.K[neg]) > 1e-15* abs(K[-1][neg]))) { + warning("K(.) is 'substantially' non-monotone for K() / diff(K) =", + immediate.=TRUE) + print(cbind(K = K[-1][Neg], diff.K = d.K[Neg])) + } + } + stopifnot(is.numeric(K), length(K) == length(u01), 0 <= K, K <= 1) + } + + ## K for d = 2 + cat("\n(7) values of K for d = 2 at u01:\n") + CT[["K"]] <- system.time(K. <- pK(u01, cop, d = 2)) + check.K.u01( print(K.) ) + cat("check if K(0) = 0 and K(1) = 1: ") + stopifnot(pK(0, cop, d = 2)==0, + pK(1, cop, d = 2)==1) + cat0("[Ok]") + ## K for d = 10 + cat("\nvalues of K for d = 10 at u01:\n") + CT[["K.10"]] <- system.time(K. <- pK(u01, cop, d = 10)) + check.K.u01( print(K.) ) + cat("check if K(0) = 0 and K(1) = 1: ") + stopifnot(pK(0, cop, d = 10)==0, + pK(1, cop, d = 10)==1) + cat0("[Ok]") + if(K.MC.do) { + ## K for d = 10 and MC + cat("\nvalues of K for d = 10 and MC at u01:\n") + CT[["K.MC"]] <- system.time(K. <- pK(u01, cop, d = 10, n.MC = 1000)) + check.K.u01( print(K.) ) + cat("check if K(0)=0 and K(1)=1: ") + stopifnot(pK(0, cop, d = 10, n.MC = 1000)==0, + pK(1, cop, d = 10, n.MC = 1000)==1) + cat0("[Ok]") + } else + stopifnot(pK(0, cop, d = 10, n.MC = 2)==0, + pK(1, cop, d = 10, n.MC = 2)==1) + + ### (8) tau, iTau + + cat("\n(8) tau at thetavec:\n") + CT[["tau"]] <- system.time(ta <- cop@tau(thetavec)) + print(ta) + if(iTau.do) { + CT[["tauI"]] <- system.time(ta.I <- cop@iTau(ta)) + cat0("check if iTau(tau(thetavec))==thetavec: ", + all.equal(ta.I, thetavec)) + } + + ### (9) lambdaL, lambdaLInv + cat("\n(9) lambdaL at thetavec:\n") + lambdaLvec <- rep(as.double(lambdaLvec), length.out= nt) + lambdaUvec <- rep(as.double(lambdaUvec), length.out= nt) + CT[["lambdaL"]] <- system.time(lT <- cop@lambdaL(thetavec)) + CT[["lT.I"]] <- system.time(lT.I <- cop@lambdaLInv(lT)) + print(lT) + cat0("check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: ", + all.equal(lT.I, lambdaLvec)) + + ### (10) lambdaU, lambdaUInv + + cat("\n(10) lambdaU at thetavec:\n") + CT[["lambdaU"]] <- system.time(uT <- cop@lambdaU(thetavec)) + CT[["uT.I"]] <- system.time(uT.I <- cop@lambdaUInv(uT)) + print(uT) + cat0("check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: ", + all.equal(uT.I, lambdaUvec)) + + ### (11) dDiag + + cat("\n(11) dDiag at u01 for d=10:\n") + CT[["dDiag"]] <- system.time( dDiag. <- cop@dDiag(u01, theta=theta0, d=10)) + print(dDiag.) + stopifnot(is.numeric(dDiag.), all(dDiag. > 0)) + cat0("[Ok]") + + class(CT) <- "proc_time_list" + CT + } > > ##' print() method for the tstCop() results > print.proc_time_list <- function (x, ...) { + stopifnot(is.list(x), !is.null(nx <- names(x))) + cat("proc.time()s not all(. == 0): user system elapsed\n") + ## 2 4 6 8 0 2 4 6 8 0 2 4 6 89|1 3 |1 3 56|1 3 5 7 + ## 1 2 2 + for(nm in nx) + if(!all(x[[nm]] == 0, na.rm=TRUE)) { + ## use 'Time ..' as that works with 'R CMD Rdiff' + m <- 1000*x[[nm]] + cat(sprintf("Time [ms] for %13s :%5.0f %6.0f %7.0f\n", + ## 2 4 6 8 0 2 4 6 8 0| (20 + (13-4)) = 29 + nm, m[1], m[2], m[3])) + ## cat(nm,":\n"); print(x[[nm]], ...) + } + invisible(x) + } > showProc.time() Time (user system elapsed): 0.1 0 0.09 > > > ### copAMH ##################################################################### > > myAMH <- setTheta(copAMH, 0.7135001) > thetavec <- if(doExtras) c(0.1,0.3,0.5,0.7,0.9) else (1:3)/4 > set.seed(1) > tstCop(myAMH, 0.9429679, thetavec = thetavec) (1) copula family: AMH, theta0 = 0.7135 (2) values of psi at x: [1] 0.1429082757 0.0114248806 0.0011742837 0.0001234456 0.0000130075 check if identical(numeric(0), cop@iPsi(numeric(0), theta = theta0)): [Ok] check if cop@iPsi(0, theta = theta0) == Inf: [Ok] values of iPsi at u01: [1] 2.290664008 1.667234728 1.326942110 1.100443711 0.934955245 0.807145091 [7] 0.704687301 0.620307473 0.549372886 0.488763912 0.436288065 0.390351357 [13] 0.349762544 0.313610761 0.281185933 0.251925275 0.225376280 0.201170507 [19] 0.179004601 0.158626312 0.139824007 0.122418699 0.106257904 0.091210850 [25] 0.077164698 0.064021539 0.051695973 0.040113158 0.029207207 0.018919878 [31] 0.009199487 check if iPsi(psi(x))==x: TRUE check if psi(iPsi(u01))==u01: TRUE values of absdPsi with degree=10 at x: [1] 5.942950e+04 1.150526e+00 5.286052e-03 1.629375e-04 1.343933e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10) = 0 and the class of absdPsi(0,theta,degree=10): [Ok] values of absdPsi with degree=10 and MC at x: [1] 6.107456e+04 1.111224e+00 5.017170e-03 1.679631e-04 1.428779e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10,n.MC=1000) = 0 and the class of absdPsi(0,theta,degree=10,n.MC=1000): [Ok] values of absdiPsi at u01: [1] 29.6894186 13.8450203 8.6476536 6.1008145 4.6072154 3.6356601 [7] 2.9592833 2.4651783 2.0909796 1.7995212 1.5673340 1.3789055 [13] 1.2235953 1.0938781 0.9842930 0.8907887 0.8103020 0.7404780 [19] 0.6794798 0.6258553 0.5784436 0.5363068 0.4986795 0.4649322 [25] 0.4345425 0.4070742 0.3821604 0.3594905 0.3388004 0.3198638 [31] 0.3024863 check the class of absdiPsi(0,theta): [Ok] (3) parameter interval: 'interval' object [0, 1) theta1=0.9429679 nesting condition for theta0 and theta1 fulfilled: TRUE (4) 32 generated V0's: Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 1.000 2.000 3.844 4.500 14.000 values of dV0 at x': [1] 0.28649990 0.14585205 0.05297789 0.02697011 0.01373001 32 generated V01's: Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 4.50 11.50 17.72 26.25 69.00 values of dV01 at x': [1] 0.19906499 0.12769957 0.06561170 0.04208970 0.02700041 (5) values of cCopula(cbind(v,rev(v)), copula = cop) for v=u01: [1] 0.98943192 0.97555788 0.95829145 0.93758791 0.91344989 0.88593218 [7] 0.85514541 0.82125827 0.78449792 0.74514854 0.70354767 0.66008044 [13] 0.61517180 0.56927695 0.52287034 0.47643373 0.43044385 0.38536018 [19] 0.34161360 0.29959618 0.25965280 0.22207476 0.18709558 0.15488915 [25] 0.12556993 0.09919522 0.07576905 0.05524737 0.03754423 0.02253848 [31] 0.01008068 (6) check dCopula(*, log = TRUE) for u being a random (20x2)-matrix: [1] 0.057458992 0.007526097 -0.132093022 -0.281914129 -0.113486227 [6] -0.078639444 -0.514951239 -0.174515266 0.485884744 -0.106753078 [11] 0.386237282 -0.773142338 -1.020017913 0.357706750 -0.832024946 [16] 0.245168331 0.148011908 0.019298458 0.036164630 -0.932109004 [Ok] check at (0,0.5) and (1,0.5): [Ok] check dCopula(*, log = TRUE) for u being a random (20x20)-matrix: [1] -1.371756517 0.455069293 -3.625293357 -2.984592929 -4.827701594 [6] -2.118241601 -0.815362979 -1.036508247 -2.050972794 0.637533081 [11] -7.457395917 -1.350907810 -3.639305458 -2.738469551 -3.885599206 [16] 0.008515402 -3.643444468 -2.538137003 -0.798263354 -5.392475003 [Ok] check dCopula(*, log = TRUE) and MC for u being a random (20x20)-matrix: [1] -1.3160304 0.4069546 -3.5824883 -2.9364556 -4.8144941 -2.0754980 [7] -0.7575656 -1.0035307 -1.9934062 0.6680134 -7.5159141 -1.2961862 [13] -3.5916895 -2.6806642 -3.8360436 0.0589474 -3.5924578 -2.4805061 [19] -0.7464196 -5.3804553 [Ok] (7) values of K for d = 2 at u01: [1] 0.1084042 0.1829213 0.2471953 0.3053765 0.3591828 0.4095079 0.4568777 [8] 0.5016278 0.5439847 0.5841078 0.6221132 0.6580878 0.6920982 0.7241963 [15] 0.7544230 0.7828115 0.8093886 0.8341765 0.8571936 0.8784553 0.8979745 [22] 0.9157624 0.9318285 0.9461810 0.9588269 0.9697724 0.9790230 0.9865834 [29] 0.9924577 0.9966498 0.9991629 check if K(0) = 0 and K(1) = 1: [Ok] values of K for d = 10 at u01: [1] 0.6991605 0.8118068 0.8766878 0.9180775 0.9454747 0.9638678 0.9762553 [8] 0.9845716 0.9901130 0.9937657 0.9961411 0.9976608 0.9986151 0.9992016 [15] 0.9995532 0.9997583 0.9998742 0.9999373 0.9999703 0.9999867 0.9999945 [22] 0.9999979 0.9999993 0.9999998 0.9999999 1.0000000 1.0000000 1.0000000 [29] 1.0000000 1.0000000 1.0000000 check if K(0) = 0 and K(1) = 1: [Ok] (8) tau at thetavec: [1] 0.05942577 0.12876479 0.21379958 check if iTau(tau(thetavec))==thetavec: Mean relative difference: 2.838892e-05 (9) lambdaL at thetavec: [1] 0 0 0 check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: TRUE (10) lambdaU at thetavec: [1] 0 0 0 check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: TRUE (11) dDiag at u01 for d=10: [1] 9.582940e-09 2.278818e-06 4.277475e-05 2.906406e-04 1.148571e-03 [6] 3.254728e-03 7.386392e-03 1.433067e-02 2.478167e-02 3.929087e-02 [11] 5.826981e-02 8.203064e-02 1.108479e-01 1.450287e-01 1.849856e-01 [16] 2.313096e-01 2.848498e-01 3.468068e-01 4.188524e-01 5.032964e-01 [21] 6.033287e-01 7.233848e-01 8.697129e-01 1.051283e+00 1.281288e+00 [26] 1.579737e+00 1.978109e+00 2.528231e+00 3.320337e+00 4.523051e+00 [31] 6.481942e+00 [Ok] proc.time()s not all(. == 0): user system elapsed > showProc.time() Time (user system elapsed): 4.7 0.03 4.73 > > > ### copClayton ################################################################# > > myClayton <- setTheta(copClayton, 0.5) > thetavec <- c(0.5,1,2,5,10) > tstCop(myClayton, 2, thetavec, lambdaL = thetavec, lambdaU = NA) (1) copula family: Clayton, theta0 = 0.5 (2) values of psi at x: [1] 0.250000000 0.055363322 0.023668639 0.013061224 0.008264463 check if identical(numeric(0), cop@iPsi(numeric(0), theta = theta0)): [Ok] check if cop@iPsi(0, theta = theta0) == Inf: [Ok] values of iPsi at u01: [1] 4.65685425 3.00000000 2.26598632 1.82842712 1.52982213 1.30940108 [7] 1.13808994 1.00000000 0.88561808 0.78885438 0.70560573 0.63299316 [13] 0.56892908 0.51185789 0.46059349 0.41421356 0.37198868 0.33333333 [19] 0.29777137 0.26491106 0.23442680 0.20604538 0.17953565 0.15470054 [25] 0.13137085 0.10940039 0.08866211 0.06904497 0.05045146 0.03279556 [31] 0.01600102 check if iPsi(psi(x))==x: TRUE check if psi(iPsi(u01))==u01: TRUE values of absdPsi with degree=10 at x: [1] 9.745313e+03 1.149446e+00 7.017710e-03 1.981797e-04 1.271872e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10) = 0 and the class of absdPsi(0,theta,degree=10): [Ok] values of absdPsi with degree=10 and MC at x: [1] 9.716890e+03 1.149827e+00 6.841307e-03 2.025161e-04 1.342252e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10,n.MC=1000) = 0 and the class of absdPsi(0,theta,degree=10,n.MC=1000): [Ok] values of absdiPsi at u01: [1] 90.5096680 32.0000000 17.4185937 11.3137085 8.0954308 6.1584029 [7] 4.8870627 4.0000000 3.3522099 2.8621670 2.4808811 2.1773242 [13] 1.9309896 1.7278376 1.5579664 1.4142136 1.2912835 1.1851852 [19] 1.0928601 1.0119289 0.9405157 0.8771239 0.8205465 0.7698004 [25] 0.7240773 0.6827079 0.6451331 0.6108828 0.5795594 0.5508243 [31] 0.5243876 check the class of absdiPsi(0,theta): [Ok] (3) parameter interval: 'interval' object [0, Inf) theta1=2 nesting condition for theta0 and theta1 fulfilled: TRUE (4) 32 generated V0's: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.2694 1.1795 1.7857 1.9399 2.5498 6.0306 values of dV0 at x': [1] 0 0 0 0 0 32 generated V01's: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.0006681 0.0437420 0.5715796 0.6303525 0.7445805 4.2959946 values of dV01 at x': [1] 9.583385e-02 3.982071e-03 9.160650e-05 8.953843e-06 9.397679e-07 (5) values of cCopula(cbind(v,rev(v)), copula = cop) for v=u01: [1] 0.991561961 0.975801217 0.955053106 0.930201843 0.901820431 0.870356988 [7] 0.836199414 0.799703923 0.761209226 0.721043943 0.679530420 0.636986488 [13] 0.593725966 0.550058363 0.506288089 0.462713347 0.419624867 0.377304572 [19] 0.336024288 0.296044567 0.257613720 0.220967149 0.186327112 0.153903092 [25] 0.123893039 0.096485994 0.071866986 0.050226178 0.031776964 0.016796719 [31] 0.005744832 (6) check dCopula(*, log = TRUE) for u being a random (20x2)-matrix: [1] -0.13000598 -0.05715175 0.29304632 0.04117006 0.10277691 -0.45063791 [7] -0.05975298 0.12228646 -0.21684696 -0.18346948 -0.51725069 -0.12879465 [13] 0.07265274 0.09454809 0.05922776 0.10245875 -0.09087793 -0.00536712 [19] 0.11734222 0.13463918 [Ok] check at (0,0.5) and (1,0.5): [Ok] check dCopula(*, log = TRUE) for u being a random (20x20)-matrix: [1] -1.6077732 1.6980805 -1.2860673 -5.4629959 -6.0926099 -5.2652269 [7] 0.5358584 0.3355620 -2.9778706 -4.1016395 -4.6107174 -0.5256412 [13] -8.0453079 -6.5253387 -2.7055502 -5.2016623 -0.5566701 -3.4536002 [19] -2.0583889 -12.3972089 [Ok] check dCopula(*, log = TRUE) and MC for u being a random (20x20)-matrix: [1] -1.6788298 1.8057500 -1.3540724 -5.5333077 -5.9857996 -5.3051720 [7] 0.5029999 0.2760961 -3.0500764 -4.0646729 -4.4755737 -0.4727136 [13] -8.1030296 -6.5960767 -2.7736911 -5.2707964 -0.6265374 -3.5212752 [19] -2.1191137 -12.3798548 [Ok] (7) values of K for d = 2 at u01: [1] 0.08270146 0.15625000 0.22384008 0.28661165 0.34522353 0.40012024 [7] 0.45162811 0.50000000 0.54543933 0.58811438 0.62816740 0.66572067 [13] 0.70088083 0.73374190 0.76438763 0.79289322 0.81932670 0.84375000 [19] 0.86621980 0.88678823 0.90550348 0.92241023 0.93755008 0.95096189 [25] 0.96268207 0.97274479 0.98118227 0.98802489 0.99330141 0.99703906 [31] 0.99926372 check if K(0) = 0 and K(1) = 1: [Ok] values of K for d = 10 at u01: [1] 0.6043546 0.8029027 0.8950084 0.9421964 0.9676135 0.9816992 0.9896335 [8] 0.9941406 0.9967077 0.9981672 0.9989923 0.9994544 0.9997101 0.9998493 [15] 0.9999237 0.9999625 0.9999822 0.9999919 0.9999965 0.9999986 0.9999994 [22] 0.9999998 0.9999999 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 [29] 1.0000000 1.0000000 1.0000000 check if K(0) = 0 and K(1) = 1: [Ok] (8) tau at thetavec: [1] 0.2000000 0.3333333 0.5000000 0.7142857 0.8333333 check if iTau(tau(thetavec))==thetavec: TRUE (9) lambdaL at thetavec: [1] 0.2500000 0.5000000 0.7071068 0.8705506 0.9330330 check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: TRUE (10) lambdaU at thetavec: [1] 0 0 0 0 0 check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: TRUE (11) dDiag at u01 for d=10: [1] 0.01681765 0.02148300 0.02630307 0.03155188 0.03739808 0.04399409 [7] 0.05150167 0.06010518 0.07002211 0.08151417 0.09490046 0.11057408 [13] 0.12902376 0.15086247 0.17686583 0.20802446 0.24561610 0.29130633 [19] 0.34729138 0.41650372 0.50291294 0.61197470 0.75131422 0.93179101 [25] 1.16920060 1.48707363 1.92143313 2.52918840 3.40360165 4.70427207 [31] 6.71888252 [Ok] proc.time()s not all(. == 0): user system elapsed Time [ms] for dV01 : 10 0 20 > showProc.time() Time (user system elapsed): 4.66 0 4.66 > > ### copFrank ################################################################### > > myFrank <- setTheta(copFrank, 1.860884) > thetavec <- c(0.5,1,2,5,10) > set.seed(11) > tstCop(myFrank, 5.736283, thetavec) (1) copula family: Frank, theta0 = 1.86088 (2) values of psi at x: [1] 1.999171e-01 1.789017e-02 1.857775e-03 1.955055e-04 2.060278e-05 check if identical(numeric(0), cop@iPsi(numeric(0), theta = theta0)): [Ok] check if cop@iPsi(0, theta = theta0) == Inf: [Ok] values of iPsi at u01: [1] 2.70456759 2.04007407 1.66298100 1.40338960 1.20805568 1.05326313 [7] 0.92636137 0.81979941 0.72870702 0.64975918 0.58058456 0.51943262 [13] 0.46497428 0.41617674 0.37222160 0.33244941 0.29632111 0.26339046 [19] 0.23328400 0.20568606 0.18032753 0.15697717 0.13543491 0.11552652 [25] 0.09709940 0.08001921 0.06416711 0.04943749 0.03573616 0.02297875 [31] 0.01108947 check if iPsi(psi(x))==x: TRUE check if psi(iPsi(u01))==u01: TRUE values of absdPsi with degree=10 at x: [1] 4.089829e+04 8.929826e-01 5.586828e-03 2.323855e-04 2.100739e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10) = 0 and the class of absdPsi(0,theta,degree=10): [Ok] values of absdPsi with degree=10 and MC at x: [1] 4.096787e+04 8.834065e-01 5.672373e-03 2.293936e-04 2.049573e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10,n.MC=1000) = 0 and the class of absdPsi(0,theta,degree=10,n.MC=1000): [Ok] values of absdiPsi at u01: [1] 31.0785754 15.0875898 9.7632648 7.1055973 5.5145843 4.4568895 [7] 3.7039385 3.1414426 2.7059069 2.3592332 2.0771763 1.8435717 [13] 1.6472280 1.4801504 1.3364755 1.2118052 1.1027755 1.0067697 [19] 0.9217215 0.8459773 0.7781983 0.7172887 0.6623430 0.6126063 [25] 0.5674444 0.5263203 0.4887765 0.4544205 0.4229141 0.3939638 [31] 0.3673140 check the class of absdiPsi(0,theta): [Ok] (3) parameter interval: 'interval' object [0, Inf) theta1=5.736283 nesting condition for theta0 and theta1 fulfilled: TRUE (4) 32 generated V0's: Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 1.00 2.00 2.75 3.00 11.00 values of dV0 at x': [1] 0.453797724 0.107870897 0.032480233 0.017371752 0.009910529 32 generated V01's: Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 5.00 18.00 45.69 37.00 504.00 values of dV01 at x': [1] 0.38291599 0.07177938 0.02723906 0.01831041 0.01346045 (5) values of cCopula(cbind(v,rev(v)), copula = cop) for v=u01: [1] 0.98831917 0.97455270 0.95848158 0.93989161 0.91858286 0.89438127 [7] 0.86715206 0.83681451 0.80335693 0.76685077 0.72746192 0.68545756 [13] 0.64120684 0.59517408 0.54790405 0.50000000 0.45209595 0.40482592 [19] 0.35879316 0.31454244 0.27253808 0.23314923 0.19664307 0.16318549 [25] 0.13284794 0.10561873 0.08141714 0.06010839 0.04151842 0.02544730 [31] 0.01168083 (6) check dCopula(*, log = TRUE) for u being a random (20x2)-matrix: [1] -0.243036993 -0.415935379 -0.430864457 -0.055210284 -0.241565824 [6] -0.321304911 0.394487923 0.246448881 0.050966876 0.296279885 [11] 0.226874442 0.052940606 0.119935814 -0.596877372 0.112339966 [16] -0.709688707 -0.117845974 0.007842286 0.348835324 0.050478456 [Ok] check at (0,0.5) and (1,0.5): [Ok] check dCopula(*, log = TRUE) for u being a random (20x20)-matrix: [1] -2.4582183 -2.1719236 -1.6184754 -1.6047069 -0.9766182 -2.2701370 [7] 0.3786114 -2.2765230 -3.0882136 -2.6359217 -2.5487671 -1.1082366 [13] 1.1965984 1.9098068 -2.2099957 0.0186907 -2.7691931 -5.1827301 [19] -5.4947497 -1.8520876 [Ok] check dCopula(*, log = TRUE) and MC for u being a random (20x20)-matrix: [1] -2.42670519 -2.14049096 -1.58719298 -1.57644623 -0.95629958 -2.24348311 [7] 0.40356194 -2.24950940 -3.05676829 -2.60565548 -2.51792484 -1.09061400 [13] 1.22098694 1.93419264 -2.17867526 0.04375092 -2.74326532 -5.15373139 [19] -5.46517218 -1.82663199 [Ok] (7) values of K for d = 2 at u01: [1] 0.1182735 0.1977154 0.2640804 0.3225048 0.3753156 0.4238225 0.4688517 [8] 0.5109627 0.5505523 0.5879112 0.6232566 0.6567534 0.6885268 0.7186719 [15] 0.7472599 0.7743423 0.7999548 0.8241194 0.8468460 0.8681342 0.8879744 [22] 0.9063480 0.9232285 0.9385820 0.9523670 0.9645352 0.9750311 0.9837924 [29] 0.9907498 0.9958271 0.9989407 check if K(0) = 0 and K(1) = 1: [Ok] values of K for d = 10 at u01: [1] 0.7438465 0.8194291 0.8650654 0.8969999 0.9207681 0.9390109 0.9532327 [8] 0.9643939 0.9731594 0.9800177 0.9853443 0.9894375 0.9925401 0.9948526 [15] 0.9965418 0.9977469 0.9985829 0.9991444 0.9995074 0.9997318 0.9998634 [22] 0.9999357 0.9999726 0.9999896 0.9999966 0.9999991 0.9999998 1.0000000 [29] 1.0000000 1.0000000 1.0000000 check if K(0) = 0 and K(1) = 1: [Ok] (8) tau at thetavec: [1] 0.05541725 0.11001854 0.21389457 0.45670096 0.66577739 check if iTau(tau(thetavec))==thetavec: Mean relative difference: 4.060136e-07 (9) lambdaL at thetavec: [1] 0 0 0 0 0 check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: TRUE (10) lambdaU at thetavec: [1] 0 0 0 0 0 check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: TRUE (11) dDiag at u01 for d=10: [1] 2.532419e-10 4.520187e-06 3.663656e-05 1.641525e-04 5.302230e-04 [6] 1.389057e-03 3.141776e-03 6.366590e-03 1.183613e-02 2.051516e-02 [11] 3.353480e-02 5.214431e-02 7.764821e-02 1.113418e-01 1.544610e-01 [16] 2.081620e-01 2.735409e-01 3.516994e-01 4.438583e-01 5.515240e-01 [21] 6.767159e-01 8.222836e-01 9.923656e-01 1.193097e+00 1.433774e+00 [26] 1.728902e+00 2.102065e+00 2.593877e+00 3.280149e+00 4.319735e+00 [31] 6.109924e+00 [Ok] proc.time()s not all(. == 0): user system elapsed Time [ms] for absdPsi : 10 0 20 > > ## with a slightly more extensive test: > tau.th <- c(0.055417, 0.11002, 0.21389, 0.4567, 0.66578) > tau.F <- myFrank@tau(thetavec) > stopifnot(all.equal(tau.th, tau.F, tolerance = 0.0001), + all.equal(.9999, copFrank@tau(copFrank@iTau(0.9999))), + all.equal(myFrank@iTau(tau.F, tol = 1e-14), thetavec, tolerance=1e-11)) > showProc.time() Time (user system elapsed): 4.65 0 4.66 > > > ### copGumbel ################################################################## > > myGumbel <- setTheta(copGumbel, 1.25) > thetavec <- c(1,2,4,6,10) > (tG <- tstCop(myGumbel, 2, thetavec, lambdaL = NA, lambdaU = thetavec)) (1) copula family: Gumbel, theta0 = 1.25 (2) values of psi at x: [1] 0.367879441 0.076728831 0.020020049 0.005824610 0.001818809 check if identical(numeric(0), cop@iPsi(numeric(0), theta = theta0)): [Ok] check if cop@iPsi(0, theta = theta0) == Inf: [Ok] values of iPsi at u01: [1] 4.72872798 3.57772384 2.93613603 2.49708540 2.16675406 1.90408650 [7] 1.68749405 1.50424758 1.34622588 1.20793909 1.08550806 0.97609427 [13] 0.87756110 0.78826247 0.70690551 0.63245820 0.56408561 0.50110472 [19] 0.44295178 0.38915849 0.33933434 0.29315353 0.25034555 0.21068858 [25] 0.17400574 0.14016420 0.10907837 0.08071964 0.05513976 0.03252923 [31] 0.01340162 check if iPsi(psi(x))==x: TRUE check if psi(iPsi(u01))==u01: TRUE values of absdPsi with degree=10 at x: [1] 1.364542e+04 3.854817e-01 3.992157e-03 1.990155e-04 1.990823e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10) = 0 and the class of absdPsi(0,theta,degree=10): [Ok] values of absdPsi with degree=10 and MC at x: [1] 1.733441e+04 3.962929e-01 3.948077e-03 1.964560e-04 1.984815e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10,n.MC=1000) = 0 and the class of absdPsi(0,theta,degree=10,n.MC=1000): [Ok] values of absdiPsi at u01: [1] 54.5768992 25.8078222 16.5384183 12.0084424 9.3379579 7.5830876 [7] 6.3446899 5.4254263 4.7167306 4.1540240 3.6965273 3.3172416 [13] 2.9975892 2.7243702 2.4879463 2.2811108 2.0983604 1.9354109 [19] 1.7888655 1.6559808 1.5344978 1.4225133 1.3183786 1.2206101 [25] 1.1278016 1.0385181 0.9511389 0.8635705 0.7726005 0.6720375 [31] 0.5446654 check the class of absdiPsi(0,theta): [Ok] (3) parameter interval: 'interval' object [1, Inf) theta1=2 nesting condition for theta0 and theta1 fulfilled: TRUE (4) 32 generated V0's: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.3876 0.6391 0.9360 1.7803 1.7125 12.5269 values of dV0 at x': [1] 0.545626959 0.040690238 0.009349937 0.005235462 0.003369051 32 generated V01's: Min. 1st Qu. Median Mean 3rd Qu. Max. 0.04989 0.44900 1.22789 11.94095 4.43535 169.06411 values of dV01 at x': [1] 0.31100276 0.05246948 0.01630657 0.01004648 0.00690673 (5) values of cCopula(cbind(v,rev(v)), copula = cop) for v=u01: [1] 0.991613815 0.978280348 0.961541145 0.941796129 0.919234785 0.893971954 [7] 0.866096805 0.835697613 0.802877271 0.767764279 0.730520788 0.691348076 [13] 0.650489363 0.608229719 0.564892903 0.520835205 0.476436678 0.432090477 [19] 0.388191273 0.345123839 0.303252870 0.262914924 0.224413109 0.188014889 [25] 0.153953205 0.122431137 0.093630749 0.067727965 0.044918876 0.025474488 [31] 0.009896087 (6) check dCopula(*, log = TRUE) for u being a random (20x2)-matrix: [1] -0.419911441 0.239525751 -0.350253541 -0.240728003 -0.027099865 [6] 0.098240732 0.132792954 0.031455940 0.153252837 0.006595209 [11] 0.078647618 -0.034281283 0.084905924 0.272722792 -0.694730362 [16] -0.785492437 -0.357968954 -0.528276657 0.078233006 0.041626032 [Ok] check at (0,0.5) and (1,0.5): [Ok] check dCopula(*, log = TRUE) for u being a random (20x20)-matrix: [1] -2.39479668 -0.65194753 -2.80524656 -1.24603172 -1.40223396 -0.03769853 [7] -2.38535723 -0.03372653 -3.78012003 -0.15012613 -3.22033566 -0.52001117 [13] -1.19381742 -6.24458406 0.08410564 -2.46252483 1.16222274 -5.26999006 [19] 0.18438729 0.14840582 [Ok] check dCopula(*, log = TRUE) and MC for u being a random (20x20)-matrix: [1] -2.40278797 -0.65865088 -2.79158487 -1.25108501 -1.36683422 -0.03526235 [7] -2.39323689 -0.03009888 -3.78309869 -0.13548057 -3.18197961 -0.48609744 [13] -1.14826299 -6.20483179 0.09140940 -2.46874652 1.15512849 -5.27507410 [19] 0.22074981 0.19453530 [Ok] (7) values of K for d = 2 at u01: [1] 0.1178934 0.2011294 0.2712843 0.3329442 0.3882872 0.4385965 0.4847195 [8] 0.5272589 0.5666650 0.6032877 0.6374062 0.6692488 0.6990056 0.7268375 [15] 0.7528821 0.7772589 0.8000721 0.8214139 0.8413660 0.8600018 0.8773871 [22] 0.8935814 0.9086390 0.9226092 0.9355375 0.9474656 0.9584318 0.9684720 [29] 0.9776191 0.9859039 0.9933552 check if K(0) = 0 and K(1) = 1: [Ok] values of K for d = 10 at u01: [1] 0.8014447 0.8545899 0.8821793 0.9002562 0.9134618 0.9237453 0.9320984 [8] 0.9390902 0.9450752 0.9502884 0.9548931 0.9590068 0.9627167 0.9660895 [15] 0.9691768 0.9720198 0.9746513 0.9770983 0.9793830 0.9815239 0.9835368 [22] 0.9854348 0.9872293 0.9889303 0.9905461 0.9920842 0.9935512 0.9949528 [29] 0.9962942 0.9975799 0.9988139 check if K(0) = 0 and K(1) = 1: [Ok] (8) tau at thetavec: [1] 0.0000000 0.5000000 0.7500000 0.8333333 0.9000000 check if iTau(tau(thetavec))==thetavec: TRUE (9) lambdaL at thetavec: [1] 0 0 0 0 0 check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: TRUE (10) lambdaU at thetavec: [1] 0.0000000 0.5857864 0.8107929 0.8775380 0.9282265 check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: TRUE (11) dDiag at u01 for d=10: [1] 6.431254e-08 2.550566e-06 2.195867e-05 1.011527e-04 3.307716e-04 [6] 8.708570e-04 1.974266e-03 4.011606e-03 7.497498e-03 1.311804e-02 [11] 2.175938e-02 3.453725e-02 5.282745e-02 7.829727e-02 1.129376e-01 [16] 1.590960e-01 2.195101e-01 2.973426e-01 3.962158e-01 5.202475e-01 [21] 6.740872e-01 8.629536e-01 1.092671e+00 1.369710e+00 1.701221e+00 [26] 2.095080e+00 2.559922e+00 3.105185e+00 3.741151e+00 4.478983e+00 [31] 5.330772e+00 [Ok] proc.time()s not all(. == 0): user system elapsed Time [ms] for dV0 : 20 0 10 Time [ms] for dV01 : 20 0 10 Time [ms] for dCopula.20 : 10 0 20 Time [ms] for K.10 : 10 0 20 > u <- seq(0,1, length=32 + 1)[-c(1,32+1)] > u <- as.matrix(expand.grid(u,u)) > myGumbel@dacopula(u, theta=1.25) [1] 1.9359759 1.7558596 1.6352755 1.5413828 1.4630181 1.3949197 1.3341502 [8] 1.2788863 1.2279092 1.1803574 1.1355948 1.0931343 1.0525908 1.0136506 [15] 0.9760511 0.9395659 0.9039941 0.8691511 0.8348618 0.8009528 0.7672449 [22] 0.7335443 0.6996304 0.6652389 0.6300342 0.5935647 0.5551778 0.5138487 [29] 0.4677753 0.4131956 0.3392743 1.7558596 1.6340359 1.5480078 1.4784069 [36] 1.4184956 1.3650462 1.3162336 1.2709110 1.2283033 1.1878562 1.1491552 [43] 1.1118783 1.0757665 1.0406046 1.0062078 0.9724124 0.9390679 0.9060305 [50] 0.8731580 0.8403029 0.8073061 0.7739877 0.7401350 0.7054848 0.6696944 [57] 0.6322928 0.5925905 0.5494944 0.5010715 0.4432754 0.3644418 1.6352755 [64] 1.5480078 1.4836083 1.4298391 1.3823658 1.3390866 1.2988004 1.2607459 [71] 1.2244035 1.1893975 1.1554437 1.1223177 1.0898359 1.0578423 1.0261992 [78] 0.9947806 0.9634669 0.9321391 0.9006746 0.8689416 0.8367924 0.8040545 [85] 0.7705186 0.7359197 0.6999068 0.6619919 0.6214542 0.5771448 0.5270253 [92] 0.4668225 0.3842144 1.5413828 1.4784069 1.4298391 1.3880188 1.3501784 [99] 1.3149572 1.2815690 1.2495100 1.2184320 1.1880800 1.1582576 1.1288068 [106] 1.0995948 1.0705056 1.0414335 1.0122779 0.9829395 0.9533156 0.9232964 [113] 0.8927596 0.8615640 0.8295406 0.7964799 0.7621129 0.7260788 0.6878722 [120] 0.6467422 0.6014879 0.5499746 0.4877236 0.4018221 1.4630181 1.4184956 [127] 1.3823658 1.3501784 1.3202752 1.2918245 1.2643369 1.2374931 1.2110682 [134] 1.1848934 1.1588352 1.1327828 1.1066394 1.0803172 1.0537321 1.0268010 [141] 0.9994381 0.9715512 0.9430380 0.9137814 0.8836429 0.8524536 0.8200013 [148] 0.7860099 0.7501075 0.7117695 0.6702144 0.6241897 0.5714672 0.5073711 [155] 0.4184258 1.3949197 1.3650462 1.3390866 1.3149572 1.2918245 1.2692508 [162] 1.2469684 1.2247956 1.2025988 1.1802732 1.1577314 1.1348965 1.1116977 [169] 1.0880663 1.0639334 1.0392272 1.0138699 0.9877753 0.9608446 0.9329623 [176] 0.9039889 0.8737528 0.8420360 0.8085533 0.7729188 0.7345867 0.6927431 [183] 0.6460824 0.5922822 0.5264710 0.4346182 1.3341502 1.3162336 1.2988004 [190] 1.2815690 1.2643369 1.2469684 1.2293655 1.2114515 1.1931622 1.1744408 [197] 1.1552337 1.1354889 1.1151534 1.0941719 1.0724850 1.0500270 1.0267246 [204] 1.0024931 0.9772340 0.9508298 0.9231373 0.8939787 0.8631275 0.8302866 [211] 0.7950532 0.7568577 0.7148508 0.6676711 0.6128986 0.5454621 0.4507715 [218] 1.2788863 1.2709110 1.2607459 1.2495100 1.2374931 1.2247956 1.2114515 [225] 1.1974656 1.1828270 1.1675149 1.1515010 1.1347508 1.1172238 1.0988727 [232] 1.0796432 1.0594725 1.0382873 1.0160020 0.9925146 0.9677026 0.9414161 [239] 0.9134682 0.8836207 0.8515613 0.8168672 0.7789415 0.7368958 0.6893083 [246] 0.6336564 0.5646616 0.4671586 1.2279092 1.2283033 1.2244035 1.2184320 [253] 1.2110682 1.2025988 1.1931622 1.1828270 1.1716233 1.1595578 1.1466208 [260] 1.1327900 1.1180327 1.1023061 1.0855575 1.0677236 1.0487287 1.0284826 [267] 1.0068773 0.9837824 0.9590381 0.9324454 0.9037510 0.8726241 0.8386180 [274] 0.8011044 0.7591508 0.7112697 0.6548280 0.5843286 0.4840066 1.1803574 [281] 1.1878562 1.1893975 1.1880800 1.1848934 1.1802732 1.1744408 1.1675149 [288] 1.1595578 1.1505971 1.1406366 1.1296625 1.1176463 1.1045468 1.0903105 [295] 1.0748708 1.0581470 1.0400421 1.0204394 0.9991976 0.9761443 0.9510651 [302] 0.9236881 0.8936600 0.8605059 0.8235610 0.7818431 0.7337916 0.6766533 [309] 0.6046970 0.5015246 1.1355948 1.1491552 1.1554437 1.1582576 1.1588352 [316] 1.1577314 1.1552337 1.1515010 1.1466208 1.1406366 1.1335623 1.1253892 [323] 1.1160913 1.1056268 1.0939390 1.0809563 1.0665911 1.0507374 1.0332674 [330] 1.0140268 0.9928271 0.9694355 0.9435579 0.9148140 0.8826960 0.8464966 [337] 0.8051767 0.7570935 0.6993616 0.6259959 0.5199214 1.0931343 1.1118783 [344] 1.1223177 1.1288068 1.1327828 1.1348965 1.1354889 1.1347508 1.1327900 [351] 1.1296625 1.1253892 1.1199655 1.1133656 1.1055457 1.0964452 1.0859860 [358] 1.0740721 1.0605868 1.0453894 1.0283101 1.0091421 0.9876305 0.9634553 [365] 0.9362046 0.9053318 0.8700804 0.8293458 0.7813921 0.7231862 0.6484641 [372] 0.5394192 1.0525908 1.0757665 1.0898359 1.0995948 1.1066394 1.1116977 [379] 1.1151534 1.1172238 1.1180327 1.1176463 1.1160913 1.1133656 1.1094431 [386] 1.1042777 1.0978037 1.0899362 1.0805699 1.0695765 1.0568011 1.0420563 [393] 1.0251146 1.0056962 0.9834508 0.9579302 0.9285425 0.8944743 0.8545451 [400] 0.8069124 0.7483760 0.6723624 0.5602651 1.0136506 1.0406046 1.0578423 [407] 1.0705056 1.0803172 1.0880663 1.0941719 1.0988727 1.1023061 1.1045468 [414] 1.1056268 1.1055457 1.1042777 1.1017738 1.0979640 1.0927565 1.0860366 [421] 1.0776641 1.0674689 1.0552451 1.0407424 1.0236532 1.0035933 0.9800728 [428] 0.9524478 0.9198384 0.8809766 0.8338965 0.7752073 0.6979871 0.5827449 [435] 0.9760511 1.0062078 1.0261992 1.0414335 1.0537321 1.0639334 1.0724850 [442] 1.0796432 1.0855575 1.0903105 1.0939390 1.0964452 1.0978037 1.0979640 [449] 1.0968524 1.0943712 1.0903971 1.0847783 1.0773299 1.0678270 1.0559952 [456] 1.0414969 1.0239105 1.0026992 0.9771596 0.9463352 0.9088561 0.8626128 [463] 0.8039948 0.7256833 0.6071975 0.9395659 0.9724124 0.9947806 1.0122779 [470] 1.0268010 1.0392272 1.0500270 1.0594725 1.0677236 1.0748708 1.0809563 [477] 1.0859860 1.0899362 1.0927565 1.0943712 1.0946784 1.0935480 1.0908183 [484] 1.0862905 1.0797215 1.0708126 1.0591950 1.0444070 1.0258597 1.0027829 [491] 0.9741320 0.9384184 0.8933655 0.8351060 0.7558635 0.6340352 0.9039941 [498] 0.9390679 0.9634669 0.9829395 0.9994381 1.0138699 1.0267246 1.0382873 [505] 1.0487287 1.0581470 1.0665911 1.0740721 1.0805699 1.0860366 1.0903971 [512] 1.0935480 1.0953552 1.0956499 1.0942221 1.0908123 1.0850991 1.0766825 [519] 1.0650591 1.0495843 1.0294134 1.0034017 0.9699232 0.9265047 0.8689771 [526] 0.7890294 0.6637696 0.8691511 0.9060305 0.9321391 0.9533156 0.9715512 [533] 0.9877753 1.0024931 1.0160020 1.0284826 1.0400421 1.0507374 1.0605868 [540] 1.0695765 1.0776641 1.0847783 1.0908183 1.0956499 1.0991014 1.1009558 [547] 1.1009412 1.0987168 1.0938540 1.0858079 1.0738761 1.0571333 1.0343225 [554] 1.0036590 0.9624389 0.9061350 0.8258031 0.6970482 0.8348618 0.8731580 [561] 0.9006746 0.9232964 0.9430380 0.9608446 0.9772340 0.9925146 1.0068773 [568] 1.0204394 1.0332674 1.0453894 1.0568011 1.0674689 1.0773299 1.0862905 [575] 1.0942221 1.1009558 1.1062738 1.1098983 1.1114752 1.1105523 1.1065473 [582] 1.0987009 1.0860032 1.0670739 1.0399476 1.0016494 0.9472243 0.8669687 [589] 0.7347066 0.8009528 0.8403029 0.8689416 0.8927596 0.9137814 0.9329623 [596] 0.9508298 0.9677026 0.9837824 0.9991976 1.0140268 1.0283101 1.0420563 [603] 1.0552451 1.0678270 1.0797215 1.0908123 1.1009412 1.1098983 1.1174090 [610] 1.1231157 1.1265522 1.1271072 1.1239700 1.1160475 1.1018285 1.0791452 [617] 1.0447067 0.9930441 0.9135320 0.7778456 0.7672449 0.8073061 0.8367924 [624] 0.8615640 0.8836429 0.9039889 0.9231373 0.9414161 0.9590381 0.9761443 [631] 0.9928271 1.0091421 1.0251146 1.0407424 1.0559952 1.0708126 1.0850991 [638] 1.0987168 1.1114752 1.1231157 1.1332914 1.1415372 1.1472281 1.1495155 [645] 1.1472324 1.1387365 1.1216403 1.0922891 1.0445975 0.9668057 0.8279475 [652] 0.7335443 0.7739877 0.8040545 0.8295406 0.8524536 0.8737528 0.8939787 [659] 0.9134682 0.9324454 0.9510651 0.9694355 0.9876305 1.0056962 1.0236532 [666] 1.0414969 1.0591950 1.0766825 1.0938540 1.1105523 1.1265522 1.1415372 [673] 1.1550662 1.1665246 1.1750516 1.1794283 1.1778957 1.1678412 1.1452019 [680] 1.1031566 1.0285367 0.8870595 0.6996304 0.7401350 0.7705186 0.7964799 [687] 0.8200013 0.8420360 0.8631275 0.8836207 0.9037510 0.9236881 0.9435579 [694] 0.9634553 0.9834508 1.0035933 1.0239105 1.0444070 1.0650591 1.0858079 [701] 1.1065473 1.1271072 1.1472281 1.1665246 1.1844304 1.2001145 1.2123498 [708] 1.2192993 1.2181454 1.2043912 1.1703511 1.1010968 0.9580888 0.6652389 [715] 0.7054848 0.7359197 0.7621129 0.7860099 0.8085533 0.8302866 0.8515613 [722] 0.8726241 0.8936600 0.9148140 0.9362046 0.9579302 0.9800728 1.0026992 [729] 1.0258597 1.0495843 1.0738761 1.0987009 1.1239700 1.1495155 1.1750516 [736] 1.2001145 1.2239687 1.2454556 1.2627410 1.2728719 1.2709385 1.2482816 [743] 1.1877791 1.0453040 0.6300342 0.6696944 0.6999068 0.7260788 0.7501075 [750] 0.7729188 0.7950532 0.8168672 0.8386180 0.8605059 0.8826960 0.9053318 [757] 0.9285425 0.9524478 0.9771596 1.0027829 1.0294134 1.0571333 1.0860032 [764] 1.1160475 1.1472324 1.1794283 1.2123498 1.2454556 1.2777803 1.3076424 [771] 1.3321153 1.3460047 1.3396548 1.2932646 1.1552216 0.5935647 0.6322928 [778] 0.6619919 0.6878722 0.7117695 0.7345867 0.7568577 0.7789415 0.8011044 [785] 0.8235610 0.8464966 0.8700804 0.8944743 0.9198384 0.9463352 0.9741320 [792] 1.0034017 1.0343225 1.0670739 1.1018285 1.1387365 1.1778957 1.2192993 [799] 1.2627410 1.3076424 1.3527340 1.3954407 1.4306384 1.4479044 1.4243713 [806] 1.2982711 0.5551778 0.5925905 0.6214542 0.6467422 0.6702144 0.6927431 [813] 0.7148508 0.7368958 0.7591508 0.7818431 0.8051767 0.8293458 0.8545451 [820] 0.8809766 0.9088561 0.9384184 0.9699232 1.0036590 1.0399476 1.0791452 [827] 1.1216403 1.1678412 1.2181454 1.2728719 1.3321153 1.3954407 1.4612309 [834] 1.5252395 1.5771493 1.5912441 1.4921389 0.5138487 0.5494944 0.5771448 [841] 0.6014879 0.6241897 0.6460824 0.6676711 0.6893083 0.7112697 0.7337916 [848] 0.7570935 0.7813921 0.8069124 0.8338965 0.8626128 0.8933655 0.9265047 [855] 0.9624389 1.0016494 1.0447067 1.0922891 1.1452019 1.2043912 1.2709385 [862] 1.3460047 1.4306384 1.5252395 1.6281128 1.7314522 1.8091008 1.7690582 [869] 0.4677753 0.5010715 0.5270253 0.5499746 0.5714672 0.5922822 0.6128986 [876] 0.6336564 0.6548280 0.6766533 0.6993616 0.7231862 0.7483760 0.7752073 [883] 0.8039948 0.8351060 0.8689771 0.9061350 0.9472243 0.9930441 1.0445975 [890] 1.1031566 1.1703511 1.2482816 1.3396548 1.4479044 1.5771493 1.7314522 [897] 1.9113886 2.0997511 2.1932656 0.4131956 0.4432754 0.4668225 0.4877236 [904] 0.5073711 0.5264710 0.5454621 0.5646616 0.5843286 0.6046970 0.6259959 [911] 0.6484641 0.6723624 0.6979871 0.7256833 0.7558635 0.7890294 0.8258031 [918] 0.8669687 0.9135320 0.9668057 1.0285367 1.1010968 1.1877791 1.2932646 [925] 1.4243713 1.5912441 1.8091008 2.0997511 2.4851439 2.9063793 0.3392743 [932] 0.3644418 0.3842144 0.4018221 0.4184258 0.4346182 0.4507715 0.4671586 [939] 0.4840066 0.5015246 0.5199214 0.5394192 0.5602651 0.5827449 0.6071975 [946] 0.6340352 0.6637696 0.6970482 0.7347066 0.7778456 0.8279475 0.8870595 [953] 0.9580888 1.0453040 1.1552216 1.2982711 1.4921389 1.7690582 2.1932656 [960] 2.9063793 4.2199074 > showProc.time() Time (user system elapsed): 3.63 0.02 3.65 > > > ### copJoe ##################################################################### > > myJoe <- setTheta(copJoe, 1.25) > thetavec <- c(1.1,2,4,6,10) > set.seed(111) > tstCop(myJoe, 2, thetavec, lambdaL = NA, lambdaU = thetavec) (1) copula family: Joe, theta0 = 1.25 (2) values of psi at x: [1] 3.071489e-01 3.114155e-02 3.270755e-03 3.446089e-04 3.632011e-05 check if identical(numeric(0), cop@iPsi(numeric(0), theta = theta0)): [Ok] check if cop@iPsi(0, theta = theta0) == Inf: [Ok] values of iPsi at u01: [1] 3.24653731 2.55741487 2.15605800 1.87257167 1.65371570 1.47577825 [7] 1.32611329 1.19717478 1.08409783 0.98356325 0.89320613 0.81128299 [13] 0.73647261 0.66775099 0.60430961 0.54550025 0.49079661 0.43976719 [19] 0.39205557 0.34736606 0.30545319 0.26611401 0.22918280 0.19452778 [25] 0.16204986 0.13168389 0.10340353 0.07723256 0.05326968 0.03174870 [31] 0.01322609 check if iPsi(psi(x))==x: TRUE check if psi(iPsi(u01))==u01: TRUE values of absdPsi with degree=10 at x: [1] 1.187561e+04 3.175536e-01 4.771870e-03 3.599450e-04 3.648897e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10) = 0 and the class of absdPsi(0,theta,degree=10): [Ok] values of absdPsi with degree=10 and MC at x: [1] 1.050962e+04 2.928855e-01 4.786327e-03 3.659237e-04 3.712451e-05 check if all values are nonnegative check absdPsi(Inf,theta,degree=10,n.MC=1000) = 0 and the class of absdPsi(0,theta,degree=10,n.MC=1000): [Ok] values of absdiPsi at u01: [1] 31.8725039 15.8698930 10.5338253 7.8642909 6.2612791 5.1914444 [7] 4.4262017 3.8512502 3.4030802 3.0435769 2.7484753 2.5015878 [13] 2.2916955 2.1107715 1.9529138 1.8136783 1.6896455 1.5781308 [19] 1.4769845 1.3844491 1.2990544 1.2195353 1.1447640 1.0736862 [25] 1.0052507 0.9383177 0.8715105 0.8029328 0.7295206 0.6451613 [31] 0.5325575 check the class of absdiPsi(0,theta): [Ok] (3) parameter interval: 'interval' object [1, Inf) theta1=2 nesting condition for theta0 and theta1 fulfilled: TRUE (4) 32 generated V0's: Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 1.000 1.000 2.062 2.000 14.000 values of dV0 at x': [1] 0.800000000 0.032000000 0.007884800 0.004539392 0.002977841 32 generated V01's: Min. 1st Qu. Median Mean 3rd Qu. Max. 1.00 1.00 2.00 5.00 3.25 32.00 values of dV01 at x': [1] 0.625000000 0.053710938 0.015696287 0.009604341 0.006591312 (5) values of cCopula(cbind(v,rev(v)), copula = cop) for v=u01: [1] 0.98675604 0.96824211 0.94684145 0.92319502 0.89762841 0.87034455 [7] 0.84148825 0.81117526 0.77950804 0.74658520 0.71250782 0.67738351 [13] 0.64132909 0.60447203 0.56695083 0.52891451 0.49052134 0.45193709 [19] 0.41333288 0.37488311 0.33676358 0.29915031 0.26221914 0.22614674 [25] 0.19111345 0.15730877 0.12494133 0.09425688 0.06557355 0.03936210 [31] 0.01648770 (6) check dCopula(*, log = TRUE) for u being a random (20x2)-matrix: [1] -0.087170867 -0.356962044 -0.019532034 0.027849194 -0.005604182 [6] 0.075085025 0.086059567 0.118102669 0.294022133 0.061139784 [11] 0.059035208 0.099386283 -0.074025045 -0.126190880 -1.179690222 [16] 0.019201639 -0.005266644 0.048533283 -0.323104991 -0.097485846 [Ok] check at (0,0.5) and (1,0.5): [Ok] check dCopula(*, log = TRUE) for u being a random (20x20)-matrix: [1] 0.1674861 -0.9956733 -1.8857311 1.9883374 -0.4397215 -1.4050733 [7] 0.6207870 0.2708725 0.1737061 -0.9132191 -0.5931607 1.0625634 [13] 0.0480284 -0.2182888 -1.9926353 -1.2343581 -1.3374726 -0.5670577 [19] -0.5561139 0.0219886 [Ok] check dCopula(*, log = TRUE) and MC for u being a random (20x20)-matrix: [1] 0.18229947 -0.98282499 -1.87523517 2.00322536 -0.42564185 -1.45725612 [7] 0.63565692 0.28575796 0.18063440 -0.90379770 -0.65329582 1.07745199 [13] 0.06290346 -0.20340237 -1.97783879 -1.25222578 -1.32262383 -0.55326846 [19] -0.62856870 0.03669320 [Ok] (7) values of K for d = 2 at u01: [1] 0.1331101 0.2236488 0.2984295 0.3631107 0.4203679 0.4717712 0.5183553 [8] 0.5608535 0.5998137 0.6356603 0.6687324 0.6993072 0.7276158 0.7538540 [15] 0.7781900 0.8007701 0.8217231 0.8411633 0.8591933 0.8759056 0.8913850 [22] 0.9057094 0.9189509 0.9311775 0.9424534 0.9528404 0.9623986 0.9711881 [29] 0.9792701 0.9867105 0.9935850 check if K(0) = 0 and K(1) = 1: [Ok] values of K for d = 10 at u01: [1] 0.8878125 0.9112200 0.9242306 0.9332805 0.9402354 0.9458939 0.9506713 [8] 0.9548112 0.9584686 0.9617485 0.9647250 0.9674525 0.9699722 0.9723157 [15] 0.9745084 0.9765704 0.9785182 0.9803656 0.9821239 0.9838029 0.9854109 [22] 0.9869552 0.9884419 0.9898766 0.9912643 0.9926093 0.9939158 0.9951876 [29] 0.9964281 0.9976410 0.9988300 check if K(0) = 0 and K(1) = 1: [Ok] (8) tau at thetavec: [1] 0.05439832 0.35506593 0.61370564 0.72259094 0.82204394 (9) lambdaL at thetavec: [1] 0 0 0 0 0 check if lambdaLInv(lambdaL(thetavec))==lambdaLvec: TRUE (10) lambdaU at thetavec: [1] 0.1221382 0.5857864 0.8107929 0.8775380 0.9282265 check if lambdaUInv(lambdaU(thetavec))==lambdaUvec: TRUE (11) dDiag at u01 for d=10: [1] 2.027560e-12 9.930070e-10 3.647835e-08 4.637395e-07 3.294254e-06 [6] 1.618657e-05 6.164311e-05 1.947279e-04 5.330731e-04 1.302921e-03 [11] 2.904346e-03 5.997559e-03 1.161026e-02 2.126351e-02 3.711045e-02 [16] 6.207902e-02 1.000062e-01 1.557485e-01 2.352502e-01 3.455500e-01 [21] 4.947042e-01 6.916033e-01 9.456570e-01 1.266317e+00 1.662395e+00 [26] 2.141113e+00 2.706771e+00 3.358814e+00 4.088786e+00 4.874769e+00 [31] 5.667751e+00 [Ok] proc.time()s not all(. == 0): user system elapsed > showProc.time() Time (user system elapsed): 3.54 0.02 3.56 > > > ### Regression tests ------------------------------------ > > chkPsi <- function(copula, t = c(0, 2^c(-1000,-500, -200,-10*(10:0)), 2:3, 2^(2:40),Inf)) { + stopifnot(is(copula, "Copula")) + if(is.unsorted(t)) t <- sort(t) + psf <- psi(copula, t) + ## and also an equidistant t --> to check convexity + ps.eq <- psi(copula, t. <- seq(0, 20, length=1+2^7)) + stopifnot(is.finite(psf), 0 <= psf, psf <= 1, + psf[1] == 1, diff(psf) <= 0, + is.na (pN <- psi(copula, c(NA, NaN))), + is.nan(pN[2]), + 0 <= ps.eq, ps.eq <= 1, diff(ps.eq) <= 0, + ## convexity (in light of finite accuracy arithmetic): + diff(ps.eq, diff=2) >= - 4*.Machine$double.eps *ps.eq[-(1:2)] + ) + + ## for plotting: + it <- sort.list(tt <- c(t,t.)) + invisible(list(x=tt[it], y= c(psf, ps.eq)[it])) + } > > ### Negative tau (and dim = 2): > > taus <- c(-1,0,1); names(taus) <- paste0("tau=",taus) > taus tau=-1 tau=0 tau=1 -1 0 1 > > ## Frank: -------------------------------------------------------- > vapply(taus, function(tau) iTau(frankCopula(), tau), 1.) tau=-1 tau=0 tau=1 -1.81e+16 0.00e+00 7.21e+16 > ## tau=-1 tau=0 tau=1 > ## -1.81e+16 0.00e+00 7.21e+16 > ## ~= - Inf 0 + Inf > > r <- chkPsi(frankCopula(-2)) > plot(r, type="o") > plot(r, type="o", log="xy") Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log) : 2 x values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 32 y values <= 0 omitted from logarithmic plot > chkPsi(frankCopula( -800))# failed before 2014-06 > chkPsi(frankCopula(-2000))# (ditto) > chkPsi(frankCopula(-1e10))# (ditto) > showProc.time() Time (user system elapsed): 0.05 0 0.05 > > > ## Clayton: ------------------------------------------------------ > > vapply(taus, function(tau) iTau(claytonCopula(), tau), 1.) tau=-1 tau=0 tau=1 -1 0 Inf > ## tau=-1 tau=0 tau=1 > ## -1 0 Inf > > stopifnot(all.equal(-2/3, iTau(claytonCopula(), -1/2))) > > tools::assertError(chkPsi(claytonCopula(-1.1))) # par. out of bound > chkPsi(claytonCopula(-1)) ## all failed before 2014-05 > chkPsi(claytonCopula(-.5)) > chkPsi(claytonCopula(-1/8)) > chkPsi(claytonCopula(-2^-10)) > showProc.time() Time (user system elapsed): 0 0 0 > > ## AMH: > tAMH <- c((5 - 8*log(2))/ 3, -1/8, 0, 1/8, 1/3) > (th.t <- vapply(tAMH, function(tau) iTau(amhCopula(), tau), 1.)) [1] -1.000000e+00 -6.463595e-01 1.368215e-13 4.875224e-01 1.000000e+00 > stopifnot(-1 <= th.t, th.t <= 1, + all.equal(th.t[c(1,3,5)], c(-1,0,1))) > ## rho: --> ../vignettes/rhoAMH-dilog.Rnw > > ## cCopula() for all three "negative" tau families: > ## -------- -------------- > cCneg <- function(tau, u1 = (1:8)/8) { + stopifnot(length(tau) == 1, is.finite(tau), -1 <= tau, tau <= 1) + u <- cbind(u1, .5) + rbind(A = cCopula(u, amhCopula(iTau( amhCopula(), tau)))[,2], + C = cCopula(u, claytonCopula(iTau(claytonCopula(), tau)))[,2], + F = cCopula(u, frankCopula(iTau( frankCopula(), tau)))[,2]) + } > > ## AMH and Frank "failed" because cop(AMH|Frank) @ absdPsi(*, log=TRUE) gave NaN > (cACF <- cCneg(tau = -0.18)) [,1] [,2] [,3] [,4] [,5] [,6] [,7] A 0.3640882 0.3976526 0.4360816 0.4803628 0.5317478 0.5918415 0.6627290 C 0.3625358 0.4569607 0.5081785 0.5424644 0.5678378 0.5877676 0.6040540 F 0.3488709 0.3974733 0.4481862 0.5000000 0.5518138 0.6025267 0.6511291 [,8] A 0.7471590 C 0.6177456 F 0.6967794 > stopifnot(is.finite(cACF), !apply(cACF, 1, is.unsorted), + 0.348 <= cACF, cACF <= 0.748, + ## *are* somewhat similar as they have same tau: + all.equal(cACF["A",], cACF["F",], tol = 0.035) + , + all.equal(cACF["C",], cACF["F",], tol = 0.079) + ) > > ## FIXME: u1 = 0 still gives NaN, and for Clayton even others > u1. <- c(0, 1e-100, 1e-20, 1e-10, 1e-5, 1e-4, 1e-3, .01) > cCneg(-0.18, u1 = u1.) [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] A NaN 0.3346008 0.3346008 0.3346008 0.3346030 0.3346229 0.3348223 0.33682553 C NaN NaN NaN NaN NaN NaN NaN 0.03284787 F NaN 0.3032206 0.3032206 0.3032206 0.3032241 0.3032558 0.3035723 0.30674776 Warning message: In log1p(-t) : NaNs produced > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > ###---- Large Tau Random Numbers ------------------------------------- > > taus <- c(.80, .85, .90, .95, .98, .99, .993, .995, .996, .997, .998, .999) > namT <- paste0("tau=", formatC(taus)) > archCops <- list(C = claytonCopula, + F = frankCopula, + G = gumbelCopula, + ## A = amhCopula, ## max tau = 1/3 = 0.33333 + J = joeCopula) > thC <- lapply(archCops, function(Cop) setNames(iTau(Cop(), taus), namT)) > simplify2array(thC) C F G J tau=0.8 8.00000 18.19154 5.000000 8.767707 tau=0.85 11.33333 24.90541 6.666667 12.086468 tau=0.9 18.00000 38.28121 10.000000 18.738669 tau=0.95 38.00000 78.31978 20.000000 38.724328 tau=0.98 98.00000 198.34131 50.000000 98.715792 tau=0.99 198.00000 398.34825 100.000000 198.712959 tau=0.993 283.71429 569.77887 142.857143 284.426396 tau=0.995 398.00000 798.35167 200.000000 398.711545 tau=0.996 498.00000 998.35235 250.000000 498.711262 tau=0.997 664.66667 1331.68636 333.333333 665.377646 tau=0.998 998.00000 1998.35371 500.000000 998.710697 tau=0.999 1998.00000 3998.35439 1000.000000 1998.710414 > showProc.time() Time (user system elapsed): 0.06 0 0.07 > > Cops <- lapply(names(thC), function(nm) lapply(thC[[nm]], function(th) archCops[[nm]](th, dim=3))) > uC <- lapply(setNames(,names(thC)), function(nm) + lapply(thC[[nm]], function(th) rCopula(n = 100, archCops[[nm]](th, dim=3)))) > > (aU <- simplify2array(uC)) C F G J tau=0.8 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.85 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.9 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.95 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.98 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.99 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.993 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.995 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.996 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.997 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.998 numeric,300 numeric,300 numeric,300 numeric,300 tau=0.999 numeric,300 numeric,300 numeric,300 numeric,300 > mima <- t(sapply(aU, range)) > stopifnot(!vapply(aU, anyNA, NA), # no NA's + 0 <= mima[,1], mima[,1] <= mima[,2], mima[,2] <= 1) > showProc.time() Time (user system elapsed): 0.16 0.03 0.18 > > proc.time() user system elapsed 23.26 0.42 23.68