R Under development (unstable) (2025-07-22 r88445 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 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. > #### A "fast" and simplified version of demo(VaRsuperadd) > #### see ../demo/VaRsuperadd.R > #### ~~~~~~~~~~~~~~~~~~~~~ > > require(simsalapar) Loading required package: simsalapar > source(system.file("xtraR/assertErr-etc.R", package="simsalapar", mustWork=TRUE)) > > ## Must be fast, rather than "interesting": > n.obs <- 16 > n.alpha <- 8 > doExtras <- FALSE > > demo(VaRsuperadd) # crucial part! demo(VaRsuperadd) ---- ~~~~~~~~~~~ > ## Copyright Marius Hofert > > ### Goal: Simulation study for determining superadditivity of VaR > ### Submit via (for example): > ### bsub -N -W 01:00 -n 48 -R "select[model==Opteron8380]" -R "span[ptile=16]" \ > ### mpirun -n 1 R CMD BATCH Brutus_demo.R > > > ### 0) Setup ################################################################### > > usr <- Sys.info()[["user"]] # better portability than Sys.getenv("USER") > if(exists("n.obs") && is.numeric(n.obs) && + exists("n.alpha") && is.numeric(n.alpha)) { + cat(">>> Using n.obs=", n.obs, ", n.alpha=", n.alpha, " <<<\n", sep="") + } else { + ## user specific settings + n.obs <- 100 # default + n.alpha <- 32 + switch(usr, + "hofertj" = { # ETH Brutus + setwd("/cluster/home/math/hofertj/VaRsuperadd") + .libPaths("~/R/library") + n.obs <- 1e4 + n.alpha <- 128 + }, + "mhofert" = { + n.obs <- 1e4 # n=1e4 ~ 4min (standard notebook) + n.alpha <- 128 + }, + "maechler" = {} + ) + } >>> Using n.obs=16, n.alpha=8 <<< > if(n.alpha > n.obs) + stop("number of quantiles, n.alpha, must be smaller than the sample size, n.obs") > if(!exists("doExtras")) doExtras <- simsalapar:::doExtras() > doExtras [1] FALSE > ## load packages > require(simsalapar) > ## list of variables > varList <- + varlist( + ## sample size + n = list(value = n.obs), + ## dimensions, and weights (vector) for each d + d = list(type="grid", value = c(4, 20, 100)), + ## copula family names + family = list(type="grid", expr = quote(C), + value = c("normal", "t", "Clayton", "Gumbel")), # t = t_4 + ## dependencies by Kendall's tau + tau = list(type="grid", value = c(0.2, 0.5, 0.8)), + ## margins + qmargin = list(type="inner", expr = quote(F[j]), + value = c(norm = qnorm, + t4 = function(p) qt(p, df=4), + Par2 = function(p) (1-p)^(-1/2))), # Pareto(2) + ## VaR confidence levels + alpha = list(type="inner", value = 0:n.alpha/n.alpha)) > ### 1) Functions ############################################################### > > ##' @title Function to Compute F_{X_1+..+X_d}(d*F_1^-(\alpha)) > ##' @author Marius Hofert > doOne <- function(n, d, family, tau, qmargin, alpha) + { + ## checks (and load required packages here for parallel computing later on) + ## stopifnot(require(copula)) + ## Note: This require() affects run time due to loading copula once on each + ## node, but this can easily be identified by *robust* analysis of all + ## run times (or is negligible due to much larger run times of the jobs) + ## or, as another alternative, can be avoided by using initExpr (see below) + + cop <- switch(family, + "normal" = + ellipCopula("normal", param=iTau(ellipCopula("normal"), tau=tau), + dim=d), + "t" = + ellipCopula("t", param=iTau(ellipCopula("t"), tau=tau), dim=d), + "Clayton" = + onacopulaL("Clayton", list(th=iTau(archmCopula("clayton"), tau), + seq_len(d))), + "Gumbel" = + onacopulaL("Gumbel", list(th=iTau(archmCopula("gumbel"), tau), + seq_len(d))), + stop("unsupported 'family'")) + U <- rCopula(n, copula=cop) + + ## compute F_{X_1+..+X_d}(d*F_1^-(\alpha)) for all confidence levels alpha + ## => VaR_alpha superadditive <=> F_{X_1+..+X_d}(d*F_1^-(\alpha)) - alpha < 0 + t(sapply(qmargin, function(FUN) ecdf(rowSums(FUN(U)))(d*FUN(alpha)) - alpha)) + ## note: t() is important here, since, otherwise, the order of the variables + ## ---- would not be correct (=> check should reveal this) + } > ### 2) Main #################################################################### > > ## check doOne > > ## manually > require(copula) # for the following call of doOne() Loading required package: copula > nonGr <- get.nonGrids(varList)$nonGrids > dd <- doOne(n= min(nonGr$n, 100), d=4, family="Clayton", tau=0.5, + qmargin=nonGr$qmargin, alpha=nonGr$alpha) > stopifnot(dim(dd) == with(nonGr, c(length(qmargin), length(alpha)))) > ## with doCheck > doCheck(doOne, varList) j = 18: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 27: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 12: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 32: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 7: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 25: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 8: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 16: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 1: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 4: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 6: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 10: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 35: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 28: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 14: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 2: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 24: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 22: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 23: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 11: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 34: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 30: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 26: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 13: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 5: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 19: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 17: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 33: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 31: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 9: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 29: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 15: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 36: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 21: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 3: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL j = 20: -> length(r) = 27 dim(r):[1] 3 9 str(dimnames(r)):List of 2 $ : chr [1:3] "norm" "t4" "Par2" $ : NULL > ## if simsalapar no longer *depends* on parallel: > makeCluster <- parallel::makeCluster > ## to be CRAN check compatible and as called from tests with explicit doExtras <- FALSE: > nc <- simsalapar:::nCores4test(); nc <- if(doExtras) nc else min(nc, 2) > nc [1] 2 > nc.win <- if(.Platform$OS.type=="windows") 1 else nc # otherwise win-builder fails > ## computation > system.time(res <- + doClusterApply(varList, cluster=makeCluster(nc, type="PSOCK"), + sfile = if(n.obs > 1000) "VaR_superadd.rds" else NULL, + doOne=doOne, monitor = printInfo[["gfile"]], + ## load copula once on each worker, to not affect run time + initExpr = require("copula"), + timer=mkTimer(gcFirst=TRUE) )) user system elapsed 0.70 0.37 6.22 > if(doExtras) { + ## doMclapply() ------------------------------------------------------------ + ## "init expression" for mclapply can happen here locally (see above) + print(system.time(res2 <- doMclapply(varList, cores=nc.win, doOne=doOne))) + stopifnot( doRes.equal(res, res2) ) + + ## doRmpi() ---------------------------------------------------------------- + ## => passing an init expression does not (easily) work + + ## doForeach() ------------------------------------------------------------- + print(system.time(res3 <- { + doForeach(varList, cluster=makeCluster(nc, type="PSOCK"), + doOne=doOne, extraPkgs = "copula", + timer=mkTimer(gcFirst=TRUE) ) + })) + stopifnot( doRes.equal(res, res3) ) + } > ### 3) Analysis ################################################################ > > ## extract results > val <- getArray(res) # array of values > err <- getArray(res, "error") # array of error indicators > warn <- getArray(res, "warning") # array of warning indicators > time <- getArray(res, "time") # array of user times in ms > ## warnings, errors > if(any(err > 0)) + ftable(100* err, col.vars="tau") # percentage of errors > if(any(warn > 0)) + ftable(100*warn, col.vars="tau") # percentage of warnings > ## run time > ftable(time, row.vars=c("family", "d"), col.vars="tau") tau 0.2 0.5 0.8 family d normal 4 20 0 20 20 20 0 10 100 20 20 20 t 4 0 0 20 20 0 0 0 100 20 20 20 Clayton 4 0 0 0 20 0 0 0 100 50 20 20 Gumbel 4 30 0 0 20 20 0 0 100 0 20 0 > ## add 'tau==' (just nicer) > str(val) num [1:3, 1:9, 1:3, 1:4, 1:3] 0 0 0 -0.125 -0.125 ... - attr(*, "dimnames")=List of 5 ..$ qmargin: chr [1:3] "norm" "t4" "Par2" ..$ alpha : NULL ..$ d : chr [1:3] "4" "20" "100" ..$ family : chr [1:4] "normal" "t" "Clayton" "Gumbel" ..$ tau : chr [1:3] "0.2" "0.5" "0.8" > dimnames(val)[["tau"]] <- paste0("tau==", dimnames(val)[["tau"]]) > ## plot of VaR estimates > ## plotting to pdf if not interactive graphics (=> R CMD BATCH VaRsuperadd.R) > doPDF <- if(usr=="mhofert") TRUE else !dev.interactive(orNone=TRUE) > if(!doPDF) par(ask=TRUE) > names(varList[["qmargin"]][["value"]]) # "norm" , "t4", "Par2" [1] "norm" "t4" "Par2" > for(m in names(varList[["qmargin"]][["value"]])) { + if(doPDF) pdf(paste0("VaR_superadd_", m, ".pdf"), width=6, height=6) + mayplot(val[qmargin=m,,,,], varList, row.vars="family", col.vars="tau", + xvar="alpha", ylim=if(n.obs > 1000) "local" else "global", + panel.first = function(...) abline(h=0, col="gray40"), # gray line + panel.last = function(x,y, col, ...) { + ## For demo: write the 'panel.last' string, dependent on (x,y) + rx <- range(x); dx <- diff(rx) + ry <- range(y); dy <- diff(ry) + text(rx[1]+.7*dx, ry[1]+.2*dy, "< 0: superadd.", col=col) + }) + if(doPDF) dev.off() + } > > ## The only difference to doOne() is the missing t(.) at the very end > do1.wrong <- function(n, d, family, tau, qmargin, alpha) + { + cop <- switch(family, + "normal" = + ellipCopula("normal", param=iTau(ellipCopula("normal"), tau=tau), + dim=d), + "t" = + ellipCopula("t", param=iTau(ellipCopula("t"), tau=tau), dim=d), + "Clayton" = + onacopulaL("Clayton", list(th=iTau(archmCopula("clayton"), tau), + seq_len(d))), + "Gumbel" = + onacopulaL("Gumbel", list(th=iTau(archmCopula("gumbel"), tau), + seq_len(d))), + stop("unsupported 'family'")) + U <- rCopula(n, copula=cop) + sapply(qmargin, function(FUN) ecdf(rowSums(FUN(U)))(d*FUN(alpha)) - alpha) + } > > assertError( doCheck(do1.wrong, varList) ) j = 25: -> length(r) = 27 dim(r):[1] 9 3 str(dimnames(r)):List of 2 $ : NULL $ : chr [1:3] "norm" "t4" "Par2" > > proc.time() user system elapsed 3.34 0.82 9.31