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(TGforecasts)
> #### see ../demo/TGforecasts.R
> ####     ~~~~~~~~~~~~~~~~~~~~~
> 
> require(simsalapar)
Loading required package: simsalapar
> ## also  require(fGarch)
> 
> 
> ### Variable list ##############################################################
> 
> ## Create varlist (default type is 'frozen'; default expressions are generated)
> ## note: no n.sim here
> vList <-
+     varlist(forecaster = list(type="grid", expr = quote(italic(forecaster)),
+                               value = c("statistician", "optimist", "pessimist")),
+ 	    scoringfun = list(type="inner", expr = quote(S),
+ 			      value = list(SE = function(x, y) mean((x-y)^2),
+ 					   AE = function(x, y) mean(abs(x-y)))),
+             ## GARCH(1,1): X_t = sigma_t * eps_t,
+             ## sigma_t^2 = alpha_0 + alpha_1*X_{t-1}^2 + beta_1*sigma_{t-1}^2
+             alpha0     = list(expr = quote(alpha[0]), value = 0.05),
+             alpha1     = list(expr = quote(alpha[1]), value = 0.2),
+             beta1      = list(expr = quote(beta[1]), value = 0.75),
+             ## constant prediction of optimist:
+             pred.optimist  = list(expr = "Optimist", value = 5),
+             ## constant prediction of pessimist:
+             pred.pessimist = list(expr = "Pessimist", value = 0.05),
+             ## n
+             n = list(value = 64),
+             ## simulation method
+             method = list(value = "fGarch"))
> 
> ### doOne() including ingredient functions #####################################
> ##  -------
> ##  instead of 'local', use a construction like
> ##  doOne <- function(x, <nonGrids>, simGARCH11, predictTG)
> 
> doOne <- local({
+ 
+     ## function for simulating a squared GARCH(1,1)
+     simGARCH11 <- function(n, a0, a1, b1, method=c("fGarch", "Gneiting"))
+     {
+         method <- match.arg(method)
+         switch(method,
+                "fGarch" = {
+                    ## specify the GARCH(1,1) process
+ 		   spec <- fGarch::garchSpec(model =
+ 					     list(omega = a0, alpha = a1, beta = b1))
+ 		   ## simulate; extended=TRUE => data.frame
+ 		   ##		with columns garch, sigma, eps  time series
+                    fGarch::garchSim(spec, n=n, extended=TRUE)
+                },
+                "Gneiting" = { # Gneiting's code (with minor adjustments)
+                    alpha <- a1
+                    beta <- b1
+                    gamma <- a0
+                    z <- rep(0, n+1)
+                    y <- rep(0, n+1)
+                    sigmasq <- rep(1, n+1)
+                    z[1] <- rnorm(1)
+                    for (t in 2:(n+1)) {
+                        sigmasq[t] <- alpha*z[t-1]^2 + beta*sigmasq[t-1] + gamma
+                        z[t] <- rnorm(1, sd=sqrt(sigmasq[t]))
+                    }
+                    ts(cbind(garch=z[-1], sigma=sqrt(sigmasq[-1])))
+                },
+                stop("wrong 'method'"))
+     }
+ 
+     ## prediction function
+     predictTG <- function(x, n, forecaster, pred.opt, pred.pes) {
+         stopifnot(is.character(forecaster))
+         switch(forecaster,
+                "optimist"  = rep(pred.opt, n),
+                "pessimist" = rep(pred.pes, n),
+                ## statistician predicts hat{X}_t = E[Y_t | sigma_t^2] = sigma_t^2
+                "statistician" = x[,"sigma"]^2,
+                stop("forecaster not supported"))
+     }
+ 
+     ##' Statistic (use it like this to include simulation in time measurement)
+ 
+     ##' @title Function to Compute the Simulation Results for One Grid Line
+     ##' @param x one-line data.frame containing a combination of grid variables
+     ##' @param nonGrids values of non-"grid"-variables
+     ##' @return return value of doCallWE()
+     ##' @author Marius Hofert
+     function(n, alpha0, alpha1, beta1, forecaster,
+              pred.optimist, pred.pessimist, method, scoringfun)
+     {
+         ## simulate squared GARCH(1,1)
+         X <- simGARCH11(n, a0=alpha0, a1=alpha1, b1=beta1, method=method)
+         Y <- X[,"garch"]^2              # see (2) in Gneiting (2011)
+         ## predict
+         pred <- predictTG(X, n=n, forecaster=forecaster,
+                           pred.opt=pred.optimist, pred.pes=pred.pessimist)
+ 
+         ## basic check
+         stopifnot(length(pred)==length(Y))
+ 
+         ## compute scoring functions (for all 'inner' variables simultaneously)
+         c(SE = scoringfun[["SE"]](pred, y=Y),
+           AE = scoringfun[["AE"]](pred, y=Y))
+     }
+ })
> 
> ## build grid and non-grid variables
> stopifnot(dim(pGrid <- mkGrid(vList)) == c(3,1), # only the forecaster here
+ 	  get.nonGrids(vList)$n.sim == 1) # the "automatic n.sim" if there is none
> 
> tex.vl <- toLatex(vList)
> 
> ### Computation ################################################################
> 
> S.T <- system.time
> 
> ## 'reproducing' Table 4 in Gneiting (2011)
> 
> ## sequentially
> S.T(res <- doLapply(vList, doOne=doOne, monitor=TRUE))
07:11:41 on CRANWIN3: forecaster=statistician; value:  num [1:2] 0.516 0.528
07:11:41 on CRANWIN3: forecaster=    optimist; value:  num [1:2] 20.55 4.47
07:11:41 on CRANWIN3: forecaster=   pessimist; value:  num [1:2] 0.792 0.499
   user  system elapsed 
   1.70    0.16    1.86 
> 
> ## our own monitoring function:
> myMoni <- function(i.sim, j, pGrid, res4, n.sim) {
+     cat(SysI(),": ",formG(pGrid, j), "; SE=",
+         format(res4$value[["SE"]], digits=9, width=15), "\n", sep = "")
+ }
> # Trick, so we can use the printInfo utility functions:
> environment(myMoni) <- environment(printInfo[["default"]])
> 
> 
> ## in parallel
> 
> ## due to 'R CMD check --as-cran' allowing only <= 2 cores
> ## note: if doExtras, the check with '--as-cran' fails
> (doExtras <- simsalapar:::doExtras())
[1] FALSE
> (nc <- if(.Platform$OS.type == "windows") {
+            1 # otherwise win-builder fails
+        } else min(simsalapar:::nCores4test(), 2)) # otherwise CRAN fails
[1] 1
> 
> ## if simsalapar no longer *depends* on parallel:
> makeCluster <- parallel::makeCluster
> 
> S.T(resM.<- doMclapply    (vList, cores = nc, doOne=doOne, monitor=TRUE))
07:11:41 on CRANWIN3: forecaster=statistician; value:  num [1:2] 0.516 0.528
07:11:41 on CRANWIN3: forecaster=    optimist; value:  num [1:2] 20.55 4.47
07:11:41 on CRANWIN3: forecaster=   pessimist; value:  num [1:2] 0.792 0.499
   user  system elapsed 
   0.22    0.00    0.22 
> S.T(resM <- doMclapply    (vList, cores = nc, doOne=doOne, monitor=myMoni))
07:11:41 on CRANWIN3: forecaster=statistician; SE=    0.516378922
07:11:41 on CRANWIN3: forecaster=    optimist; SE=     20.5543582
07:11:41 on CRANWIN3: forecaster=   pessimist; SE=     0.79204769
   user  system elapsed 
   0.06    0.00    0.06 
> S.T(resF <- doForeach     (vList, cluster=makeCluster(nc, type="PSOCK"),
+                            doOne=doOne, monitor=printInfo[["fileEach"]]))
Loading required package: foreach
Loading required package: doParallel
Loading required package: iterators
Loading required package: parallel
getDoParWorkers(): 1
   user  system elapsed 
   0.08    0.03    2.43 
> ## For this problem, these are much slower on a typical 4-core desktop:
> S.T(resC <- doClusterApply(vList, cluster=makeCluster(nc, type="PSOCK"),
+                            doOne=doOne, monitor=printInfo[["gfile"]]))
   user  system elapsed 
   0.12    0.00    2.50 
> ## Note: Rmpi has a bug under openmpi-2.x (see https://stat.ethz.ch/pipermail/r-sig-hpc/2017-September/002063.html)
> if(FALSE) {
+     (hasRmpi <- .Platform$OS.type!="windows" # <- as Rmpi produces errors on the winbuilder
+         && require("Rmpi"))
+ }
> hasRmpi <- FALSE # for now...
> if(hasRmpi)
+     S.T(resR <- doRmpi(vList, nslaves=nc, doOne=doOne, monitor=TRUE))
> if(doExtras) { # use all available cores
+   print(S.T(resC.<- doClusterApply(vList, cluster=makeCluster(nc, type="PSOCK"),
+                                    doOne=doOne, monitor=myMoni)))
+   if(hasRmpi)
+       print(S.T(resR.<- doRmpi(vList, nslaves=nc, doOne=doOne, monitor=myMoni)))
+ }
> 
> ## check that we have the same results
> stopifnot(doRes.equal(resM, res),
+ 	  doRes.equal(resM, resM.),
+ 	  doRes.equal(resC, res),
+           if(hasRmpi) doRes.equal(resR, res) else TRUE,
+ 	  if(doExtras)
+ 	  doRes.equal(resC, resC.) &&
+ 	  (if(hasRmpi) doRes.equal(resR, resR.) else TRUE) else TRUE,
+ 	  doRes.equal(resF, res))
> 
> 
> ### 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="forecaster") # percentage of errors
> if(any(warn > 0))
+ ftable(100*warn, col.vars="forecaster") # percentage of warnings
> 
> ## 'reproduce' Table 4 of Gneiting (2011)
> tab1 <- apply(val, c("forecaster", "scoringfun"), mean)
> print(apply(tab1, 2, format, digits=3, scipen=-1), quote=FALSE)
              scoringfun
forecaster     SE     AE   
  statistician  0.516 0.528
  optimist     20.554 4.471
  pessimist     0.792 0.499
> 
> names(attributes(tab1)[["dim"]]) <- NULL # so that check below works
> ## not ok names(dim(tab1)) <- NULL # so that check below works
> 
> ## *without* mkAL (as check):
> res0 <- doLapply(vList, doAL=FALSE, doOne=doOne, monitor=FALSE)
> tab0 <- t(sapply(res0, `[[`, "value")) # values for the 3 forecasters: statistician, optimist, pessimist (see monitoring output)
> dimnames(tab0)[[1]] <- c("statistician", "optimist", "pessimist")
> names(dimnames(tab0)) <- c("forecaster", "scoringfun")
> stopifnot(identical(tab0, tab1)) # check
> 
> 
> ### Testing other demo()s  etc
> 
> if(doExtras) {
+     .checking <- TRUE # to run a smaller demo:
+     demo(robcovMCD)
+ }
> 
> proc.time()
   user  system elapsed 
   3.17    0.28    8.14