R Under development (unstable) (2024-08-15 r87022 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 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 . > > ### Checks for rstable __and__ rCopula() more generally > > require(copula) Loading required package: copula > > source(system.file("Rsource", "utils.R", package="copula", mustWork=TRUE)) Loading required package: tools > source(system.file("Rsource", "cops.R", package="copula", mustWork=TRUE)) > ## --> copcl, copObs, copBnds, excl.2 , copO.2, copBnd.2 > > set.seed(101) > X <- rstable1(1e4, alpha=.9999, beta=1, gamma= .25, delta=1) > summary(X) Min. 1st Qu. Median Mean 3rd Qu. Max. 1592 1592 1593 1594 1593 2092 > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## 1592 1592 1593 1594 1593 2092 > stopifnot(all(X > 1500)) > X <- rstable1(1e4, alpha=.999999, beta=1, gamma= .25, delta=1) > stopifnot(150000 < X, X < 170000) > showProc.time() Time (user system elapsed): 0.06 0 0.06 > > r1R <- copula:::rstable1R > r1C <- copula:::rstable1C > set.seed(123) > > ## speed is *very* similar (!): > showSys.time(Z <- r1R(10000, .2, beta=1, gamma= 45)) Time user system elapsed Time 0 0 0 > showSys.time(Z. <- r1C(10000, .2, beta=1, gamma= 45)) Time user system elapsed Time 0 0 0 > > ks.test(Z, Z.) # p-value ~ 0.50 they "are the same" Asymptotic two-sample Kolmogorov-Smirnov test data: Z and Z. D = 0.0081, p-value = 0.8982 alternative hypothesis: two-sided > > if(!dev.interactive(orNone=TRUE)) pdf("rstable-ex.pdf") > qqplot(Z, Z., log = "xy") # looks nice > acf(log(Z)) # "nice" > acf(log(Z.)) # ditto > > showProc.time() Time (user system elapsed): 0.52 0 0.51 > > ### Ensure basic properties of rCopula() --------------------- > > pkg <- "package:copula" > ## all "acopula"s : > aCops <- sapply(ls(pkg, pattern="^cop[A-Z]"), + get, pkg, simplify=FALSE) > stopifnot(sapply(aCops, is, "acopula")) > > taus <- 1/ c(16, 8, 5, 4) # < 1/3 which is max{ AMH } > > aCt2 <- unlist(lapply(taus, function(TAU) { + lapply(aCops, function(COP) { + th <- iTau(COP, TAU); onacopulaL(COP, list(th, 1:2))})})) > > ## copcl (etc) from ../inst/Rsource/cops.R > stopList <- c("khoudrajiCopula", #"khoudrajiBivCopula", "khoudrajiExplicitCopula", + "indepCopula", "rotCopula", "lowfhCopula", "upfhCopula", "empCopula", "moCopula") > copcl <- copcl[is.na(match(copcl, stopList))] > str(cfn <- sapply(copcl, get, pkg, simplify=FALSE)) List of 13 $ fgmCopula :function (param = NA_real_, dim = 2L) $ plackettCopula :function (param = NA_real_) $ normalCopula :function (param = NA_real_, dim = 2L, dispstr = "ex") $ tCopula :function (param = NA_real_, dim = 2L, dispstr = "ex", df = 4, df.fixed = FALSE, df.min = 0.01) $ claytonCopula :function (param = NA_real_, dim = 2L, use.indepC = c("message", "TRUE", "FALSE")) $ frankCopula :function (param = NA_real_, dim = 2L, use.indepC = c("message", "TRUE", "FALSE")) $ amhCopula :function (param = NA_real_, dim = 2L, use.indepC = c("message", "TRUE", "FALSE")) $ joeCopula :function (param = NA_real_, dim = 2L, use.indepC = c("message", "TRUE", "FALSE")) $ gumbelCopula :function (param = NA_real_, dim = 2L, use.indepC = c("message", "TRUE", "FALSE")) $ galambosCopula :function (param = NA_real_) $ huslerReissCopula:function (param = NA_real_) $ tawnCopula :function (param = NA_real_) $ tevCopula :function (param = NA_real_, df = 4, df.fixed = FALSE) > str(th.25 <- lapply(cfn, function(F) iTau(F(), 0.25))) List of 13 $ fgmCopula : num 1 $ plackettCopula : num 3.14 $ normalCopula : num 0.383 $ tCopula : num 0.383 $ claytonCopula : num 0.667 $ frankCopula : num 2.37 $ amhCopula : num 0.838 $ joeCopula : num 1.6 $ gumbelCopula : num 1.33 $ galambosCopula : num 0.598 $ huslerReissCopula: num 0.987 $ tawnCopula : num 0.651 $ tevCopula : num 0.605 Warning message: In iTauFgmCopula(copula, tau) : For the FGM copula, tau must be in [-2/9, 2/9]. Replacing too small (large) values by lower (upper) bound. > Ct2 <- unlist(lapply(taus, function(TAU) { + lapply(cfn, function(cFn) { th <- iTau(cFn(), TAU); cFn(th)})})) Warning message: In iTauFgmCopula(copula, tau) : For the FGM copula, tau must be in [-2/9, 2/9]. Replacing too small (large) values by lower (upper) bound. > > ## A list of "fully specified" (parameter, dimension) copula models : > length(cops <- c(aCt2, Ct2)) [1] 72 > ## 72 of them ... > > set.seed(17) > Uc <- lapply(cops, rCopula, n = 1024) > > ## rPosStable [coverage] > gC <- cops[[match("gumbelCopula", names(cops))]] # the 1st > set.seed(17); Uc <- c(Uc, list(gumbelC.x = rCopula(gC, n = nrow(Uc[[1]])))) > options('copula:rstable1' = "rPosStable") > set.seed(17); Uc <- c(Uc, list(gumbelC.rP = rCopula(gC, n = nrow(Uc[[1]])))) > options('copula:rstable1' = NULL) # check that things *are* reproducible > set.seed(17); stopifnot(all.equal(Uc[["gumbelC.x"]], rCopula(gC, n = nrow(Uc[[1]])))) > > ## NOTE: The *names* are not unique ==> uses integer indices > > ## 1) Check the tau's > tau.c <- t(sapply(seq_along(cops), function(i) { + U <- Uc[[i]] + c(tauC = tau(cops[[i]]), + tauHat = cor(U[,1], U[,2], method="kendall")) + })) > relE <- apply(tau.c, 1, function(r) (1 - r[[2]]/r[[1]])) # relative "errors" > stopifnot(-0.41 < relE, relE < 0.67, median(abs(relE)) < 0.09) > round(cbind(tau.c, relE), 3) tauC tauHat relE [1,] 0.062 0.037 0.402 [2,] 0.062 0.047 0.243 [3,] 0.062 0.082 -0.313 [4,] 0.062 0.042 0.325 [5,] 0.062 0.060 0.046 [6,] 0.125 0.121 0.034 [7,] 0.125 0.135 -0.078 [8,] 0.125 0.121 0.034 [9,] 0.125 0.132 -0.059 [10,] 0.125 0.096 0.231 [11,] 0.200 0.187 0.064 [12,] 0.200 0.219 -0.093 [13,] 0.200 0.202 -0.010 [14,] 0.200 0.227 -0.137 [15,] 0.200 0.243 -0.216 [16,] 0.250 0.222 0.113 [17,] 0.250 0.263 -0.054 [18,] 0.250 0.244 0.024 [19,] 0.250 0.247 0.013 [20,] 0.250 0.251 -0.005 [21,] 0.062 0.071 -0.142 [22,] 0.062 0.021 0.668 [23,] 0.062 0.080 -0.284 [24,] 0.062 0.041 0.341 [25,] 0.062 0.056 0.110 [26,] 0.062 0.055 0.127 [27,] 0.063 0.059 0.059 [28,] 0.063 0.083 -0.335 [29,] 0.062 0.030 0.516 [30,] 0.062 0.088 -0.409 [31,] 0.062 0.082 -0.314 [32,] 0.062 0.056 0.109 [33,] 0.063 0.056 0.107 [34,] 0.125 0.120 0.038 [35,] 0.125 0.132 -0.057 [36,] 0.125 0.126 -0.008 [37,] 0.125 0.133 -0.062 [38,] 0.125 0.126 -0.009 [39,] 0.125 0.102 0.185 [40,] 0.125 0.137 -0.095 [41,] 0.125 0.106 0.151 [42,] 0.125 0.105 0.160 [43,] 0.125 0.165 -0.320 [44,] 0.125 0.094 0.250 [45,] 0.125 0.124 0.010 [46,] 0.125 0.110 0.123 [47,] 0.200 0.175 0.124 [48,] 0.200 0.246 -0.231 [49,] 0.200 0.209 -0.043 [50,] 0.200 0.231 -0.157 [51,] 0.200 0.188 0.060 [52,] 0.200 0.184 0.080 [53,] 0.200 0.112 0.442 [54,] 0.200 0.206 -0.032 [55,] 0.200 0.203 -0.017 [56,] 0.200 0.232 -0.162 [57,] 0.200 0.200 -0.001 [58,] 0.200 0.174 0.132 [59,] 0.200 0.203 -0.015 [60,] 0.222 0.239 -0.074 [61,] 0.250 0.257 -0.029 [62,] 0.250 0.260 -0.039 [63,] 0.250 0.218 0.127 [64,] 0.250 0.247 0.012 [65,] 0.250 0.258 -0.034 [66,] 0.250 0.269 -0.076 [67,] 0.250 0.239 0.045 [68,] 0.250 0.223 0.108 [69,] 0.250 0.249 0.003 [70,] 0.250 0.259 -0.034 [71,] 0.250 0.242 0.030 [72,] 0.250 0.245 0.018 > ## we could look at cases which were "bad" .. > plot(tau.c[,"tauC"], relE) # error is larger for small true taus > > > ## 2) Check the margins: all must be uniform > ad.test <- ADGofTest::ad.test > pp <- sapply(Uc, function(U2) + apply(U2, 2, + function(u) ad.test(u)$p.value)) > > ## The P-values should simply be uniform in [0,1]: > hh <- hist(c(pp), breaks = (0:20)/20)## should "look" uniform > summary(hh$counts) Min. 1st Qu. Median Mean 3rd Qu. Max. 4.0 6.0 7.0 7.4 9.0 11.0 > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## 4.00 5.75 7.00 7.20 9.00 11.00 > > ## and hence their test should typically *not* be significant > (ad.pp <- ad.test(c(pp))) Anderson-Darling GoF Test data: c(pp) AD = 0.64863, p-value = 0.603 alternative hypothesis: NA > stopifnot(hh$counts < 15, + ad.pp$p.value > 0.10) > > proc.time() user system elapsed 3.43 0.14 3.56