Merge pull request #8 from Techtonique/dov070

Do version 0.7.0

Author: thierrymoudiki

DESCRIPTION

@@ -1,13 +1,14 @@
 Package: ahead
 Type: Package
 Title: Time Series Forecasting
-Version: 0.6.2
-Date: 2023-08-08
+Version: 0.7.0
+Date: 2023-08-23
 Author: T. Moudiki
 Maintainer: T. Moudiki <thierry.moudiki@gmail.com>
 Description: Univariate and multivariate time series forecasting.
 License: BSD_3_clause Clear + file LICENSE
-Imports: Rcpp (>= 1.0.6), forecast, vars, randtoolbox, fGarch
+Imports: Rcpp (>= 1.0.6) 
+Depends: forecast, vars, randtoolbox, fGarch, VineCopula
 LinkingTo: Rcpp
 RoxygenNote: 7.2.2
 Encoding: UTF-8

NAMESPACE

@@ -1,13 +1,29 @@
 useDynLib(ahead, .registration=TRUE)
 S3method(plot, mtsforecast)
+S3method(plot, pairs)
 S3method(predict, ridge)
 importFrom(Rcpp, evalCpp)
-importFrom("stats", "as.ts", "complete.cases", "end", "frequency",
-               "median", "predict", "predict.lm", "qnorm", "quantile",
-               "residuals", "rnorm", "rt", "sd", "start", "ts", "tsp",
-               "tsp<-")
+importFrom(VineCopula, C2RVine)
+importFrom(VineCopula, D2RVine)
+importFrom(VineCopula, RVineSim)
+importFrom(VineCopula, RVineStructureSelect)
+importFrom(forecast, ets)
+importFrom(forecast, auto.arima)
+importFrom(forecast, thetaf)
+importFrom(vars, VAR)
+importFrom(randtoolbox, sobol)
+importFrom(randtoolbox, halton)
+importFrom(fGarch, garchFit)
+importFrom("grDevices", "colorRampPalette")
+  importFrom("graphics", "abline", "par", "points", "polygon")
+  importFrom("stats", "t.test")
+importFrom("stats", "as.ts", "complete.cases", "ecdf", "end", "frequency",
+           "lag", "median", "predict", "predict.lm", "qnorm", "quantile",
+           "residuals", "rnorm", "rt", "sd", "setNames", "start",
+           "ts", "tsp", "tsp<-")
 importFrom("utils", "tail")
 importFrom("graphics", "lines")
+importFrom("graphics", "matplot")
 importFrom("stats", "is.ts", "runif", "time")
 export("armagarchf")
 export("basicf")
@@ -17,3 +33,5 @@ export("loessf")
 export("ridge2f")
 export("ridge")
 export("varf")
+export("getreturns")
+export("getsimulations")

NEWS.md

@@ -1,3 +1,14 @@
+# version 0.7.0
+
+- moving block bootstrap in `ridge2f`, `basicf` and `loessf`, in addition to circular block bootstrap from 0.6.2
+- adjust R-Vine copulas on residuals for `ridge2f` simulation
+- new plots for simulations see (new) vignettes
+- split conformal prediction intervals (**very very experimental** and basic right now, too conservative)
+- `Depends` and selective `Imports` (beneficial to Python and rpy2 for installation time?)
+- `getsimulations` extracts simulations from a given time series (from `ridge2f` and `basicf`)
+- `getreturns` extracts returns/log-returns from multivariate time series
+- `splitts` splits time series using a proportion of data
+
 # version 0.6.2
 
 - Add Block Bootstrap to `ridge2f`

R/RcppExports.R

@@ -5,10 +5,6 @@ embedc <- function(x, lags) {
     .Call(`_ahead_embedc`, x, lags)
 }
 
-rcpp_hello_world <- function() {
-    .Call(`_ahead_rcpp_hello_world`)
-}
-
 forecast_innovs_loop_cpp <- function(eps, rts, eps_prev, omega, alpha, beta, df, h) {
     .Call(`_ahead_forecast_innovs_loop_cpp`, eps, rts, eps_prev, omega, alpha, beta, df, h)
 }

R/basic.R

@@ -6,8 +6,9 @@
 #' @param h Forecasting horizon
 #' @param level Confidence level for prediction intervals
 #' @param method forecasting method, either "mean", "median", or random walk ("rw")
-#' @param type_pi type of prediction interval currently, "gaussian" or "bootstrap"
-#' @param block_length length of block for circular block bootstrap (\code{type_pi == 'blockbootstrap'})
+#' @param type_pi type of prediction interval currently, "gaussian", "bootstrap",
+#' "blockbootstrap" or "movingblockbootstrap"
+#' @param block_length length of block for (circular) "blockbootstrap" or "movingblockbootstrap"
 #' @param seed reproducibility seed for \code{type_pi == 'bootstrap'}
 #' @param B Number of bootstrap replications for \code{type_pi == 'bootstrap'}
 #'
@@ -52,11 +53,24 @@
 #' plot(res5, "Quotes")
 #' plot(res5, "TV.advert")
 #'
+#'
+#' # moving block bootstrap
+#' res6 <- ahead::basicf(fpp::insurance, h=10,
+#'                       type_pi = "movingblockbootstrap", B=10,
+#'                       block_length = 4, method = "rw")
+#'
+#' par(mfrow=c(1, 2))
+#' plot(res6, "Quotes")
+#' plot(res6, "TV.advert")
+#'
 basicf <- function(y,
                    h = 5,
                    level = 95,
                    method = c("mean", "median", "rw"),
-                   type_pi = c("gaussian", "bootstrap", "blockbootstrap"),
+                   type_pi = c("gaussian",
+                               "bootstrap",
+                               "blockbootstrap",
+                               "movingblockbootstrap"),
                    block_length = NULL,
                    seed = 1,
                    B = 100)
@@ -90,6 +104,7 @@ basicf <- function(y,
   fits <- ts(t(replicate(n = n_inputs, expr = point_forecast)),
              start = start_x, frequency = freq_x)
   resids <- y - fits
+  colnames(resids) <- colnames(y)
 
   # out-of-sample
   fcast <- ts(t(replicate(n = h, expr = point_forecast)),
@@ -117,24 +132,45 @@ basicf <- function(y,
     return(structure(out, class = "mtsforecast"))
   }
 
-  if (type_pi %in% c("bootstrap", "blockbootstrap"))
+  if(type_pi %in% c("blockbootstrap", "movingblockbootstrap"))
+  {
+    if (nrow(y) <= 2 * freq_x)
+      freq_x <- 1L
+
+    if (is.null(block_length)) {
+      block_length <- ifelse(freq_x > 1, 2 * freq_x, min(8, floor(nrow(y) / 2)))
+    }
+  }
+
+  if (type_pi %in% c("bootstrap", "blockbootstrap", "movingblockbootstrap"))
   {
     sims <- vector("list", length = B)
     for (i in 1:B)
     {
       # sampling from the residuals
       set.seed(seed + i*100 + nchar(method))
 
-      idx <- switch(type_pi,
-                    bootstrap = sample.int(n = nrow(resids),
-                        size = h, replace = TRUE),
-                    blockbootstrap = mbb(
-                      r = resids,
-                      n = h,
-                      b = block_length,
-                      return_indices = TRUE,
-                      seed = seed + i*100 + nchar(method)
-                    )
+      idx <- switch(
+        type_pi,
+        bootstrap = sample.int(
+          n = nrow(resids),
+          size = h,
+          replace = TRUE
+        ),
+        blockbootstrap = mbb(
+          r = resids,
+          n = h,
+          b = block_length,
+          return_indices = TRUE,
+          seed = seed + i * 100 + nchar(method)
+        ),
+        movingblockbootstrap = mbb2(
+          r = resids,
+          n = h,
+          b = block_length,
+          return_indices = TRUE,
+          seed = seed + i * 100 + nchar(method)
+        )
       )
 
       sims[[i]] <- ts(as.matrix(fcast) + as.matrix(resids[idx, ]),

R/getreturns.R

@@ -0,0 +1,37 @@
+#' Calculate returns or log-returns for multivariate time series
+#'
+#' @param x Multivariate time series
+#' @param type Type of return: basic return ("basic") or log-return ("log")
+#'
+#' @return The returns
+#' @export
+#'
+#' @examples
+#'
+#' returns <- getreturns(EuStockMarkets)
+#' log_returns <- getreturns(EuStockMarkets,
+#'                           type = "log")
+#'
+#' par(mfrow=c(1, 3))
+#' matplot(EuStockMarkets, type = 'l', main = "Closing Prices of \n European stocks (1991-1998)",
+#' xlab = "time")
+#' matplot(returns, type = 'l', main = "Returns", xlab = "time")
+#' matplot(log_returns, type = 'l', main = "Log-returns", xlab = "time")
+#'
+getreturns <- function(x, type = c("basic", "log"))
+{
+  # ?rlang::abort #?
+  stopifnot(is.ts(x))
+  names_x <- colnames(x)
+  type <- match.arg(type)
+
+  if (identical(type, "basic"))
+  {
+    res <- diff(x)/lag(x)
+    colnames(res) <- names_x
+    return (res)
+  }
+  res <- diff(log(x))
+  colnames(res) <- names_x
+  return (res)
+}

R/getsimulations.R

@@ -0,0 +1,42 @@
+#' Obtain simulations (when relevant) from a selected time series
+#'
+#' @param obj result from ridge2f (multivariate time series forecast with simulations)
+#' @param selected_series name of the time series selected
+#' @param transpose return a transposed time series
+#'
+#' @export
+#'
+#' @examples
+#'
+#' require(fpp)
+#'
+#' obj <- ahead::ridge2f(fpp::insurance, h = 7,
+#'                       type_pi = "bootstrap", B = 5)
+#' print(getsimulations(obj, selected_series = "TV.advert"))
+#' print(getsimulations(obj, selected_series = "Quotes"))
+#' print(getsimulations(obj, selected_series = "TV.advert", transpose = TRUE))
+#' print(getsimulations(obj, selected_series = "Quotes", transpose = TRUE))
+#'
+getsimulations <- function(obj, selected_series, transpose = FALSE)
+{
+  n_sims <- length(obj$sims)
+  start_preds <- stats::start(obj$mean)
+  frequency_preds <- stats::frequency(obj$mean)
+  selected_series_sims <- sapply(1:n_sims,
+                                 FUN = function(i) obj$sims[[i]][, selected_series])
+
+  temp <- ts(selected_series_sims,
+             start = start_preds,
+             frequency = frequency_preds)
+
+  if (transpose)
+  {
+    temp <- t(temp)
+    colnames(temp) <- paste0("date", 1:nrow(obj$mean))
+    return(list(series = temp,
+                name = selected_series))
+  } else {
+    return(list(series = temp,
+                name = selected_series))
+  }
+}

R/loess.R

@@ -7,7 +7,7 @@
 #' @param level Confidence level for prediction intervals
 #' @param span the parameter which controls the degree of smoothing
 #' @param degree the degree of the polynomials to be used, normally 1 or 2. (Degree 0 is also allowed, but see \code{stats::loess})
-#' @param type_pi Type of prediction interval currently (independent) "bootstrap" or (circular) "blockbootstrap"
+#' @param type_pi Type of prediction interval currently (independent) "bootstrap", (circular) "blockbootstrap", or "movingblockbootstrap"
 #' @param b block length for circular block bootstrap
 #' @param B number of bootstrap replications
 #' @param type_aggregation Type of aggregation, ONLY for bootstrapping; either "mean" or "median"
@@ -28,20 +28,29 @@
 #'
 #' @examples
 #'
+#' par(mfrow = c(3, 1))
+#'
 #' plot(loessf(Nile, h=20, level=95, B=10))
 #'
+#' plot(loessf(Nile, h=20, level=95, B=10,
+#'      type_pi = "blockbootstrap"))
+#'
+#' plot(loessf(Nile, h=20, level=95, B=10,
+#'      type_pi = "movingblockbootstrap"))
+#'
 loessf <- function(y,
                    h = 5,
                    level = 95,
                    span = 0.75,
                    degree = 2,
-                   type_pi = c("bootstrap", "blockbootstrap"),
+                   type_pi = c("bootstrap",
+                               "blockbootstrap",
+                               "movingblockbootstrap"),
                    b = NULL,
                    B = 250,
                    type_aggregation = c("mean", "median"),
                    seed = 123)
 {
-  # adapted from forecast:::bld.mbb.bootstrap ---
   freq_y <- frequency(y)
   if (length(y) <= 2 * freq_y)
     freq_y <- 1L
@@ -100,6 +109,7 @@ loessf <- function(y,
     idx <- sample.int(n = length(resids),
                       size = h,
                       replace = TRUE)
+
     bootstrapped_residuals <- ts(switch(
       type_pi,
       bootstrap = matrix(resids, ncol = 1L)[idx,],
@@ -111,9 +121,17 @@ loessf <- function(y,
           seed = 100 * i + 3,
           return_indices =
             FALSE
-        )
-      )
-    ),
+        )),
+        movingblockbootstrap = drop(
+          mbb2(
+            r = matrix(resids, ncol = 1L),
+            n = h,
+            b = b,
+            seed = 100 * i + 3,
+            return_indices =
+              FALSE
+          ))
+      ),
     start = start_preds,
     frequency = freq_y)