FoReco 0.2

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: FoReco
 Title: Point Forecast Reconciliation
-Version: 0.1.2
+Version: 0.2.0
 Authors@R: 
     c(person(given = "Daniele",
              family = "Girolimetto",
@@ -10,16 +10,17 @@ Authors@R:
       person(given = "Tommaso",
              family = "Di Fonzo",
              role = c("aut", "fnd")))
-Description: Provides classical (bottom-up), optimal and heuristic combination 
-    forecast reconciliation procedures for cross-sectional, temporal, and 
-    cross-temporal linearly constrained time series.
+Description:Classical (bottom-up and top-down), optimal and heuristic combination forecast 
+    reconciliation procedures for cross-sectional, temporal, and cross-temporal linearly 
+    constrained time series.
 License: GPL-3
 URL: https://github.com/daniGiro/FoReco,
     https://danigiro.github.io/FoReco/
 BugReports: https://github.com/daniGiro/FoReco/issues
 Depends: R (>= 2.10), Matrix, osqp, stats
-Imports: cli, corpcor, methods, pracma
+Imports: cli, corpcor, methods, mathjaxr
 Suggests: knitr, rmarkdown
+RdMacros: mathjaxr
 VignetteBuilder: knitr
 Encoding: UTF-8
 LazyData: true

LICENSE.md

@@ -553,7 +553,7 @@ and each file should have at least the “copyright” line and a pointer to
 where the full notice is found.
 
     <one line to give the program's name and a brief idea of what it does.>
-    Copyright (C) 2020 Tommaso Di Fonzo; Daniele Girolimetto
+    Copyright (C) 2020 Daniele Girolimetto; Tommaso Di Fonzo
 
     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
@@ -573,7 +573,7 @@ Also add information on how to contact you by electronic and paper mail.
 If the program does terminal interaction, make it output a short notice like this
 when it starts in an interactive mode:
 
-    FoReco Copyright (C) 2020 Tommaso Di Fonzo; Daniele Girolimetto
+    FoReco Copyright (C) 2020 Daniele Girolimetto; Tommaso Di Fonzo
     This program comes with ABSOLUTELY NO WARRANTY; for details type 'show w'.
     This is free software, and you are welcome to redistribute it
     under certain conditions; type 'show c' for details.

NAMESPACE

@@ -15,12 +15,17 @@ export(ctf_tools)
 export(hts_tools)
 export(htsrec)
 export(iterec)
+export(lccrec)
+export(oct_bounds)
 export(octrec)
 export(score_index)
 export(shrink_estim)
+export(srref)
 export(tcsrec)
+export(tdrec)
 export(thf_tools)
 export(thfrec)
+export(ut2c)
 import(Matrix)
 import(methods)
 import(osqp)

NEWS.md

@@ -1,22 +1,31 @@
-# FoReco 0.1.2
-**Note**: *Documentation for this version is still under development*
+# FoReco 0.2.0
 
 ##### Major changes
-* Add the possibility to use part of factors of *m* (max. order of temporal aggregation)
-* Add the possibility for `htsrec()`, `thfrec()` and `octrec()` to introduce a list of *h* covariance matrices in the parameters `W` and `Omega`, where *h* stands for the forecast horizon (note that for `thfrec()` and `octrec()` this is the forecast horizon of the entire cycle)
+* It's possible to use a subset of factors of *m* (max. order of temporal aggregation);
+* Added the possibility for `htsrec()`, `thfrec()` and `octrec()` to introduce a list of *h* covariance matrices in the parameters `W` and `Omega`, where *h* stands for the forecast horizon (note that for `thfrec()` and `octrec()` this is the forecast horizon of the entire cycle);
+* Param `Sstruc` is no more avaible in `octrec()` and `ctf_tools()`. FoReco uses a fast algorithm to compute **Scheck**, so no external input is needed;
+* Modified output of `ctf_tools()` (added `Ccheck`, `Htcheck`, `Scheck`, removed `Cstruc`, `Sstruc`), `hts_tools()` (added `C`) and `thf_tools()` (added `m`);
+* Added two new not negative reconciliation techniques ("**KAnn**" and "**sntz**") with a new parameter (`nn_type`) in `htsrec()`, `thfrec()` and `octrec()`;
+* Added the top-down reconciliation function `tdrec()`;
+* Added the level conditional forecast reconciliation (with and without not-negative constraints) for genuine hierarchical/grouped time series `levrec()` (cross-sectional, temporal and cross-temporal).
 
 ##### Minor changes
-* Now in `octrec()` it is also possible to introduce the **Ω** covariance matrix variant through the `Omega` parameter and not only the **W** variant with the`W` parameter
-* Renew `tcsrec()`, `cstrec()` and `iterec()`. Also in the iterec function the `maxit` parameter has been replaced by `itmax`, however for the moment `maxit` is still supported
+* Now in `octrec()` it is also possible to introduce the **Ω** covariance matrix variant through the `Omega` parameter and not only the **W** variant with the `W` parameter;
+* Updated `tcsrec()`, `cstrec()` and `iterec()`. In the iterec function the `maxit` parameter has been replaced by `itmax`, however for the moment `maxit` is still supported;
+* Now FoReco removes null rows from the cross-sectional aggregation matrix **C** and it warns the user if the balanced version of an unbalanced hierarchy is considering duplicated variables;
+* Redesigned the console output and added a new convergence norm as *default* for `iterec()` (`norm` parameter).
 
 ##### Experimental
-* Add the possibility to introduce constraints through the `bounds` param in `htsrec()`, `thfrec()` and `octrec()`.
-* Add a (hidden) function `oct_bounds()` to organize the bounds on a singular dimension (i.e. only cross-sectional or only temporal) in a cross-temporal framework
+* Add the possibility to introduce constraints through the `bounds` param in `htsrec()`, `thfrec()` and `octrec()`;
+* Add a function `oct_bounds()` to organize the bounds on a specific dimension (i.e. only cross-sectional or only temporal) in a cross-temporal framework;
+* Added `ut2c()` and `srref()` to develop a cross-sectional structural representation starting from a zero constraints kernel matrix;
+* Added in `score_index()` the calculation of multiple forecast horizons index (like 1:6) and multiple cross-sectional levels for a forecasting experiment.
+
 
 # FoReco 0.1.1
-This is a small release focusing on fixing some bugs and the documentation
+Minore release, fixing some bugs and the documentation
 
-* Fixed a bug in `iterec()` when calculating incoherence
+* Fixed a bug in `iterec()` when calculating the incoherence
 * Fixed documentation 
 * Changed the contact mail (now it's daniele.girolimetto@phd.unipd.it)
 * Corrected the second section of the vignette "`Average relative accuracy indices`"

R/Cmatrix.R

@@ -1,9 +1,10 @@
 #' Cross-sectional (contemporaneous) aggregation matrix
 #'
-#'
-#' This function allows the user to easily build the (\code{na x nb}) cross-sectional
-#' (contemporaneous) matrix mapping the \code{nb} bottom level series into the \code{na} higher level
-#' ones. (Experimental version)
+#' @description
+#' \loadmathjax
+#' This function allows the user to easily build the (\mjseqn{n_a \times n_b})
+#' cross-sectional (contemporaneous) matrix mapping the \mjseqn{n_b} bottom
+#' level series into the \mjseqn{n_a} higher level ones. (\emph{Experimental version})
 #'
 #' @param formula  Specification of the hierarchical structure: grouped hierarchies are specified
 #' using \code{~ g1 * g2} and nested hierarchies are specified using \code{~ parent / child}.
@@ -16,6 +17,9 @@
 #'
 #' @return A (\code{na x nb}) matrix.
 #'
+#' @keywords utilities
+#' @family utilities
+#'
 #' @examples
 #' ## Balanced hierarchy
 #' #         T

R/Dmat.R

@@ -1,17 +1,48 @@
 # @title Maps a vectorized matrix of the reconciled forecasts into a differently organized vector
 #
-# @description This function returns the [\code{hn(ks+m) x hn(ks+m)}]
-# permutation matrix transforming vec(Y') into vec(Y^*) (see notaH.pdf)
+# @description
+# \loadmathjax
+# This function returns the [\mjseqn{hn(k^\ast+m) \times hn(k^\ast+m)}]
+# permutation matrix transforming \mjseqn{\mbox{vec}(\mathbf{Y}')} into
+# \mjseqn{\mbox{vec}(\mathbf{Y}_h')}:
+# \mjsdeqn{\mathbf{D}_h \mbox{vec}(\mathbf{Y}') = \mbox{vec}(\mathbf{Y}_h')}
+# where \mjseqn{\mathbf{Y}_h'} is a
+# \mjseqn{h \times n(k^\ast+m)} with column ordered as [lowest_freq' ...  highest_freq']
+# repeat for \mjseqn{n} series.
 #
 # @param h forecast horizon for the lowest frequency (most temporally agregated)
 # time series;
-# @param kset set of p factors of \code{m}, from \code{m} to 1;
-# @param n number of variables (\code{n} = \code{n_a} + \code{n_b}).
+# @param m Highest available sampling frequency per seasonal cycle (max. order
+# of temporal aggregation, \mjseqn{m}), or a subset of \mjseqn{p} factors
+# of \mjseqn{m}.
+# @param n number of the cross-sectional variables (\code{n} = \code{n_a} + \code{n_b}).
 #
-# @return a sparse [\code{hn(ks+m) x hn(ks+m)}] matrix
+# @return A sparse [\code{hn(ks+m) x hn(ks+m)}] matrix.
 #
-Dmat <- function(h, kset, n){
-  m <- max(kset)
+# @examples
+# data(FoReco_data)
+# # An example in the temporal frameworks
+# # top ts residuals ([h*lowest_freq' ...  h*highest_freq']')
+# topres <- FoReco_data$res[1, ]
+# D1 <- FoReco:::Dmat(m = c(1,2,3,4,6,12), n = 1, h = 14)
+# Yh_t <- matrix(D1%*%topres, nrow = 14, byrow = TRUE)
+# # check
+# all(as.vector(D1%*%topres) == as.vector(t(Yh_t)))
+#
+# # An example in the cross-temporal frameworks
+# Dn <- FoReco:::Dmat(m = c(1,2,3,4,6,12), n = 8, h = 14)
+# Yh_ct <- matrix(Dn%*%as.vector(t(FoReco_data$res)), nrow = 14, byrow = TRUE)
+# # check
+# all(Dn%*%as.vector(t(FoReco_data$res)) == as.vector(t(Yh_ct)))
+#
+Dmat <- function(h, m, n){
+  if(length(m)==1){
+    kset <- 1
+  }else{
+    kset <- sort(m, decreasing = TRUE)
+    m <- max(m)
+  }
+
   I <- .sparseDiagonal(h * sum(m/kset) * n)
   i <- rep(rep(rep(1:h, length(kset)), rep(m/kset, each = h)), n)
   D <- I[order(i), ]

R/FoReco-hts.R

@@ -49,28 +49,33 @@
 #' }
 #' colnames(res) <- colnames(data)
 #'
+#' ## Comparisons
 #' # ols
+#' # two commands in hts...
 #' Y_hts_forecast <- forecast(htseg1, method = "comb", fmethod = "arima", weights = "ols")
 #' Y_hts_ols <- combinef(BASE, nodes = get_nodes(htseg1), keep = "all")
-#' sum((allts(Y_hts_forecast) - Y_hts_ols) > 1e-10)
+#' # ...with the same results:
+#' sum(abs(allts(Y_hts_forecast) - Y_hts_ols) > 1e-10)
+#'
 #' Y_FoReco_ols <- htsrec(BASE, C = C, comb = "ols")$recf
-#' sum((Y_hts_ols - Y_FoReco_ols) > 1e-10)
+#' sum(abs(Y_hts_ols - Y_FoReco_ols) > 1e-10)
 #'
 #' # struc
 #' w <- 1 / apply(smatrix(htseg1), 1, sum)
 #' Y_hts_struc <- combinef(BASE, nodes = get_nodes(htseg1), weights = w, keep = "all")
 #' Y_FoReco_struc <- htsrec(BASE, C = C, comb = "struc")$recf
-#' sum((Y_hts_struc - Y_FoReco_struc) > 1e-10)
+#' sum(abs(Y_hts_struc - Y_FoReco_struc) > 1e-10)
 #'
 #' # shr
-#' Y_hts_shr <- MinT(BASE, nodes = get_nodes(htseg1), keep = "all", covariance = "shr", residual = res)
+#' Y_hts_shr <- MinT(BASE, nodes = get_nodes(htseg1), keep = "all",
+#'                   covariance = "shr", residual = res)
 #' Y_FoReco_shr <- htsrec(BASE, C = C, comb = "shr", res = res)$recf
-#' sum((Y_hts_shr - Y_FoReco_shr) > 1e-10)
+#' sum(abs(Y_hts_shr - Y_FoReco_shr) > 1e-10)
 #'
-#' # sam - hts error "MinT needs covariance matrix to be positive definite"
-#' #       The covariance matrix is non-singular, but its condition number is very low,
-#' #       and hts considers it as non-invertible
-#' Y_hts_sam <- MinT(BASE, nodes = get_nodes(htseg1), keep = "all", covariance = "sam", residual = res)
+#' # sam - hts error "MinT needs covariance matrix to be positive definite."
+#' #       The covariance matrix is ill-conditioned, hts considers it as non-invertible
+#' Y_hts_sam <- MinT(BASE, nodes = get_nodes(htseg1), keep = "all",
+#'                   covariance = "sam", residual = res)
 #' Y_FoReco_sam <- htsrec(BASE, C = C, comb = "sam", res = res)$recf
 #' # sum((Y_hts_sam-Y_FoReco_sam)>1e-10)
 #'
@@ -81,7 +86,8 @@
 #' na <- n-nb
 #' C <- smatrix(htseg2)[1:na, ]
 #'
-#' ## In FoReco, forecasts must be obtained externally
+#' ## Computation of the base forecasts
+#' # using the auto.arima() function of the package forecast (loaded by hts)
 #' # List containing the base forecasts
 #' # Forecast horizon: 10
 #' base <- list()
@@ -103,11 +109,9 @@
 #' }
 #' colnames(res) <- colnames(data)
 #'
-#' ## Comparison
+#' ## Comparisons
 #' # ols
-#' Y_hts_forecast <- forecast(htseg2, method = "comb", fmethod = "arima", weights = "ols")
 #' Y_hts_ols <- combinef(BASE, nodes = get_nodes(htseg2), keep = "all")
-#' sum((allts(Y_hts_forecast) - Y_hts_ols) > 1e-10)
 #' Y_FoReco_ols <- htsrec(BASE, C = C, comb = "ols")$recf
 #' sum(abs(Y_hts_ols - Y_FoReco_ols) > 1e-10)
 #'
@@ -123,3 +127,4 @@
 #' sum(abs(Y_hts_shr- Y_FoReco_shr) > 1e-10)
 #' }
 NULL
+

R/FoReco-package.R

@@ -1,20 +1,23 @@
 #' FoReco: point forecast reconciliation
 #'
-#' An R package offering classical (bottom-up), optimal and heuristic combination
+#' An R package offering classical (bottom-up and top-down), and modern (optimal and heuristic combination)
 #' forecast reconciliation procedures for cross-sectional, temporal, and cross-temporal
 #' linearly constrained time series.
 #'
 #' @author Tommaso Di Fonzo and Daniele Girolimetto, Department of Statistical Sciences, University of Padua (Italy).
 #'
 #' @details
 #' The \code{FoReco} package is designed for point forecast reconciliation, a
-#' post-forecasting process aimed to improve the quality of the base
+#' post-forecasting process aimed to improve the accuracy of the base
 #' forecasts for a system of linearly constrained (e.g. hierarchical/grouped) time series.
 #' The main functions are:
 #'
 #' \describe{
 #'   \item{\code{\link[FoReco]{htsrec}():}}{cross-sectional (contemporaneous) forecast reconciliation.}
 #'   \item{\code{\link[FoReco]{thfrec}():}}{forecast reconciliation for a single time series through temporal hierarchies.}
+#'   \item{\code{\link[FoReco]{lccrec}():}}{level conditional forecast reconciliation for genuine hierarchical/grouped time series.}
+#'   \item{\code{\link[FoReco]{tdrec}():}}{top-down (cross-sectional, temporal, cross-temporal) forecast reconciliation for genuine hierarchical/grouped time series.}
+#'   \item{\code{\link[FoReco]{ctbu}():}}{bottom-up cross-temporal forecast reconciliation.}
 #'   \item{\code{\link[FoReco]{tcsrec}():}}{heuristic first-temporal-then-cross-sectional cross-temporal forecast reconciliation.}
 #'   \item{\code{\link[FoReco]{cstrec}():}}{heuristic first-cross-sectional-then-temporal cross-temporal forecast reconciliation.}
 #'   \item{\code{\link[FoReco]{iterec}():}}{heuristic iterative cross-temporal forecast reconciliation.}
@@ -26,6 +29,8 @@
 #' Optimal Combination Method and Heuristic Alternatives, Department of Statistical
 #' Sciences, University of Padua, \href{https://arxiv.org/abs/2006.08570}{arXiv:2006.08570}.
 #'
+#' Di Fonzo, T., Girolimetto, D. (2021), \emph{Forecast combination based forecast reconciliation: insights and extensions} (mimeo).
+#'
 #' @docType package
 #' @name FoReco-package
 #' @keywords package
@@ -36,4 +41,3 @@ NULL
     right = paste0("FoReco ", utils::packageVersion("FoReco"))
   ))
 }
-

R/FoReco-thief.R

@@ -37,6 +37,7 @@
 #' }
 #'
 #' # OLS
+#' # two commands in thief...
 #' obj_thief <- thief(dataset, m = 52, h = 2 * 52, comb = "ols", usemodel = "arima")
 #' obj_thief <- tsaggregates(obj_thief$mean)
 #' y_thief <- NULL
@@ -48,7 +49,9 @@
 #' for (i in 6:1) {
 #'   y_thief_ols <- c(y_thief_ols, obj_thief_ols[[i]]$mean)
 #' }
+#' # ...with the same results:
 #' sum(abs(y_thief_ols - y_thief) > 1e-10)
+#'
 #' y_FoReco_ols <- thfrec(base_vec, 52, comb = "ols")$recf
 #' sum(abs(y_FoReco_ols - y_thief_ols) > 1e-10)
 #'

R/FoReco2ts.R

@@ -1,18 +1,19 @@
 #' Reconciled forecasts matrix/vector to time-series class
 #'
-#' Function to transform the matrix/vector of reconciled forecasts (\code{recf} from
-#' \link[FoReco]{htsrec}, \link[FoReco]{thfrec}, \link[FoReco]{octrec},
-#' \link[FoReco]{tcsrec}, \link[FoReco]{cstrec}, \link[FoReco]{iterec},
-#' \link[FoReco]{ctbu}) into a list of time series objects.
+#' @description
+#' \loadmathjax
+#' Function to transform the matrix/vector of reconciled forecasts
+#' (\code{recf} from \link[FoReco]{htsrec}, \link[FoReco]{thfrec}, \link[FoReco]{tdrec},
+#' \link[FoReco]{octrec}, \link[FoReco]{lccrec}, \link[FoReco]{tcsrec}, \link[FoReco]{cstrec},
+#' \link[FoReco]{iterec}, \link[FoReco]{ctbu}) into a list of time series objects.
 #'
-#' @param recf (\code{h(k* + m) x 1}) reconciled forecasts vector from \link[FoReco]{thfrec},
-#' (\code{h x n}) reconciled forecasts matrix from \link[FoReco]{htsrec} or
-#' (\code{n x h(k* + m)}) reconciled forecasts matrix from \link[FoReco]{octrec},
+#' @param recf (\mjseqn{h(k^\ast + m) \times 1}) reconciled forecasts vector from \link[FoReco]{thfrec},
+#' (\mjseqn{h \times n}) reconciled forecasts matrix from \link[FoReco]{htsrec} or
+#' (\mjseqn{n \times h(k^\ast + m)}) reconciled forecasts matrix from \link[FoReco]{octrec},
 #' \link[FoReco]{tcsrec}, \link[FoReco]{cstrec}, \link[FoReco]{iterec},
 #' \link[FoReco]{ctbu}.
 #' @param ... optional arguments to \link[stats]{ts} (i.e. starting date);
-#' frequency is only required for the cross-sectional system.
-#' .
+#' frequency is required only for the cross-sectional case.
 #'
 #' @return
 #' A list of class \code{"ts"} objects
@@ -39,12 +40,13 @@
 #' # Temporal framework
 #' # top ts base forecasts ([lowest_freq' ...  highest_freq']')
 #' topbase <- FoReco_data$base[1, ]
-#'  # top ts residuals ([lowest_freq' ...  highest_freq']')
+#' # top ts residuals ([lowest_freq' ...  highest_freq']')
 #' topres <- FoReco_data$res[1, ]
 #' thf_recf <- thfrec(topbase, m = 12, comb = "acov", res = topres)$recf
 #' obj_thf <- FoReco2ts(recf = thf_recf, start = c(15, 1))
 #'
 #' @keywords utilities
+#' @family utilities
 #'
 #' @export
 FoReco2ts <- function(recf, ...) {