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 . > > require(copula) Loading required package: copula > > if(getRversion() < "2.15") + paste0 <- function(...) paste(..., sep="") > > source(system.file("Rsource", "utils.R", package="copula", mustWork=TRUE)) Loading required package: tools > ##-> setPar(), showProc.time() etc > ## All non-virtual copula classes: > source(system.file("Rsource", "cops.R", package="copula", mustWork=TRUE)) > ## --> copcl, copcl., copObs, copBnds, excl.2 , copO.2, copBnd.2 > > options(width = 125)# -> nicer table printing > showProc.time() Time (user system elapsed): 0.05 0 0.04 > > copcl ## the classes (incl. 'indepCopula') [1] "moCopula" "fgmCopula" "plackettCopula" "normalCopula" "tCopula" "claytonCopula" [7] "frankCopula" "amhCopula" "joeCopula" "gumbelCopula" "galambosCopula" "huslerReissCopula" [13] "tawnCopula" "tevCopula" > str(copObs, max.level=1)# copula objects (w/o 'indepCopula') List of 13 $ fgmCopula :Formal class 'fgmCopula' [package "copula"] with 8 slots $ plackettCopula :Formal class 'plackettCopula' [package "copula"] with 7 slots $ normalCopula :Formal class 'normalCopula' [package "copula"] with 8 slots $ tCopula :Formal class 'tCopula' [package "copula"] with 9 slots $ claytonCopula :Formal class 'claytonCopula' [package "copula"] with 7 slots $ frankCopula :Formal class 'frankCopula' [package "copula"] with 7 slots $ amhCopula :Formal class 'amhCopula' [package "copula"] with 7 slots $ joeCopula :Formal class 'joeCopula' [package "copula"] with 7 slots $ gumbelCopula :Formal class 'gumbelCopula' [package "copula"] with 7 slots $ galambosCopula :Formal class 'galambosCopula' [package "copula"] with 7 slots $ huslerReissCopula:Formal class 'huslerReissCopula' [package "copula"] with 7 slots $ tawnCopula :Formal class 'tawnCopula' [package "copula"] with 7 slots $ tevCopula :Formal class 'tevCopula' [package "copula"] with 7 slots > ## and their parameter bounds: > t(copBnds) min max fgmCopula -1 1 plackettCopula 0 Inf normalCopula -1 1 tCopula -1 1 claytonCopula -1 Inf frankCopula -Inf Inf amhCopula -1 1 joeCopula 1 Inf gumbelCopula 1 Inf galambosCopula 0 Inf huslerReissCopula 0 Inf tawnCopula 0 1 tevCopula -1 1 > > ###-------- tau & and inverse --------------------------------------------------- > > ## currently fails: --- FIXME?: should AMH also warn like the others? > tau.s <- c(-.999, -.1, 0, (1:3)/10, .5, .999) > ### give different warnings , but "work" : { .5 , even 1/3, gives error for AMH FIXME} > > ## > tau(tevCopula(0)) # 0.05804811 [1] 0.05804811 > ## restricted tau-range works better > tau.s <- c( -.1, 0, 0.05805, (1:2)/9, 0.3) > names(tau.s) <- paste0("tau=", sub("0[.]", ".", formatC(tau.s))) > tTau <- sapply(tau.s, function(tau) vapply(copObs, iTau, numeric(1), tau = tau)) Warning messages: 1: In .local(copula, tau, ...) : For the Joe copula, tau must be >= 0. Replacing negative values by 0. 2: In iTauGumbelCopula(copula, tau) : For the Gumbel copula, tau must be >= 0. Replacing negative values by 0. 3: In iTauGalambosCopula(copula, tau) : For the Galambos copula, tau must be >= 0. Replacing negative values by 0. 4: In iTauHuslerReissCopula(copula, tau) : For the Husler-Reiss copula, tau must be >= 0. Replacing negative values by 0. 5: In iTauTawnCopula(copula, tau) : For the Tawn copula, tau must be in [0, 0.4183992]. Replacing too small (large) values by lower (upper) bound. 6: In iTauTevCopula(copula, tau) : For the t-ev copula, tau must be >= 0. Replacing negative values by 0. 7: 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. > ## -> 7 warnings: (Joe, Gumbel, Galambos, Husler-Reiss, Tawn, t-ev, FGM) > tTau tau=-.1 tau=0 tau=.05805 tau=.1111 tau=.2222 tau=.3 fgmCopula -0.4500000 0.000000e+00 0.261225000 0.5000000 1.0000000 1.0000000 plackettCopula 0.6395724 1.001153e+00 1.297034794 1.6430737 2.7545621 3.9867469 normalCopula -0.1564345 0.000000e+00 0.091058418 0.1736482 0.3420201 0.4539905 tCopula -0.1564345 0.000000e+00 0.091058418 0.1736482 0.3420201 0.4539905 claytonCopula -0.1818182 0.000000e+00 0.123254950 0.2500000 0.5714286 0.8571429 frankCopula -0.9073675 0.000000e+00 0.523881124 1.0101317 2.0843872 2.9174344 amhCopula -0.5030297 1.368215e-13 0.244596257 0.4404231 0.7715254 0.9429734 joeCopula 1.0000000 1.000000e+00 1.107182397 1.2190610 1.5088681 1.7721048 gumbelCopula 1.0000000 1.000000e+00 1.061627475 1.1250000 1.2857143 1.4285714 galambosCopula 0.0000000 0.000000e+00 0.269205006 0.3581411 0.5458864 0.6992740 huslerReissCopula 0.0000000 0.000000e+00 0.562768150 0.6830113 0.9228086 1.1110641 tawnCopula 0.0000000 0.000000e+00 0.168274582 0.3124090 0.5876047 0.7613025 tevCopula 0.0000000 0.000000e+00 0.000739897 0.2597801 0.5545006 0.6818544 > > stopifnot(rep(copBnds["min",],ncol(tTau)) <= tTau + 1e-7, + tTau <= rep(copBnds["max",],ncol(tTau)), + tTau["joeCopula", "tau=-.1"] == 1, + ## theta and tau are comonotone : + apply(tTau, 1, diff) >= -1e9) > > showProc.time() Time (user system elapsed): 0.06 0 0.07 > > tautau <- t(sapply(names(copObs), function(cNam) + sapply(tTau[cNam,], + function(th) tau(setPar(copObs[[cNam]], th))))) > tautau tau=-.1 tau=0 tau=.05805 tau=.1111 tau=.2222 tau=.3 fgmCopula -1.000000e-01 0.000000e+00 0.05805000 0.1111111 0.2222222 0.2222222 plackettCopula -9.999994e-02 5.136148e-04 0.05804998 0.1111111 0.2222221 0.2999999 normalCopula -1.000000e-01 0.000000e+00 0.05805000 0.1111111 0.2222222 0.3000000 tCopula -1.000000e-01 0.000000e+00 0.05805000 0.1111111 0.2222222 0.3000000 claytonCopula -1.000000e-01 0.000000e+00 0.05805000 0.1111111 0.2222222 0.3000000 frankCopula -1.000000e-01 0.000000e+00 0.05805000 0.1111111 0.2222222 0.3000000 amhCopula -1.000000e-01 3.040477e-14 0.05805000 0.1111111 0.2222222 0.3000000 joeCopula 4.440892e-16 -1.978007e-10 0.05805000 0.1111111 0.2222222 0.3000000 gumbelCopula 0.000000e+00 0.000000e+00 0.05805000 0.1111111 0.2222222 0.3000000 galambosCopula 5.545057e-69 5.545057e-69 0.05804706 0.1111090 0.2222222 0.3000004 huslerReissCopula -6.111270e-113 -6.111270e-113 0.05804589 0.1111088 0.2222220 0.3000004 tawnCopula NaN NaN 0.05805000 0.1111111 0.2222222 0.3000000 tevCopula 5.804811e-02 5.804811e-02 0.05805055 0.1111113 0.2222222 0.3000000 > showProc.time() Time (user system elapsed): 0.03 0 0.03 > > xctTau <- matrix(tau.s, nrow = nrow(tautau), ncol=length(tau.s), + byrow=TRUE) > ## The absolute errors > errTau <- tautau-xctTau > print.table(round(10000*errTau), zero.print=".", na.print="_NA_") tau=-.1 tau=0 tau=.05805 tau=.1111 tau=.2222 tau=.3 fgmCopula . . . . . -778 plackettCopula . 5 . . . . normalCopula . . . . . . tCopula . . . . . . claytonCopula . . . . . . frankCopula . . . . . . amhCopula . . . . . . joeCopula 1000 . . . . . gumbelCopula 1000 . . . . . galambosCopula 1000 . . . . . huslerReissCopula 1000 . . . . . tawnCopula _NA_ _NA_ . . . . tevCopula 1580 580 . . . . > ## has two NaN .. ok, for now: > errTau["tawnCopula", 1:2] <- 0 > ## These families do not support tau < 0 > errTau[c("gumbelCopula", "joeCopula", + "galambosCopula", "huslerReissCopula","tevCopula"), + "tau=-.1"] <- 0 > ## the tevCopula cannot get a tau = 0 (for now) __FIXME?__ > errTau["tevCopula", 2] <- 0 > ## "fgmCopula" has tau in [-2/9, 2/9] : > errTau["fgmCopula", "tau=.3"] <- 0 > stopifnot(max(abs(errTau)) <= 0.00052)# ok for IJ-taus > ## show the blown up deviations: > print.table(round(1e7*errTau, 4), zero.print=".", na.print="_NA_") tau=-.1 tau=0 tau=.05805 tau=.1111 tau=.2222 tau=.3 fgmCopula . . . . . . plackettCopula 0.5935 5136.1485 -0.1921 -0.3289 -1.4855 -0.8832 normalCopula . . . . . . tCopula . . . . . . claytonCopula . . . . . . frankCopula . . . 0.0011 . -0.0001 amhCopula . . . -0.0060 0.0837 -0.0001 joeCopula . -0.0020 0.0048 0.0006 -0.0437 0.0009 gumbelCopula . . . . . . galambosCopula . . -29.4444 -21.4249 -0.2935 3.5491 huslerReissCopula . . -41.1256 -23.1996 -2.2515 4.2730 tawnCopula . . -0.0001 -0.0016 -0.0356 -0.0046 tevCopula . . 5.4781 2.2798 -0.3300 -0.2011 > > ## now remove the current worst: > errT <- errTau; errT["plackettCopula", "tau=0"] <- 0 > 10e6 * head(sort(abs(errT), decreasing=TRUE)) [1] 41.125574 29.444436 23.199626 21.424885 5.478115 4.273005 > stopifnot(max(abs(errT)) <= 5e6) > > showProc.time() Time (user system elapsed): 0 0 0 > > > ###-------- rho & and inverse --------------------------------------------------- > > ## Give different warnings, but "work" (not using AMH, Joe, t copula > ## as iRho() does not work for them) > rho.s <- c(-.999, -.1, 0, (1:3)/9, .5, .9, .999) > names(rho.s) <- paste0("rho=", sub("0[.]", ".", formatC(rho.s))) > tRho <- sapply(rho.s, function(rho) + sapply(copO.2, iRho, rho = rho)) There were 16 warnings (use warnings() to see them) > tRho rho=-.999 rho=-.1 rho=0 rho=.1111 rho=.2222 rho=.3333 rho=.5 rho=.9 rho=.999 fgmCopula -1.000000e+00 -0.3000000 0.000000e+00 0.3333333 0.6666667 1.0000000 1.0000000 1.0000000 1.000000e+00 plackettCopula 5.699805e-05 0.7401479 1.000000e+00 1.3973489 1.9676413 2.8169605 5.1156609 66.1154147 1.754744e+04 normalCopula -9.990930e-01 -0.1046719 0.000000e+00 0.1162897 0.2321858 0.3472964 0.5176381 0.9079810 9.990930e-01 claytonCopula -9.989032e-01 -0.1257142 NA 0.1604531 0.3508275 0.5861439 1.0759811 5.5660650 NA frankCopula -1.390113e+02 -0.6029018 2.965981e-15 0.6706543 1.3661619 2.1164969 3.4459877 12.2614867 1.390113e+02 gumbelCopula 1.001106e+00 1.0011059 1.001106e+00 1.0818155 1.1761467 1.2933287 1.5442144 3.7355389 3.823601e+01 galambosCopula 0.000000e+00 0.0000000 0.000000e+00 0.2975800 0.4214715 0.5554504 0.8133918 3.0076496 Inf huslerReissCopula 0.000000e+00 0.0000000 0.000000e+00 0.6017263 0.7651298 0.9336457 1.2451931 3.7282213 Inf tawnCopula 0.000000e+00 0.0000000 0.000000e+00 0.2150646 0.4167466 0.6063103 0.8703185 1.0000000 1.000000e+00 tevCopula 0.000000e+00 0.0000000 0.000000e+00 0.1023452 0.3863229 0.5661123 0.7459998 0.9716375 9.997916e-01 > warnings()## 16 warnings [2014-05] Warning messages: 1: In iRhoFgmCopula(copula, rho) : For the FGM copula, rho must be in [-1/3, 1/3]. Replacing too small (large) values by lower (upper) bound. 2: In iRhoGumbelCopula(copula, rho) : For the Gumbel copula, rho must be >= 0. Replacing negative values by 0. 3: In iRhoGalambosCopula(copula, rho) : For the Galambos copula, rho must be >= 0. Replacing negative values by 0. 4: In iRhoHuslerReissCopula(copula, rho) : For the Husler-Reiss copula, rho must be >= 0. Replacing negative values by 0. 5: In iRhoTawnCopula(copula, rho) : For the Tawn copula, rho must be in [0, 0.58743682]. Replacing too small (large) values by lower (upper) bound. 6: In iRhoTevCopula(copula, rho) : For the t-ev copula, rho must be >= 0. Replacing negative values by 0. 7: In iRhoGumbelCopula(copula, rho) : For the Gumbel copula, rho must be >= 0. Replacing negative values by 0. 8: In iRhoGalambosCopula(copula, rho) : For the Galambos copula, rho must be >= 0. Replacing negative values by 0. 9: In iRhoHuslerReissCopula(copula, rho) : For the Husler-Reiss copula, rho must be >= 0. Replacing negative values by 0. 10: In iRhoTawnCopula(copula, rho) : For the Tawn copula, rho must be in [0, 0.58743682]. Replacing too small (large) values by lower (upper) bound. 11: In iRhoTevCopula(copula, rho) : For the t-ev copula, rho must be >= 0. Replacing negative values by 0. 12: In iRhoFgmCopula(copula, rho) : For the FGM copula, rho must be in [-1/3, 1/3]. Replacing too small (large) values by lower (upper) bound. 13: In iRhoFgmCopula(copula, rho) : For the FGM copula, rho must be in [-1/3, 1/3]. Replacing too small (large) values by lower (upper) bound. 14: In iRhoTawnCopula(copula, rho) : For the Tawn copula, rho must be in [0, 0.58743682]. Replacing too small (large) values by lower (upper) bound. 15: In iRhoFgmCopula(copula, rho) : For the FGM copula, rho must be in [-1/3, 1/3]. Replacing too small (large) values by lower (upper) bound. 16: In iRhoTawnCopula(copula, rho) : For the Tawn copula, rho must be in [0, 0.58743682]. Replacing too small (large) values by lower (upper) bound. > ## and from now on, show them as they happen: > options(warn = 1) > > ##--> oops! clayton [rho=0] is NA __FIXME__ > tRho["claytonCopula", "rho=0"] <- 0 > ## and it has NA also for .999 {but that maybe consider ok}: > tRho["claytonCopula", "rho=.999"] <- 10^100 > > stopifnot(rep(copBnd.2["min",],ncol(tRho)) <= tRho, + tRho <= rep(copBnd.2["max",],ncol(tRho)), + ## theta and rho are comonotone : + apply(tRho, 1, diff) >= 0) > > rhorho <- t(sapply(names(copO.2), function(cNam) + sapply(tRho[cNam,], + function(th) rho(setPar(copO.2[[cNam]], th))))) > rhorho rho=-.999 rho=-.1 rho=0 rho=.1111 rho=.2222 rho=.3333 rho=.5 rho=.9 rho=.999 fgmCopula -3.333333e-01 -1.000000e-01 0.000000e+00 0.1111111 0.2222222 0.3333333 0.3333333 0.3333333 0.3333333 plackettCopula -9.989998e-01 -9.999999e-02 1.490116e-08 0.1111111 0.2222222 0.3333333 0.5000000 0.9000000 0.9990000 normalCopula -9.990000e-01 -1.000000e-01 0.000000e+00 0.1111111 0.2222222 0.3333333 0.5000000 0.9000000 0.9990000 claytonCopula -9.990000e-01 -9.999988e-02 -2.581600e-05 0.1111124 0.2222239 0.3333382 0.5000016 0.9000001 0.9953108 frankCopula -9.990000e-01 -1.000000e-01 4.943302e-16 0.1111111 0.2222222 0.3333333 0.5000000 0.9000000 0.9990000 gumbelCopula 2.382468e-05 2.382468e-05 2.382468e-05 0.1111487 0.2222907 0.3334172 0.5000239 0.9000001 0.9990000 galambosCopula -1.312275e-74 -1.312275e-74 -1.312275e-74 0.1111109 0.2222216 0.3333334 0.5000008 0.9000001 0.9965914 huslerReissCopula -5.919435e-129 -5.919435e-129 -5.919435e-129 0.1111084 0.2222217 0.3333335 0.5000006 0.9000005 0.9957213 tawnCopula 0.000000e+00 0.000000e+00 0.000000e+00 0.1111111 0.2222222 0.3333333 0.5000000 0.5874368 0.5874368 tevCopula 8.404841e-02 8.404841e-02 8.404841e-02 0.1111110 0.2222221 0.3333333 0.5000000 0.8999971 0.9989899 > showProc.time() Time (user system elapsed): 0.29 0 0.28 > > xctRho <- matrix(rho.s, nrow = nrow(rhorho), ncol=length(rho.s), + byrow=TRUE) > ## The absolute errors > errRho <- rhorho-xctRho > round(10000*errRho) rho=-.999 rho=-.1 rho=0 rho=.1111 rho=.2222 rho=.3333 rho=.5 rho=.9 rho=.999 fgmCopula 6657 0 0 0 0 0 -1667 -5667 -6657 plackettCopula 0 0 0 0 0 0 0 0 0 normalCopula 0 0 0 0 0 0 0 0 0 claytonCopula 0 0 0 0 0 0 0 0 -37 frankCopula 0 0 0 0 0 0 0 0 0 gumbelCopula 9990 1000 0 0 1 1 0 0 0 galambosCopula 9990 1000 0 0 0 0 0 0 -24 huslerReissCopula 9990 1000 0 0 0 0 0 0 -33 tawnCopula 9990 1000 0 0 0 0 0 -3126 -4116 tevCopula 10830 1840 840 0 0 0 0 0 0 > ## the tevCopula cannot get a rho <= 0 (for now) __FIXME?__ > errRho["tevCopula", 1:3] <- 0 > ## These three families do not support rho < 0 (currently): > errRho[c("gumbelCopula", "galambosCopula", "huslerReissCopula", "tawnCopula"), + c("rho=-.1","rho=-.999")] <- 0 > errRho["tawnCopula", rho.s >= 0.9] <- 0 > ## "fgmCopula" has rho in [-1/3, 1/3] : > errRho["fgmCopula", abs(rho.s) > 1/3] <- 0 > > stopifnot(max(abs(errRho)) <= 0.00369, + max(abs(errRho[,rho.s <= 0.9])) <= 0.0002)# ok for IJ-rhos > > showProc.time() Time (user system elapsed): 0 0 0 > > > > proc.time() user system elapsed 1.39 0.21 1.56