From 9ace22dd39e20a736600dfafcf29a89f89f90aca Mon Sep 17 00:00:00 2001 From: dani Date: Fri, 21 May 2021 11:30:17 +0200 Subject: [PATCH] FoReco 0.2 --- DESCRIPTION | 11 +- LICENSE.md | 4 +- NAMESPACE | 5 + NEWS.md | 29 +- R/Cmatrix.R | 12 +- R/Dmat.R | 45 +- R/FoReco-hts.R | 31 +- R/FoReco-package.R | 10 +- R/FoReco-thief.R | 3 + R/FoReco2ts.R | 22 +- R/FoReco_data.R | 16 +- R/commat.R | 20 +- R/cstrec.R | 106 +- R/ctbu.R | 50 +- R/ctf_tools.R | 160 +-- R/hts_tools.R | 109 +- R/htsrec.R | 359 ++++-- R/iterec.R | 240 ++-- R/lccrec.R | 1044 +++++++++++++++++ R/oct_bounds.R | 151 ++- R/octrec.R | 401 ++++--- R/reco.R | 671 +++++++---- R/score_index.R | 172 ++- R/shrink_estim.R | 1 + R/srref.R | 108 ++ R/tcsrec.R | 116 +- R/tdrec.R | 170 +++ R/thf_tools.R | 11 +- R/thfrec.R | 245 ++-- R/ut2c.R | 169 +++ README.Rmd | 15 +- README.md | 33 +- _pkgdown.yml | 6 + cran-comments.md | 11 + docs/404.html | 14 +- docs/CODE_OF_CONDUCT.html | 14 +- docs/LICENSE.html | 18 +- docs/articles/FoReco_package.html | 180 +-- .../anchor-sections-1.0/anchor-sections.css | 4 + .../anchor-sections-1.0/anchor-sections.js | 33 + .../header-attrs-2.5/header-attrs.js | 12 + .../header-attrs-2.7/header-attrs.js | 12 + docs/articles/accuracy_indices.html | 208 +++- .../anchor-sections-1.0/anchor-sections.css | 4 + .../anchor-sections-1.0/anchor-sections.js | 33 + .../header-attrs-2.5/header-attrs.js | 12 + .../header-attrs-2.7/header-attrs.js | 12 + docs/articles/index.html | 14 +- docs/authors.html | 16 +- docs/index.html | 42 +- docs/man/figures/ite.png | Bin 73324 -> 69762 bytes docs/news/index.html | 46 +- docs/pkgdown.yml | 4 +- docs/reference/Cmatrix.html | 43 +- docs/reference/Dmat.html | 278 +++++ docs/reference/FoReco-hts.html | 44 +- docs/reference/FoReco-package.html | 24 +- docs/reference/FoReco-thief.html | 17 +- docs/reference/FoReco2ts.html | 58 +- docs/reference/FoReco_data.html | 41 +- docs/reference/commat.html | 52 +- docs/reference/cstrec.html | 120 +- docs/reference/ctbu.html | 81 +- docs/reference/ctf_tools.html | 102 +- docs/reference/figures/ite.png | Bin 73324 -> 69762 bytes docs/reference/hts_tools.html | 68 +- docs/reference/htsrec.html | 269 ++++- docs/reference/index.html | 51 +- docs/reference/iterec-1.png | Bin 52944 -> 49826 bytes docs/reference/iterec.html | 136 ++- docs/reference/lccrec.html | 538 +++++++++ docs/reference/levrec.html | 420 +++++++ docs/reference/oct_bounds.html | 301 +++++ docs/reference/octrec.html | 301 +++-- docs/reference/score_index.html | 85 +- docs/reference/shrink_estim.html | 27 +- docs/reference/srref.html | 276 +++++ docs/reference/tcsrec.html | 124 +- docs/reference/tdrec.html | 338 ++++++ docs/reference/thf_tools.html | 29 +- docs/reference/thfrec.html | 224 +++- docs/reference/ut2c.html | 352 ++++++ docs/sitemap.xml | 15 + inst/CITATION | 2 +- man/Cmatrix.Rd | 22 +- man/FoReco-hts.Rd | 30 +- man/FoReco-package.Rd | 9 +- man/FoReco-thief.Rd | 3 + man/FoReco2ts.Rd | 34 +- man/FoReco_data.Rd | 16 +- man/commat.Rd | 32 +- man/cstrec.Rd | 90 +- man/ctbu.Rd | 54 +- man/ctf_tools.Rd | 78 +- man/figures/ite.png | Bin 73324 -> 69762 bytes man/hts_tools.Rd | 49 +- man/htsrec.Rd | 266 ++++- man/iterec.Rd | 116 +- man/lccrec.Rd | 296 +++++ man/oct_bounds.Rd | 74 ++ man/octrec.Rd | 278 +++-- man/score_index.Rd | 54 +- man/shrink_estim.Rd | 14 + man/srref.Rd | 50 + man/tcsrec.Rd | 92 +- man/tdrec.Rd | 104 ++ man/thf_tools.Rd | 16 + man/thfrec.Rd | 214 +++- man/ut2c.Rd | 123 ++ vignettes/FoReco_package.Rmd | 41 +- vignettes/accuracy_indices.Rmd | 41 + 111 files changed, 9588 insertions(+), 2258 deletions(-) mode change 100755 => 100644 NAMESPACE create mode 100644 R/lccrec.R create mode 100644 R/srref.R create mode 100644 R/tdrec.R create mode 100644 R/ut2c.R create mode 100644 cran-comments.md create mode 100644 docs/articles/FoReco_package_files/anchor-sections-1.0/anchor-sections.css create mode 100644 docs/articles/FoReco_package_files/anchor-sections-1.0/anchor-sections.js create mode 100644 docs/articles/FoReco_package_files/header-attrs-2.5/header-attrs.js create mode 100644 docs/articles/FoReco_package_files/header-attrs-2.7/header-attrs.js create mode 100644 docs/articles/accuracy_indices_files/anchor-sections-1.0/anchor-sections.css create mode 100644 docs/articles/accuracy_indices_files/anchor-sections-1.0/anchor-sections.js create mode 100644 docs/articles/accuracy_indices_files/header-attrs-2.5/header-attrs.js create mode 100644 docs/articles/accuracy_indices_files/header-attrs-2.7/header-attrs.js create mode 100644 docs/reference/Dmat.html create mode 100644 docs/reference/lccrec.html create mode 100644 docs/reference/levrec.html create mode 100644 docs/reference/oct_bounds.html create mode 100644 docs/reference/srref.html create mode 100644 docs/reference/tdrec.html create mode 100644 docs/reference/ut2c.html create mode 100644 man/lccrec.Rd create mode 100644 man/oct_bounds.Rd create mode 100644 man/srref.Rd create mode 100644 man/tdrec.Rd create mode 100644 man/ut2c.Rd diff --git a/DESCRIPTION b/DESCRIPTION index be5e8bb..6022df1 100755 --- a/DESCRIPTION +++ b/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 diff --git a/LICENSE.md b/LICENSE.md index 65bd6cd..f8435af 100755 --- a/LICENSE.md +++ b/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. - 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. diff --git a/NAMESPACE b/NAMESPACE old mode 100755 new mode 100644 index 0392e65..75f4760 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index 6546dea..4b3e01e 100755 --- a/NEWS.md +++ b/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`" diff --git a/R/Cmatrix.R b/R/Cmatrix.R index 4454c6b..8056a96 100755 --- a/R/Cmatrix.R +++ b/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 diff --git a/R/Dmat.R b/R/Dmat.R index c12b406..814dce1 100755 --- a/R/Dmat.R +++ b/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), ] diff --git a/R/FoReco-hts.R b/R/FoReco-hts.R index e49a9c7..4ddde2c 100755 --- a/R/FoReco-hts.R +++ b/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 + diff --git a/R/FoReco-package.R b/R/FoReco-package.R index 803a842..fad5e22 100755 --- a/R/FoReco-package.R +++ b/R/FoReco-package.R @@ -1,6 +1,6 @@ #' 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. #' @@ -8,13 +8,16 @@ #' #' @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")) )) } - diff --git a/R/FoReco-thief.R b/R/FoReco-thief.R index 2825c39..99557ba 100755 --- a/R/FoReco-thief.R +++ b/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) #' diff --git a/R/FoReco2ts.R b/R/FoReco2ts.R index 7fcfb8e..1eaad64 100755 --- a/R/FoReco2ts.R +++ b/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, ...) { diff --git a/R/FoReco_data.R b/R/FoReco_data.R index 708125b..12495f3 100755 --- a/R/FoReco_data.R +++ b/R/FoReco_data.R @@ -1,10 +1,11 @@ -#' Forecast reconciliation for a simulated hierarchical time series +#' Forecast reconciliation for a simulated linearly constrained, genuine hierarchical multiple time series #' #' A two-level hierarchy with \code{n = 8} monthly time series. In the cross-sectional framework, -#' at any time it is \code{Tot = A + B + C}, \code{A = AA + AB} and \code{B = BA + BB} (\code{nb = 5} at the -#' bottom level). For monthly data, the observations are aggregated to annual (\code{k = 12}), +#' at any time it is \code{Tot = A + B + C}, \code{A = AA + AB} and \code{B = BA + BB} +#' (the bottom time series being \code{AA}, \code{AB}, \code{BA}, \code{BB}, and \code{C}, it is \code{nb = 5}). +#' The monthly observations are aggregated to their annual (\code{k = 12}), #' semi-annual (\code{k = 6}), four-monthly (\code{k = 4}), quarterly (\code{k = 3}), and -#' bi-monthly (\code{k = 2}) observations. The monthly bottom time series are simulated +#' bi-monthly (\code{k = 2}) counterparts. The monthly bottom time series are simulated #' from five different SARIMA models (see #' \href{https://danigiro.github.io/FoReco/articles/FoReco_package.html}{\code{Using the `FoReco` package}}). #' There are 180 (15 years) monthly observations: the first 168 values (14 years) are used as @@ -40,7 +41,7 @@ #' mbase <- t(FoReco_data$base[, id]) #' # residuals #' id <- which(simplify2array(strsplit(colnames(FoReco_data$res), -#' split = "_"))[1, ] == "k1") +#' split = "_"))[1, ] == paste("k", K[i], sep="")) #' mres <- t(FoReco_data$res[, id]) #' hts_recf[[i]] <- htsrec(mbase, C = FoReco_data$C, comb = "shr", #' res = mres, keep = "recf") @@ -48,6 +49,7 @@ #' names(hts_recf) <- paste("k", K, sep="") #' #' # Forecast reconciliation through temporal hierarchies for all time series +#' # comb = "acov" #' n <- NROW(FoReco_data$base) #' thf_recf <- matrix(NA, n, NCOL(FoReco_data$base)) #' dimnames(thf_recf) <- dimnames(FoReco_data$base) @@ -61,22 +63,26 @@ #' } #' #' # Iterative cross-temporal reconciliation +#' # Each iteration: t-acov + cs-shr #' ite_recf <- iterec(FoReco_data$base, note=FALSE, #' m = 12, C = FoReco_data$C, #' thf_comb = "acov", hts_comb = "shr", #' res = FoReco_data$res, start_rec = "thf")$recf #' #' # Heuristic first-cross-sectional-then-temporal cross-temporal reconciliation +#' # cs-shr + t-acov #' cst_recf <- cstrec(FoReco_data$base, m = 12, C = FoReco_data$C, #' thf_comb = "acov", hts_comb = "shr", #' res = FoReco_data$res)$recf #' #' # Heuristic first-temporal-then-cross-sectional cross-temporal reconciliation +#' # t-acov + cs-shr #' tcs_recf <- tcsrec(FoReco_data$base, m = 12, C = FoReco_data$C, #' thf_comb = "acov", hts_comb = "shr", #' res = FoReco_data$res)$recf #' #' # Optimal cross-temporal reconciliation +#' # comb = "bdshr" #' oct_recf <- octrec(FoReco_data$base, m = 12, C = FoReco_data$C, #' comb = "bdshr", res = FoReco_data$res, keep = "recf") #' } diff --git a/R/commat.R b/R/commat.R index 3a48048..ba58aa2 100755 --- a/R/commat.R +++ b/R/commat.R @@ -1,18 +1,22 @@ #' @title Commutation matrix #' -#' @description This function returns the (\code{(r c) x (r c)}) commutation matrix \code{P} such that -#' \deqn{P vec(Y) = vec(Y'),} -#' where \code{Y} is a (\code{r x c}) matrix. +#' @description This function returns the (\mjseqn{(r c) \times (r c)}) +#' commutation matrix \mjseqn{\textbf{P}} such that +#' \mjsdeqn{\textbf{P} \mbox{vec}(\textbf{Y}) = \mbox{vec}(\textbf{Y}'),} +#' where \mjseqn{\textbf{Y}} is a (\mjseqn{r \times c}) matrix. #' -#' @param r Number of rows of \code{Y}. -#' @param c Number of columns of \code{Y}. +#' @param r Number of rows of \mjseqn{\textbf{Y}}. +#' @param c Number of columns of \mjseqn{\textbf{Y}}. #' -#' @return A sparse (\code{(r c) x (r c)}) matrix. +#' @return A sparse (\mjseqn{(r c) \times (r c)}) matrix, \mjseqn{\textbf{P}}. #' -#' @references Magnus, J.R., Neudecker, H. (2019), Matrix Differential Calculus with Applications -#' in Statistics and Econometrics, third edition, New York, Wiley, pp. 54-55. +#' @references Magnus, J.R., Neudecker, H. (2019), Matrix Differential Calculus +#' with Applications in Statistics and Econometrics, third edition, New York, +#' Wiley, pp. 54-55. #' #' @keywords utilities +#' @family utilities +#' #' @examples #' Y <- matrix(rnorm(30), 5, 6) #' P <- commat(5, 6) diff --git a/R/cstrec.R b/R/cstrec.R index 3a5c5b1..c41cc08 100755 --- a/R/cstrec.R +++ b/R/cstrec.R @@ -1,35 +1,51 @@ #' @title Heuristic first-cross-sectional-then-temporal cross-temporal forecast reconciliation #' #' @description -#' The order of application of the two reconciliation steps proposed by Kourentzes and Athanasopoulos (2019), -#' implemented in the function \code{\link[FoReco]{tcsrec}}, may be inverted. The function -#' \code{\link[FoReco]{cstrec}} performs cross-sectional reconciliation (\code{\link[FoReco]{htsrec}}) -#' first, then temporal reconciliation (\code{\link[FoReco]{thfrec}}), and finally applies the -#' average of the projection matrices obtained in the second step to the one dimensional reconciled -#' values obtained in the first step. +#' \loadmathjax +#' Cross-temporal forecast reconciliation according to the heuristic procedure by +#' Kourentzes and Athanasopoulos (2019), where the order of application of the +#' two reconciliation steps (temporal-first-then-cross-sectional, as in the function +#' \code{\link[FoReco]{tcsrec}()}), is inverted. The function +#' \code{\link[FoReco]{cstrec}()} performs cross-sectional reconciliation +#' (\code{\link[FoReco]{htsrec}()}) first, then temporal reconciliation +#' (\code{\link[FoReco]{thfrec}()}), and finally applies the average of the +#' projection matrices obtained in the second step to the one dimensional +#' reconciled values obtained in the first step. #' -#' @param basef (\code{n x h(k* + m)}) matrix of base forecasts to be reconciled; -#' \code{n} is the total number of variables, \code{m} is the highest frequency, -#' \code{k*} is the sum of (\code{p-1}) factors of \code{m}, excluding \code{m}, -#' and \code{h} is the forecast horizon. Each row identifies, a time series, and the forecasts -#' are ordered as [lowest_freq' ... highest_freq']'. -#' @param hts_comb,thf_comb Type of covariance matrix (respectively (\code{n x n}) and -#' (\code{(k* + m) x (k* + m)})) to be used in the cross-sectional and temporal reconciliation, -#' see more in \code{comb} param of \code{\link[FoReco]{htsrec}} and \code{\link[FoReco]{thfrec}}. -#' @param res (\code{n x N(k* + m)}) matrix containing the residuals at all the -#' temporal frequencies, ordered [lowest_freq' ... highest_freq']' (columns) for -#' each variable (row), needed to estimate the covariance matrix when \code{hts_comb =} -#' \code{\{"wls",} \code{"shr",} \code{"sam"\}} and/or \code{hts_comb =} \code{\{"wlsv",} -#' \code{"wlsh",} \code{"acov",} \code{"strar1",} \code{"sar1",} \code{"har1",} -#' \code{"shr",} \code{"sam"\}}. The rows must be in the same order as \code{basef}. -#' @param ... any other options useful for \code{\link[FoReco]{htsrec}} and -#' \code{\link[FoReco]{thfrec}}, e.g. \code{m}, \code{C} (or \code{Ut} and \code{nb}), -#' \code{nn} (for non negativity reconciliation only at first step), ... +#' @param basef (\mjseqn{n \times h(k^\ast+m)}) matrix of base forecasts to be +#' reconciled, \mjseqn{\widehat{\mathbf{Y}}}; \mjseqn{n} is the total number of variables, +#' \mjseqn{m} is the highest time frequency, \mjseqn{k^\ast} is the sum of (a +#' subset of) (\mjseqn{p-1}) factors of \mjseqn{m}, excluding \mjseqn{m}, and +#' \mjseqn{h} is the forecast horizon for the lowest frequency time series. +#' Each row identifies a time series, and the forecasts are ordered as +#' [lowest_freq' ... highest_freq']'. +#' @param hts_comb,thf_comb Type of covariance matrix (respectively +#' (\mjseqn{n \times n}) and (\mjseqn{(k^\ast + m) \times (k^\ast + m)})) to +#' be used in the cross-sectional and temporal reconciliation, see more in +#' \code{comb} param of \code{\link[FoReco]{htsrec}()} and +#' \code{\link[FoReco]{thfrec}()}. +#' @param res (\mjseqn{n \times N(k^\ast + m)}) matrix containing the residuals +#' at all the temporal frequencies ordered as [lowest_freq' ... highest_freq']' +#' (columns) for each variable (row), needed to estimate the covariance matrix +#' when \code{hts_comb =} \code{\{"wls",} \code{"shr",} \code{"sam"\}} and/or +#' \code{hts_comb =} \code{\{"wlsv",} \code{"wlsh",} \code{"acov",} +#' \code{"strar1",} \code{"sar1",} \code{"har1",} \code{"shr",} \code{"sam"\}}. +#' The rows must be in the same order as \code{basef}. +#' @param ... any other options useful for \code{\link[FoReco]{htsrec}()} and +#' \code{\link[FoReco]{thfrec}()}, e.g. \code{m}, \code{C} (or \code{Ut} and +#' \code{nb}), \code{nn} (for non-negative reconciliation only at the first step), +#' \code{mse}, \code{corpcor}, \code{type}, \code{sol}, \code{settings}, +#' \code{W}, \code{Omega},... +#' +#' @details +#' \strong{Warning}, +#' the two-step heuristic reconciliation allows non negativity constraints only in the first step. +#' This means that it is not guaranteed the non-negativity of the final reconciled values. #' #' @return #' The function returns a list with two elements: -#' \item{\code{recf}}{(\code{n x h(k* + m)}) reconciled forecasts matrix.} -#' \item{\code{M}}{Projection matrix (projection approach).} +#' \item{\code{recf}}{(\mjseqn{n \times h(k^\ast + m)}) reconciled forecasts matrix, \mjseqn{\widetilde{\textbf{Y}}}.} +#' \item{\code{M}}{Matrix which transforms the uni-dimensional reconciled forecasts of step 1 (projection approach) .} #' #' @references #' Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: @@ -47,23 +63,26 @@ #' Matrix Estimation and Implications for Functional Genomics, \emph{Statistical #' Applications in Genetics and Molecular Biology}, 4, 1. #' -#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2018). OSQP: -#' An Operator Splitting Solver for Quadratic Programs, \href{https://arxiv.org/abs/1711.08013}{arXiv:1711.08013}. +#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2020). OSQP: +#' An Operator Splitting Solver for Quadratic Programs, \emph{Mathematical Programming Computation}, +#' 12, 4, 637-672. #' #' Stellato, B., Banjac, G., Goulart, P., Boyd, S., Anderson, E. (2019), OSQP: #' Quadratic Programming Solver using the 'OSQP' Library, R package version 0.6.0.3 #' (October 10, 2019), \href{https://CRAN.R-project.org/package=osqp}{https://CRAN.R-project.org/package=osqp}. #' -#' @keywords reconciliation heuristic +#' @keywords heuristic +#' @family reconciliation procedures +#' #' @examples #' data(FoReco_data) -#' obj <- cstrec(FoReco_data$base, m = 12, C = FoReco_data$C, thf_comb = "acov", -#' hts_comb = "shr", res = FoReco_data$res) +#' obj <- cstrec(FoReco_data$base, m = 12, C = FoReco_data$C, +#' hts_comb = "shr", thf_comb = "acov", res = FoReco_data$res) #' -#' @usage cstrec(basef, thf_comb, hts_comb, res, ...) +#' @usage cstrec(basef, hts_comb, thf_comb, res, ...) #' #' @export -cstrec <- function(basef, thf_comb, hts_comb, res, ...) { +cstrec <- function(basef, hts_comb, thf_comb, res, ...) { arg_input <- list(...) @@ -84,7 +103,7 @@ cstrec <- function(basef, thf_comb, hts_comb, res, ...) { tools <- thf_tools(m) kset <- tools$kset - m <- max(kset) + m <-tools$m kt <- tools$kt # Step 1 @@ -118,41 +137,32 @@ cstrec <- function(basef, thf_comb, hts_comb, res, ...) { # Step 2 arg_thf <- names(as.list(args(thfrec))) - arg_thf <- arg_thf[!(arg_thf %in% c("basef", "keep", "res", "", "comb", "m", "nn", "bounds"))] + arg_thf <- arg_thf[!(arg_thf %in% c("basef", "keep", "res", "", "comb", "m", "nn", "bounds", "nn_type"))] if (missing(res)) { M <- lapply(split(Y1, row(Y1)), function(x) { obj <- do.call("thfrec", c(list(basef = x, m = kset, comb = thf_comb), arg_input[which(names(arg_input) %in% arg_thf)])) - out <- obj[["M"]] - if(is.null(out)){ - return(obj[["G"]]) - }else{ - return(out) - } + return(obj[["M"]]) }) } else { M <- mapply(function(Y, X) { obj <- do.call("thfrec", c(list(basef = Y, m = kset, comb = thf_comb, res = X), arg_input[which(names(arg_input) %in% arg_thf)])) - out <- obj[["M"]] - if(is.null(out)){ - return(obj[["G"]]) - }else{ - return(out) - } + return(obj[["M"]]) }, Y = split(Y1, row(Y1)), X = split(res, row(res)), SIMPLIFY = FALSE) } ## Step 3: Cross-Temporal reconciled forecasts with heuristic - meanM <- apply(simplify2array(lapply(M, as.matrix)), c(1, 2), sum) / length(M) + #meanM <- apply(simplify2array(lapply(M, as.matrix)), c(1, 2), sum) / length(M) + meanM <- Reduce("+", M)/length(M) if (h == 1) { Y3 <- Y1 %*% t(meanM) } else { h_pos <- rep(rep(1:h, length(kset)), rep((m/kset), each = h)) - D <- Dmat(h = h, kset = kset, n = 1) + D <- Dmat(h = h, m = kset, n = 1) Y3 <- list() for (i in 1:h) { diff --git a/R/ctbu.R b/R/ctbu.R index 72a6834..5591aa9 100755 --- a/R/ctbu.R +++ b/R/ctbu.R @@ -1,20 +1,44 @@ -#' @title Bottom-up Cross-temporal forecast reconciliation +#' @title Bottom-up cross-temporal forecast reconciliation #' -#' @description Cross temporal reconciled forecasts for all series at any temporal -#' aggregation level can be easily computed by appropriate summation of the high-frequency -#' bottom base forecasts \eqn{\hat{\textbf{b}}_i, i = 1,...,n_b, \;}{} according to a +#' @description +#' \loadmathjax +#' Cross temporal reconciled forecasts for all series at any temporal +#' aggregation level are computed by appropriate summation of the high-frequency +#' bottom base forecasts \mjseqn{\widehat{\mathbf{b}}_i, i = 1,...,n_b}, according to a #' bottom-up procedure like what is currently done in both the cross-sectional #' and temporal frameworks. #' -#' @param Bmat (\code{nb x (h m)}) matrix of high-frequency bottom time series base forecasts. +#' @param Bmat (\mjseqn{n_b \times h m}) matrix of high-frequency bottom time +#' series base forecasts (\mjseqn{\widehat{\mathbf{B}}^{[1]}}). #' \code{h} is the forecast horizon for the lowest frequency (most temporally aggregated) #' time series. -#' @param C (\code{na x nb}) cross-sectional (contemporaneous) matrix mapping the bottom -#' level series into the higher level ones. -#' @param m Highest available sampling frequency per seasonal cycle (max. order of temporal aggregation). +#' @param C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones. +#' @param m Highest available sampling frequency per seasonal cycle (max. order +#' of temporal aggregation, \mjseqn{m}), or a subset of the \mjseqn{p} factors +#' of \mjseqn{m}. #' -#' @return The function returns a (\code{n x (h (k* + m))}) matrix of -#' bottom-up cross-temporally reconciled forecasts. +#' @details +#' Denoting by \mjseqn{\ddot{\mathbf{Y}}} the (\mjseqn{n \times h (k^\ast + m)}) matrix containing +#' the bottom-up cross temporal reconciled forecasts, it is: +#' \mjsdeqn{\ddot{\mathbf{Y}} = \left[\begin{array}{cc} +#' \mathbf{C}\widehat{\mathbf{B}}^{[1]}\mathbf{K}_1' & \mathbf{C}\widehat{\mathbf{B}}^{[1]} \cr +#' \widehat{\mathbf{B}}^{[1]} \mathbf{K}_1' & \widehat{\mathbf{B}}^{[1]} +#' \end{array}\right],} +#' where \mjseqn{\mathbf{C}} is the cross-sectional (contemporaneous) aggregation matrix, +#' \mjseqn{\mathbf{K}_1} is the temporal aggregation matrix with \mjseqn{h=1}, and +#' \mjseqn{\widehat{\mathbf{B}}^{[1]}} is the matrix containing the high-frequency bottom +#' time series base forecasts. This expression is equivalent to +#' \mjseqn{\mbox{vec}(\ddot{\mathbf{Y}}') = \widetilde{\mathbf{S}} +#' \mbox{vec}(\widehat{\mathbf{Y}}')} for \mjseqn{h = 1}, where +#' \mjseqn{\widetilde{\mathbf{S}}} is the cross-temporal summing matrix for +#' \mjseqn{\mbox{vec}(\widehat{\mathbf{Y}}')}, and \mjseqn{\widehat{\mathbf{Y}}} +#' is the (\mjseqn{n \times h (k^\ast + m)}) matrix containing all the base forecasts +#' at any temporal aggregation order. +#' +#' +#' @return The function returns a (\mjseqn{n \times h (k^\ast + m)}) matrix of +#' bottom-up cross-temporally reconciled forecasts, \mjseqn{\ddot{\mathbf{Y}}}. #' #' @references #' Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: @@ -23,13 +47,15 @@ #' #' @examples #' data(FoReco_data) +#' # monthly base forecasts #' id <- which(simplify2array(strsplit(colnames(FoReco_data$base), #' split = "_"))[1, ] == "k1") #' hfbts <- FoReco_data$base[-c(1:3), id] #' obj <- ctbu(Bmat = hfbts, m = 12, C = FoReco_data$C) #' rownames(obj) <- rownames(FoReco_data$base) #' -#' @keywords reconciliation bottom-up +#' @keywords bottom-up +#' @family reconciliation procedures #' @import Matrix #' @export ctbu <- function(Bmat, m, C) { @@ -38,7 +64,7 @@ ctbu <- function(Bmat, m, C) { tools <- thf_tools(m) kset <- tools$kset - m <- max(kset) + m <- tools$m if (NCOL(Bmat) %% m != 0) { stop("The number of Bmat's columns doesn't match with m*h", call. = FALSE) diff --git a/R/ctf_tools.R b/R/ctf_tools.R index 1a91d8f..7b216bb 100755 --- a/R/ctf_tools.R +++ b/R/ctf_tools.R @@ -1,110 +1,120 @@ #' Cross-temporal reconciliation tools #' -#' Some useful tools for the cross-temporal forecast reconciliation of cross-sectionally and -#' temporally constrained time series +#' @description +#' \loadmathjax +#' Some useful tools for the cross-temporal forecast reconciliation of a linearly constrained +#' (hierarchical/grouped) multiple time series. #' -#' @param C (\code{na x nb}) cross-sectional (contemporaneous) matrix mapping the bottom -#' level series into the higher level ones. -#' @param m Highest available sampling frequency per seasonal cycle (max. order of -#' temporal aggregation). +#' @param m Highest available sampling frequency per seasonal cycle (max. order +#' of temporal aggregation, \mjseqn{m}), or a subset of the \mjseqn{p} factors +#' of \mjseqn{m}. #' @param h Forecast horizon for the lowest frequency (most temporally aggregated) time #' series (\emph{default} is \code{1}). +#' @param C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones. #' @param Ut Zero constraints cross-sectional (contemporaneous) kernel matrix -#' \eqn{(\textbf{U}'\textbf{Y} = \mathbf{0})}{} spanning the null space valid for the reconciled -#' forecasts. It can be used instead of parameter \code{C}, but in this case \code{nb} (n = na + nb) is needed. If -#' the hierarchy admits a structural representation, \code{Ut} has dimension (\code{na x n}). -#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} is not used. -#' @param Sstruc If \code{Sstruc = TRUE} the function returns also the structural representation matrix of -#' a cross-temporal system (\emph{default} is \code{FALSE}). -#' +#' \mjseqn{(\textbf{U}'\textbf{y} = \mathbf{0})} spanning the null space valid +#' for the reconciled forecasts. It can be used instead of parameter +#' \code{C}, but \code{nb} (\mjseqn{n = n_a + n_b}) is needed if +#' \mjseqn{\textbf{U}' \neq [\textbf{I} \ -\textbf{C}]}{}. If the hierarchy +#' admits a structural representation, \mjseqn{\textbf{U}'} has dimension +#' (\mjseqn{n_a \times n}). +#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} +#' and \code{Ut} are not used. +#' @param sparse Option to return sparse object (\emph{default} is \code{TRUE}). #' #' @return #' \strong{ctf} list with: -#' \item{\code{Ht}}{Full row-rank cross-temporal zero-valued constraints (kernel) -#' matrix\eqn{,\; \textbf{H}'\textbf{y} = \mathbf{0}}{}.} -#' \item{\code{Htbreve}}{Complete, not full row-rank cross-temporal zero-valued -#' constraints (kernel) matrix.} -#' \item{\code{Htstruc}}{Zero constraints full row-rank cross-temporal kernel matrix -#' (structural representation) \eqn{,\; \check{\textbf{H}}'}{}.} -#' \item{\code{Sstruc}}{Cross-temporal summing matrix (structural -#' representation)\eqn{,\; \check{\textbf{S}}}{}.} +#' \item{\code{Ht}}{Full row-rank cross-temporal zero constraints (kernel) +#' matrix coherent with \mjseqn{\mathbf{y} = \mbox{vec}(\mathbf{Y}')}: \mjseqn{\mathbf{H}'\mathbf{y} = \mathbf{0}}.} +#' \item{\code{Hbrevet}}{Complete, not full row-rank cross-temporal zero +#' constraints (kernel) matrix coherent with \mjseqn{\mathbf{y} = \mbox{vec}(\mathbf{Y}')}: +#' \mjseqn{\breve{\mathbf{H}}'\mathbf{y} = \mathbf{0}}.} +#' \item{\code{Hcheckt}}{Full row-rank cross-temporal zero constraints (kernel) matrix coherent with +#' \mjseqn{\check{\mathbf{y}}} (structural representation): +#' \mjseqn{\check{\mathbf{H}}' \check{\mathbf{y}} = \mathbf{0}}.} +#' \item{\code{Ccheck}}{Cross-temporal aggregation matrix \mjseqn{\check{\mathbf{C}}} +#' coherent with \mjseqn{\check{\mathbf{y}}} (structural representation).} +#' \item{\code{Scheck}}{Cross-temporal summing matrix \mjseqn{\check{\mathbf{S}}} +#' coherent with \mjseqn{\check{\mathbf{y}}} (structural representation).} +#' \item{\code{Stilde}}{Cross-temporal summing matrix \mjseqn{\widetilde{\mathbf{S}}} +#' coherent with \mjseqn{\mathbf{y} = \mbox{vec}(\mathbf{Y}')}.} #' -#' \strong{hts} list from \code{\link{hts_tools}} +#' \strong{hts} list from \code{\link{hts_tools}} . #' -#' \strong{thf} list from \code{\link{thf_tools}} +#' \strong{thf} list from \code{\link{thf_tools}} . #' #' @examples #' # One level hierarchy (na = 1, nb = 2) with quarterly data -#' obj <- ctf_tools(C = matrix(c(1, 1), 1), m = 4, Sstruc = TRUE) +#' obj <- ctf_tools(C = matrix(c(1, 1), 1), m = 4) +#' #' @keywords utilities +#' @family utilities +#' #' @import Matrix #' #' @export -#' -ctf_tools <- function(C, m, h = 1, Ut, nb, Sstruc = FALSE) { +ctf_tools <- function(C, m, h = 1, Ut, nb, sparse = TRUE) { + # Using Ut or C if (missing(C)) { - n <- ncol(Ut) - na <- n - nb + if (missing(Ut)) { + stop("Please, give C or Ut.", call. = FALSE) + } else if(missing(nb)){ + hts <- hts_tools(Ut = Ut, h = 1, sparse = TRUE) + } else { + hts <- hts_tools(Ut = Ut, nb = nb, h = 1, sparse = TRUE) + } } else { - nb <- ncol(C) - na <- nrow(C) - n <- nb + na - S <- rbind(C, Diagonal(nb)) - Ut <- cbind(Diagonal(na), -C) + hts <- hts_tools(C = C, h = 1, sparse = TRUE) } + n <- hts$n + na <- hts$na + nb <- hts$nb + C <- hts$C + S <- hts$S + Ut <- hts$Ut + tmp <- thf_tools(m, h = h, sparse = TRUE) - m <- max(tmp$kset) + m <- tmp$m Zt <- tmp$Zt + K <- tmp$K + R <- tmp$R ks <- tmp$ks kt <- tmp$kt - Htbreve <- rbind(kronecker(Ut, Diagonal(h * kt)), kronecker(Diagonal(n), Zt)) - D <- commat(h * kt, n) - Us <- cbind(Matrix(0, h * NROW(Ut) * m, n * h * ks), kronecker(Diagonal(h * m), Ut)) %*% D - Ht <- rbind(Us, kronecker(Diagonal(n), Zt)) + Hbrevet <- rbind(kronecker(Ut, .sparseDiagonal(h * kt)), kronecker(.sparseDiagonal(n), Zt)) + P <- commat(h * kt, n) + Us <- cbind(Matrix(0, h * NROW(Ut) * m, n * h * ks), kronecker(.sparseDiagonal(h * m), Ut)) %*% P + Ht <- rbind(Us, kronecker(.sparseDiagonal(n), Zt)) out1 <- list() - out1$Ht <- Ht - out1$Htbreve <- Htbreve - - if (Sstruc) { - if (missing(C)) { - stop("The argument C is not specified, but it's necessary for the structural representation of a cross-temporal system") - } - Qtilde <- commat(nb, kt) %*% bdiag(commat(ks, nb), commat(m, nb)) - Q <- bdiag(Diagonal(na * kt), Qtilde) - - Htstruc <- Matrix(pracma::rref(as.matrix(Ht %*% Q)), sparse = TRUE) - r <- nrow(Htstruc) - c <- ncol(Htstruc) - Cstruc <- -Htstruc[, (r + 1):c] - Sstruc <- rbind(Cstruc, Diagonal(m * nb)) - - # Trivial checks on matrix Hstruc - # 1. The left square block must be an identity matrix - # 2. The last m*nb elements of the first row mast be all equal to -1 - flag <- 0 - Htcheck1 <- all(abs(Htstruc[, 1:r] - Diagonal(r)) < 1e-6) - Htcheck2 <- all(abs(Htstruc[1, (r + 1):c] + rep(1, m * nb)) < 1e-6) - if (!Htcheck1 | !Htcheck2) { - flag <- 1 - warning("Failed checks on matrix Htstruc. Check the results!", call. = FALSE) - } - - out1$Htstruc <- Htstruc - out1$Sstruc <- Sstruc + if(!sparse){ + out1$Ht <- as.matrix(Ht) + out1$Hbrevet <- as.matrix(Hbrevet) + }else{ + out1$Ht <- Ht + out1$Hbrevet <- Hbrevet } + if(!is.null(C)){ + Ccheck <- rbind(kronecker(C, R),kronecker(.sparseDiagonal(nb), K)) + Hcheckt <- cbind(.sparseDiagonal(h*(na*m + n*ks)), -Ccheck) + Scheck <- rbind(Ccheck, .sparseDiagonal(nb*m*h)) + Stilde <- kronecker(S, R) - out2 <- list() - if (!missing(C)) { - out2$S <- kronecker(S, Diagonal(h)) + if(!sparse){ + out1$Hcheckt <- as.matrix(Hcheckt) + out1$Ccheck <- as.matrix(Ccheck) + out1$Scheck <- as.matrix(Scheck) + out1$Stilde <- as.matrix(Stilde) + }else{ + out1$Hcheckt <- Hcheckt + out1$Ccheck <- Ccheck + out1$Scheck <- Scheck + out1$Stilde <- Stilde + } } - out2$Ut <- kronecker(Ut, Diagonal(h)) - out2$n <- n - out2$na <- na - out2$nb <- nb - return(list(ctf = out1, hts = out2, thf = tmp)) + return(list(ctf = out1, hts = hts, thf = tmp)) } diff --git a/R/hts_tools.R b/R/hts_tools.R index 447a9d7..922ca49 100755 --- a/R/hts_tools.R +++ b/R/hts_tools.R @@ -1,60 +1,124 @@ #' Cross-sectional reconciliation tools #' -#' Some useful tools for the cross-sectional reconciliation of linearly and hierarchically constrained time series +#' @description +#' \loadmathjax +#' Some useful tools for the cross-sectional forecast reconciliation of a +#' linearly constrained (e.g., hierarchical/grouped) multiple time series. #' -#' @param C (\code{na x nb}) cross-sectional (contemporaneous) matrix mapping the bottom -#' level series into the higher level ones. +#' @param C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones. +#' @param Ut Zero constraints cross-sectional (contemporaneous) kernel matrix +#' \mjseqn{(\mathbf{U}'\mathbf{y} = \mathbf{0})} spanning the null space valid +#' for the reconciled forecasts. It can be used instead of parameter +#' \code{C}, but \code{nb} is needed if +#' \mjseqn{\mathbf{U}' \neq [\mathbf{I} \ -\mathbf{C}]}. If the hierarchy +#' admits a structural representation, \mjseqn{\mathbf{U}'} has dimension +#' (\mjseqn{n_a \times n}). +#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} +#' and \code{Ut} are not used. #' @param h Forecast horizon (\emph{default} is \code{1}). -#' @param Ut (\code{na x n}) zero constraints cross-sectional (contemporaneous) kernel -#' matrix \eqn{\textbf{U}'\textbf{Y} = \mathbf{0}_{\left[n_a \times (k^*+m)\right]}}{} -#' spanning the null space valid for the reconciled forecasts. It can be used instead of -#' parameter \code{C}, but needs \code{nb} (n = na + nb). -#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} is not used. -#' @param sparse Option to return sparse object (\emph{default} is \code{TRUE}). +#' @param sparse Option to return sparse matrices (\emph{default} is \code{TRUE}). #' #' @return A list of five elements: -#' \item{\code{S}}{(\code{n x nb}) cross-sectional (contemporaneous) summing matrix.} -#' \item{\code{Ut}}{(\code{na x n}) zero constraints cross-sectional (contemporaneous) -#' kernel matrix.} -#' \item{\code{n}}{Number of variables \code{na + nb}.} +#' \item{\code{C}}{(\mjseqn{n \times n_b}) cross-sectional (contemporaneous) aggregation matrix.} +#' \item{\code{S}}{(\mjseqn{n \times n_b}) cross-sectional (contemporaneous) summing matrix, +#' \mjseqn{\mathbf{S} = \left[\begin{array}{c} \mathbf{C} \cr \mathbf{I}_{n_b}\end{array}\right].}} +#' \item{\code{Ut}}{(\mjseqn{n_a \times n}) zero constraints cross-sectional (contemporaneous) +#' kernel matrix. If the hierarchy admits a structural representation \mjseqn{\mathbf{U}' = [\mathbf{I} \ -\mathbf{C}]}} +#' \item{\code{n}}{Number of variables \mjseqn{n_a + n_b}.} #' \item{\code{na}}{Number of upper level variables.} #' \item{\code{nb}}{Number of bottom level variables.} #' #' @examples #' # One level hierarchy (na = 1, nb = 2) #' obj <- hts_tools(C = matrix(c(1, 1), 1), sparse = FALSE) +#' #' @keywords utilities +#' @family utilities +#' #' @import Matrix #' #' @export #' hts_tools <- function(C, h = 1, Ut, nb, sparse = TRUE) { if (missing(C)) { - if(missing(Ut) | missing(nb)){ - stop("the argument C (or Ut AND nb) is not specified", call. = FALSE) + if (missing(Ut)) { + stop("Please, give C or Ut.", call. = FALSE) + } else { + if (!(is.matrix(Ut) | is(Ut, "Matrix"))){ + stop("Ut must be a matrix.", call. = FALSE) + } + + if(!is(Ut, "Matrix")){ + Ut <- Matrix(Ut, sparse = TRUE) + } + + if(all(.sparseDiagonal(NROW(Ut))==Ut[,1:NROW(Ut)])){ + C <- -Ut[, -c(1:NROW(Ut))] + if(any(rowSums(abs(C)) == 0)){ + message("Removed a zeros row in C matrix") + C <- C[rowSums(abs(C)) != 0,] + Ut <- cbind(.sparseDiagonal(NROW(C)), -C) + } + n <- NCOL(Ut) + na <- NROW(Ut) + nb <- n - na + S <- rbind(C, .sparseDiagonal(nb)) + } else if(missing(nb)){ + stop("Ut is not in form [I -C], give also nb", call. = FALSE) + }else{ + n <- NCOL(Ut) + na <- n-nb + C <- NULL + S <- NULL + } } - n <- ncol(Ut) - na <- n - nb } else { if (!(is.matrix(C) | is(C, "Matrix"))){ stop("C must be a matrix.", call. = FALSE) } + if(!is(C, "Matrix")){ + C <- Matrix(C, sparse = TRUE) + } + + if(any(rowSums(abs(C)) == 0)){ + message("Removed a zeros row in C matrix") + C <- C[rowSums(abs(C)) != 0,] + } nb <- ncol(C) na <- nrow(C) n <- nb + na - S <- rbind(C, Diagonal(nb)) - Ut <- cbind(Diagonal(na), -C) + S <- rbind(C, .sparseDiagonal(nb)) + Ut <- cbind(.sparseDiagonal(na), -C) + } + + if (n <= nb) { + stop("n <= nb, total number of TS is less (or equal) than the number of bottom TS.", call. = FALSE) } out2 <- list() - if (!missing(C)) { - out2$S <- kronecker(S, Diagonal(h)) + if (!is.null(C)) { + if(any(rowSums(abs(C)) == 1)){ + message(paste0("There is only one non-zero value in ", sum(rowSums(abs(C)) == 1), + " row(s) of C.", + "\nRemember that Foreco can also work with unbalanced hierarchies (recommended).")) + } + + if(h>1){ + out2$C <- kronecker(C, .sparseDiagonal(h)) + out2$S <- kronecker(S, .sparseDiagonal(h)) + }else{ + out2$S <- S + out2$C <- C + } + if (!sparse) { + out2$C <- as.matrix(out2$C) out2$S <- as.matrix(out2$S) } } - out2$Ut <- kronecker(Ut, Diagonal(h)) + out2$Ut <- kronecker(Ut, .sparseDiagonal(h)) if (!sparse) { out2$Ut <- as.matrix(out2$Ut) @@ -63,6 +127,5 @@ hts_tools <- function(C, h = 1, Ut, nb, sparse = TRUE) { out2$n <- n out2$na <- na out2$nb <- nb - return(out2) } diff --git a/R/htsrec.R b/R/htsrec.R index 1ed6618..69255c6 100755 --- a/R/htsrec.R +++ b/R/htsrec.R @@ -1,20 +1,21 @@ #' @title Cross-sectional (contemporaneous) forecast reconciliation #' #' @description -#' Cross-sectional (contemporaneous) forecast reconciliation of hierarchical and -#' grouped time series. The reconciled forecasts are calculated either through a +#' Cross-sectional (contemporaneous) forecast reconciliation of a linearly constrained +#' (e.g., hierarchical/grouped) multiple time series. +#' The reconciled forecasts are calculated either through a #' projection approach (Byron, 1978, see also van Erven and Cugliari, 2015, and #' Wickramasuriya et al., 2019), or the equivalent structural approach by Hyndman #' et al. (2011). Moreover, the classic bottom-up approach is available. #' #' @usage htsrec(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FALSE, -#' type = "M", sol = "direct", nn = FALSE, keep = "list", -#' settings = osqpSettings(), bounds = NULL, W = NULL) +#' type = "M", sol = "direct", keep = "list", nn = FALSE, +#' nn_type = "osqp", settings = list(), bounds = NULL, W = NULL) #' -#' @param basef (\code{h x n}) matrix of base forecasts to be reconciled; -#' \code{h} is the forecast horizon and \code{n} is the total number of time series. +#' @param basef (\mjseqn{h \times n}) matrix of base forecasts to be reconciled; +#' \mjseqn{h} is the forecast horizon and \mjseqn{n} is the total number of time series. #' @param comb Type of the reconciliation. Except for Bottom-up, each -#' option corresponds to a specified (\code{n x n}) covariance matrix: +#' option corresponds to a specific (\mjseqn{n \times n}) covariance matrix: #' \itemize{ #' \item \bold{bu} (Bottom-up); #' \item \bold{ols} (Identity); @@ -24,70 +25,186 @@ #' \item \bold{sam} (Sample covariance matrix - MinT-sam) - uses res; #' \item \bold{w} use your personal matrix W in param \code{W}. #' } -#' @param C (\code{na x nb}) cross-sectional (contemporaneous) matrix mapping the bottom -#' level series into the higher level ones. -#' @param res (\code{N x n}) in-sample residuals matrix needed when \code{comb =} +#' @param C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones. +#' @param Ut Zero constraints cross-sectional (contemporaneous) kernel matrix +#' \mjseqn{(\mathbf{U}'\mathbf{y} = \mathbf{0})} spanning the null space valid +#' for the reconciled forecasts. It can be used instead of parameter +#' \code{C}, but \code{nb} (\mjseqn{n = n_a + n_b}) is needed if +#' \mjseqn{\mathbf{U}' \neq [\mathbf{I} \ -\mathbf{C}]}{}. If the hierarchy +#' admits a structural representation, \mjseqn{\mathbf{U}'} has dimension +#' (\mjseqn{n_a \times n}). +#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} +#' and \code{Ut} are not used. +#' @param res (\mjseqn{N \times n}) in-sample residuals matrix needed when \code{comb =} #' \code{\{"wls",} \code{"shr",} \code{"sam"\}}. The columns must be in #' the same order as \code{basef}. -#' @param Ut Zero constraints cross-sectional (contemporaneous) kernel matrix -#' \eqn{(\textbf{U}'\textbf{Y} = \mathbf{0})}{} spanning the null space valid for the reconciled -#' forecasts. It can be used instead of parameter \code{C}, but in this case \code{nb} (n = na + nb) is needed. If -#' the hierarchy admits a structural representation, \code{Ut} has dimension (\code{na x n}). -#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} is not used. #' @param W This option permits to directly enter the covariance matrix: #' \enumerate{ -#' \item \code{W} must be a p.d. (\code{n x n}) matrix or a list of \code{h} matrix (one for each forecast horizon); +#' \item \code{W} must be a p.d. (\mjseqn{n \times n}) matrix or a list of \mjseqn{h} matrix (one for each forecast horizon); #' \item if \code{comb} is different from "\code{w}", \code{W} is not used. #' } -#' If you want to reconcile h step #' @param mse Logical value: \code{TRUE} (\emph{default}) calculates the -#' covariance matrix of the in-sample residuals (when necessary) according to the original -#' \pkg{hts} and \pkg{thief} formulation: no mean correction, T as denominator. +#' covariance matrix of the in-sample residuals (when necessary) according to +#' the original \pkg{hts} and \pkg{thief} formulation: no mean correction, +#' T as denominator. #' @param corpcor Logical value: \code{TRUE} if \pkg{corpcor} (\enc{Schäfer}{Schafer} et #' al., 2017) must be used to shrink the sample covariance matrix according to -#' \enc{Schäfer}{Schafer} and Strimmer (2005), otherwise the function uses the same -#' implementation as package \pkg{hts}. +#' \enc{Schäfer}{Schafer} and Strimmer (2005), otherwise the function uses the +#' same implementation as package \pkg{hts}. #' @param type Approach used to compute the reconciled forecasts: \code{"M"} for #' the projection approach with matrix M (\emph{default}), or \code{"S"} for the #' structural approach with summing matrix S. -#' @param keep Return a list object of the reconciled forecasts at all levels. -#' @param sol Solution technique for the reconciliation problem: either \code{"direct"} (\emph{default}) for the direct -#' solution or \code{"osqp"} for the numerical solution (solving a linearly constrained quadratic -#' program using \code{\link[osqp]{solve_osqp}}). -#' @param nn Logical value: \code{TRUE} if non-negative reconciled forecasts are wished. -#' @param settings Settings for \pkg{osqp} (object \code{\link[osqp]{osqpSettings}}). The default options -#' are: \code{verbose = FALSE}, \code{eps_abs = 1e-5}, \code{eps_rel = 1e-5}, -#' \code{polish_refine_iter = 100} and \code{polish = TRUE}. For details, see the -#' \href{https://osqp.org/}{\pkg{osqp} documentation} (Stellato et al., 2019). -#' @param bounds (\code{n x 2}) matrix with bounds on the variables: the first column is the lower bound, -#' and the second column is the upper bound +#' @param keep Return a list object of the reconciled forecasts at all levels +#' (if keep = "list") or only the reconciled forecasts matrix (if keep = "recf"). +#' @param sol Solution technique for the reconciliation problem: either +#' \code{"direct"} (\emph{default}) for the closed-form matrix solution, or +#' \code{"osqp"} for the numerical solution (solving a linearly constrained +#' quadratic program using \code{\link[osqp]{solve_osqp}}). +#' @param nn Logical value: \code{TRUE} if non-negative reconciled forecasts +#' are wished. +#' @param nn_type "osqp" (default), "KAnn" (only \code{type == "M"}) or "sntz". +#' @param settings Settings for \pkg{osqp} (object \code{\link[osqp]{osqpSettings}}). +#' The default options are: \code{verbose = FALSE}, \code{eps_abs = 1e-5}, +#' \code{eps_rel = 1e-5}, \code{polish_refine_iter = 100} and \code{polish = TRUE}. +#' For details, see the \href{https://osqp.org/}{\pkg{osqp} documentation} +#' (Stellato et al., 2019). +#' @param bounds (\mjseqn{n \times 2}) matrix containing the cross-sectional bounds: +#' the first column is the lower bound, and the second column is the upper bound. #' #' @details -#' In case of non-negativity constraints, there are two ways: -#' \enumerate{ -#' \item \code{sol = "direct"} and \code{nn = TRUE}: the base forecasts -#' will be reconciled at first without non-negativity constraints, then, if negative reconciled -#' values are present, the \code{"osqp"} solver is used. -#' \item \code{sol = "osqp"} and \code{nn = TRUE}: the base forecasts will be -#' reconciled through the \code{"osqp"} solver. +#' \loadmathjax +#' Let \mjseqn{\mathbf{y}} be a (\mjseqn{n \times 1}) vector of target point +#' forecasts which are wished to satisfy the system of linearly independent +#' constraints \mjsdeqn{\mathbf{U}'\mathbf{y} = \mathbf{0}_{(r \times 1)},} +#' where \mjseqn{\mathbf{U}'} is a (\mjseqn{r \times n}) matrix, with +#' rank\mjseqn{(\mathbf{U}') = r \leq n}, and \mjseqn{\mathbf{0}_{(r \times 1)}} +#' is a (\mjseqn{r \times 1}) null vector. Let \mjseqn{\widehat{\mathbf{y}}} +#' be a (\mjseqn{n \times 1}) vector of unbiased point forecasts, not +#' fulfilling the linear constraints (i.e., \mjseqn{\mathbf{U}'\widehat{\mathbf{y}} +#' \ne \mathbf{0}}). +#' +#' We consider a regression-based reconciliation method assuming +#' that \mjseqn{\widehat{\mathbf{y}}} is related to \mjseqn{\mathbf{y}} by +#' \mjsdeqn{\widehat{\mathbf{y}} = \mathbf{y} + \mathbf{\varepsilon},} +#' where \mjseqn{\mathbf{\varepsilon}} is a (\mjseqn{n \times 1}) vector +#' of zero mean disturbances, with known p.d. covariance matrix +#' \mjseqn{\mathbf{W}}. The reconciled forecasts \mjseqn{\widetilde{\mathbf{y}}} +#' are found by minimizing the generalized least squares (GLS) objective +#' function \mjseqn{\left(\widehat{\mathbf{y}} - \mathbf{y}\right)'\mathbf{W}^{-1} +#' \left(\widehat{\mathbf{y}} - \mathbf{y}\right)} constrained by +#' \mjseqn{\mathbf{U}'\mathbf{y} = \mathbf{0}_{(r \times 1)}}: +#' +#' \mjsdeqn{\widetilde{\mathbf{y}} = \mbox{argmin}_\mathbf{y} \left(\mathbf{y} - +#' \widehat{\mathbf{y}} \right)' \mathbf{W}^{-1} \left(\mathbf{y} - +#' \widehat{\mathbf{y}} \right), \quad \mbox{s.t. } \mathbf{U}'\mathbf{y} = +#' \mathbf{0}.} +#' +#' The solution is given by +#' \mjsdeqn{\widetilde{\mathbf{y}}= \widehat{\mathbf{y}} - \mathbf{W}\mathbf{U} +#' \left(\mathbf{U}’\mathbf{WU}\right)^{-1}\mathbf{U}'\widehat{\mathbf{y}}= +#' \mathbf{M}\widehat{\mathbf{y}},} +#' where \mjseqn{\mathbf{M} = \mathbf{I}_n - \mathbf{W}\mathbf{U}\left( +#' \mathbf{U}’\mathbf{WU}\right)^{-1}\mathbf{U}’} +#' is a (\mjseqn{n \times n}) projection matrix. This solution is used by +#' \code{\link[FoReco]{htsrec}} when \code{type = "M"}. +#' +#' Denoting with \mjseqn{\mathbf{d}_{\widehat{\mathbf{y}}} = \mathbf{0} - +#' \mathbf{U}'\widehat{\mathbf{y}}} the (\mjseqn{r \times 1}) vector containing +#' the \emph{coherency} errors of the base forecasts, we can re-state the solution as +#' \mjsdeqn{\widetilde{\mathbf{y}}= \widehat{\mathbf{y}} + \mathbf{WU} \left( +#' \mathbf{U}'\mathbf{WU}\right)^{-1}\mathbf{d}_{\widehat{y}},} +#' which makes it clear that the reconciliation formula simply adjusts the +#' vector \mjseqn{\widehat{\mathbf{y}}} with a linear +#' combination -- according to the smoothing matrix +#' \mjseqn{\mathbf{L} = \mathbf{WU} \left(\mathbf{U}’\mathbf{WU}\right)^{-1}} -- +#' of the coherency errors of the base forecasts. +#' +#' In addition, if the error term \mjseqn{\mathbf{\varepsilon}} is gaussian, the reconciliation +#' error \mjseqn{\widetilde{\varepsilon} = \widetilde{\mathbf{y}} - \mathbf{y}} is +#' a zero-mean gaussian vector with covariance matrix +#' \mjsdeqn{E\left( \widetilde{\mathbf{y}} - \mathbf{y}\right) \left( +#' \widetilde{\mathbf{y}} - \mathbf{y}\right)' = \mathbf{W} - \mathbf{WU} +#' \left(\mathbf{U}'\mathbf{WU}\right)^{-1}\mathbf{U}' = \mathbf{MW}.} +#' +#' Hyndman et al. (2011, see also Wickramasuriya et al., 2019) propose an +#' equivalent, alternative formulation as for the reconciled estimates, obtained +#' by GLS estimation of the model +#' \mjsdeqn{\widehat{\mathbf{y}} = \mathbf{S}\mathbf{\beta} + \mathbf{\varepsilon},} +#' where \mjseqn{\mathbf{S}} is the \emph{structural summation matrix} describing +#' the aggregation relationships operating on \mjseqn{\mathbf{y}}, and +#' \mjseqn{\mathbf{\beta}} is a subset of \mjseqn{\mathbf{y}} containing the +#' target forecasts of the bottom level series, such that +#' \mjseqn{\mathbf{y} = \mathbf{S}\mathbf{\beta}}. Since the hypotheses on +#' \mjseqn{\mathbf{\varepsilon}} remain unchanged, +#' \mjsdeqn{\widetilde{\mathbf{\beta}} = \left(\mathbf{S}'\mathbf{W}^{-1}\mathbf{S} +#' \right)^{-1}\mathbf{S}'\mathbf{W}^{-1}\widehat{\mathbf{y}}} +#' is the best linear unbiased estimate of \mjseqn{\mathbf{\beta}}, and +#' the whole reconciled forecasts vector is given by +#' \mjsdeqn{\widetilde{\mathbf{y}} = \mathbf{S}\widetilde{\mathbf{\beta}} = \mathbf{SG} +#' \widehat{\mathbf{y}},} +#' where \mjseqn{\mathbf{G} = \left(\mathbf{S}'\mathbf{W}^{-1} +#' \mathbf{S}\right)^{-1}\mathbf{S}'\mathbf{W}^{-1}}, and \mjseqn{\mathbf{M}=\mathbf{SG}}. +#' This solution is used by \code{\link[FoReco]{htsrec}} when \code{type = "S"}. +#' +#' \strong{Bounds on the reconciled forecasts} +#' +#' The user may impose bounds on the reconciled forecasts. +#' The parameter \code{bounds} permits to consider lower (\mjseqn{\mathbf{a}}) and +#' upper (\mjseqn{\mathbf{b}}) bounds like \mjseqn{\mathbf{a} \leq +#' \widetilde{\mathbf{y}} \leq \mathbf{b}} such that: +#' \mjsdeqn{ \begin{array}{c} +#' a_1 \leq \widetilde{y}_1 \leq b_1 \cr +#' \dots \cr +#' a_n \leq \widetilde{y}_n \leq b_n \cr +#' \end{array} \Rightarrow +#' \mbox{bounds} = [\mathbf{a} \; \mathbf{b}] = +#' \left[\begin{array}{cc} +#' a_1 & b_1 \cr +#' \vdots & \vdots \cr +#' a_n & b_n \cr +#' \end{array}\right],} +#' where \mjseqn{a_i \in [- \infty, + \infty]} and \mjseqn{b_i \in [- \infty, + \infty]}. +#' If \mjseqn{y_i} is unbounded, the i-th row of \code{bounds} would be equal +#' to \code{c(-Inf, +Inf)}. +#' Notice that if the \code{bounds} parameter is used, \code{sol = "osqp"} must be used. +#' This is not true in the case of non-negativity constraints: +#' \itemize{ +#' \item \code{sol = "direct"}: first the base forecasts +#' are reconciled without non-negativity constraints, then, if negative reconciled +#' values are present, the \code{"osqp"} solver is used; +#' \item \code{sol = "osqp"}: the base forecasts are +#' reconciled using the \code{"osqp"} solver. +#' } +#' In this case it is not necessary to build a matrix containing +#' the bounds, and it is sufficient to set \code{nn = "TRUE"}. +#' +#' Non-negative reconciled forecasts may be obtained by setting \code{nn_type} alternatively as: +#' \itemize{ +#' \item \code{nn_type = "KAnn"} (Kourentzes and Athanasopoulos, 2021) +#' \item \code{nn_type = "sntz"} ("set-negative-to-zero") +#' \item \code{nn_type = "osqp"} (Stellato et al., 2020) #' } #' #' @return #' If the parameter \code{keep} is equal to \code{"recf"}, then the function -#' returns only the (\code{h x n}) reconciled forecasts matrix, otherwise (\code{keep="all"}) +#' returns only the (\mjseqn{h \times n}) reconciled forecasts matrix, otherwise (\code{keep="all"}) #' it returns a list that mainly depends on what type of representation (\code{type}) -#' and methodology (\code{sol}) have been used: -#' \item{\code{recf}}{(\code{h x n}) reconciled forecasts matrix.} -#' \item{\code{W}}{Covariance matrix used for forecast reconciliation.} +#' and solution technique (\code{sol}) have been used: +#' \item{\code{recf}}{(\mjseqn{h \times n}) reconciled forecasts matrix, \mjseqn{\widetilde{\mathbf{Y}}}.} +#' \item{\code{W}}{Covariance matrix used for forecast reconciliation, \mjseqn{\mathbf{W}}.} #' \item{\code{nn_check}}{Number of negative values (if zero, there are no values below zero).} #' \item{\code{rec_check}}{Logical value: has the hierarchy been respected?} -#' \item{\code{M} (\code{type="M"} and \code{type="direct"})}{Projection matrix (projection approach)} -#' \item{\code{G} (\code{type="S"} and \code{type="direct"})}{Projection matrix (structural approach).} -#' \item{\code{S} (\code{type="S"} and \code{type="direct"})}{Cross-sectional summing matrix, \eqn{\textbf{S}}{S}.} +#' \item{\code{varf} (\code{type="direct"})}{(\mjseqn{n \times 1}) reconciled forecasts variance vector for \mjseqn{h=1}, \mjseqn{\mbox{diag}(\mathbf{MW}}).} +#' \item{\code{M} (\code{type="direct"})}{Projection matrix, \mjseqn{\mathbf{M}} (projection approach).} +#' \item{\code{G} (\code{type="S"} and \code{type="direct"})}{Projection matrix, \mjseqn{\mathbf{G}} (structural approach, \mjseqn{\mathbf{M}=\mathbf{S}\mathbf{G}}).} +#' \item{\code{S} (\code{type="S"} and \code{type="direct"})}{Cross-sectional summing matrix, \mjseqn{\mathbf{S}}.} #' \item{\code{info} (\code{type="osqp"})}{matrix with information in columns -#' for each forecast horizon \code{h} (rows): run time (\code{run_time}) number of iteration, -#' norm of primal residual (\code{pri_res}), status of osqp's solution (\code{status}) and -#' polish's status (\code{status_polish}).} +#' for each forecast horizon \mjseqn{h} (rows): run time (\code{run_time}), +#' number of iteration (\code{iter}), norm of primal residual (\code{pri_res}), +#' status of osqp's solution (\code{status}) and polish's status +#' (\code{status_polish}). It will also be returned with \code{nn = TRUE} if +#' a solver (see \code{nn_type}) will be used.} #' #' Only if \code{comb = "bu"}, the function returns \code{recf}, \code{S} and \code{M}. #' @@ -99,10 +216,19 @@ #' 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., Marini, M. (2011), Simultaneous and two-step reconciliation of +#' systems of time series: methodological and practical issues, +#' \emph{Journal of the Royal Statistical Society. Series C (Applied Statistics)}, +#' 60, 2, 143-164 +#' #' Hyndman, R.J., Ahmed, R.A., Athanasopoulos, G., Shang, H.L.(2011), #' Optimal combination forecasts for hierarchical time series, #' \emph{Computational Statistics & Data Analysis}, 55, 9, 2579-2589. #' +#' Kourentzes, N., Athanasopoulos, G. (2021), +#' Elucidate structure in intermittent demand series, +#' \emph{European Journal of Operational Research}, 288, 1, pp. 141–152. +#' #' \enc{Schäfer}{Schafer}, J.L., Opgen-Rhein, R., Zuber, V., Ahdesmaki, M., #' Duarte Silva, A.P., Strimmer, K. (2017), \emph{Package `corpcor'}, R #' package version 1.6.9 (April 1, 2017), \href{https://CRAN.R-project.org/package=corpcor}{https://CRAN.R-project.org/package= corpcor}. @@ -111,11 +237,12 @@ #' Matrix Estimation and Implications for Functional Genomics, \emph{Statistical #' Applications in Genetics and Molecular Biology}, 4, 1. #' -#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2018). OSQP: -#' An Operator Splitting Solver for Quadratic Programs, \href{https://arxiv.org/abs/1711.08013}{arXiv:1711.08013}. +#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2020). OSQP: +#' An Operator Splitting Solver for Quadratic Programs, \emph{Mathematical Programming Computation}, +#' 12, 4, 637-672. #' #' Stellato, B., Banjac, G., Goulart, P., Boyd, S., Anderson, E. (2019), OSQP: -#' Quadratic Programming Solver using the 'OSQP' Library, R package version 0.6.0.3 +#' Quadratic Programming Solver using the `OSQP' Library, R package version 0.6.0.3 #' (October 10, 2019), \href{https://CRAN.R-project.org/package=osqp}{https://CRAN.R-project.org/package=osqp}. #' #' van Erven, T., Cugliari, J. (2015), Game-theoretically Optimal Reconciliation @@ -127,26 +254,36 @@ #' reconciliation for hierarchical and grouped time series through trace minimization, #' \emph{Journal of the American Statistical Association}, 114, 526, 804-819. #' -#' @keywords reconciliation +#' @keywords bottom-up +#' @family reconciliation procedures +#' #' @examples #' data(FoReco_data) #' # monthly base forecasts -#' id <- which(simplify2array(strsplit(colnames(FoReco_data$base), -#' split = "_"))[1, ] == "k1") +#' id <- which(simplify2array(strsplit(colnames(FoReco_data$base),split = "_"))[1, ] == "k1") #' mbase <- t(FoReco_data$base[, id]) #' # monthly residuals -#' id <- which(simplify2array(strsplit(colnames(FoReco_data$res), -#' split = "_"))[1, ] == "k1") +#' id <- which(simplify2array(strsplit(colnames(FoReco_data$res),split = "_"))[1, ] == "k1") #' mres <- t(FoReco_data$res[, id]) #' obj <- htsrec(mbase, C = FoReco_data$C, comb = "shr", res = mres) #' +#' # FoReco is able to work also with covariance matrix that are not equal +#' # across all the forecast horizon. For example, we can consider the +#' # normalized squared differences (see Di Fonzo and Marini, 2011) where +#' # Wh = diag(|yh|): +#' Wh <- lapply(split(mbase, row(mbase)), function(x) diag(abs(x))) +#' +#' # Now we can introduce the list of the covariance matrix in htsrec throught +#' # the parameter "W" and setting comb = "w". +#' objh <- htsrec(mbase, C = FoReco_data$C, W = Wh, comb = "w") +#' #' @export #' #' @import Matrix osqp methods #' htsrec <- function(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FALSE, - type = "M", sol = "direct", nn = FALSE, keep = "list", - settings = osqpSettings(), bounds = NULL, W = NULL) { + type = "M", sol = "direct", keep = "list", nn = FALSE, + nn_type = "osqp", settings = list(), bounds = NULL, W = NULL){ if(missing(comb)){ stop("The argument comb is not specified.", call. = FALSE) @@ -159,13 +296,14 @@ htsrec <- function(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FALSE, #' @export htsrec.default <- function(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FALSE, - type = "M", sol = "direct", nn = FALSE, keep = "list", - settings = osqpSettings(), bounds = NULL, W) { + type = "M", sol = "direct", keep = "list", nn = FALSE, + nn_type = "osqp", settings = list(), bounds = NULL, W) { comb <- match.arg(comb, c("bu", "ols", "struc", "wls", "shr", "sam", "w")) type <- match.arg(type, c("M", "S")) keep <- match.arg(keep, c("list", "recf")) + nn_type <- match.arg(nn_type, c("osqp", "KAnn", "fbpp", "sntz")) # base forecasts condition if (missing(basef)) { @@ -177,57 +315,48 @@ htsrec.default <- function(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FA n <- NCOL(basef) - # Using Ut AND nb or C + # Using Ut or C if (missing(C)) { - if (missing(Ut) | missing(nb)) { - stop("Please, give C (or Ut AND nb).", call. = FALSE) - } - - if (comb == "bu" | comb == "struc") { - stop("bu and struc use C matrix.", call. = FALSE) + if (missing(Ut)) { + stop("Please, give C or Ut.", call. = FALSE) + } else if(missing(nb)){ + hts <- hts_tools(Ut = Ut, h = 1, sparse = TRUE) + } else { + hts <- hts_tools(Ut = Ut, nb = nb, h = 1, sparse = TRUE) } + } else { + hts <- hts_tools(C = C, h = 1, sparse = TRUE) + } - if (type == "S") { - stop("Type = S use C matrix.", call. = FALSE) - } + n <- hts$n + na <- hts$na + nb <- hts$nb + C <- hts$C + S <- hts$S + Ut <- hts$Ut - if (n < 3) { - stop("The hierarchy must have, at least, three series.", call. = FALSE) + if(nn){ + if(nn_type == "fbpp" | nn_type == "KAnn"){ + type = "M" } + } - if (NCOL(Ut) != n) { - stop("Incorrect dimension of Ut or basef (they don't have same columns).", call. = FALSE) + if(is.null(S)){ + if (type == "S") { + stop("Type = S needs C matrix. \nPlease consider using ut2c() to get the right C when working with Ut.", call. = FALSE) } - if (n <= nb) { - stop("n <= nb, total number of TS is less (or equal) than the number of bottom TS.", call. = FALSE) + if(nn_type == "sntz"){ + stop("nn_type = sntz needs C matrix.\nPlease consider using ut2c() to get the right C when working with Ut.", call. = FALSE) } - na <- n - nb - } else { - if (!(is.matrix(C) | is(C, "Matrix"))) stop("C must be a matrix.", call. = FALSE) - - C <- Matrix(C, sparse = TRUE) - - nb <- NCOL(C) - na <- NROW(C) - - if (comb != "bu" & nb != n) { - if (n < 3) { - stop("The hierarchy must have, at least, three series.", call. = FALSE) - } - - if (n <= nb) { - stop("n <= nb, total number of TS is less (or equal) than the number of bottom TS.", call. = FALSE) - } - - if (na + nb != n) { - stop("na + nb != n, matrix C or basef has incorrect dimension.", call. = FALSE) - } + if (comb == "bu" | comb == "struc") { + stop("Param comb equal to bu or struc needs C matrix.\nPlease consider using ut2c() to get the right C when working with Ut.", call. = FALSE) } + } - S <- rbind(C, .sparseDiagonal(nb)) - Ut <- cbind(.sparseDiagonal(na), -C) + if (NCOL(basef) != n) { + stop("Incorrect dimension of Ut or basef (they don't have same columns).", call. = FALSE) } if (any(comb == c("wls", "shr", "sam"))) { @@ -254,19 +383,29 @@ htsrec.default <- function(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FA switch(comb, bu = { - if (n == nb) { - outf <- basef %*% t(S) - } else { - outf <- basef[, (na + 1):n] %*% t(S) + if (n != nb) { + basef_bu <- basef[, (na + 1):n] + names_bu <- if (is.null(colnames(basef))) paste("serie", 1:n, sep = "") else colnames(basef) } + if(nn){ + basef_bu <- basef_bu * (basef_bu > 0) + } + + outf <- basef_bu %*% t(S) + rownames(outf) <- paste("h", 1:NROW(outf), sep = "") - colnames(outf) <- if (is.null(colnames(basef))) paste("serie", 1:n, sep = "") else colnames(basef) + if(length(names_bu) == NCOL(outf)){ + colnames(outf) <- names_bu + }else{ + colnames(outf) <- paste("serie", 1:n, sep = "") + } + if (keep == "list") { - return(list(recf = outf, S = S, M = S %*% cbind(matrix(0, nb, na), diag(nb)))) + return(list(recf = as.matrix(outf), S = S, M = S %*% cbind(matrix(0, nb, na), diag(nb)))) } else { - return(outf) + return(as.matrix(outf)) } }, ols = @@ -296,12 +435,12 @@ htsrec.default <- function(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FA if (type == "S") { rec_sol <- recoS( basef = basef, W = W, S = S, sol = sol, nn = nn, keep = keep, - settings = settings, b_pos = b_pos, bounds = bounds + settings = settings, b_pos = b_pos, bounds = bounds, nn_type = nn_type ) } else { rec_sol <- recoM( - basef = basef, W = W, H = Ut, sol = sol, nn = nn, keep = keep, - settings = settings, b_pos = b_pos, bounds = bounds + basef = basef, W = W, Ht = Ut, sol = sol, nn = nn, keep = keep, S = S, + settings = settings, b_pos = b_pos, bounds = bounds, nn_type = nn_type ) } diff --git a/R/iterec.R b/R/iterec.R index 8f0bde9..750aa43 100755 --- a/R/iterec.R +++ b/R/iterec.R @@ -1,48 +1,62 @@ #' @title Iterative heuristic cross-temporal forecast reconciliation #' #' @description +#' \loadmathjax #' Iterative procedure which produces cross-temporally reconciled #' forecasts by alternating forecast reconciliation along one single dimension #' (either cross-sectional or temporal) at each iteration step. \strong{Each iteration} -#' consists in the first two steps of the heuristic KA procedure, so the forecasts are -#' reconciled by alternating cross-sectional (contemporaneous) reconciliation, and -#' reconciliation through temporal hierarchies in a cyclic fashion. -#' The choice of the dimension along with the first reconciliation step in each +#' consists in the first two steps of the heuristic procedure by Kourentzes and Athanasopoulos (2019), +#' so the forecasts are reconciled by alternating cross-sectional (contemporaneous) reconciliation, +#' and reconciliation through temporal hierarchies in a cyclic fashion. +#' The choice of the dimension along which the first reconciliation step in each #' iteration is performed is up to the user (param \code{start_rec}), and there is #' no particular reason why one should perform the temporal reconciliation first, and #' the cross-sectional reconciliation then. #' The iterative procedure allows the user to get non-negative reconciled forecasts. #' -#' @param basef (\code{n x h(k* + m)}) matrix of base forecasts to be reconciled; -#' \code{n} is the total number of variables, \code{m} is the highest frequency, -#' \code{k*} is the sum of (\code{p-1}) factors of \code{m}, excluding \code{m}, -#' and \code{h} is the forecast horizon. Each row identifies, a time series, and the forecasts -#' are ordered as [lowest_freq' ... highest_freq']'. -#' @param hts_comb,thf_comb Type of covariance matrix (respectively (\code{n x n}) and -#' (\code{(k* + m) x (k* + m)})) to be used in the cross-sectional and temporal reconciliation, -#' see more in \code{comb} param of \code{\link[FoReco]{htsrec}} and \code{\link[FoReco]{thfrec}}. -#' @param res (\code{n x N(k* + m)}) matrix containing the residuals at all the -#' temporal frequencies ordered [lowest_freq' ... highest_freq']' (columns) for -#' each variable (row), needed to estimate the covariance matrix when \code{hts_comb =} -#' \code{\{"wls",} \code{"shr",} \code{"sam"\}} and/or \code{hts_comb =} \code{\{"wlsv",} -#' \code{"wlsh",} \code{"acov",} \code{"strar1",} \code{"sar1",} \code{"har1",} -#' \code{"shr",} \code{"sam"\}}. The row must be in the same order as \code{basef}. +#' @param basef (\mjseqn{n \times h(k^\ast+m)}) matrix of base forecasts to be +#' reconciled, \mjseqn{\widehat{\mathbf{Y}}}; \mjseqn{n} is the total number of variables, +#' \mjseqn{m} is the highest time frequency, \mjseqn{k^\ast} is the sum of (a +#' subset of) (\mjseqn{p-1}) factors of \mjseqn{m}, excluding \mjseqn{m}, and +#' \mjseqn{h} is the forecast horizon for the lowest frequency time series. +#' Each row identifies a time series, and the forecasts are ordered as +#' [lowest_freq' ... highest_freq']'. +#' @param hts_comb,thf_comb Type of covariance matrix (respectively +#' (\mjseqn{n \times n}) and (\mjseqn{(k^\ast + m) \times (k^\ast + m)})) to +#' be used in the cross-sectional and temporal reconciliation, see more in +#' \code{comb} param of \code{\link[FoReco]{htsrec}()} and +#' \code{\link[FoReco]{thfrec}()}. +#' @param res (\mjseqn{n \times N(k^\ast + m)}) matrix containing the residuals +#' at all the temporal frequencies ordered [lowest_freq' ... highest_freq']' +#' (columns) for each variable (row), needed to estimate the covariance matrix +#' when \code{hts_comb =} \code{\{"wls",} \code{"shr",} \code{"sam"\}} and/or +#' \code{hts_comb =} \code{\{"wlsv",} \code{"wlsh",} \code{"acov",} +#' \code{"strar1",} \code{"sar1",} \code{"har1",} \code{"shr",} \code{"sam"\}}. +#' The row must be in the same order as \code{basef}. #' @param tol Convergence tolerance (\code{1e-5}, \emph{default}). -#' @param itmax Max number of iteration (\code{100}, \emph{default}) (old version \code{maxit}). +#' @param itmax Max number of iteration (\code{100}, \emph{default}) +#' (old version \code{maxit}). #' @param start_rec Dimension along with the first reconciliation step in each -#' iteration is performed: it start from temporal reconciliation with "\code{thf}" (\emph{default}), -#' from cross-sectional with "\code{hts}" and it does both reconciliation with "\code{auto}". -#' @param note If \code{note = TRUE} (\emph{default}) the function writes some notes to the console, -#' otherwise no note is produced (also no plot). -#' @param plot Some useful plots: \code{"mti"} (\emph{default}) marginal trend inconsistencies, -#' \code{"pat"} step by step inconsistency pattern for each iteration, \code{"distf"} distance -#' forecasts iteration i and i-1, \code{"all"} all the plots. -#' @param ... any other options useful for \code{\link[FoReco]{htsrec}} and -#' \code{\link[FoReco]{thfrec}}, e.g. \code{m}, \code{C} (or \code{Ut} and \code{nb}), -#' \code{nn} (for non negativity reconciliation only at first step), ... +#' iteration is performed: it start from temporal reconciliation with +#' "\code{thf}" (\emph{default}), from cross-sectional with "\code{hts}" and +#' it does both reconciliation with "\code{auto}". +#' @param norm Norm used to calculate the temporal and the cross-sectional +#' incoherence. There are two alternatives: "\code{inf}" (\mjseqn{\max\{|x_1|, +#' |x_2|,\dots\}}, \emph{default}) or "\code{one}" (\mjseqn{\sum |x_i|}). +#' @param note If \code{note = TRUE} (\emph{default}) the function writes some +#' notes to the console, otherwise no note is produced (also no plot). +#' @param plot Some useful plots: \code{"mti"} (\emph{default}) marginal trend +#' inconsistencies, \code{"pat"} step by step inconsistency pattern for each +#' iteration, \code{"distf"} distance forecasts iteration i and i-1, +#' \code{"all"} all the plots. +#' @param ... any other options useful for \code{\link[FoReco]{htsrec}()} and +#' \code{\link[FoReco]{thfrec}()}, e.g. \code{m}, \code{C} (or \code{Ut} and +#' \code{nb}), \code{nn} (for non negativity reconciliation only at first step), +#' \code{mse}, \code{corpcor}, \code{type}, \code{sol}, \code{settings}, +#' \code{W}, \code{Omega},... #' #' @details -#' The above procedure can be seen as an extension of the well known iterative proportional +#' This reconciliation procedure can be seen as an extension of the well known iterative proportional #' fitting procedure (Deming and Stephan, 1940, Johnston and Pattie, 1993), also known as #' RAS method (Miller and Blair, 2009), to adjust the internal cell values of a #' two-dimensional matrix until they sum to some predetermined row and column @@ -63,7 +77,7 @@ #' } #' #' @return \code{iterec} returns a list with: -#' \item{\code{recf}}{(\code{n x h(k* + m)}) matrix of reconciled forecasts.} +#' \item{\code{recf}}{(\mjseqn{n \times h(k^\ast + m)}) reconciled forecasts matrix, \mjseqn{\widetilde{\textbf{Y}}}.} #' \item{\code{d_cs}}{Cross-sectional incoherence at each iteration.} #' \item{\code{d_te}}{Temporal incoherence at each iteration.} #' \item{\code{start_rec}}{Starting coherence dimension (thf or hts).} @@ -93,38 +107,50 @@ #' #' \enc{Schäfer}{Schafer}, J.L., Opgen-Rhein, R., Zuber, V., Ahdesmaki, M., #' Duarte Silva, A.P., Strimmer, K. (2017), \emph{Package `corpcor'}, R -#' package version 1.6.9 (April 1, 2017), \href{https://CRAN.R-project.org/package=corpcor}{https://CRAN.R-project.org/package= corpcor}. +#' package version 1.6.9 (April 1, 2017), +#' \href{https://CRAN.R-project.org/package=corpcor}{https://CRAN.R-project.org/package= corpcor}. #' #' \enc{Schäfer}{Schafer}, J.L., Strimmer, K. (2005), A Shrinkage Approach to Large-Scale Covariance #' Matrix Estimation and Implications for Functional Genomics, \emph{Statistical #' Applications in Genetics and Molecular Biology}, 4, 1. #' -#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2018). OSQP: -#' An Operator Splitting Solver for Quadratic Programs, \href{https://arxiv.org/abs/1711.08013}{arXiv:1711.08013}. +#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2020). OSQP: +#' An Operator Splitting Solver for Quadratic Programs, \emph{Mathematical Programming Computation}, +#' 12, 4, 637-672. #' #' Stellato, B., Banjac, G., Goulart, P., Boyd, S., Anderson, E. (2019), OSQP: -#' Quadratic Programming Solver using the 'OSQP' Library, R package version 0.6.0.3 +#' Quadratic Programming Solver using the `OSQP' Library, R package version 0.6.0.3 #' (October 10, 2019), \href{https://CRAN.R-project.org/package=osqp}{https://CRAN.R-project.org/package=osqp}. #' -#' @keywords reconciliation heuristic +#' @keywords heuristic +#' @family reconciliation procedures +#' #' @examples #' \donttest{ #' data(FoReco_data) -#' obj <- iterec(FoReco_data$base, note=TRUE, +#' obj <- iterec(FoReco_data$base, note = FALSE, #' m = 12, C = FoReco_data$C, thf_comb = "acov", #' hts_comb = "shr", res = FoReco_data$res, start_rec = "thf") #' } #' #' @usage iterec(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, -#' start_rec = "thf", note = TRUE, plot = "mti", ...) +#' start_rec = "thf", norm = "inf", note = TRUE, plot = "mti", ...) #' #' @export iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, - start_rec = "thf", note = TRUE, plot = "mti", ...) { + start_rec = "thf", norm = "inf", note = TRUE, plot = "mti", ...) { arg_input <- list(...) start_rec <- match.arg(start_rec, c("thf", "hts", "auto")) + norm <- match.arg(norm, c("one", "inf")) + + if(norm == "one"){ + ite_norm <- function(x) sum(x) + }else{ + ite_norm <- function(x) max(x) + } + if(any(names(arg_input)=="maxit")){ maxit <- arg_input$maxit }else{ @@ -152,13 +178,13 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, stop("the argument hts_comb is not specified", call. = FALSE) } - if(any(names(arg_input)=="C")){ + if(any(names(arg_input)=="Ut") & any(names(arg_input)=="nb")){ + Ut <- arg_input$Ut + r_u <- NROW(Ut) + }else if(any(names(arg_input)=="C")){ C <- arg_input$C r_u <- NROW(C) Ut <- cbind(diag(1, r_u), -C) - }else if(any(names(arg_input)=="Ut") & any(names(arg_input)=="nb")){ - Ut <- arg_input$Ut - r_u <- NROW(Ut) }else{ stop("Please, give C (or Ut AND nb)", call. = FALSE) } @@ -167,7 +193,7 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, # Tools for cross-sectional reconciliation tools <- thf_tools(m) kset <- tools$kset - m <- max(kset) + m <- tools$m kt <- tools$kt ks <- tools$ks @@ -178,13 +204,13 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, # Incoherence vector d_cs_thf <- rep(NA, maxit + 1) d_te_thf <- rep(NA, maxit + 1) - d_cs_thf[1] <- sum(abs(Ut %*% basef)) - d_te_thf[1] <- sum(abs(Zt %*% t(basef))) + d_cs_thf[1] <- ite_norm(abs(Ut %*% basef)) + d_te_thf[1] <- ite_norm(abs(Zt %*% t(basef))) d_cs_hts <- rep(NA, maxit + 1) d_te_hts <- rep(NA, maxit + 1) - d_cs_hts[1] <- sum(abs(Ut %*% basef)) - d_te_hts[1] <- sum(abs(Zt %*% t(basef))) + d_cs_hts[1] <- ite_norm(abs(Ut %*% basef)) + d_te_hts[1] <- ite_norm(abs(Zt %*% t(basef))) # Check: reconciliation needed? if (d_cs_thf[1] < tol & d_te_thf[1] < tol) { @@ -239,32 +265,33 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, # Time starting start_time <- Sys.time() if (note == TRUE) { - message("---------------------------------------") + message(paste0(rep("-", 80), collapse = "")) } if (start_thf) { if (note == TRUE) { - message("Iter # (Cross-sectional incoherence) (Temporal incoherence)", sep = "") - message("thf: ", formatC(0, width = 3, format = "d", flag = "0"), " (", d_cs_thf[1], ") (", - d_te_thf[1], ")", - sep = "" - ) + message("Iter ", formatC("#", width = 7), + " | ",format("Cross-sec. incoherence", width = 30, justify = "centre"), + " | ",format("Temporal incoherence", width = 30, justify = "centre"), " |", sep = "") + message("thf: ", formatC(0, width = 7, format = "d"), " | ", + console_incoherence(x = d_cs_thf[1], y = d_cs_thf[1]), " | ", + console_incoherence(x = d_te_thf[1], y = d_te_thf[1]), " |", sep = "") } for (i in 1:maxit) { # Step1 Y1_thf <- step_thf(basef = Y2_thf, res = res, arg_input = arg_input, thf_comb = thf_comb) - d_cs_thf[i + 1] <- sum(abs(Ut %*% Y1_thf)) - check <- sum(abs(Zt %*% t(Y1_thf))) + d_cs_thf[i + 1] <- ite_norm(abs(Ut %*% Y1_thf)) + check <- ite_norm(abs(Zt %*% t(Y1_thf))) d_thf <- rbind(d_thf, c(i, d_cs_thf[i + 1], check)) if (check > tol) { Y2_thf <- Y1_thf flag_thf <- -2 warning(paste("thf: Program (starting from thf) stops at iteration number ", i, - "! \nthf: Temporal reconciliation does not work. (", check, ")", - sep = "" - ), call. = FALSE) + "! \nthf: Temporal reconciliation does not work. (", + console_incoherence(x = check, y = d_te_thf[1], m = 8), ")", + sep = ""), call. = FALSE) break } else if (i > 1) { if (d_cs_thf[(i + 1)] > d_cs_thf[i] & d_cs_thf[i] > d_cs_thf[(i - 1)] & flag_thf < 2) { @@ -278,25 +305,24 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, Y2_thf <- step_hts(basef = Y1_thf, kset = kset, h = h, res = res, arg_input = arg_input, hts_comb = hts_comb) - d_te_thf[i + 1] <- sum(abs(Zt %*% t(Y2_thf))) - check <- sum(abs(Ut %*% Y2_thf)) + d_te_thf[i + 1] <- ite_norm(abs(Zt %*% t(Y2_thf))) + check <- ite_norm(abs(Ut %*% Y2_thf)) d_thf <- rbind(d_thf, c(i, check, d_te_thf[i + 1])) if (check > tol) { flag_thf <- -2 - warning(paste("thf: Program (starting from thf) stops at iteration number ", i, - "! \nthf: Cross-sectional reconciliation does not work.(", check, ")", - sep = "" - ), call. = FALSE) + warning(paste("thf: Program (starting from thf) stops at iteration number ", + i, "! \nthf: Cross-sectional reconciliation does not work. (", + console_incoherence(x = check, y = d_cs_thf[1], m = 8), ")", + sep = ""), call. = FALSE) break } else if (d_te_thf[i + 1] < tol) { if (flag_thf == -1) flag_thf <- 0 if (note == TRUE) { - message("thf: Convergence (starting from thf) achieved at iteration number ", i, - "! \nthf: Temporal incoherence ", d_te_thf[i + 1], " < ", tol, - " tolerance", - sep = "" - ) + message("thf: Convergence (starting from thf) achieved at iteration number ", + i, "! \nthf: Temporal incoherence ", + console_incoherence(x = d_te_thf[i + 1], y = d_te_thf[1], m = 8), + " < ", tol, " tolerance", sep = "") } rel <- ((Y2_thf - Yi_thf) / abs(Yi_thf)) dist_thf[i, ] <- c( @@ -323,23 +349,22 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, ) Yi_thf <- Y2_thf if (note == TRUE) { - message("thf: ", formatC(i, width = 3, format = "d", flag = "0"), " (", d_cs_thf[i + 1], ") (", - d_te_thf[i + 1], ")", - sep = "" - ) + message("thf: ", formatC(i, width = 7, format = "d"), " | ", + console_incoherence(x = d_cs_thf[i + 1], y = d_cs_thf[1]), " | ", + console_incoherence(x = d_te_thf[i + 1], y = d_te_thf[1]), " |", + sep = "") } } if ((flag_thf == -1 | i == maxit) & start_rec != "hts") { flag_thf <- -1 - warning("thf: Convergence NOT achieved (starting from thf)! Maximum number of iterations reached (", - maxit, ")", - sep = "", call. = FALSE - ) + warning("thf: Convergence NOT achieved (starting from thf)!", + "\nthf: Maximum number of iterations reached (", + maxit, ")", sep = "", call. = FALSE) } dist_thf <- dist_thf[1:i, , drop = FALSE] if (note == TRUE) { - message("---------------------------------------") + message(paste0(rep("-", 80), collapse = "")) } } @@ -347,11 +372,13 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, if (start_hts) { if (note == TRUE) { - message("Iter # (Cross-sectional incoherence) (Temporal incoherence)", sep = "") - message("hts: ", formatC(0, width = 3, format = "d", flag = "0"), " (", d_cs_hts[1], ") (", - d_te_hts[1], ")", - sep = "" - ) + message("Iter ", formatC("#", width = 7), + " | ",format("Cross-sec. incoherence", width = 30, justify = "centre"), + " | ",format("Temporal incoherence", width = 30, justify = "centre"), " |", sep = "") + message("hts: ", formatC(0, width = 7, format = "d"), " | ", + console_incoherence(x = d_cs_hts[1], y = d_cs_hts[1]), " | ", + console_incoherence(x = d_te_hts[1], y = d_te_hts[1]), " |", + sep = "") } for (j in 1:maxit) { @@ -359,16 +386,16 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, Y1_hts <- step_hts(basef = Y2_hts, kset = kset, h = h, res = res, arg_input = arg_input, hts_comb = hts_comb) - d_te_hts[j + 1] <- sum(abs(Zt %*% t(Y1_hts))) - check <- sum(abs(Ut %*% Y1_hts)) + d_te_hts[j + 1] <- ite_norm(abs(Zt %*% t(Y1_hts))) + check <- ite_norm(abs(Ut %*% Y1_hts)) d_hts <- rbind(d_hts, c(j, check, d_te_hts[j + 1])) if (check > tol) { flag_hts <- -2 warning(paste("hts: Program (starting from hts) stops at iteration number ", j, - "! \nhts: Cross-sectional reconciliation does not work.(", check, ")", - sep = "" - ), call. = FALSE) + "! \nhts: Cross-sectional reconciliation does not work. (", + console_incoherence(x = check, y = d_cs_hts[1]), ")", + sep = ""), call. = FALSE) break } else if (j > 1) { if (d_te_hts[(j + 1)] > d_te_hts[j] & d_te_hts[j] > d_te_hts[(j - 1)] & flag_hts < 2) { @@ -381,25 +408,24 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, # Step 2 Y2_hts <- step_thf(basef = Y1_hts, res = res, arg_input = arg_input, thf_comb = thf_comb) - d_cs_hts[j + 1] <- sum(abs(Ut %*% Y2_hts)) - check <- sum(abs(Zt %*% t(Y2_hts))) + d_cs_hts[j + 1] <- ite_norm(abs(Ut %*% Y2_hts)) + check <- ite_norm(abs(Zt %*% t(Y2_hts))) d_hts <- rbind(d_hts, c(j, d_cs_hts[j + 1], check)) if (check > tol) { flag_hts <- -2 - warning(paste("hts: Program (starting from hts) stops at iteration number ", j, - "! \nhts: Temporal reconciliation does not work.(", check, ")", - sep = "" - ), call. = FALSE) + warning(paste("hts: Program (starting from hts) stops at iteration number ", + j, "! \nhts: Temporal reconciliation does not work. (", + console_incoherence(x = check, y = d_te_hts[1], m = 8), ")", + sep = ""), call. = FALSE) break } else if (d_cs_hts[j + 1] < tol) { if (flag_hts == -1) flag_hts <- 0 if (note == TRUE) { message("hts: Convergence (starting from hts) achieved at iteration number ", j, - "! \nhts: Cross-sectional incoherence ", d_cs_hts[j + 1], " < ", tol, - " tolerance", - sep = "" - ) + "! \nhts: Cross-sectional incoherence ", + console_incoherence(x = d_cs_hts[j + 1], y = d_cs_hts[1], m = 8), + " < ", tol, " tolerance", sep = "") } rel <- ((Y2_hts - Yi_hts) / abs(Yi_hts)) dist_hts[j, ] <- c( @@ -419,8 +445,9 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, } if (note == TRUE) { - message("hts: ", formatC(j, width = 3, format = "d", flag = "0"), " (", d_cs_hts[j + 1], ") (", - d_te_hts[j + 1], ")", + message("hts: ", formatC(j, width = 7, format = "d"), " | ", + console_incoherence(x = d_cs_hts[j + 1], y = d_cs_hts[1]), " | ", + console_incoherence(x = d_te_hts[j + 1], y = d_te_hts[1]), " |", sep = "" ) } @@ -437,14 +464,15 @@ iterec <- function(basef, thf_comb, hts_comb, res, itmax = 100, tol = 1e-5, if ((flag_hts == -1 | j == maxit) & start_rec != "thf") { flag_hts <- -1 - warning("hts: Convergence NOT achieved (starting from hts)! Maximum number of iterations reached (", + warning("hts: Convergence NOT achieved (starting from hts)!", + "\nhts: Maximum number of iterations reached (", maxit, ")", sep = "", call. = FALSE ) } dist_hts <- dist_hts[1:j, , drop = FALSE] if (note == TRUE) { - message("---------------------------------------") + message(paste0(rep("-", 80), collapse = "")) } } @@ -864,3 +892,9 @@ step_thf <- function(basef, res, arg_input, thf_comb) { dimnames(Y1) <- NULL return(Y1) } + + +console_incoherence <- function(x, y, m = 30){ + formatC(x, width = m, + format = ifelse(x<1,"e", "f"), digits = 2) +} diff --git a/R/lccrec.R b/R/lccrec.R new file mode 100644 index 0000000..9dd1f80 --- /dev/null +++ b/R/lccrec.R @@ -0,0 +1,1044 @@ +#' @title Level conditional coherent forecast reconciliation for genuine hierarchical/grouped time series +#' +#' @description +#' \loadmathjax +#' Forecast reconciliation procedure built on and extending the original +#' proposal by Hollyman et al. (2021). Level conditional coherent +#' reconciled forecasts may be computed in cross-sectional, temporal, and +#' cross-temporal frameworks. The reconciled forecasts are conditional to (i.e., +#' constrained by) the base forecasts of a specific upper level of the hierarchy +#' (exogenous constraints). The linear constraints linking the variables may be +#' dealt with endogenously as well (Di Fonzo and Girolimetto, 2021). +#' \emph{Combined Conditional Coherent} (CCC) +#' forecasts may be calculated as simple averages of LCC and bottom-up +#' reconciled forecasts, with either endogenous or exogenous constraints. +#' +#' @usage lccrec(basef, m, C, nl, weights, bnaive = NULL, const = "exogenous", +#' CCC = TRUE, nn = FALSE, nn_type = "osqp", settings = list(), ...) +#' +#' @param basef matrix/vector of base forecasts to be reconciled: +#' (\mjseqn{h \times n}) matrix in the cross-sectional framework; +#' (\mjseqn{h(k^\ast + m) \times 1}) vector in the temporal framework; +#' (\mjseqn{n \times h(k^\ast+m)}) matrix in the cross-temporal framework. +#' \mjseqn{n} is the total number of variables, \mjseqn{m} is the highest time +#' frequency, \mjseqn{k^\ast} is the sum of (a subset of) (\mjseqn{p-1}) factors +#' of \mjseqn{m}, excluding \mjseqn{m}, and \mjseqn{h} is the forecast horizon. +#' @param C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones (or a list +#' of matrices forming \mjseqn{\textbf{C} = [\textbf{C}_1' \; \textbf{C}_2' \; ... \; \textbf{C}_L']'}, +#' \mjseqn{1, ..., L} being the number of the cross-sectional upper levels. +#' @param nl (\mjseqn{L \times 1}) vector containing the number of time series +#' in each level of the hierarchy (\code{nl[1] = 1}). +#' @param m Highest available sampling frequency per seasonal cycle (max. order +#' of temporal aggregation, \mjseqn{m}), or a subset of the \mjseqn{p} factors +#' of \mjseqn{m}. +#' @param weights covariance matrix or a vector +#' @param bnaive matrix/vector of naive bts base forecasts +#' (e.g., seasonal averages, as in Hollyman et al., 2021): +#' (\mjseqn{h \times n_b}) matrix in the cross-sectional framework; +#' (\mjseqn{h m \times 1}) vector in the temporal framework; +#' (\mjseqn{n_b \times h m}) matrix in the cross-temporal framework. +#' Ignore it, if only basef are to be used as base forecasts. +#' @param const \strong{exo}genous (\emph{default}) or \strong{endo}genous +#' constraints +#' @param CCC Option to return Combined Conditional Coherent reconciled +#' forecasts (\emph{default} is TRUE). +#' @param nn Logical value: \code{TRUE} if non-negative reconciled forecasts +#' are wished. +#' @param nn_type Non-negative method: "osqp" (\emph{default}) or "sntz" +#' (\emph{set-negative-to-zero}, only if \code{CCC = TRUE}) with exogenous +#' constraints (\code{const = "exo"}); "osqp" (\emph{default}), "KAnn" +#' (only \code{type == "M"}) or "sntz" with endogenous constraints +#' (\code{const = "endo"}). +#' @param settings Settings for \pkg{osqp} (object \code{\link[osqp]{osqpSettings}}). +#' The default options are: \code{verbose = FALSE}, \code{eps_abs = 1e-5}, +#' \code{eps_rel = 1e-5}, \code{polish_refine_iter = 100} and \code{polish = TRUE}. +#' For details, see the \href{https://osqp.org/}{\pkg{osqp} documentation} +#' (Stellato et al., 2019). +#' @param ... Additional functional arguments passed to \link[FoReco]{htsrec} of +#' FoReco. +#' +#' @details +#' \strong{Cross-sectional hierarchies} +#' +#' To be as simple as possible, we fix the forecast horizon equal to 1. +#' Let the base forecasts be a vector \mjseqn{\widehat{\mathbf{y}} = +#' \left[\widehat{\mathbf{a}}' \; \widehat{\mathbf{b}}'\right]'}, where +#' vector \mjseqn{\widehat{\mathbf{a}}} consists of the sub-vectors forming each +#' of the upper \mjseqn{L} levels of the hierarchy/grouping: +#' \mjsdeqn{\widehat{\mathbf{a}} = \left[\begin{array}{c} +#' \widehat{a}_1 \cr \widehat{\mathbf{a}}_2 \cr \vdots \cr \widehat{\mathbf{a}}_L +#' \end{array}\right],} +#' where \mjseqn{\widehat{\mathbf{a}}_l}, \mjseqn{l=1,\ldots, L}, has dimension +#' \mjseqn{(n_l \times 1)} and \mjseqn{\sum_{l=1}^{L} n_l = n_a}. Denote +#' \mjseqn{\mathbf{C}_l} the \mjseqn{(n_l \times n_b)} matrix mapping the +#' bts into the level-\mjseqn{l} uts, then the aggregation matrix +#' \mjseqn{\mathbf{C}} may be written as +#' \mjsdeqn{\mathbf{C} = \left[\begin{array}{c} +#' \mathbf{C}_1 \cr\mathbf{C}_2 \cr \vdots \cr \mathbf{C}_L +#' \end{array}\right],} +#' where the generic matrix \mjseqn{\mathbf{C}_L} is (\mjseqn{n_L \times n_b}), +#' \mjseqn{l=1, \ldots, L}. Given a generic level \mjseqn{l}, the reconciled +#' forecasts coherent with the base forecasts of that level are the solution to +#' a quadratic minimization problem with linear +#' exogenous constraints (\code{const = "exo"}): +#' \mjsdeqn{\begin{array}{c}\widetilde{\mathbf{b}}_{l} = \arg\min_{\mathbf{b}} +#' \left(\mathbf{b} - \widehat{\mathbf{b}}\right)'\mathbf{W}_b^{-1} +#' \left(\mathbf{b} - \widehat{\mathbf{b}}\right) \quad \mbox{ s.t. } +#' \mathbf{C}_l\mathbf{b} = \widehat{\mathbf{a}}_l, \quad l=1, \ldots, L ,\cr +#' \Downarrow \cr +#' \widetilde{\mathbf{b}}_{l} = \widehat{\mathbf{b}} + +#' \mathbf{W}_b\mathbf{C}_l'\left(\mathbf{C}_l\mathbf{W}_b\mathbf{C}_l' +#' \right)^{-1}\left(\widehat{\mathbf{a}}_l -\mathbf{C}_l\widehat{\mathbf{b}} +#' \right), \quad l=1,\ldots,L,\end{array}} +#' where \mjseqn{\mathbf{W}_b} is a (\mjseqn{n_b \times n_b}) p.d. matrix +#' (in Hollyman et al., 2021, \mjseqn{\mathbf{W}_b} is diagonal). +#' If endogenous constraints (\code{const = "endo"}) are considered, +#' denote \mjseqn{\widehat{\mathbf{y}}_l = +#' \left[\widehat{\mathbf{a}}_l' \; \widehat{\mathbf{b}}'\right]'} and +#' \mjseqn{\mathbf{U}_l' = \left[\mathbf{I}_{n_l}' \; \mathbf{C}_l'\right]'}, +#' then +#' \mjsdeqn{\begin{array}{c}\widetilde{\mathbf{y}}_{l} = \arg\min_{\mathbf{y}_l} +#' \left(\mathbf{y}_l - \widehat{\mathbf{y}}_l\right)'\mathbf{W}_l^{-1} +#' \left(\mathbf{y}_l - \widehat{\mathbf{y}}_l\right) \quad \mbox{ s.t. } +#' \mathbf{U}_l'\mathbf{y}_l = \mathbf{0}, \quad l=1, \ldots, L ,\cr +#' \Downarrow \cr +#' \widetilde{\mathbf{y}}_{l} = \left(\mathbf{I}_{n_b+n_l} - +#' \mathbf{W}_l\mathbf{U}_l\left(\mathbf{U}_l'\mathbf{W}_l +#' \mathbf{U}_l\right)^{-1}\mathbf{U}_l'\right)\widehat{\mathbf{y}}_{l}, +#' \quad l=1,...,L, +#' \end{array}} +#' where \mjseqn{\mathbf{W}_l} is a (\mjseqn{n_l + n_b \times n_l + n_b}) +#' p.d. matrix. +#' Combined Conditional Coherent (CCC) forecasts are obtained by taking +#' the simple average of the \mjseqn{L} level conditional, and of the bottom-up +#' reconciled forecasts, respectively (Di Fonzo and Girolimetto, 2021): +#' \mjsdeqn{\widetilde{\mathbf{y}}_{CCC}=\frac{1}{L+1}\sum_{l=1}^{L+1} \mathbf{S}\widetilde{\mathbf{b}}_l,} +#' with \mjsdeqn{\widetilde{\mathbf{b}}_{L+1} = \widehat{\mathbf{b}}.} +#' +#' Non-negative reconciled forecasts may be obtained by setting \code{nn_type} +#' alternatively as: +#' \itemize{ +#' \item to apply non-negative constraints to each level: +#' \itemize{ +#' \item \code{nn_type = "KAnn"} (only \code{const = "endo"}) +#' \item \code{nn_type = "osqp"} +#' } +#' \item to apply non-negative constraints only to the CCC forecasts: +#' \itemize{ +#' \item \code{nn_type = "sntz"} ("set-negative-to-zero") +#' } +#' } +#' +#' \strong{Temporal hierarchies} +#' +#' The extension to the case of \strong{temporal hierarchies} is quite simple. +#' Using the same notation as in \code{\link[FoReco]{thfrec}()}, the +#' following `equivalences' hold between the symbols used +#' for the level conditional cross-sectional reconciliation and the ones +#' used in temporal reconciliation: +#' \mjsdeqn{L \equiv p-1, \quad (n_a, n_b, n) \equiv (k^*, m, k^*+m),} +#' and +#' \mjsdeqn{\mathbf{C} \equiv \mathbf{K} , \; +#' \mathbf{C}_1 \equiv \mathbf{K}_1 = \mathbf{1}_m', \; +#' \mathbf{C}_2 \equiv \mathbf{K}_2 = \mathbf{I}_{\frac{m}{k_{p-1}}},\; \ldots, \; +#' \mathbf{C}_L \equiv \mathbf{K}_{p-1} = \mathbf{I}_{\frac{m}{k_{2}}} \otimes +#' \mathbf{1}_{k_{2}}',\; \mathbf{S} \equiv \mathbf{R}.} +#' +#' The description of the \strong{cross-temporal extension} is currently under progress. +#' +#' @return The function returns the level reconciled forecasts +#' in case of an elementary hierarchy with one level. +#' Otherwise it returns a list with +#' \item{\code{recf}}{the first level reconciled forecasts, if \code{CCC = FALSE}, or the CCC reconciled forecasts.} +#' \item{\code{levrecf}}{list of level conditional reconciled forecasts (+ BU).} +#' \item{\code{info} (\code{type="osqp"})}{matrix with some useful indicators (columns) +#' for each forecast horizon \mjseqn{h} (rows): run time (\code{run_time}), number of iteration, +#' norm of primal residual (\code{pri_res}), status of osqp's solution (\code{status}) and +#' polish's status (\code{status_polish}).} +#' +#' @references +#' Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: +#' 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). +#' +#' Hollyman, R., Petropoulos, F., Tipping, M.E. (2021), Understanding Forecast Reconciliation, +#' \emph{European Journal of Operational Research} (in press). +#' +#' @family reconciliation procedures +#' +#' @examples +#' data(FoReco_data) +#' ### LCC and CCC CROSS-SECTIONAL FORECAST RECONCILIATION +#' # Cross sectional aggregation matrix +#' C <- rbind(FoReco_data$C, c(0,0,0,0,1)) +#' # monthly base forecasts +#' id <- which(simplify2array(strsplit(colnames(FoReco_data$base), split = "_"))[1, ] == "k1") +#' mbase <- t(FoReco_data$base[, id])[,c("T", "A","B","C","AA","AB","BA","BB","C")] +#' # residuals +#' id <- which(simplify2array(strsplit(colnames(FoReco_data$res), split = "_"))[1, ] == "k1") +#' mres <- t(FoReco_data$res[, id])[,c("T", "A","B","C","AA","AB","BA","BB","C")] +#' # covariance matrix of all the base forecasts, computed using the in-sample residuals +#' Wres <- cov(mres) +#' # covariance matrix of the bts base forecasts, computed using the in-sample residuals +#' Wb <- Wres[c("AA","AB", "BA", "BB", "C"), +#' c("AA","AB", "BA", "BB", "C")] +#' # bts seasonal averages +#' obs_1 <- FoReco_data$obs$k1 +#' bts_sm <- apply(obs_1, 2, function(x) tapply(x[1:168],rep(1:12, 14), mean)) +#' bts_sm <- bts_sm[,c("AA", "AB", "BA", "BB", "C")] +#' +#' ## EXOGENOUS CONSTRAINTS AND DIAGONAL COVARIANCE MATRIX (Hollyman et al., 2021) +#' # Forecast reconciliation for an elementary hierarchy: +#' # 1 top-level series + 5 bottom-level series (Level 2 base forecasts are not considered). +#' # The input is given by the base forecasts of the top and bottom levels series, +#' # along with a vector of positive weights for the bts base forecasts +#' exo_EHd <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")], +#' weights = diag(Wb), bnaive = bts_sm) +#' +#' # Level conditional reconciled forecasts +#' # recf/L1: Level 1 reconciled forecasts for the whole hierarchy +#' # L2: Middle-Out reconciled forecasts hinging on Level 2 conditional reconciled forecasts +#' # L3: Bottom-Up reconciled forecasts +#' exo_LCd <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = diag(Wb), +#' CCC = FALSE, bnaive = bts_sm) +#' +#' # Combined Conditional Coherent (CCC) reconciled forecasts +#' # recf: CCC reconciled forecasts matrix +#' # L1: Level 1 conditional reconciled forecasts for the whole hierarchy +#' # L2: Middle-Out reconciled forecasts hinging on Level 2 conditional reconciled forecasts +#' # L3: Bottom-Up reconciled forecasts +#' exo_CCCd <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = diag(Wb)) +#' +#' ## EXOGENOUS CONSTRAINTS AND FULL COVARIANCE MATRIX +#' # Simply substitute weights=diag(Wb) with weights=Wb +#' exo_EHf <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")], weights = Wb, bnaive = bts_sm) +#' exo_LCf <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = Wb, CCC = FALSE, bnaive = bts_sm) +#' exo_CCCf <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = Wb, bnaive = bts_sm) +#' +#' ## ENDOGENOUS CONSTRAINTS AND DIAGONAL COVARIANCE MATRIX +#' # parameters of function htsrec(), like "type" and "nn_type" may be used +#' +#' # Forecast reconciliation for an elementary hierarchy with endogenous constraints +#' W1 <- Wres[c("T","AA","AB", "BA", "BB", "C"), +#' c("T","AA","AB", "BA", "BB", "C")] +#' endo_EHd <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")], +#' weights = diag(W1), const = "endo", CCC = FALSE) +#' +#' # Level conditional reconciled forecasts with endogenous constraints +#' endo_LCd <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = diag(Wres), +#' const = "endo", CCC = FALSE) +#' +#' # Combined Conditional Coherent (CCC) reconciled forecasts with endogenous constraints +#' endo_CCCd <- lccrec(basef = mbase, C = C, nl = c(1,3), +#' weights = diag(Wres), const = "endo") +#' +#' ## ENDOGENOUS CONSTRAINTS AND FULL COVARIANCE MATRIX +#' # Simply substitute weights=diag(Wres) with weights=Wres +#' endo_EHf <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")], +#' weights = W1, +#' const = "endo") +#' endo_LCf <- lccrec(basef = mbase, C = C, nl = c(1,3), +#' weights = Wres, const = "endo", CCC = FALSE) +#' endo_CCCf <- lccrec(basef = mbase-40, C = C, nl = c(1,3), +#' weights = Wres, const = "endo") +#' +#' ### LCC and CCC TEMPORAL FORECAST RECONCILIATION +#' # top ts base forecasts ([lowest_freq' ... highest_freq']') +#' topbase <- FoReco_data$base[1, ] +#' # top ts residuals ([lowest_freq' ... highest_freq']') +#' topres <- FoReco_data$res[1, ] +#' Om_bt <- cov(matrix(topres[which(simplify2array(strsplit(names(topres), +#' "_"))[1,]=="k1")], ncol = 12, byrow = TRUE)) +#' t_exo_LCd <- lccrec(basef = topbase, m = 12, weights = diag(Om_bt), CCC = FALSE) +#' t_exo_CCCd <- lccrec(basef = topbase, m = 12, weights = diag(Om_bt)) +#' +#' ### LCC and CCC CROSS-TEMPORAL FORECAST RECONCILIATION +#' idr <- which(simplify2array(strsplit(colnames(FoReco_data$res), "_"))[1,]=="k1") +#' bres <- FoReco_data$res[-c(1:3), idr] +#' bres <- lapply(1:5, function(x) matrix(bres[x,], nrow=14, byrow = TRUE)) +#' bres <- do.call(cbind, bres) +#' ctbase <- FoReco_data$base[c("T", "A","B","C","AA","AB","BA","BB","C"),] +#' ct_exo_LCf <- lccrec(basef = ctbase, m = 12, C = C, nl = c(1,3), +#' weights = diag(cov(bres)), CCC = FALSE) +#' ct_exo_CCCf <- lccrec(basef = ctbase, m = 12, C = C, nl = c(1,3), +#' weights = diag(cov(bres)), CCC = TRUE) +#' +#' @export +#' @import Matrix +#' +lccrec <- function(basef, m, C, nl, weights, bnaive = NULL, + const = "exogenous", CCC = TRUE, nn = FALSE, + nn_type = "osqp", settings = list(), ...){ + + const <- match.arg(const, c("exogenous", "endogenous")) + + if(missing(C) | missing(m)){ + if(missing(m)){ + # lccrec cross-sectional + if(missing(C)){ + C <- NULL + } + out <- lccrec_hts(basef = basef, C = C, nl = nl, weights = weights, + bnaive = bnaive, nn = nn, settings = settings, + const = const, CCC = CCC, nn_type = nn_type, ...) + }else{ + # lccrec temporal + out <- lccrec_thf(basef = basef, m = m, weights = weights, + bnaive = bnaive, nn = nn, settings = settings, + const = const, CCC = CCC, nn_type = nn_type, ...) + } + }else{ + # lccrec cross-temporal + out <- lccrec_ctf(basef = basef, m = m, C = C, nl = nl, weights = weights, + bnaive = bnaive, nn = nn, settings = settings, + const = const, CCC = CCC, nn_type = nn_type, ...) + } + return(out) +} + +lccrec_thf <- function(basef, m, weights, bnaive = NULL, + const = "exogenous", CCC = TRUE, nn = FALSE, + nn_type = "osqp", settings = list(), ...){ + # m condition + if (missing(m)) { + stop("The argument m is not specified", call. = FALSE) + } + tools <- thf_tools(m) + kset <- tools$kset + m <- tools$m + p <- tools$p + kt <- tools$kt + ks <- tools$ks + + # matrix + K <- tools$K + R <- tools$R + Zt <- tools$Zt + + # base forecasts condition + if (missing(basef)) { + stop("The argument basef is not specified", call. = FALSE) + } + + if (NCOL(basef) != 1) { + stop("basef must be a vector", call. = FALSE) + } + + # Base Forecasts matrix + h <- length(basef) / kt + Dh <- Dmat(h = h, m = kset, n = 1) + basef <- matrix(Dh %*% basef, nrow = h, byrow = TRUE) + + idr <- rep(1:length(kset[-1]), rev(kset[-1])) + Kl <- lapply(unique(idr), function(i) K[idr==i, , drop = FALSE]) + + check_bil <- any(sapply(Kl, function(x) any(colSums(x)!=1))) + if(check_bil){ + warning("The hierarchy is not balanced. \nThis could produce results that are inconsistent with the reconciliation procedure.", call. = FALSE) + } + + idc <- rep(1:length(kset), rev(kset)) + lal <- lapply(unique(idc), function(i) basef[,idc==i, drop = FALSE]) + b <- lal[[max(idc)]] + lal <- lal[-max(idc)] + + if(!is.null(bnaive)){ + if(is.vector(bnaive)){ + if(length(bnaive) != m*h){ + stop("bnaive must be a vector (mh x 1)", call. = FALSE) + } + }else{ + stop("bnaive must be a vector (mh x 1)", call. = FALSE) + } + + Db <- Dmat(h = h, m = m, n = 1) + bm <- matrix(Db%*%bnaive, nrow = h, byrow = TRUE) + }else{ + bm <- b + } + + if(is.vector(weights)){ + weights <- Diagonal(x = weights) + } + + nn_type <- match.arg(nn_type, c("osqp", "fbpp", "KAnn", "sntz")) + if(nn_type == "sntz"){ + nn <- FALSE + } + + if(const=="exogenous"){ + if(NCOL(weights) != m | NROW(weights) != m){ + stop("weights must be a m x m matrix", call. = FALSE) + }else if(!is(weights, "Matrix")){ + weights <- Matrix(weights, sparse = TRUE) + } + out <- Map(function(al, cl){ + recoLEV(al = al, Wb = weights, bl = bm, cl = cl, nn = nn, + settings = settings) + }, al = lal, cl = Kl) + + bl <- lapply(out, function(x){ + extract_data(x = x, name = "recf") + }) + + info <- lapply(out, function(x){ + extract_data(x = x, name = "info") + }) + names(info) <- paste0("id_k_", kset[-length(kset)]) + info <- info[!is.na(info)] + + uts0 <- do.call(c, lapply(out, function(x){ + extract_data(x = x, name = "uts0") + })) + if(any(uts0)){ + message("Some aggregation orders (greater than m) have base forecasts less then 0,", + " then they have been set equal to zero") + } + }else{ + if(NCOL(weights) != kt | NROW(weights) != kt){ + stop("weights must be a (k* + m) x (k* + m) matrix", call. = FALSE) + }else if(!is(weights, "Matrix")){ + weights <- Matrix(weights, sparse = TRUE) + } + + id_nl <- rep(1:length(kset[-1]), rev(kset[-1])) + Ol_id <- lapply(1:length(kset[-1]), function(x) c(which(id_nl == x), (ks+1):(m+ks))) + out <- Map(function(al, cl, w){ + suppressMessages(htsrec(basef = cbind(al, bm), + C = cl, + comb = "w", + W = weights[w,w,drop = FALSE], + nn = nn, + settings = settings, + nn_type = nn_type, ...)) + }, al = lal, cl = Kl, w = Ol_id) + + bl <- lapply(out, function(x){ + out <- extract_data(x = x, name = "recf") + out <- out[,-(1:(NCOL(out)-m)), drop = FALSE] + }) + + info <- lapply(out, function(x){ + extract_data(x = x, name = "info") + }) + names(info) <- paste0("id_k_", kset[-length(kset)]) + info <- info[!is.na(info)] + } + + if(nn){ + b[b<0] <- 0 + } + bl[[length(bl)+1]] <- b + + yl <- lapply(bl, function(x){ + out <- as.matrix(x%*%t(R)) + return(out) + }) + + # return condition + names(yl) <- paste0("id_k_", kset) + if(CCC){ + out <- list() + if(nn_type == "sntz"){ + b_sntz <- Reduce("+", bl)/length(bl) + b_sntz <- b_sntz*(b_sntz>0) + ccc_recf <- as.matrix(b_sntz%*%t(R)) + }else{ + ccc_recf <- Reduce("+", yl)/length(yl) + } + + out$recf <- stats::setNames(as.vector(t(Dh) %*% as.vector(t(ccc_recf))), paste0("k", rep(kset, h * (m/kset)), "h", + do.call("c", as.list(sapply( + (m/kset) * h, + function(x) seq(1:x) + ))))) + out$levrecf <- lapply(yl, function(x){ + y <- as.vector(t(Dh) %*% as.vector(t(x))) + y <- stats::setNames(y, paste0("k", rep(kset, h * (m/kset)), "h", + do.call("c", as.list(sapply( + (m/kset) * h, + function(x) seq(1:x) + ))))) + return(y) + }) + if(length(info)>0){ + out$info <- info + } + return(out) + }else{ + if(length(info)>0){ + lev <- lapply(yl, function(x){ + y <- as.vector(t(Dh) %*% as.vector(t(x))) + y <- stats::setNames(y, paste0("k", rep(kset, h * (m/kset)), "h", + do.call("c", as.list(sapply( + (m/kset) * h, + function(x) seq(1:x) + ))))) + return(y) + }) + out <- list() + out$recf <- lev[[1]] + out$levrecf <- lev + out$info <- info + }else if(nn_type == "sntz"){ + yl0 <- lapply(bl, function(x){ + out <- as.matrix((x*(x>0))%*%t(R)) + return(out) + }) + lev <- lapply(yl0, function(x){ + y <- as.vector(t(Dh) %*% as.vector(t(x))) + y <- stats::setNames(y, paste0("k", rep(kset, h * (m/kset)), "h", + do.call("c", as.list(sapply( + (m/kset) * h, + function(x) seq(1:x) + ))))) + return(y) + }) + out <- list() + out$recf <- lev[[1]] + out$levrecf <- lev + }else{ + lev <- lapply(yl, function(x){ + y <- as.vector(t(Dh) %*% as.vector(t(x))) + y <- stats::setNames(y, paste0("k", rep(kset, h * (m/kset)), "h", + do.call("c", as.list(sapply( + (m/kset) * h, + function(x) seq(1:x) + ))))) + return(y) + }) + out <- list() + out$recf <- lev[[1]] + out$levrecf <- lev + } + return(out) + } +} + +lccrec_hts <- function(basef, C, nl, weights, bnaive = NULL, + const = "exogenous", CCC = TRUE, nn = FALSE, + nn_type = "osqp", settings = list(), ...){ + + if(is.vector(basef)){ + basef <- t(basef) + } + + h <- NROW(basef) + n <- NCOL(basef) + + if(is.null(C)){ + utd <- TRUE + C <- list(Matrix(1, ncol = n-1)) + nl <- 1 + }else{ + utd <- FALSE + } + + if(is.list(C)){ + nb <- NCOL(C[[1]]) + nl <- unlist(Map("NROW", C)) + na <- sum(nl) + + if(n != nb + na){ + stop("columns numb of basef != nb + na!", call. = FALSE) + } + Cl <- lapply(C, function(x) Matrix(x, sparse = TRUE)) + C <- do.call("rbind", Cl) + S <- rbind(C, Diagonal(nb)) + }else{ + nb <- NCOL(C) + na <- NROW(C) + + if(n != nb + na){ + stop("columns numb of basef != nb + na!", call. = FALSE) + } + + if(missing(nl)){ + warning("nl? There is only one level of upper time series?", call. = FALSE) + nl <- na + }else if(sum(nl) != na){ + stop("Please, provide a valid nl vector s.t. sum(nl) == na", call. = FALSE) + }else if(nl[1] != 1){ + stop("Please, provide a valid nl vector s.t. nl[1] == 1", call. = FALSE) + } + + + idr <- rep(1:length(nl), nl) + if(!is(C, "Matrix")){ + C <- Matrix(C, sparse = TRUE) + } + + Cl <- lapply(unique(idr), function(i) C[idr==i, , drop = FALSE]) + S <- rbind(C, Diagonal(nb)) + } + + check_bil <- any(sapply(Cl, function(x) any(colSums(x)!=1))) + if(check_bil){ + warning("The hierarchy is not balanced. \nThis could produce results that are inconsistent with the reconciliation procedure.", call. = FALSE) + } + + idc <- rep(1:(length(nl)+1), c(nl, nb)) + lal <- lapply(unique(idc), function(i) basef[,idc==i, drop = FALSE]) + b <- lal[[max(idc)]] + lal <- lal[-max(idc)] + + if(!is.null(bnaive)){ + if(is.vector(bnaive)){ + bnaive <- rbind(bnaive) + } + if(is.matrix(bnaive)){ + if(any(dim(b) != dim(bnaive))){ + bnaive <- t(bnaive) + } + + if(any(dim(b) != dim(bnaive))){ + stop("bnaive must be a matrix (h x nb)") + } + }else{ + stop("bnaive must be a matrix (h x nb)") + } + bm <- bnaive + }else{ + bm <- b + } + + if(is.vector(weights)){ + weights <- .sparseDiagonal(x = weights) + } + + nn_type <- match.arg(nn_type, c("osqp", "fbpp", "KAnn", "sntz")) + if(nn_type == "sntz"){ + nn <- FALSE + } + + if(const=="exogenous"){ + if(NCOL(weights) != nb | NROW(weights) != nb){ + stop("weights must be a nb x nb matrix", call. = FALSE) + }else if(!is(weights, "Matrix")){ + weights <- Matrix(weights, sparse = TRUE) + } + out <- Map(function(al, cl){ + recoLEV(al = al, Wb = weights, bl = bm, cl = cl, nn = nn, + settings = settings) + }, al = lal, cl = Cl) + + bl <- lapply(out, function(x){ + extract_data(x = x, name = "recf") + }) + + info <- lapply(out, function(x){ + extract_data(x = x, name = "info") + }) + names(info) <- paste0("id_l_", 1:length(bl)) + info <- info[!is.na(info)] + + uts0 <- do.call(c, lapply(out, function(x){ + extract_data(x = x, name = "uts0") + })) + if(any(uts0)){ + message("Some upper times series have base forecasts less then 0,", + " then they have been set equal to zero") + } + + }else{ + if(NCOL(weights) != n | NROW(weights) != n){ + stop("weights must be a n x n matrix", call. = FALSE) + }else if(!is(weights, "Matrix")){ + weights <- Matrix(weights, sparse = TRUE) + } + + id_nl <- rep(1:length(nl), nl) + Wl_id <- lapply(1:length(nl), function(x) c(which(id_nl == x), (na+1):(nb+na))) + out <- Map(function(al, cl, w){ + suppressMessages(htsrec(basef = cbind(al, bm), + C = cl, + comb = "w", + W = weights[w,w,drop = FALSE], + nn = nn, + settings = settings, + nn_type = nn_type, ...)) + }, al = lal, cl = Cl, w = Wl_id) + + bl <- lapply(out, function(x){ + out <- extract_data(x = x, name = "recf") + out <- out[,-(1:(NCOL(out)-nb)), drop = FALSE] + }) + + info <- lapply(out, function(x){ + extract_data(x = x, name = "info") + }) + names(info) <- paste0("id_l_", 1:length(bl)) + info <- info[!is.na(info)] + } + + if(!utd){ + if(nn){ + b[b<0] <- 0 + } + bl[[length(bl)+1]] <- b + } + + yl <- lapply(bl, function(x){ + out <- as.matrix(x%*%t(S)) + rownames(out) <- paste0("h", 1:h) + colnames(out) <- if(is.null(colnames(basef))) paste("serie", 1:n, sep = "") else colnames(basef) + return(out) + }) + + # return condition + if(length(yl)==1){ + if(length(info)>0){ + return(list(recf = yl[[1]], + info = info)) + }else if(nn_type == "sntz"){ + b_sntz <- bl[[1]]*(bl[[1]]>0) + y_sntz <- as.matrix(b_sntz%*%t(S)) + rownames(y_sntz) <- paste0("h", 1:h) + colnames(y_sntz) <- if(is.null(colnames(basef))) paste("serie", 1:n, sep = "") else colnames(basef) + return(y_sntz) + }else{ + return(yl[[1]]) + } + }else{ + names(yl) <- paste0("id_l_", 1:length(yl)) + if(CCC){ + out <- list() + if(nn_type == "sntz"){ + b_sntz <- Reduce("+", bl)/length(bl) + b_sntz <- b_sntz*(b_sntz>0) + out$recf <- as.matrix(b_sntz%*%t(S)) + rownames(out$recf) <- paste0("h", 1:h) + colnames(out$recf) <- if(is.null(colnames(basef))) paste("serie", 1:n, sep = "") else colnames(basef) + }else{ + out$recf <- Reduce("+", yl)/length(yl) + } + + out$levrecf <- yl + if(length(info)>0){ + out$info <- info + } + return(out) + }else{ + out <- list() + if(length(info)>0){ + out$recf <- yl[[1]] + out$levrecf <- yl + out$info <- info + }else if(nn_type == "sntz"){ + yl0 <- lapply(bl, function(x){ + out <- as.matrix((x*(x>0))%*%t(S)) + rownames(out) <- paste0("h", 1:h) + colnames(out) <- if(is.null(colnames(basef))) paste("serie", 1:n, sep = "") else colnames(basef) + return(out) + }) + out$recf <- yl0[[1]] + out$levrecf <- yl0 + }else{ + out$recf <- yl[[1]] + out$levrecf <- yl + } + return(out) + } + } +} + +lccrec_ctf <- function(basef, C, nl, m, weights, bnaive = NULL, + const = "exogenous", CCC = TRUE, nn = FALSE, + nn_type = "osqp", settings = list(), ...){ + + ctf <- suppressMessages(ctf_tools(C = C, m = m, h = 1, sparse = TRUE)) + S <- ctf$ctf$Stilde + + n <- ctf$hts$n + na <- ctf$hts$na + nb <- ctf$hts$nb + + kset <- ctf$thf$kset + m <- ctf$thf$m + p <- ctf$thf$p + kt <- ctf$thf$kt + ks <- ctf$thf$ks + + if (NROW(basef) != n) { + stop("Incorrect dimension of Ut or basef (they don't have same columns).", call. = FALSE) + } + + # Base forecast + if (NCOL(basef) %% kt != 0) { + stop("basef has a number of row not in line with the frequency of the series.", call. = FALSE) + } + + h <- NCOL(basef) / kt + Dh <- Dmat(h = h, m = kset, n = n) + Ybase <- matrix(Dh %*% as.vector(t(basef)), nrow = h, byrow = TRUE) + + if(missing(nl)){ + warning("nl? There is only one level of upper time series?", call. = FALSE) + nl <- na + }else if(sum(nl) != na){ + stop("Please, provide a valid nl vector s.t. sum(nl) == na", call. = FALSE) + }else if(nl[1] != 1){ + stop("Please, provide a valid nl vector s.t. nl[1] == 1", call. = FALSE) + } + + nk_nb <- c(nl, nb) + nLk <- kronecker(nk_nb, m/kset) + MLk <- matrix(1:length(nLk), ncol = length(kset), byrow = TRUE) + id <- as.vector(sapply(rep(split(MLk, 1:NROW(MLk)), nk_nb), + function(x) rep(x, m/kset))) + lal <- lapply(unique(id), function(i) Ybase[,id==i, drop = FALSE]) + b <- lal[[max(id)]] + lal <- lal[-max(id)] + + Cl <- lapply(1:(max(id)-1), function(x) S[id==x,, drop = FALSE]) + check_bil <- any(sapply(Cl, function(x) any(colSums(x)!=1))) + if(check_bil){ + warning("The hierarchy is not balanced. \nThis could produce results that are inconsistent with the reconciliation procedure.", call. = FALSE) + } + + names_lev <- t(sapply(unique(id), function(x) which(MLk == x, arr.ind = T))) + names_lev[,2] <- kset[names_lev[,2]] + + if(!is.null(bnaive)){ + if(is.vector(bnaive)){ + bnaive <- rbind(bnaive) + } + if(is.matrix(bnaive)){ + if(any(c(nb, m*h) != dim(bnaive))){ + bnaive <- t(bnaive) + } + + if(any(c(nb, m*h) != dim(bnaive))){ + stop("bnaive must be a matrix (nb x mh)", call. = FALSE) + } + }else{ + stop("bnaive must be a matrix (nb x mh)", call. = FALSE) + } + Db <- Dmat(h = h, m = m, n = nb) + bm <- matrix(Db%*%as.vector(t(bnaive)), nrow = h, byrow = TRUE) + }else{ + bm <- b + } + + if(is.vector(weights)){ + weights <- Diagonal(x = weights) + } + + nn_type <- match.arg(nn_type, c("osqp", "fbpp", "KAnn", "sntz")) + if(nn_type == "sntz"){ + nn <- FALSE + } + + if(const=="exogenous"){ + if(NCOL(weights) != nb*m | NROW(weights) != nb*m){ + stop("weights must be a nb*m x nb*m matrix", call. = FALSE) + }else if(!is(weights, "Matrix")){ + weights <- Matrix(weights, sparse = TRUE) + } + out <- Map(function(al, cl){ + recoLEV(al = al, Wb = weights, bl = bm, cl = cl, nn = nn, + settings = settings) + }, al = lal, cl = Cl) + + bl <- lapply(out, function(x){ + extract_data(x = x, name = "recf") + }) + + info <- lapply(out, function(x){ + extract_data(x = x, name = "info") + }) + names(info) <- paste0("id_l_", names_lev[-NROW(names_lev),1], "_k_", names_lev[-NROW(names_lev),2]) + info <- info[!is.na(info)] + + uts0 <- do.call(c, lapply(out, function(x){ + extract_data(x = x, name = "uts0") + })) + if(any(uts0)){ + message("Some times series (except for the high frequency bottom time ", + "series) have base forecasts less then 0,", + " then they have been set equal to zero") + } + }else{ + if(NCOL(weights) != kt*n | NROW(weights) != kt*n){ + stop("weights must be a kt*n x kt*n matrix", call. = FALSE) + }else if(!is(weights, "Matrix")){ + weights <- Matrix(weights, sparse = TRUE) + } + + Wl_id <- lapply(1:(max(id)-1), function(x) + c(which(id == x), which(id == max(id)))) + out <- Map(function(al, cl, w){ + suppressMessages(htsrec(basef = cbind(al, bm), + C = cl, + comb = "w", + W = weights[w,w,drop = FALSE], + nn = nn, + settings = settings, + nn_type = nn_type, ...)) + }, al = lal, cl = Cl, w = Wl_id) + + bl <- lapply(out, function(x){ + out <- extract_data(x = x, name = "recf") + out <- out[,-(1:(NCOL(out)-NCOL(b))), drop = FALSE] + }) + + info <- lapply(out, function(x){ + extract_data(x = x, name = "info") + }) + names(info) <- paste0("id_l_", names_lev[-NROW(names_lev),1], "_k_", names_lev[-NROW(names_lev),2]) + info <- info[!is.na(info)] + } + + if(nn){ + b[b<0] <- 0 + } + bl[[length(bl)+1]] <- b + + yl <- lapply(bl, function(x){ + out <- as.matrix(x%*%t(S)) + return(out) + }) + + names(yl) <- paste0("id_l_", names_lev[,1], "_k_", names_lev[,2]) + + if(CCC){ + out <- list() + if(nn_type == "sntz"){ + b_sntz <- Reduce("+", bl)/length(bl) + b_sntz <- b_sntz*(b_sntz>0) + out$recf <- as.matrix(b_sntz%*%t(S)) + }else{ + out$recf <- Reduce("+", yl)/length(yl) + } + out$recf <- v2m_oct(Dh = Dh, y = out$recf, n = n, + nam = rownames(basef), m = m, + kset = kset, h = h) + out$levrecf <- lapply(yl, v2m_oct, Dh = Dh, n = n, + nam = rownames(basef), m = m, + kset = kset, h = h) + if(length(info)>0){ + out$info <- info + } + return(out) + }else{ + out <- list() + if(length(info)>0){ + yl <- lapply(yl, v2m_oct, Dh = Dh, n = n, + nam = rownames(basef), m = m, + kset = kset, h = h) + out$recf <- yl[[1]] + out$levrecf <- yl + out$info <- info + }else if(nn_type == "sntz"){ + yl <- lapply(bl, function(x){ + out <- as.matrix((x*(x>0))%*%t(S)) + return(out) + }) + yl <- lapply(yl, v2m_oct, Dh = Dh, n = n, + nam = rownames(basef), m = m, + kset = kset, h = h) + out$recf <- yl[[1]] + out$levrecf <- yl + }else{ + yl <- lapply(yl, v2m_oct, Dh = Dh, n = n, + nam = rownames(basef), m = m, + kset = kset, h = h) + out$recf <- yl[[1]] + out$levrecf <- yl + } + return(out) + } +} + +v2m_oct <- function(y, Dh, n, nam, m, kset, h){ + out <- matrix(t(Dh) %*% as.vector(t(y)), nrow = n, byrow = TRUE) + rownames(out) <- if (is.null(nam)) paste("serie", 1:n, sep = "") else nam + colnames(out) <- paste("k", rep(kset, h * (m/kset)), "h", + do.call("c", as.list(sapply( + (m/kset) * h, + function(x) seq(1:x) + ))), + sep = "" + ) + return(out) +} + +# Generic lccrec +recoLEV <- function(al, Wb, bl, cl, nn = FALSE, settings) { + out <- list() + rh <- cl%*%Wb%*%t(cl) + lh <- t(al) - cl%*%t(bl) + out$recf <- t(bl) + Wb%*%t(cl)%*%solveLin(rh, lh, verb = FALSE) + out$recf <- t(out$recf) + out$uts0 <- FALSE + if(nn & any(out$recf < (-sqrt(.Machine$double.eps)))){ + + if(isDiagonal(Wb)){ + P <- .sparseDiagonal(x = diag(Wb)^(-1)) + }else{ + R <- chol(Wb) + P <- chol2inv(R) + } + id <- which(rowSums(out$recf<(-sqrt(.Machine$double.eps)))!=0) + + if(any(al[id,]<0)){ + out$uts0 <- TRUE + al[id,] <- al[id,]*(al[id,]>0) + } + + rec <- Map(function(ah, bh){ + lev_osqp(a = ah, b = bh, P = P, C = cl, settings = settings) + }, ah = split(al[id,,drop = FALSE], id), bh = split(bl[id,,drop = FALSE], id)) + + rec <- do.call("rbind",rec) + out$recf[id,] <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", + "status", "status_polish") + rownames(out$info) <- id + + out$recf[-id,,drop=FALSE] <- out$recf[-id,,drop=FALSE] * (out$recf[-id,,drop=FALSE] > 0) + }else if(nn){ + out$recf <- out$recf * (out$recf > 0) + } + + return(out) +} + +# lccrec osqp +lev_osqp <- function(a, b, P, C, settings) { + nb <- NCOL(C) + + l <- c(a, rep(0, nb)) + u <- c(a, rep(Inf, nb)) + A <- rbind(C, .sparseDiagonal(nb)) + + q <- (-1) * t(P) %*% as.vector(b) + + if(length(settings)==0){ + settings = osqpSettings(verbose = FALSE, + eps_abs = 1e-5, + eps_rel = 1e-5, + polish_refine_iter = 100, + polish=TRUE) + } + + rec <- solve_osqp(P, q, A, l, u, settings) + + out <- list() + out$recf <- rec$x + + if(rec$info$status_val != 1){ + warning(paste("OSQP flag", rec$info$status_val, "OSQP pri_res", rec$info$pri_res), call. = FALSE) + } + + out$recf[which(out$recf < 0)] <- 0 + + out$info <- c(rec$info$obj_val, rec$info$run_time, rec$info$iter, rec$info$pri_res, + rec$info$status_val, rec$info$status_polish) + + return(out) +} diff --git a/R/oct_bounds.R b/R/oct_bounds.R index 654c625..9b5dc99 100644 --- a/R/oct_bounds.R +++ b/R/oct_bounds.R @@ -1,75 +1,110 @@ -# Title -# -# @param hts_bounds -# @param thf_bounds -# @param m -# @param C -# @param Ut -# @param nb -# -# @return -# @export -# -# @examples -# hts_bounds <- matrix(c(rnorm(3), rnorm(3, 100)), nrow = 3) -# hts_bounds[1,2] <- Inf -# hts_bounds[2,1] <- -Inf -# thf_bounds <- matrix(c(rnorm(7), rnorm(7, 100)), nrow = 7) -# obj <- oct_bounds(hts_bounds = hts_bounds, m = c(4,1), C = matrix(1, ncol = 2)) -oct_bounds <- function(hts_bounds, thf_bounds, m, C, Ut, nb){ +#' Optimal cross-temporal bounds +#' +#' @description +#' \loadmathjax +#' Function to export the constraints designed for the cross-sectional and/or +#' temporal reconciled forecasts +#' +#' @param hts_bounds (\mjseqn{n \times 2}) matrix with cross-sectional bounds: +#' the first column is the lower bound, and the second column is the upper bound. +#' @param thf_bounds (\mjseqn{(k^\ast + m) \times 2}) matrix with temporal bounds: +#' the first column is the lower bound, and the second column is the upper bound. +#' @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 C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones. +#' @param Ut Zero constraints cross-sectional (contemporaneous) kernel matrix +#' \mjseqn{(\textbf{U}'\textbf{y} = \mathbf{0})} spanning the null space valid +#' for the reconciled forecasts. It can be used instead of parameter +#' \code{C}, but \code{nb} (\mjseqn{n = n_a + n_b}) is needed if +#' \mjseqn{\textbf{U}' \neq [\textbf{I} \ -\textbf{C}]}{}. If the hierarchy +#' admits a structural representation, \mjseqn{\textbf{U}'} has dimension +#' (\mjseqn{n_a \times n}). +#' +#' @return A matrix with the cross-temporal bounds. +#' +#' @keywords utilities +#' @family utilities +#' +#' @examples +#' data(FoReco_data) +#' # monthly base forecasts +#' mbase <- t(FoReco_data$base[, which(simplify2array(strsplit( +#' colnames(FoReco_data$base), split = "_"))[1, ] == "k1")]) +#' #' # monthly residuals +#' mres <- t(FoReco_data$res[, which(simplify2array(strsplit( +#' colnames(FoReco_data$res),split = "_"))[1, ] == "k1")]) +#' +#' # For example, in FoReco_data we want that BA > 78, and C > 50 +#' cs_bound <- matrix(c(rep(-Inf, 5), 78, -Inf, 50, rep(+Inf, 8)), ncol = 2) +#' ## Cross-sectional reconciliation +#' csobj <- htsrec(mbase, C = FoReco_data$C, comb = "shr", res = mres, bounds = cs_bound) +#' +#' # Extension of the constraints to the cross-temporal case +#' ct_bound <- oct_bounds(hts_bounds = cs_bound, m = 12) +#' ## Cross-temporal reconciliation +#' obj <- octrec(FoReco_data$base, m = 12, C = FoReco_data$C, comb = "bdshr", +#' res = FoReco_data$res, bounds = ct_bound) +#' +#' @export +oct_bounds <- function(hts_bounds, thf_bounds, m, C, Ut){ if(missing(thf_bounds) & missing(hts_bounds)){ stop("ciao") - }else if(missing(thf_bounds)){ - hts_nrow <- NROW(hts_bounds) - thf_nrow <- thf_tools(m = m)$kt + }else if(missing(thf_bounds) | missing(hts_bounds)){ + print("c") + if(missing(thf_bounds)){ + hts_nrow <- NROW(hts_bounds) + thf_nrow <- thf_tools(m = m)$kt - if(!missing(C)){ - n <- hts_tools(C = C)$n - if(hts_nrow != n){ - stop("hts_bounds must be a matrix (", n, " x 2)", call. = FALSE) + if(!missing(C)){ + n <- NROW(C)+NCOL(C) + if(hts_nrow != n){ + stop("hts_bounds must be a matrix (", n, " x 2)", call. = FALSE) + } + }else if(!missing(Ut)){ + n <- NCOL(Ut) + if(hts_nrow != n){ + stop("hts_bounds must be a matrix (", n, " x 2)", call. = FALSE) + } } - }else if(!(missing(Ut) & missing(nb))){ - n <- hts_tools(Ut = Ut, nb = nb)$n - if(hts_nrow != n){ - stop("hts_bounds must be a matrix (", n, " x 2)", call. = FALSE) - } - } - bhts <- t(simplify2array(rep(split(hts_bounds, 1:hts_nrow), each = thf_nrow))) - dimnames(bhts) <- NULL - colnames(bhts) <- c("lower", "upper") - return(bhts) - }else if(missing(hts_bounds)){ - if(!missing(C)){ - hts_nrow <- hts_tools(C = C)$n - }else if(!(missing(Ut) & missing(nb))){ - hts_nrow <- hts_tools(Ut = Ut, nb = nb)$n - }else{ - stop("the argument C (or Ut AND nb) is not specified", call. = FALSE) - } + bhts <- t(simplify2array(rep(split(hts_bounds, 1:hts_nrow), each = thf_nrow))) + dimnames(bhts) <- NULL + colnames(bhts) <- c("lower", "upper") + return(bhts) + }else if(missing(hts_bounds)){ + if(!missing(C)){ + hts_nrow <- NROW(C)+NCOL(C) + }else if(!missing(Ut)){ + hts_nrow <- NCOL(Ut) + }else{ + stop("the argument C or Ut is not specified", call. = FALSE) + } - thf_nrow <- NROW(thf_bounds) + thf_nrow <- NROW(thf_bounds) - if(!missing(m)){ - kt <- thf_tools(m = m)$kt - if(thf_nrow != kt){ - stop("thf_bounds must be a matrix (", kt, " x 2)", call. = FALSE) + if(!missing(m)){ + kt <- thf_tools(m = m)$kt + if(thf_nrow != kt){ + stop("thf_bounds must be a matrix (", kt, " x 2)", call. = FALSE) + } } - } - bthf <- t(simplify2array(rep(split(thf_bounds, 1:thf_nrow), hts_nrow))) - dimnames(bthf) <- NULL - colnames(bthf) <- c("lower", "upper") - return(bthf) + bthf <- t(simplify2array(rep(split(thf_bounds, 1:thf_nrow), hts_nrow))) + dimnames(bthf) <- NULL + colnames(bthf) <- c("lower", "upper") + return(bthf) + } }else{ kt <- thf_tools(m = m)$kt if(!missing(C)){ - n <- hts_tools(C = C)$n - }else if(!(missing(Ut) & missing(nb))){ - n <- hts_tools(Ut = Ut, nb = nb)$n + n <- NROW(C)+NCOL(C) + }else if(!missing(Ut)){ + n <- NCOL(Ut) }else{ - stop("the argument C (or Ut AND nb) is not specified", call. = FALSE) + stop("the argument C or Ut is not specified", call. = FALSE) } thf_nrow <- NROW(thf_bounds) hts_nrow <- NROW(hts_bounds) diff --git a/R/octrec.R b/R/octrec.R index d133e90..c0af4cc 100755 --- a/R/octrec.R +++ b/R/octrec.R @@ -1,29 +1,35 @@ #' @title Optimal combination cross-temporal forecast reconciliation #' #' @description -#' Optimal (in least squares sense) combination cross-temporal forecast reconciliation. -#' The reconciled forecasts are calculated either through a projection approach -#' (Byron, 1978), or the equivalent structural approach by Hyndman et al. (2011). +#' \loadmathjax +#' Optimal (in least squares sense) combination cross-temporal forecast +#' reconciliation. The reconciled forecasts are calculated either through a +#' projection approach (Byron, 1978), or the equivalent structural approach +#' by Hyndman et al. (2011). #' -#' @usage octrec(basef, m, C, comb, res, Ut, nb, Sstruc, -#' mse = TRUE, corpcor = FALSE, type = "M", sol = "direct", -#' nn = FALSE, keep = "list", settings = osqpSettings(), +#' @usage octrec(basef, m, C, comb, res, Ut, nb, mse = TRUE, +#' corpcor = FALSE, type = "M", sol = "direct", keep = "list", +#' nn = FALSE, nn_type = "osqp", settings = list(), #' bounds = NULL, W = NULL, Omega = NULL) #' -#' @param basef (\code{n x h(k* + m)}) matrix of base forecasts to be reconciled; -#' \code{n} is the total number of variables, \code{m} is the highest frequency, -#' \code{k*} is the sum of (\code{p-1}) factors of \code{m}, excluding \code{m}, -#' and \code{h} is the forecast horizon. Each row identifies, a time series, and the forecasts -#' are ordered as [lowest_freq' ... highest_freq']'. -#' @param m Highest available sampling frequency per seasonal cycle (max. order of temporal aggregation, \code{m}), -#' or a subset of the \code{p} factors of \code{m}. -#' @param comb Type of the reconciliation, it corrispond to a different covariance -#' matrix (\code{n(k* + m) x n(k* + m)}), \code{k*} is the sum of (\code{p-1}) -#' factors of \code{m} (exclude \code{m} as factors) and \code{n} is the number +#' @param basef (\mjseqn{n \times h(k^\ast+m)}) matrix of base forecasts to be +#' reconciled, \mjseqn{\widehat{\mathbf{Y}}}; \mjseqn{n} is the total number of variables, +#' \mjseqn{m} is the highest time frequency, \mjseqn{k^\ast} is the sum of (a +#' subset of) (\mjseqn{p-1}) factors of \mjseqn{m}, excluding \mjseqn{m}, and +#' \mjseqn{h} is the forecast horizon for the lowest frequency time series. +#' Each row identifies a time series, and the forecasts are ordered as +#' [lowest_freq' ... highest_freq']'. +#' @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 comb Type of the reconciliation. It corresponds to a specific +#' (\mjseqn{n(k\ast + m) \times n(k^\ast + m)}) covariance matrix, where +#' \mjseqn{k^\ast} is the sum of (a subset of) (\mjseqn{p-1}) factors of +#' \mjseqn{m} (\mjseqn{m} is not considered) and \mjseqn{n} is the number #' of variables: #' \itemize{ #' \item \bold{ols} (Identity); -#' \item \bold{struc} (Cross-temporal summing matrix, use the \code{Sstruc} param to reduce computation time); +#' \item \bold{struc} (Cross-temporal summing matrix); #' \item \bold{wlsh} (Hierarchy variances matrix); #' \item \bold{wlsv} (Series variances matrix); #' \item \bold{bdshr} (Shrunk cross-covariance matrix, cross-sectional framework); @@ -36,75 +42,189 @@ #' \item \bold{w} use your personal matrix W in param \code{W}; #' \item \bold{omega} use your personal matrix Omega in param \code{Omega}. #' } -#' @param C (\code{na x nb}) cross-sectional (contemporaneous) matrix mapping the bottom -#' level series into the higher level ones. +#' @param C (\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones. #' @param Ut Zero constraints cross-sectional (contemporaneous) kernel matrix -#' \eqn{(\textbf{U}'\textbf{Y} = \mathbf{0})}{} spanning the null space valid for the reconciled -#' forecasts. It can be used instead of parameter \code{C}, but in this case \code{nb} (n = na + nb) is needed. If -#' the hierarchy admits a structural representation, \code{Ut} has dimension (\code{na x n}). -#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} is not used. -#' @param res (\code{n x N(k* + m)}) matrix containing the residuals at all the -#' temporal frequencies ordered [lowest_freq' ... highest_freq']' (columns) for -#' each variable (row), needed to estimate the covariance matrix when \code{comb =} -#' \code{\{"sam",} \code{"wlsv",} \code{"wlsh",} \code{"acov",} \code{"Ssam",} -#' \code{"Sshr",} \code{"Sshr1",} \code{"shr"\}}. +#' \mjseqn{(\textbf{U}'\textbf{y} = \mathbf{0})} spanning the null space valid +#' for the reconciled forecasts. It can be used instead of parameter +#' \code{C}, but \code{nb} (\mjseqn{n = n_a + n_b}) is needed if +#' \mjseqn{\textbf{U}' \neq [\textbf{I} \ -\textbf{C}]}. If the hierarchy +#' admits a structural representation, \mjseqn{\textbf{U}'} has dimension +#' (\mjseqn{n_a \times n}). +#' @param nb Number of bottom time series; if \code{C} is present, \code{nb} +#' and \code{Ut} are not used. +#' @param res (\mjseqn{n \times N(k^\ast + m)}) matrix containing the residuals at +#' all the temporal frequencies ordered [lowest_freq' ... highest_freq']' +#' (columns) for each variable (row), needed to estimate the covariance matrix +#' when \code{comb =} \code{\{"sam",} \code{"wlsv",} \code{"wlsh",} +#' \code{"acov",} \code{"Ssam",} \code{"Sshr",} \code{"Sshr1",} \code{"shr"\}}. #' @param W,Omega This option permits to directly enter the covariance matrix: #' \enumerate{ -#' \item \code{W} must be a p.d. (\code{n(k* + m) x n(k* + m)}) matrix or a list -#' of \code{h} matrix (one for each forecast horizon); -#' \item \code{Omega} must be a p.d. (\code{n(k* + m) x n(k* + m)}) matrix or a list -#' of \code{h} matrix (one for each forecast horizon); -#' \item if \code{comb} is different from "\code{w}" or "\code{omega}", \code{W} or \code{Omega} is not used. +#' \item \code{W} must be a p.d. (\mjseqn{n(k^\ast + m) \times n(k^\ast + m)}) +#' matrix or a list of \mjseqn{h} matrix (one for each forecast horizon); +#' \item \code{Omega} must be a p.d. (\mjseqn{n(k^\ast + m) \times n(k^\ast + m)}) +#' matrix or a list of \code{h} matrix (one for each forecast horizon); +#' \item if \code{comb} is different from "\code{w}" or "\code{omega}", +#' \code{W} or \code{Omega} is not used. #' } -#' @param Sstruc Cross-temporal summing matrix (structural representation)\eqn{,\; \check{\textbf{S}}}{}; -#' can be obtained through the function \link[FoReco]{ctf_tools}. #' @param mse Logical value: \code{TRUE} (\emph{default}) calculates the -#' covariance matrix of the in-sample residuals (when necessary) according to the original -#' \pkg{hts} and \pkg{thief} formulation: no mean correction, T as denominator. +#' covariance matrix of the in-sample residuals (when necessary) according to +#' the original \pkg{hts} and \pkg{thief} formulation: no mean correction, +#' T as denominator. #' @param corpcor Logical value: \code{TRUE} if \pkg{corpcor} (\enc{Schäfer}{Schafer} et #' al., 2017) must be used to shrink the sample covariance matrix according to -#' \enc{Schäfer}{Schafer} and Strimmer (2005), otherwise the function uses the same -#' implementation as package \pkg{hts}. +#' \enc{Schäfer}{Schafer} and Strimmer (2005), otherwise the function uses the +#' same implementation as package \pkg{hts}. #' @param type Approach used to compute the reconciled forecasts: \code{"M"} for #' the projection approach with matrix M (\emph{default}), or \code{"S"} for the #' structural approach with summing matrix S. -#' @param keep Return a list object of the reconciled forecasts at all levels. -#' @param sol Solution technique for the reconciliation problem: either \code{"direct"} (\emph{default}) for the direct -#' solution or \code{"osqp"} for the numerical solution (solving a linearly constrained quadratic -#' program using \code{\link[osqp]{solve_osqp}}). -#' @param nn Logical value: \code{TRUE} if non-negative reconciled forecasts are wished. -#' @param settings Settings for \pkg{osqp} (object \code{\link[osqp]{osqpSettings}}). The default options -#' are: \code{verbose = FALSE}, \code{eps_abs = 1e-5}, \code{eps_rel = 1e-5}, -#' \code{polish_refine_iter = 100} and \code{polish = TRUE}. For details, see the -#' \href{https://osqp.org/}{\pkg{osqp} documentation} (Stellato et al., 2019). -#' @param bounds (\code{n(k* + m) x 2}) matrix with bounds on the variables: the first column is the lower bound, -#' and the second column is the upper bound. +#' @param keep Return a list object of the reconciled forecasts at all levels +#' (if keep = "list") or only the reconciled forecasts matrix (if keep = "recf"). +#' @param sol Solution technique for the reconciliation problem: either +#' \code{"direct"} (\emph{default}) for the closed-form matrix solution, or +#' \code{"osqp"} for the numerical solution (solving a linearly constrained +#' quadratic program using \code{\link[osqp]{solve_osqp}}). +#' @param nn Logical value: \code{TRUE} if non-negative reconciled forecasts +#' are wished. +#' @param nn_type "osqp" (default), "KAnn" (only \code{type == "M"}) or "sntz". +#' @param settings Settings for \pkg{osqp} (object \code{\link[osqp]{osqpSettings}}). +#' The default options are: \code{verbose = FALSE}, \code{eps_abs = 1e-5}, +#' \code{eps_rel = 1e-5}, \code{polish_refine_iter = 100} and \code{polish = TRUE}. +#' For details, see the \href{https://osqp.org/}{\pkg{osqp} documentation} +#' (Stellato et al., 2019). +#' @param bounds (\mjseqn{n(k^\ast + m) \times 2}) matrix of the bounds on the +#' variables: the first column is the lower bound, and the second column is the +#' upper bound. #' #' @details -#' In case of non-negativity constraints, there are two ways: -#' \enumerate{ -#' \item \code{sol = "direct"} and \code{nn = TRUE}: the base forecasts -#' will be reconciled at first without non-negativity constraints, then, if negative reconciled -#' values are present, the \code{"osqp"} solver is used. -#' \item \code{sol = "osqp"} and \code{nn = TRUE}: the base forecasts will be -#' reconciled through the \code{"osqp"} solver. +#' \loadmathjax +#' Considering contemporaneous and temporal dimensions in the +#' same framework requires to extend and adapt the notations +#' used in \link[FoReco]{htsrec} and \link[FoReco]{thfrec}. +#' To do that, we define the matrix containing the base forecasts +#' at any considered temporal frequency as +#' \mjsdeqn{ +#' \widehat{\textbf{Y}}_{n \times h(k^\ast+m))} = +#' \left[\begin{array}{ccccc} +#' \widehat{\textbf{A}}^{[m]} & \widehat{\textbf{A}}^{[k_{p-1}]} & \cdots & +#' \widehat{\textbf{A}}^{[k_2]} & \widehat{\textbf{A}}^{[1]} \cr +#' \widehat{\textbf{B}}^{[m]} & \widehat{\textbf{B}}^{[k_{p-1}]} & \cdots & +#' \widehat{\textbf{B}}^{[k_2]} & \widehat{\textbf{B}}^{[1]} +#' \end{array} +#' \right] \qquad k \in {\cal K},} +#' where \mjseqn{\cal K} is a subset of \mjseqn{p} factors of \mjseqn{m} and, +#' \mjseqn{\widehat{\textbf{B}}^{[k]}} and \mjseqn{\widehat{\textbf{A}}^{[k]}} +#' are the matrices containing the \mjseqn{k}-order temporal aggregates of the +#' bts and uts, of dimension (\mjseqn{n_b \times h m/k}) and +#' (\mjseqn{n_a \times h m/k}), respectively. +#' +#' Let us consider the multivariate regression model +#' \mjsdeqn{\widehat{\mathbf{Y}} = \mathbf{Y} + \mathbf{E} ,} +#' where the involved matrices have each dimension +#' \mjseqn{[n \times (k^\ast+m)]} and contain, respectively, the base +#' (\mjseqn{\widehat{\mathbf{Y}}}) and the target forecasts +#' (\mjseqn{\mathbf{Y}}), and the coherency errors (\mjseqn{\mathbf{E}}) for +#' the \mjseqn{n} component variables of the linearly constrained time series +#' of interest. For each variable, \mjseqn{k^\ast + m} base forecasts are +#' available, pertaining to all aggregation levels of the temporal hierarchy +#' for a complete cycle of high-frequency observation, \mjseqn{m}. Consider +#' now two vectorized versions of model, by transforming the matrices either +#' in original form: +#' \mjsdeqn{\mbox{vec}\left(\widehat{\mathbf{Y}}\right) = +#' \mbox{vec}\left(\mathbf{Y}\right) + \mathbf{\varepsilon} \; \mbox{ with } +#' \; \mathbf{\varepsilon} = \mbox{vec}\left(\mathbf{E}\right)} +#' or in transposed form: +#' \mjsdeqn{\mbox{vec}\left(\widehat{\mathbf{Y}}'\right) = +#' \mbox{vec}\left(\mathbf{Y}'\right) + \mathbf{\eta} \; \mbox{ with } +#' \; \mathbf{\eta} = \mbox{vec}\left(\mathbf{E}'\right).} +#' Denote with \mjseqn{\mathbf{P}} the \mjseqn{[n(k^\ast+m) \times n(k^\ast+m)]} +#' commutation matrix such that +#' \mjseqn{\mathbf{P}\mbox{vec}(\mathbf{Y}) = \mbox{vec}(\mathbf{Y}')}, +#' \mjseqn{\mathbf{P}\mbox{vec}(\widehat{\mathbf{Y}}) = \mbox{vec}(\widehat{\mathbf{Y}}')} +#' and \mjseqn{\mathbf{P}\mathbf{\varepsilon} = {\bf \eta}}. +#' Let \mjseqn{\mathbf{W} = \mathrm{E}[\mathbf{\varepsilon\varepsilon}']} be the +#' covariance matrix of vector \mjseqn{\mathbf{\varepsilon}}, and +#' \mjseqn{\mathbf{\Omega} = \mathrm{E}[\mathbf{\eta\eta}']} the covariance matrix of +#' vector \mjseqn{\mathbf{\eta}}. Clearly, \mjseqn{\mathbf{W}} and +#' \mjseqn{\mathbf{\Omega}} are different parameterizations of the same +#' statistical object for which the following relationships hold: +#' \mjsdeqn{\mathbf{\Omega} = \mathbf{P}\mathbf{W}\mathbf{P}', +#' \qquad \mathbf{W} = \mathbf{P}' \mathbf{\Omega}\mathbf{P} .} +#' In order to apply the general point forecast reconciliation according to the +#' projection approach (\code{type = "M"}) to a cross-temporal forecast +#' reconciliation problem, we may consider either two \emph{vec}-forms , e.g. +#' if we follow the first: +#' \mjsdeqn{ +#' \tilde{\mathbf{y}}= \hat{\mathbf{y}} - \mathbf{\Omega}\mathbf{H}\left( +#' \mathbf{H}'\mathbf{\Omega}\mathbf{H}\right)^{-1}\mathbf{H}'\hat{\mathbf{y}} = +#' {\mathbf{M}}\hat{\mathbf{y}},} +#' where \mjseqn{\widehat{\mathbf{y}} = \mbox{vec}(\widehat{\mathbf{Y}}')} is the +#' row vectorization of the base forecasts matrix \mjseqn{\widehat{\mathbf{Y}}} +#' The alternative equivalent solution (\code{type = "S"}) (following the +#' structural reconciliation approach by Hyndman et al., 2011) is +#' \mjsdeqn{\widetilde{\mathbf{y}} = \widetilde{\mathbf{S}}\left(\widetilde{\mathbf{S}}' +#' \mathbf{\Omega}^{-1}\widetilde{\mathbf{S}}\right)^{-1}\widetilde{\mathbf{S}}' +#' \mathbf{\Omega}^{-1}\widehat{\mathbf{y}} = \widetilde{\mathbf{S}}\mathbf{G}\widehat{\mathbf{y}}.} +#' where \mjseqn{\widetilde{\mathbf{S}}} is the cross-temporal summing matrix. +#' +#' \strong{Bounds on the reconciled forecasts} +#' +#' When the reconciliation uses the optimization package osqp, +#' the user may impose bounds on the reconciled forecasts. +#' The parameter \code{bounds} permits to consider lower (\mjseqn{\mathbf{a}}) and +#' upper (\mjseqn{\mathbf{b}}) bounds like \mjseqn{\mathbf{a} \leq +#' \widetilde{\mathbf{y}} \leq \mathbf{b}}, where \mjseqn{\widehat{\mathbf{y}} = +#' \mbox{vec}(\widehat{\mathbf{Y}}')}, such that: +#' \mjsdeqn{ \begin{array}{c} +#' a_1 \leq \widetilde{y}_1 \leq b_1 \cr +#' \dots \cr +#' a_{n(k^\ast + m)} \leq \widetilde{y}_{n(k^\ast + m)} \leq b_{n(k^\ast + m)} \cr +#' \end{array} \Rightarrow +#' \mbox{bounds} = [\mathbf{a} \; \mathbf{b}] = +#' \left[\begin{array}{cc} +#' a_1 & b_1 \cr +#' \vdots & \vdots \cr +#' a_{n(k^\ast + m)} & b_{n(k^\ast + m)} \cr +#' \end{array}\right],} +#' where \mjseqn{a_i \in [- \infty, + \infty]} and \mjseqn{b_i \in [- \infty, + \infty]}. +#' If \mjseqn{y_i} is unbounded, the i-th row of \code{bounds} would be equal +#' to \code{c(-Inf, +Inf)}. +#' Notice that if the \code{bounds} parameter is used, \code{sol = "osqp"} must be used. +#' This is not true in the case of non-negativity constraints: +#' \itemize{ +#' \item \code{sol = "direct"}: first the base forecasts +#' are reconciled without non-negativity constraints, then, if negative reconciled +#' values are present, the \code{"osqp"} solver is used; +#' \item \code{sol = "osqp"}: the base forecasts are +#' reconciled using the \code{"osqp"} solver. +#' } +#' In this case it is not necessary to build a matrix containing +#' the bounds, and it is sufficient to set \code{nn = "TRUE"}. +#' +#' Non-negative reconciled forecasts may be obtained by setting \code{nn_type} alternatively as: +#' \itemize{ +#' \item \code{nn_type = "KAnn"} (Kourentzes and Athanasopoulos, 2021) +#' \item \code{nn_type = "sntz"} ("set-negative-to-zero") +#' \item \code{nn_type = "osqp"} (Stellato et al., 2020) #' } #' #' @return #' If the parameter \code{keep} is equal to \code{"recf"}, then the function -#' returns only the (\code{n x h(k* + m)}) reconciled forecasts matrix, otherwise (\code{keep="all"}) -#' it returns a list that mainly depends on what type of representation (\code{type}) -#' and methodology (\code{sol}) have been used: -#' \item{\code{recf}}{(\code{n x h(k* + m)}) reconciled forecasts matrix.} -#' \item{\code{Omega}}{Covariance matrix used for reconciled forecasts (\code{vec}(\code{Y}\enc{'}{'}) representation).} -#' \item{\code{W}}{Covariance matrix used for reconciled forecasts (\code{vec(Y)} representation).} +#' returns only the (\mjseqn{n \times h(k^\ast + m)}) reconciled forecasts +#' matrix, otherwise (\code{keep="all"}) it returns a list that mainly depends +#' on what type of representation (\code{type}) and solution technique +#' (\code{sol}) have been used: +#' \item{\code{recf}}{(\mjseqn{n \times h(k^\ast + m)}) reconciled forecasts matrix, \mjseqn{\widetilde{\textbf{Y}}}.} +#' \item{\code{Omega}}{Covariance matrix used for reconciled forecasts (\mjseqn{\mbox{vec}(\widehat{\textbf{Y}}')} representation).} +#' \item{\code{W}}{Covariance matrix used for reconciled forecasts (\mjseqn{\mbox{vec}(\widehat{\textbf{Y}})} representation).} #' \item{\code{nn_check}}{Number of negative values (if zero, there are no values below zero).} -#' \item{\code{rec_check}}{Logical value: has the hierarchy been respected?} -#' \item{\code{M} (\code{type="M"} and \code{type="direct"})}{Projection matrix (projection approach).} -#' \item{\code{G} (\code{type="S"} and \code{type="direct"})}{Projection matrix (structural approach).} -#' \item{\code{S} (\code{type="S"} and \code{type="direct"})}{Cross-temporal summing matrix\eqn{, \; \textbf{Q}\check{\textbf{S}}}{} (\code{vec}(\code{Y}\enc{'}{'}) representation).} +#' \item{\code{rec_check}}{Logical value: \code{rec_check = TRUE} when the constraints have been fulfilled,} +#' \item{\code{varf} (\code{type="direct"})}{(\mjseqn{n \times (k^\ast + m)}) reconciled forecasts variance matrix for \mjseqn{h=1}, \mjseqn{\mbox{diag}(\mathbf{MW}}).} +#' \item{\code{M} (\code{type="direct"})}{Projection matrix (projection approach).} +#' \item{\code{G} (\code{type="S"} and \code{type="direct"})}{Projection matrix (structural approach, \mjseqn{\mathbf{M}=\mathbf{S}\mathbf{G}}).} +#' \item{\code{S} (\code{type="S"} and \code{type="direct"})}{Cross-temporal summing matrix (\mjseqn{\widetilde{\textbf{S}}\mbox{vec}(\widehat{\textbf{Y}}')} representation).} #' \item{\code{info} (\code{type="osqp"})}{matrix with some useful indicators (columns) -#' for each forecast horizon \code{h} (rows): run time (\code{run_time}) number of iteration, +#' for each forecast horizon \mjseqn{h} (rows): run time (\code{run_time}), number of iteration, #' norm of primal residual (\code{pri_res}), status of osqp's solution (\code{status}) and #' polish's status (\code{status_polish}).} #' @@ -124,14 +244,15 @@ #' Matrix Estimation and Implications for Functional Genomics, \emph{Statistical #' Applications in Genetics and Molecular Biology}, 4, 1. #' -#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2018). OSQP: -#' An Operator Splitting Solver for Quadratic Programs, \href{https://arxiv.org/abs/1711.08013}{arXiv:1711.08013}. +#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2020). OSQP: +#' An Operator Splitting Solver for Quadratic Programs, \emph{Mathematical Programming Computation}, +#' 12, 4, 637-672. #' #' Stellato, B., Banjac, G., Goulart, P., Boyd, S., Anderson, E. (2019), OSQP: -#' Quadratic Programming Solver using the 'OSQP' Library, R package version 0.6.0.3 +#' Quadratic Programming Solver using the `OSQP' Library, R package version 0.6.0.3 #' (October 10, 2019), \href{https://CRAN.R-project.org/package=osqp}{https://CRAN.R-project.org/package=osqp}. #' -#' @keywords reconciliation +#' @family reconciliation procedures #' @examples #' data(FoReco_data) #' obj <- octrec(FoReco_data$base, m = 12, C = FoReco_data$C, @@ -141,9 +262,9 @@ #' #' @import Matrix osqp methods #' -octrec <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, - corpcor = FALSE, type = "M", sol = "direct", nn = FALSE, keep = "list", - settings = osqpSettings(), bounds = NULL, W = NULL, Omega = NULL) { +octrec <- function(basef, m, C, comb, res, Ut, nb, mse = TRUE, corpcor = FALSE, + type = "M", sol = "direct", keep = "list", nn = FALSE, + nn_type = "osqp", settings = list(), bounds = NULL, W = NULL, Omega = NULL) { if(missing(comb)){ stop("The argument comb is not specified.", call. = FALSE) @@ -160,9 +281,11 @@ octrec <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, } #' @export -octrec.default <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, - corpcor = FALSE, type = "M", sol = "direct", nn = FALSE, keep = "list", - settings = osqpSettings(), bounds = NULL, W = NULL, Omega = NULL) { +octrec.default <- function(basef, m, C, comb, res, Ut, nb, mse = TRUE, + corpcor = FALSE, type = "M", sol = "direct", + keep = "list", nn = FALSE, nn_type = "osqp", + settings = osqpSettings(), bounds = NULL, + W = NULL, Omega = NULL) { if (missing(m)) { stop("The argument m is not specified", call. = FALSE) } @@ -174,79 +297,74 @@ octrec.default <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, type <- match.arg(type, c("M", "S")) keep <- match.arg(keep, c("list", "recf")) + nn_type <- match.arg(nn_type, c("osqp", "KAnn", "fbpp", "sntz")) + if (missing(basef)) { stop("The argument basef is not specified", call. = FALSE) } - tools <- thf_tools(m, sparse = TRUE) - kset <- tools$kset - m <- max(kset) - p <- tools$p - kt <- tools$kt - ks <- tools$ks - - # matrix - Zt <- tools$Zt - - n <- NROW(basef) - if (missing(C)) { # Using Ut AND nb or C - if (missing(Ut) | missing(nb)) { - stop("Please, give C (or Ut AND nb)", call. = FALSE) - } - - if (comb == "struc") { - stop("struc use C matrix", call. = FALSE) + if (missing(C)) { + if (missing(Ut)) { + stop("Please, give C or Ut.", call. = FALSE) + } else if(missing(nb)){ + ctf <- ctf_tools(Ut = Ut, m = m, h = 1, sparse = TRUE) + } else { + ctf <- ctf_tools(Ut = Ut, nb = nb, m = m, h = 1, sparse = TRUE) } + }else{ + ctf <- ctf_tools(C = C, m = m, h = 1, sparse = TRUE) + } - if (n < 3) { - stop("The hierarchy must have, at least, three series", call. = FALSE) - } + n <- ctf$hts$n + na <- ctf$hts$na + nb <- ctf$hts$nb + #C <- ctf$hts$C + #S <- ctf$hts$S + #Ut <- ctf$hts$Ut - if (NCOL(Ut) != n) { - stop("Incorrect dimension of Ut or basef (they don't have same columns)", call. = FALSE) - } + kset <- ctf$thf$kset + m <- ctf$thf$m + p <- ctf$thf$p + kt <- ctf$thf$kt + ks <- ctf$thf$ks - if (n <= nb) { - stop("n <= nb, total number of TS is less (or equal) than the number of bottom TS", call. = FALSE) - } + # matrix + #Zt <- tools$Zt - na <- n - nb - Ut <- Matrix(Ut, sparse = TRUE) - } else { - if (!(is.matrix(C) | is(C, "Matrix"))) stop("C must be a matrix", call. = FALSE) + if (NROW(basef) != n) { + stop("Incorrect dimension of Ut or basef (they don't have same columns).", call. = FALSE) + } - C <- Matrix(C, sparse = TRUE) + Ht <- ctf$ctf$Ht + S <- ctf$ctf$Stilde - nb <- NCOL(C) - na <- NROW(C) + # Base forecast + if (NCOL(basef) %% kt != 0) { + stop("basef has a number of row not in line with the frequency of the series", call. = FALSE) + } - if (n < 3) { - stop("The hierarchy must have, at least, three series", call. = FALSE) + if(nn){ + if(nn_type == "fbpp" | nn_type == "KAnn"){ + type = "M" } + } - if (n <= nb) { - stop("n <= nb, total number of TS is less (or equal) than the number of bottom TS", call. = FALSE) + if(is.null(S)){ + if (type == "S") { + stop("Type = S needs the cross-sectional C matrix. \nPlease consider using ut2c() to get the right C when working with Ut.", call. = FALSE) } - if (na + nb != n) { - stop("na + nb != n, matrix C or basef has incorrect dimension", call. = FALSE) + if(nn_type == "sntz"){ + stop("nn_type = sntz needs the cross-sectional C matrix.\nPlease consider using ut2c() to get the right C when working with Ut.", call. = FALSE) } - Ut <- cbind(Diagonal(na), -C) - } - - # H - P <- commat(n, kt) - Us <- cbind(Matrix(0, NROW(Ut) * m, n * ks), kronecker(Diagonal(m), Ut)) %*% t(P) - Ht <- rbind(Us, kronecker(Diagonal(n), Zt)) - - # Base forecast - if (NCOL(basef) %% kt != 0) { - stop("basef has a number of row not in line with the frequency of the series", call. = FALSE) + if (comb == "struc") { + stop("Param comb equal struc needs the cross-sectional C matrix.\nPlease consider using ut2c() to get the right C when working with Ut.", call. = FALSE) + } } h <- NCOL(basef) / kt - Dh <- Dmat(h = h, kset = kset, n = n) + Dh <- Dmat(h = h, m = kset, n = n) Ybase <- matrix(Dh %*% as.vector(t(basef)), nrow = h, byrow = TRUE) # In-sample errors @@ -271,7 +389,7 @@ octrec.default <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, N <- NCOL(res) / kt - DN <- Dmat(h = N, kset = kset, n = n) + DN <- Dmat(h = N, m = kset, n = n) E <- matrix(DN %*% as.vector(t(res)), nrow = N, byrow = TRUE) if (comb == "sam" & N < n * kt) { @@ -297,17 +415,6 @@ octrec.default <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, shr_mod <- function(x, ...) shrink_estim(x, minT = mse)[[1]] } - if (type == "S" | comb == "struc") { - Qtilde <- commat(nb, kt) %*% bdiag(commat(ks, nb), commat(m, nb)) - Q <- bdiag(Diagonal(na * kt), Qtilde) - if (missing(Sstruc)) { - Htstruc <- Matrix(pracma::rref(as.matrix(Ht %*% Q)), sparse = TRUE) - Cstruc <- -Htstruc[, (nrow(Htstruc) + 1):ncol(Htstruc)] - Sstruc <- rbind(Cstruc, Diagonal(m * nb)) - } - S <- Q %*% Sstruc - } - switch(comb, ols = { Omega <- .sparseDiagonal(n * kt) @@ -346,6 +453,7 @@ octrec.default <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, if (is.null(W)) { stop("Please, put in option W your covariance matrix", call. = FALSE) } + P <- commat(n, kt) Omega <- P %*% W %*% t(P) }, omega = { @@ -373,16 +481,17 @@ octrec.default <- function(basef, m, C, comb, res, Ut, nb, Sstruc, mse = TRUE, if (type == "S") { rec_sol <- recoS( basef = Ybase, W = Omega, S = S, sol = sol, nn = nn, keep = keep, - settings = settings, b_pos = b_pos, bounds = bounds + settings = settings, b_pos = b_pos, bounds = bounds, nn_type = nn_type ) } else { rec_sol <- recoM( - basef = Ybase, W = Omega, H = Ht, sol = sol, nn = nn, keep = keep, - settings = settings, b_pos = b_pos, bounds = bounds + basef = Ybase, W = Omega, Ht = Ht, sol = sol, nn = nn, keep = keep, S = S, + settings = settings, b_pos = b_pos, bounds = bounds, nn_type = nn_type ) } if (keep == "list") { + P <- commat(n, kt) rec_sol$nn_check <- sum(rec_sol$recf < 0) rec_sol$rec_check <- all(rec_sol$recf %*% t(Ht) < 1e-6) @@ -453,7 +562,7 @@ octrec.list <- function(basef, m, ..., W, Omega){ n <- NROW(basef) - Dh <- Dmat(h = h, kset = tools$kset, n = n) + Dh <- Dmat(h = h, m = tools$kset, n = n) Ybase <- matrix(Dh %*% as.vector(t(basef)), nrow = h, byrow = TRUE) baseh <- lapply(split(Ybase,1:h), matrix, nrow = 3, byrow = TRUE) diff --git a/R/reco.R b/R/reco.R index 015ce29..1c39d86 100755 --- a/R/reco.R +++ b/R/reco.R @@ -1,10 +1,11 @@ # basef: base forecasts (h x col) # W: covariance # H: Ht in ctr, Zt in thf e Ut in hts -recoM <- function(basef, W, H, sol = "direct", nn = FALSE, settings, b_pos = NULL, - keep = "list", bounds = NULL) { +recoM <- function(basef, W, Ht, sol = "direct", nn = FALSE, nn_type = "osqp", settings, b_pos = NULL, + keep = "list", bounds = NULL, S) { sol <- match.arg(sol, c("direct", "osqp")) + nn_type <- match.arg(nn_type, c("osqp", "fbpp", "KAnn", "sntz")) if(!is.null(bounds)){ if(!is.matrix(bounds) | NCOL(bounds) != 2 | NROW(bounds) != NCOL(basef)){ @@ -15,87 +16,151 @@ recoM <- function(basef, W, H, sol = "direct", nn = FALSE, settings, b_pos = NUL } switch(sol, - direct = { - out <- list() - c_W <- ncol(W) - - lm_dx <- methods::as(H %*% t(basef), "CsparseMatrix") - lm_sx <- Matrix::Matrix(H %*% W %*% t(H), sparse = TRUE, forceCheck = TRUE) - out$recf <- t(.sparseDiagonal(NCOL(W)) %*% t(basef) - W %*% t(H) %*% solveLin(lm_sx, lm_dx)) - - if (nn) { - if (any(out$recf < (-1e-6))) { - if (isDiagonal(W)) { - P <- .sparseDiagonal(x = diag(W)^(-1)) - } else { - R <- chol(W) - P <- chol2inv(R) - } - - rec <- apply(basef, 1, function(x) { - M_osqp( - y = x, H = H, P = P, nn = nn, bounds = bounds, - settings = settings, b_pos = b_pos - ) - }) - - out <- list() - out$recf <- do.call("rbind", extract(rec, 1)) - out$info <- do.call("rbind", extract(rec, 2)) - colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", "status", "status_polish") - - if (keep == "list") { - out$W <- as.matrix(W) - } - } else { - out$recf <- out$recf * (out$recf > 0) - - if (keep == "list") { - M <- .sparseDiagonal(c_W) - W %*% t(H) %*% solveLin(H %*% W %*% t(H), H) - out$varf <- diag(M %*% W) - out$M <- as.matrix(M) - out$W <- as.matrix(W) - } - } - } else if (keep == "list") { - M <- .sparseDiagonal(c_W) - W %*% t(H) %*% solveLin(H %*% W %*% t(H), H) - out$varf <- diag(M %*% W) - out$M <- as.matrix(M) - out$W <- as.matrix(W) - } - }, - osqp = { - if (isDiagonal(W)) { - P <- .sparseDiagonal(x = diag(W)^(-1)) - } else { - R <- chol(W) - P <- chol2inv(R) - } - - rec <- apply(basef, 1, function(x) { - M_osqp( - y = x, H = H, P = P, nn = nn, bounds = bounds, - settings = settings, b_pos = b_pos - ) - }) - - out <- list() - out$recf <- do.call("rbind", extract(rec, 1)) - out$info <- do.call("rbind", extract(rec, 2)) - colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", "status", "status_polish") - - if (keep == "list") { - out$W <- as.matrix(W) - } - } + direct = { + out <- list() + + lm_dx <- methods::as(Ht %*% t(basef), "CsparseMatrix") + lm_sx <- Matrix::Matrix(Ht %*% W %*% t(Ht), sparse = TRUE, forceCheck = TRUE) + out$recf <- t(.sparseDiagonal(NCOL(W)) %*% t(basef) - W %*% t(Ht) %*% solveLin(lm_sx, lm_dx)) + + if(nn & any(out$recf < (-sqrt(.Machine$double.eps)))){ + switch(nn_type, + osqp = { + if(isDiagonal(W)){ + P <- .sparseDiagonal(x = diag(W)^(-1)) + }else{ + R <- chol(W) + P <- chol2inv(R) + } + id <- which(rowSums(out$recf<(-sqrt(.Machine$double.eps)))!=0) + rec <- apply(basef[id,,drop = FALSE], 1, function(x){ + M_osqp(y = x, Ht = Ht, P = P, nn = nn, + bounds = bounds, settings = settings, b_pos = b_pos) + }) + rec <- do.call("rbind",rec) + out$recf[id,] <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", + "status", "status_polish") + rownames(out$info) <- id + + out$recf[-id,,drop=FALSE] <- out$recf[-id,,drop=FALSE] * (out$recf[-id,,drop=FALSE] > 0) + + if (keep == "list") { + out$W <- as.matrix(W) + } + }, + fbpp = { + param <- list(tol = 1e-6, + itmax = 10) + param[names(settings)] <- settings + + rec <- apply(out$recf, 1, function(x) { + M_fbpp(y = x, Ht = Ht, W = W, b_pos = b_pos, + param = param, keep = keep) + }) + + rec <- do.call("rbind",rec) + out$recf <- do.call("rbind", rec[,"recf"]) + out$idx <- unlist(rec[,"idx"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("run_time (s)", "iter", "status", "tol", "itmax", "lidx") + rownames(out$info) <- 1:NROW(out$recf) + out$info <- out$info[out$info[,2]!=0, , drop=FALSE] + + if(keep == "list"){ + out$varf <- do.call("rbind", rec[,"varf"]) + out$M <- unlist(rec[,"M"]) + if(length(out$M)){ + out$M <- out$M[[1]] + } + out$W <- as.matrix(W) + } + }, + KAnn = { + param <- list(tol = 1e-6, + itmax = 10) + param[names(settings)] <- settings + + rec <- apply(out$recf, 1, function(x) { + M_KAnn(y = x, Ht = Ht, W = W, b_pos = b_pos, + param = param, keep = keep) + }) + + rec <- do.call("rbind",rec) + out$recf <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("run_time (s)", "iter", "status", "tol", "itmax") + rownames(out$info) <- 1:NROW(out$recf) + out$info <- out$info[out$info[,2]!=0, , drop=FALSE] + + if(keep == "list"){ + out$varf <- do.call("rbind", rec[,"varf"]) + out$M <- unlist(rec[,"M"]) + if(length(out$M)){ + out$M <- out$M[[1]] + } + out$W <- as.matrix(W) + } + }, + sntz = { + na <- NROW(S)-NCOL(S) + bottom <- out$recf[,-c(1:na),drop = FALSE] + bottom <- bottom*(bottom > 0) + out$recf <- bottom %*% t(S) + + if (keep == "list") { + out$W <- as.matrix(W) + } + }) + }else{ + if(nn){ + out$recf <- out$recf * (out$recf > 0) + } + + if(keep == "list"){ + M <- unname(.sparseDiagonal(NCOL(W)) - W %*% t(Ht) %*% solveLin(Ht %*% W %*% t(Ht), Ht)) + out$varf <- diag(M %*% W) + out$M <- as.matrix(M) + out$W <- as.matrix(W) + } + } + }, + osqp = { + if(isDiagonal(W)){ + P <- .sparseDiagonal(x = diag(W)^(-1)) + }else{ + R <- chol(W) + P <- chol2inv(R) + } + + rec <- apply(basef, 1, function(x){ + M_osqp(y = x, Ht = Ht, P = P, nn = nn, bounds = bounds, + settings = settings, b_pos = b_pos) + }) + rec <- do.call("rbind",rec) + out <- list() + out$recf <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", + "status", "status_polish") + rownames(out$info) <- 1:NROW(out$recf) + + if(keep == "list"){ + out$W <- as.matrix(W) + } + } ) return(out) } recoS <- function(basef, W, S, sol = "direct", nn = FALSE, settings, b_pos = NULL, - keep = "list", bounds = NULL) { + keep = "list", bounds = NULL, nn_type = "osqp") { sol <- match.arg(sol, c("direct", "osqp")) + + nn_type <- match.arg(nn_type, c("osqp", "sntz")) + if(!is.null(bounds)){ if(!is.matrix(bounds) | NCOL(bounds) != 2 | NROW(bounds) != NCOL(basef)){ stop("bounds must be a matrix (", NCOL(basef), "x2)", call. = FALSE) @@ -105,148 +170,162 @@ recoS <- function(basef, W, S, sol = "direct", nn = FALSE, settings, b_pos = NUL } switch(sol, - direct = { - out <- list() - - if (isDiagonal(W)) { - Wm1 <- .sparseDiagonal(x = diag(W)^(-1)) - lm_dx1 <- methods::as(t(S) %*% Wm1 %*% t(basef), "CsparseMatrix") - lm_sx1 <- t(S) %*% Wm1 %*% S - out$recf <- t(S %*% solveLin(lm_sx1, lm_dx1)) - - if (nn) { - if (any(out$recf < (-1e-6))) { - P <- t(S) %*% Wm1 %*% S - q <- (-1) * t(Wm1 %*% S) - rec <- apply(basef, 1, function(x) { - S_osqp( - y = x, S = S, P = P, q = q, nn = nn, bounds = bounds, - settings = settings, b_pos = b_pos - ) - }) - - out <- list() - out$recf <- do.call("rbind", extract(rec, 1)) - out$info <- do.call("rbind", extract(rec, 2)) - colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", "status", "status_polish") - - if (keep == "list") { - out$W <- W - } - } else { - out$recf <- out$recf * (out$recf > 0) - if (keep == "list") { - lm_dx2 <- methods::as(t(S) %*% Wm1, "CsparseMatrix") - G <- S %*% solveLin(lm_sx1, lm_dx2) - out$varf <- diag(G %*% W) - out$G <- as.matrix(G) - out$W <- W - } - } - } else if (keep == "list") { - lm_dx2 <- methods::as(t(S) %*% Wm1, "CsparseMatrix") - G <- S %*% solveLin(lm_sx1, lm_dx2) - out$varf <- diag(G %*% W) - out$G <- as.matrix(G) - out$W <- W - } - } else { - Q <- solveLin(W, S) - - lm_dx1 <- methods::as(t(Q) %*% t(basef), "CsparseMatrix") - lm_sx1 <- t(S) %*% Q - out$recf <- t(S %*% solveLin(lm_sx1, lm_dx1)) - - if (nn) { - if (any(out$recf < (-1e-6))) { - P <- t(S) %*% Q - q <- (-1) * t(Q) - rec <- apply(basef, 1, function(x) { - S_osqp( - y = x, S = S, P = P, q = q, nn = nn, bounds = bounds, - settings = settings, b_pos = b_pos - ) - }) - - out <- list() - out$recf <- do.call("rbind", extract(rec, 1)) - out$info <- do.call("rbind", extract(rec, 2)) - colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", "status", "status_polish") - - if (keep == "list") { - out$W <- W - } - } else { - out$recf <- out$recf * (out$recf > 0) - if (keep == "list") { - G <- S %*% solveLin(lm_sx1, t(Q)) - out$varf <- diag(G %*% W) - out$G <- as.matrix(G) - out$W <- W - } - } - } else if (keep == "list") { - G <- S %*% solveLin(lm_sx1, t(Q)) - out$varf <- diag(G %*% W) - out$G <- as.matrix(G) - out$W <- W - } - } - }, - osqp = { - if (isDiagonal(W)) { - Q <- .sparseDiagonal(x = diag(W)^(-1)) - P <- t(S) %*% Q %*% S - q <- (-1) * t(Q %*% S) - } else { - Q <- solveLin(W, S) - P <- t(S) %*% Q - q <- (-1) * t(Q) - } - - rec <- apply(basef, 1, function(x) { - S_osqp( - y = x, S = S, q = q, P = P, nn = nn, bounds = bounds, - settings = settings, b_pos = b_pos - ) - }) - - out <- list() - out$recf <- do.call("rbind", extract(rec, 1)) - out$info <- do.call("rbind", extract(rec, 2)) - colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", "status", "status_polish") - - if (keep == "list") { - out$W <- W - } - } + direct = { + out <- list() + + if(isDiagonal(W)){ + Wm1 <- .sparseDiagonal(x = diag(W)^(-1)) + lm_dx1 <- methods::as(t(S) %*% Wm1 %*% t(basef), "CsparseMatrix") + lm_sx1 <- t(S) %*% Wm1 %*% S + out$recf <- t(S %*% solveLin(lm_sx1, lm_dx1)) + + if(nn & any(out$recf < (-sqrt(.Machine$double.eps)))){ + switch(nn_type, + osqp = { + P <- t(S) %*% Wm1 %*% S + q <- (-1) * t(Wm1 %*% S) + + id <- which(rowSums(out$recf<(-sqrt(.Machine$double.eps)))!=0) + rec <- apply(basef[id,,drop = FALSE], 1, function(x) { + S_osqp(y = x, S = S, P = P, q = q, nn = nn, bounds = bounds, + settings = settings, b_pos = b_pos) + }) + rec <- do.call("rbind",rec) + out$recf[id,] <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", + "status", "status_polish") + rownames(out$info) <- id + + out$recf[-id,,drop=FALSE] <- out$recf[-id,,drop=FALSE] * (out$recf[-id,,drop=FALSE] > 0) + }, + sntz = { + na <- NROW(S)-NCOL(S) + bottom <- out$recf[,-c(1:na),drop = FALSE] + bottom <- bottom*(bottom > 0) + out$recf <- bottom %*% t(S) + }) + + if (keep == "list") { + out$W <- as.matrix(W) + } + }else{ + if(nn){ + out$recf <- out$recf * (out$recf > 0) + } + if(keep == "list"){ + lm_dx2 <- methods::as(t(S) %*% Wm1, "CsparseMatrix") + G <- unname(solveLin(lm_sx1, lm_dx2)) + M <- S %*% G + out$varf <- diag(M %*% W) + out$G <- as.matrix(G) + out$M <- as.matrix(M) + out$W <- W + } + } + }else{ + Q <- solveLin(W, S) + + lm_dx1 <- methods::as(t(Q) %*% t(basef), "CsparseMatrix") + lm_sx1 <- t(S) %*% Q + out$recf <- t(S %*% solveLin(lm_sx1, lm_dx1)) + + if(nn & any(out$recf < (-sqrt(.Machine$double.eps)))){ + switch(nn_type, + osqp = { + P <- t(S) %*% Q + q <- (-1) * t(Q) + id <- which(rowSums(out$recf<(-sqrt(.Machine$double.eps)))!=0) + rec <- apply(basef[id, , drop = FALSE], 1, function(x) { + S_osqp(y = x, S = S, P = P, q = q, nn = nn, bounds = bounds, + settings = settings, b_pos = b_pos) + }) + rec <- do.call("rbind",rec) + out$recf[id,] <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", + "status", "status_polish") + rownames(out$info) <- id + }, + sntz = { + na <- NROW(S)-NCOL(S) + bottom <- out$recf[,-c(1:na),drop = FALSE] + bottom <- bottom*(bottom > 0) + out$recf <- bottom %*% t(S) + }) + + if (keep == "list") { + out$W <- as.matrix(W) + } + }else{ + if(nn){ + out$recf <- out$recf * (out$recf > 0) + } + if(keep == "list"){ + G <- unname(solveLin(lm_sx1, t(Q))) + M <- S %*% G + out$varf <- diag(M %*% W) + out$G <- as.matrix(G) + out$M <- as.matrix(M) + out$W <- W + } + } + } + }, + osqp = { + if(isDiagonal(W)){ + Q <- .sparseDiagonal(x = diag(W)^(-1)) + P <- t(S) %*% Q %*% S + q <- (-1) * t(Q %*% S) + }else{ + Q <- solveLin(W, S) + P <- t(S) %*% Q + q <- (-1) * t(Q) + } + + rec <- apply(basef, 1, function(x) { + S_osqp(y = x, S = S, q = q, P = P, nn = nn, bounds = bounds, + settings = settings, b_pos = b_pos) + }) + rec <- do.call("rbind",rec) + out <- list() + out$recf <- do.call("rbind", rec[,"recf"]) + out$info <- do.call("rbind", rec[,"info"]) + colnames(out$info) <- c("obj_val", "run_time", "iter", "pri_res", + "status", "status_polish") + rownames(out$info) <- 1:NROW(out$recf) + + if(keep == "list"){ + out$W <- W + } + } ) return(out) } # y = riga di basef -M_osqp <- function(y, P = NULL, H = NULL, nn = NULL, settings = NULL, +M_osqp <- function(y, P = NULL, Ht = NULL, nn = NULL, settings = NULL, b_pos = NULL, bounds = NULL) { - c <- ncol(H) - r <- nrow(H) + c <- ncol(Ht) + r <- nrow(Ht) # Linear constrains H = 0 l <- rep(0, r) u <- rep(0, r) - A <- H + A <- Ht q <- (-1) * t(P) %*% as.vector(y) - if (nn) { + if(nn){ # bts >= 0 - A <- rbind(A, Diagonal(c)[b_pos == 1, ]) + A <- rbind(A, .sparseDiagonal(c)[b_pos == 1, ]) l <- c(l, rep(0, sum(b_pos))) u <- c(u, rep(Inf, sum(b_pos))) } if(!is.null(bounds)){ bounds_rows <- rowSums(abs(bounds) == Inf) < 2 - A <- rbind(A, Diagonal(c)[bounds_rows, ]) + A <- rbind(A, .sparseDiagonal(c)[bounds_rows, ]) l <- c(l, bounds[bounds_rows, 1, drop = TRUE]) u <- c(u, bounds[bounds_rows, 2, drop = TRUE]) } @@ -264,25 +343,23 @@ M_osqp <- function(y, P = NULL, H = NULL, nn = NULL, settings = NULL, out <- list() out$recf <- rec$x - if (rec$info$status_val != 1 | rec$info$pri_res > 1e-6) { + if(rec$info$status_val != 1){ warning(paste("OSQP flag", rec$info$status_val, "OSQP pri_res", rec$info$pri_res), call. = FALSE) } - if (nn) { - out$recf[which(out$recf < 0 & out$recf > -1e-6)] <- 0 + if(nn){ + out$recf[which(out$recf < 0)] <- 0 } - out$info <- c( - rec$info$obj_val, rec$info$run_time, rec$info$iter, rec$info$pri_res, - rec$info$status_val, rec$info$status_polish - ) + out$info <- c(rec$info$obj_val, rec$info$run_time, rec$info$iter, rec$info$pri_res, + rec$info$status_val, rec$info$status_polish) return(out) } -extract <- function(l, num) { - do.call("rbind", l)[, num] -} +#extract <- function(l, num) { +# do.call("rbind", l)[, num] +#} # S_sol <- function(basef, W, S){ # out <- list() @@ -317,7 +394,7 @@ S_osqp <- function(y, q = NULL, P = NULL, S = NULL, nn = NULL, l <- NULL u <- NULL - if (nn) { + if(nn){ A <- .sparseDiagonal(c) l <- rep(0, sum(b_pos)) u <- rep(Inf, sum(b_pos)) @@ -331,11 +408,11 @@ S_osqp <- function(y, q = NULL, P = NULL, S = NULL, nn = NULL, } if(length(settings)==0){ - settings = osqpSettings(verbose = FALSE, - eps_abs = 1e-5, - eps_rel = 1e-5, - polish_refine_iter = 100, - polish=TRUE) + settings <- osqpSettings(verbose = FALSE, + eps_abs = 1e-5, + eps_rel = 1e-5, + polish_refine_iter = 100, + polish=TRUE) } rec <- solve_osqp(P, q, A, l, u, settings) @@ -343,36 +420,164 @@ S_osqp <- function(y, q = NULL, P = NULL, S = NULL, nn = NULL, out <- list() out$recf <- as.vector(S %*% rec$x) - if (rec$info$status_val != 1 | rec$info$pri_res > 1e-6) { + if(rec$info$status_val != 1){ warning(paste("OSQP flag", rec$info$status_val, "OSQP pri_res", rec$info$pri_res), call. = FALSE) } - if (nn) { - out$recf[which(out$recf < 0 & out$recf > -1e-6)] <- 0 + if(nn){ + out$recf[which(out$recf < 0)] <- 0 } - out$info <- c( - rec$info$obj_val, rec$info$run_time, rec$info$iter, - rec$info$pri_res, rec$info$status_val, rec$info$status_polish - ) + out$info <- c(rec$info$obj_val, rec$info$run_time, rec$info$iter, + rec$info$pri_res, rec$info$status_val, rec$info$status_polish) return(out) } -solveLin <- function(msx, mdx) { - tryCatch(solve(msx, mdx), error = function(cond) { +M_fbpp <- function(y, Ht, W, b_pos, param, keep){ + tol <- param$tol + itmax <- param$itmax + if(all(y>(-tol))){ + out <- list() + out$recf <- y*(y>0) + out$info <- c(0, 0, tol, 0, itmax, 0) + out$idx <- NA + + if(keep == "list"){ + out$M <- .sparseDiagonal(NCOL(W)) - W %*% t(Ht) %*% solveLin(Ht %*% W %*% t(Ht), Ht) + out$varf <- diag(M%*%W) + } + return(out) + } + start <- Sys.time() + bid <- Position(function(x) x == 1, b_pos) + idx <- c() + id <- list() + Htnn <- Ht + if(keep == "list"){ + #M <- list() + #M[[1]] <- .sparseDiagonal(NCOL(W)) - W %*% t(Ht) %*% solveLin(Ht %*% W %*% t(Ht), Ht) + M <- .sparseDiagonal(NCOL(W)) - W %*% t(Ht) %*% solveLin(Ht %*% W %*% t(Ht), Ht) + } + + for(i in 1:itmax){ + id[[i]] <- which(y < -tol) + idx <- unique(c(idx, id[[i]][id[[i]]>=bid])) + block <- sparseMatrix(i = 1:length(idx), + j = idx, + x = 1, + dims = c(length(idx), NCOL(Htnn))) + Htnn <- rbind(Ht, block) + lm_dx <- methods::as(Htnn %*% y, "CsparseMatrix") + lm_sx <- Matrix::Matrix(Htnn %*% W %*% t(Htnn), sparse = TRUE, forceCheck = TRUE) + y <- as.vector(.sparseDiagonal(NCOL(W)) %*% y - W %*% t(Htnn) %*% solveLin(lm_sx, lm_dx, verb = FALSE)) + + if(keep == "list"){ + M2 <- .sparseDiagonal(NCOL(W)) - W %*% t(Htnn) %*% solveLin(lm_sx, Htnn, verb = FALSE) + M <- M2%*%M + #M[[i+1]] <- M2 + } - # browser() - warning( - "An error in LU decomposition occurred, with the following message:\n", - cond$message, "\n Trying QR decomposition instead...", call. = FALSE - ) - tryCatch(solve(qr(msx), mdx), error = function(cond) { + if(all(y >= (-tol))){ + flag <- 1 + break + }else{ + flag <- -2 + } + } - # browser() + if(flag == 1){ + y <- (y > 0)*y + (y < 0) * 0 + if(i == itmax){ + flag <- 2 + } + } + + end <- Sys.time() + out <- list() + out$recf <- y + out$info <- c(difftime(end,start, units = "secs"), i, flag, tol, itmax, length(idx)) + out$idx <- id + + if(keep == "list"){ + out$M <- M + out$varf <- diag(M%*%W) + } + return(out) +} + + +M_KAnn <- function(y, Ht, W, b_pos, param, keep){ + tol <- param$tol + itmax <- param$itmax + if(all(y>(-tol))){ + out <- list() + out$recf <- y*(y>0) + out$info <- c(0, 0, tol, 0, itmax) + + if(keep == "list"){ + out$M <- .sparseDiagonal(NCOL(W)) - W %*% t(Ht) %*% solveLin(Ht %*% W %*% t(Ht), Ht) + out$varf <- diag(M%*%W) + } + return(out) + } + start <- Sys.time() + lm_sx <- Matrix::Matrix(Ht %*% W %*% t(Ht), sparse = TRUE, forceCheck = TRUE) + if(keep == "list"){ + M <- .sparseDiagonal(NCOL(W)) - W %*% t(Ht) %*% solveLin(lm_sx, Ht) + } + id <- list() + for(i in 1:itmax){ + id[[i]] <- (y > -tol) + yc <- y * id[[i]] + + lm_dx <- methods::as(Ht %*% yc, "CsparseMatrix") + y <- as.vector(.sparseDiagonal(NCOL(W)) %*% yc - W %*% t(Ht) %*% solveLin(lm_sx, lm_dx, verb = FALSE)) + + if(keep == "list"){ + M1 <- M%*%diag(id[[i]]) + M <- M1%*%M + } + + if(all(y >= (-tol))){ + flag <- 1 + break + }else{ + flag <- -2 + } + } + + if(flag == 1){ + y <- (y > 0)*y + (y < 0) * 0 + if(i == itmax){ + flag <- 2 + } + } + + end <- Sys.time() + out <- list() + out$recf <- y + out$info <- c(difftime(end,start, units = "secs"), i, flag, tol, itmax) + + if(keep == "list"){ + out$M <- M + out$varf <- diag(M%*%W) + } + return(out) +} + +solveLin <- function(msx, mdx, verb = FALSE) { + tryCatch(solve(msx, mdx), error = function(cond){ + if(verb){ warning( - "An error in QR decomposition occurred, with the following message:\n", - cond$message, "\n Trying chol decomposition instead...", call. = FALSE - ) + "An error in LU decomposition occurred, with the following message:\n", + cond$message, "\n Trying QR decomposition instead...", call. = FALSE) + } + tryCatch(solve(qr(msx), mdx), error = function(cond){ + if(verb){ + warning( + "An error in QR decomposition occurred, with the following message:\n", + cond$message, "\n Trying chol decomposition instead...", call. = FALSE) + } backsolve(chol(msx), mdx) }) }) diff --git a/R/score_index.R b/R/score_index.R index 33deb4e..d705366 100755 --- a/R/score_index.R +++ b/R/score_index.R @@ -1,24 +1,37 @@ -#' @title Measuring forecasting accuracy +#' @title Measuring accuracy in a rolling forecast experiment #' -#' @description Function to calculate the accuracy indices of the reconciled point +#' @description +#' \loadmathjax +#' Function to calculate the accuracy indices of the reconciled point #' forecasts of a cross-temporal (not only, see examples) system (more in #' \href{https://danigiro.github.io/FoReco/articles/accuracy_indices.html}{Average relative accuracy indices}). +#' (\emph{Experimental version}) #' -#' @param recf list of q (forecast origins) reconciled forecasts' matrices (\code{[n x h(k* + m)]} in the -#' cross-temporal, \code{[h x n]} in cross-sectional and vectors of length \code{[h(k* + m)]} in temporal framework). -#' @param base list of q (forecast origins) base forecasts' matrices (\code{[n x h(k* + m)]} in the -#' cross-temporal, \code{[n x h]} in cross-sectional and \code{[h(k* + m) x 1]} in temporal framework). -#' @param test list of q (forecast origins) test observations' matrices (\code{[n x h(k* + m)]} in the -#' cross-temporal, \code{[n x h]} in cross-sectional and \code{[h(k* + m) x 1]} in temporal framework). -#' @param type type of accuracy measure ("\code{mse}" Mean Square Error, "\code{rmse}" Root Mean Square Error -#' or "\code{mae}" Mean Absolute Error) -#' @param m highest frequency of the forecasted time series. +#' @param recf list of q (forecast origins) reconciled forecasts' matrices +#' (\mjseqn{[n \times h(k^\ast + m)]} in the cross-temporal case, +#' \mjseqn{[h \times n]} in the cross-sectional case, and vectors of length +#' \mjseqn{[h(k^\ast \times m)]} in the temporal framework). +#' @param base list of q (forecast origins) base forecasts' matrices +#' (\mjseqn{[n \times h(k^\ast + m)]} in the cross-temporal case, +#' \mjseqn{[h \times n]} in the cross-sectional case, and vectors of length +#' \mjseqn{[h(k^\ast \times m)]} in the temporal framework). +#' @param test list of q (forecast origins) test observations' matrices +#' (\mjseqn{[n \times h(k^\ast + m)]} in the cross-temporal case, +#' \mjseqn{[h \times n]} in the cross-sectional case, and vectors of length +#' \mjseqn{[h(k^\ast \times m)]} in the temporal framework). +#' @param type type of accuracy measure ("\code{mse}" Mean Square Error, +#' "\code{rmse}" Root Mean Square Error or "\code{mae}" Mean Absolute Error). +#' @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 nb number of bottom time series in the cross-sectional framework. -#' @param compact if TRUE return only the summary matrix. +#' @param compact if TRUE returns only the summary matrix. +#' @param nl (\mjseqn{L \times 1}) vector containing the number of time series +#' in each cross-sectional level of the hierarchy (\code{nl[1] = 1}). #' #' @return #' It returns a summary table called \code{Avg_mat} (if \code{compact} option is \code{TRUE}, -#' \emph{default}), otherwise it returns a list of four tables (more in +#' \emph{default}), otherwise it returns a list of six tables (more in #' \href{https://danigiro.github.io/FoReco/articles/accuracy_indices.html}{Average relative accuracy indices}). #' #' @references @@ -52,7 +65,6 @@ #' hts_score <- score_index(recf = mrecf, base = mbase, test = mtest, nb = 5) #' #' # Temporal framework -#' data(FoReco_data) #' # top ts base forecasts ([lowest_freq' ... highest_freq']') #' topbase <- FoReco_data$base[1, ] #' # top ts residuals ([lowest_freq' ... highest_freq']') @@ -66,9 +78,10 @@ #' } #' #' @keywords utilities +#' @family utilities #' #' @export -score_index <- function(recf, base, test, m, nb, type = "mse", compact = TRUE) { +score_index <- function(recf, base, test, m, nb, nl, type = "mse", compact = TRUE) { type <- match.arg(type, c("mse", "mae", "rmse")) @@ -127,24 +140,59 @@ score_index <- function(recf, base, test, m, nb, type = "mse", compact = TRUE) { } # Some usefull tools - kset <- rev(divisors(m)) + if(length(m)==1){ + if(m == 1){ + kset <- 1 + }else{ + thf_obj <- thf_tools(m = m) + m <- thf_obj$m + kset <- thf_obj$kset + } + }else{ + thf_obj <- thf_tools(m = m) + m <- thf_obj$m + kset <- thf_obj$kset + } kt <- sum(kset) p <- length(kset) n <- dim(Ebase)[1] na <- n - nb h <- dim(Ebase)[2] / kt q <- dim(Ebase)[3] - kpos <- rep(kset, rep(rev(kset * h))) # position of k in colums of E[,,q] + kpos <- rep(kset, rep(m/kset * h)) # position of k in colums of E[,,q] kpos <- factor(kpos, kset, ordered = TRUE) if (m == 1) { - kh <- paste("k", kpos, "h", 1:h, sep = "") + kh <- paste("h", 1:h, sep = "") + khc <- paste("h1:", 1:h, sep = "") } else { kh <- paste("k", kpos, "h", - do.call("c", as.list(sapply(rev(kset) * h, function(x) seq(1:x)))), - sep = "" + do.call("c", lapply(rev(kset) * h, function(x) seq(1:x))), + sep = "" + ) + + khc <- paste("k", kpos, "h1:", + do.call("c", lapply(rev(kset) * h, function(x) seq(1:x))), + sep = "" ) } + if(missing(nl)){ + nl <- na + }else if(sum(nl) != na){ + stop("Please, provide a valid nl vector s.t. sum(nl) == na", call. = FALSE) + } + + nlnb <- c(nl, nb) + L <- length(nlnb) + nl <- c(1,2) + levpos <- rep(1:L, nlnb) + + if(L>2){ + csname <- paste0("L", 1:L) + }else{ + csname <- c("uts", "bts") + } + IND_base <- apply(Qbase, c(1, 2), mean) IND_recf <- apply(Qrecf, c(1, 2), mean) @@ -156,8 +204,17 @@ score_index <- function(recf, base, test, m, nb, type = "mse", compact = TRUE) { RelIND_ikh <- IND_recf / IND_base colnames(RelIND_ikh) <- kh logRelIND_ikh <- log(RelIND_ikh) - logRelIND_ikh[abs(logRelIND_ikh) == Inf] <- NA + #logRelIND_ikh[abs(logRelIND_ikh) == Inf] <- NA + logRelIND_ikh[logRelIND_ikh == -Inf] <- log(1e-16) + if(any(logRelIND_ikh == Inf)){ + if(all(logRelIND_ikh == Inf)){ + logRelIND_ikh[abs(logRelIND_ikh) == Inf] <- NA + warning("Base forecasts perfect prediction",call. = FALSE) + }else{ + logRelIND_ikh[logRelIND_ikh == Inf] <- 10*max(logRelIND_ikh[logRelIND_ikh != Inf]) + } + } # 24 AvgRelIND <- exp(mean(logRelIND_ikh, na.rm = TRUE)) @@ -165,29 +222,25 @@ score_index <- function(recf, base, test, m, nb, type = "mse", compact = TRUE) { AvgRelIND_i <- exp(rowMeans(logRelIND_ikh, na.rm = TRUE)) # 27 AvgRelIND__kh <- exp(colMeans(logRelIND_ikh, na.rm = TRUE)) - AvgRelIND_akh <- exp(colMeans(logRelIND_ikh[1:na, , drop = FALSE], na.rm = TRUE)) - AvgRelIND_bkh <- exp(colMeans(logRelIND_ikh[-c(1:na), , drop = FALSE], na.rm = TRUE)) - Avg_k <- rbind(AvgRelIND__kh, AvgRelIND_akh, AvgRelIND_bkh) - rownames(Avg_k) <- c("all", "uts", "bts") + AvgRelIND_lkh <- lapply(1:L, function(x) exp(colMeans(logRelIND_ikh[levpos == x, , drop = FALSE], na.rm = TRUE))) + AvgRelIND_lkh <- do.call(rbind, AvgRelIND_lkh) + Avg_k <- rbind(AvgRelIND__kh, AvgRelIND_lkh) + rownames(Avg_k) <- c("all", csname) - # solo k + # only k AvgRelIND__k <- exp(tapply(colMeans(logRelIND_ikh, na.rm = TRUE), kpos, mean, na.rm = TRUE)) # bottom and aggregate - AvgRelIND_ak <- exp(tapply(colMeans(logRelIND_ikh[1:na, , drop = FALSE], na.rm = TRUE), - kpos, mean, na.rm = TRUE)) - AvgRelIND_bk <- exp(tapply(colMeans(logRelIND_ikh[-c(1:na), , drop = FALSE], na.rm = TRUE), - kpos, mean, na.rm = TRUE)) - AvgRelIND_a <- exp(mean(logRelIND_ikh[1:na, , drop = FALSE], na.rm = TRUE)) - AvgRelIND_b <- exp(mean(logRelIND_ikh[-c(1:na), , drop = FALSE], na.rm = TRUE)) - - Avg_matrix <- data.frame( - all = c(AvgRelIND__k, AvgRelIND), - uts = c(AvgRelIND_ak, AvgRelIND_a), bts = c(AvgRelIND_bk, AvgRelIND_b) - ) + AvgRelIND_lk <- lapply(1:L, function(x) exp(tapply(colMeans(logRelIND_ikh[levpos == x, , drop = FALSE], na.rm = TRUE), + kpos, mean, na.rm = TRUE))) + AvgRelIND_l <- lapply(1:L, function(x) exp(mean(logRelIND_ikh[levpos == x, , drop = FALSE], na.rm = TRUE))) + Avg_matrix <- cbind(c(AvgRelIND__k, AvgRelIND), + as.data.frame(t(cbind(do.call(rbind, AvgRelIND_lk), + do.call(rbind, AvgRelIND_l))))) + colnames(Avg_matrix) <- c("all", csname) rownames(Avg_matrix) <- c(unique(kset), "all") - # per ogni serie + # for each series AvgRelIND_ik <- exp(t(apply(logRelIND_ikh, 1, function(z) tapply(z, kpos, mean, na.rm = TRUE)))) all <- AvgRelIND_i @@ -197,26 +250,49 @@ score_index <- function(recf, base, test, m, nb, type = "mse", compact = TRUE) { AvgRelIND_ik <- cbind(AvgRelIND_ik, all) } + # CUMULATIVE forecast horizon + AvgRelIND__khc <- exp(do.call(c, lapply(unique(kpos), function(y){ + cummean(colMeans(logRelIND_ikh[, kpos == y, drop = FALSE], na.rm = TRUE)) + }))) + AvgRelIND_lkhc <- lapply(1:L, function(x) exp(do.call(c, lapply(unique(kpos), function(y){ + cummean(colMeans(logRelIND_ikh[levpos == x, kpos == y, drop = FALSE], na.rm = TRUE)) + })))) + AvgRelIND_lkhc <- do.call(rbind, AvgRelIND_lkhc) + Avg_kc <- rbind(AvgRelIND__khc, AvgRelIND_lkhc) + rownames(Avg_kc) <- c("all", csname) + colnames(Avg_kc) <- khc + + RelIND_ikhc <- t(apply(logRelIND_ikh, 1, + function(x) + exp(do.call(c, lapply(unique(kpos), + function(y){ + cummean(x[kpos == y]) + }))))) + colnames(RelIND_ikhc) <- khc + if (compact == TRUE) { if (nb == 1) { - return(Avg_matrix[, 1, drop = FALSE]) + return(Avg_matrix[, "all", drop = FALSE]) } else if (m == 1) { - return(Avg_matrix[2, , drop = FALSE]) + return(Avg_matrix["all", , drop = FALSE]) } else { return(Avg_matrix) } } else { if (nb == 1) { out <- list() - out$Avg_mat <- Avg_matrix[, 1, drop = FALSE] - out$Avg_k <- Avg_k[1, , drop = FALSE] + out$Avg_mat <- Avg_matrix[, "all", drop = FALSE] + out$Avg_k <- Avg_k["all", , drop = FALSE] + out$Avg_k_cum <- Avg_kc["all", , drop = FALSE] return(out) } else if (m == 1) { out <- list() - out$Avg_mat <- Avg_matrix[2, , drop = FALSE] - out$Avg_ik <- AvgRelIND_ik[, 2, drop = FALSE] + out$Avg_mat <- Avg_matrix["all", , drop = FALSE] + out$Avg_ik <- AvgRelIND_ik[, "all", drop = FALSE] out$Rel_mat <- RelIND_ikh out$Avg_k <- Avg_k + out$Rel_mat_cum <- RelIND_ikhc + out$Avg_k_cum <- Avg_kc return(out) } else { out <- list() @@ -224,7 +300,15 @@ score_index <- function(recf, base, test, m, nb, type = "mse", compact = TRUE) { out$Avg_ik <- AvgRelIND_ik out$Rel_mat <- RelIND_ikh out$Avg_k <- Avg_k + out$Rel_mat_cum <- RelIND_ikhc + out$Avg_k_cum <- Avg_kc return(out) } } } + + +cummean <- function(x){ + #x[is.na(x)] <- 0 + cumsum(x) / seq_along(x) +} diff --git a/R/shrink_estim.R b/R/shrink_estim.R index 2b36276..6a28255 100755 --- a/R/shrink_estim.R +++ b/R/shrink_estim.R @@ -23,6 +23,7 @@ #' \href{https://CRAN.R-project.org/package=hts}{https://CRAN.R-project.org/package=hts}. #' #' @keywords utilities +#' @family utilities #' #' @export shrink_estim <- function(x, minT = T) { diff --git a/R/srref.R b/R/srref.R new file mode 100644 index 0000000..7d1c024 --- /dev/null +++ b/R/srref.R @@ -0,0 +1,108 @@ +#' (Sparse) Reduced Row Echelon Form of a matrix +#' +#' Returns the reduced row echelon form of a (possibly sparse) matrix through Gauss-Jordan +#' elimination. This function is used by the \code{\link{ut2c}} function +#' to obtain a 'structural representation' of a general linearly constrained +#' multiple time series, from its zero-constraints kernel representation +#' (Di Fonzo and Girolimetto, 2020). +#' +#' @param A A numeric (possibly sparse) matrix. +#' @param tol Tolerance for checking for 0 pivot. +#' @param verbose If \code{TRUE}, print intermediate steps. +#' @param sparse Option to return sparse matrices (\emph{default} is \code{TRUE}). +#' +#' @return A (sparse) matrix with the same dimension as \mjseqn{\mathbf{A}} +#' +#' @author Originally written by John Fox and modified by Geoffrey Brent to fix a bug, +#' extended to deal with sparse matrices. +#' +#' @family utilities +#' +#' @references +#' Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: +#' Optimal Combination Method and Heuristic Alternatives, Department of Statistical +#' Sciences, University of Padua, \href{https://arxiv.org/abs/2006.08570}{arXiv:2006.08570}. +#' +#' @export +#' +#' @import Matrix +srref <- function(A, tol=sqrt(.Machine$double.eps), verbose=FALSE, sparse = TRUE){ + if(is.matrix(A) | is(A, "Matrix")){ + A <- as(A, "dgCMatrix") + }else{ + stop("A must be a matrix", call. = FALSE) + } + + n <- nrow(A) + m <- ncol(A) + + # Verbose condition + if(verbose){ + quant <- round(stats::quantile(1:min(n,m), probs = seq(0, 1, 0.25)[-1])) + if(min(n,m)<100){ + quant <- n + m + }else{ + message("Reduced Row Echelon Form (%): 0%...", appendLF = FALSE) + } + } + + x.position <- 1 # col + y.position <- 1 # row + # change loop: + while((x.position<=m) & (y.position<=n)){ + col <- as(A[,x.position], "sparseVector") + col[1:n < y.position] <- 0 + if(sum(abs(col)) < tol){ + pivot <- 0 + }else{ # find maximum pivot in current column at or below current row + which <- col@i[which.max(abs(col@x))] + pivot <- col[which] + } + if(abs(pivot) <= tol){ + x.position <- x.position+1 # check for 0 pivot + }else{ + if(which > y.position){ # exchange rows + id <- A@i + lp <- A@i + lp[id==which-1] <- y.position-1 + + lp[id==y.position-1] <- which-1 + A <- sparseMatrix(i = lp+1, p=A@p, x=A@x) + } + + A[y.position,]<-A[y.position,]/pivot # pivot + + row <-as(A[y.position,], "sparseVector") + A <- A - A[,x.position, drop = FALSE] %*% t(row) # sweep + A[y.position,]<-row # restore current row + x.position<-x.position+1 + y.position<-y.position+1 + } + + + if(verbose){ + if(min(n,m)==n){ + step <- y.position + }else{ + step <- x.position + } + printstep <- (step%%quant == 0) + if(any(printstep)){ + message(names(quant[printstep]), appendLF = FALSE) + quant <- quant[-1] + if(length(quant) > 0){ + message("...", appendLF = FALSE) + } + } + + } + } + + A <- drop0(A, tol = tol) + + if(!sparse){ + A <- as.matrix(A) + } + + return(A) +} diff --git a/R/tcsrec.R b/R/tcsrec.R index 8471990..b385168 100755 --- a/R/tcsrec.R +++ b/R/tcsrec.R @@ -1,36 +1,47 @@ #' @title Heuristic first-temporal-then-cross-sectional cross-temporal forecast reconciliation #' -#' @description The cross-temporal forecast reconciliation procedure by -#' Kourentzes and Athanasopoulos (2019) can be viewed as an ensemble forecasting -#' procedure which exploits the simple averaging of different forecasts. -#' First, for each time series the forecasts at any temporal aggregation order are -#' reconciled using temporal hierarchies (\code{\link[FoReco]{thfrec}}), then -#' time-by-time cross-sectional reconciliation is performed (\code{\link[FoReco]{htsrec}}). -#' The projection matrices obtained at this step are then averaged and used to +#' @description +#' \loadmathjax +#' The cross-temporal forecast reconciliation procedure by +#' Kourentzes and Athanasopoulos (2019) can be viewed as an ensemble +#' forecasting procedure which exploits the simple averaging of different +#' forecasts. First, for each time series the forecasts at any temporal +#' aggregation order are reconciled using temporal hierarchies +#' (\code{\link[FoReco]{thfrec}()}), then time-by-time cross-sectional +#' reconciliation is performed (\code{\link[FoReco]{htsrec}()}). The +#' projection matrices obtained at this step are then averaged and used to #' cross-sectionally reconcile the forecasts obtained at step 1, by this way #' fulfilling both cross-sectional and temporal constraints. #' -#' @param basef (\code{n x h(k* + m)}) matrix of base forecasts to be reconciled; -#' \code{n} is the total number of variables, \code{m} is the highest frequency, -#' \code{k*} is the sum of (\code{p-1}) factors of \code{m}, excluding \code{m}, -#' and \code{h} is the forecast horizon. Each row identifies, a time series, and the forecasts -#' are ordered as [lowest_freq' ... highest_freq']'. -#' @param hts_comb,thf_comb Type of covariance matrix (respectively (\code{n x n}) and -#' (\code{(k* + m) x (k* + m)})) to be used in the cross-sectional and temporal reconciliation, -#' see more in \code{comb} param of \code{\link[FoReco]{htsrec}} and \code{\link[FoReco]{thfrec}}. -#' @param res (\code{n x N(k* + m)}) matrix containing the residuals at all the -#' temporal frequencies ordered [lowest_freq' ... highest_freq']' (columns) for -#' each variable (row), needed to estimate the covariance matrix when \code{hts_comb =} -#' \code{\{"wls",} \code{"shr",} \code{"sam"\}} and/or \code{hts_comb =} \code{\{"wlsv",} -#' \code{"wlsh",} \code{"acov",} \code{"strar1",} \code{"sar1",} \code{"har1",} -#' \code{"shr",} \code{"sam"\}}. The row must be in the same order as \code{basef}. -#' @param avg If \code{avg = "KA"} (\emph{default}), the final projection matrix \code{M} is the one proposed -#' by Kourentzes and Athanasopoulos (2019), otherwise it is calculated as simple average of -#' all the involved projection matrices at step 2 of th procedure (see Di Fonzo and -#' Girolimetto, 2020). -#' @param ... any other options useful for \code{\link[FoReco]{htsrec}} and -#' \code{\link[FoReco]{thfrec}}, e.g. \code{m}, \code{C} (or \code{Ut} and \code{nb}), -#' \code{nn} (for non negativity reconciliation only at first step), ... +#' @param basef (\mjseqn{n \times h(k^\ast+m)}) matrix of base forecasts to be +#' reconciled, \mjseqn{\widehat{\mathbf{Y}}}; \mjseqn{n} is the total number of variables, +#' \mjseqn{m} is the highest time frequency, \mjseqn{k^\ast} is the sum of (a +#' subset of) (\mjseqn{p-1}) factors of \mjseqn{m}, excluding \mjseqn{m}, and +#' \mjseqn{h} is the forecast horizon for the lowest frequency time series. +#' Each row identifies a time series, and the forecasts are ordered as +#' [lowest_freq' ... highest_freq']'. +#' @param hts_comb,thf_comb Type of covariance matrix (respectively +#' (\mjseqn{n \times n}) and (\mjseqn{(k^\ast + m) \times (k^\ast + m)})) to +#' be used in the cross-sectional and temporal reconciliation, see more in +#' \code{comb} param of \code{\link[FoReco]{htsrec}()} and +#' \code{\link[FoReco]{thfrec}()}. +#' @param res (\mjseqn{n \times N(k^\ast + m)}) matrix containing the residuals +#' at all the temporal frequencies ordered [lowest_freq' ... highest_freq']' +#' (columns) for each variable (row), needed to estimate the covariance matrix +#' when \code{hts_comb =} \code{\{"wls",} \code{"shr",} \code{"sam"\}} and/or +#' \code{hts_comb =} \code{\{"wlsv",} \code{"wlsh",} \code{"acov",} +#' \code{"strar1",} \code{"sar1",} \code{"har1",} \code{"shr",} \code{"sam"\}}. +#' The row must be in the same order as \code{basef}. +#' @param avg If \code{avg = "KA"} (\emph{default}), the final projection +#' matrix \mjseqn{\mathbf{M}} is the one proposed by Kourentzes and +#' Athanasopoulos (2019), otherwise it is calculated as simple average of +#' all the involved projection matrices at step 2 of the procedure (see +#' Di Fonzo and Girolimetto, 2020). +#' @param ... any other options useful for \code{\link[FoReco]{htsrec}()} and +#' \code{\link[FoReco]{thfrec}()}, e.g. \code{m}, \code{C} (or \code{Ut} and +#' \code{nb}), \code{nn} (for non negativity reconciliation only at first step), +#' \code{mse}, \code{corpcor}, \code{type}, \code{sol}, \code{settings}, +#' \code{W}, \code{Omega},... #' #' @details #' This function performs a two-step cross-temporal forecast reconciliation using @@ -41,10 +52,11 @@ #' \strong{Warning}, #' the two-step heuristic reconciliation allows non negativity constraints only in the first step. #' This means that non-negativity is not guaranteed in the final reconciled values. +#' #' @return #' The function returns a list with two elements: -#' \item{\code{recf}}{(\code{n x h(k* + m)}) reconciled forecasts matrix.} -#' \item{\code{M}}{Projection matrix (projection approach).} +#' \item{\code{recf}}{(\mjseqn{n \times h(k^\ast + m)}) reconciled forecasts matrix, \mjseqn{\widetilde{\textbf{Y}}}.} +#' \item{\code{M}}{Matrix which transforms the uni-dimensional reconciled forecasts of step 1 (projection approach) .} #' #' @references #' Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: @@ -62,22 +74,25 @@ #' Matrix Estimation and Implications for Functional Genomics, \emph{Statistical #' Applications in Genetics and Molecular Biology}, 4, 1. #' -#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2018). OSQP: -#' An Operator Splitting Solver for Quadratic Programs, \href{https://arxiv.org/abs/1711.08013}{arXiv:1711.08013}. +#' Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S. (2020). OSQP: +#' An Operator Splitting Solver for Quadratic Programs, \emph{Mathematical Programming Computation}, +#' 12, 4, 637-672. #' #' Stellato, B., Banjac, G., Goulart, P., Boyd, S., Anderson, E. (2019), OSQP: -#' Quadratic Programming Solver using the 'OSQP' Library, R package version 0.6.0.3 +#' Quadratic Programming Solver using the `OSQP' Library, R package version 0.6.0.3 #' (October 10, 2019), \href{https://CRAN.R-project.org/package=osqp}{https://CRAN.R-project.org/package=osqp}. #' -#' @keywords reconciliation heuristic +#' @keywords heuristic +#' @family reconciliation procedures +#' #' @examples #' data(FoReco_data) -#' obj <- tcsrec(FoReco_data$base, thf_comb = "acov", hts_comb = "shr", -#' res = FoReco_data$res, m = 12, C = FoReco_data$C) -#' -#' @export +#' obj <- tcsrec(FoReco_data$base, m = 12, C = FoReco_data$C, +#' thf_comb = "acov", hts_comb = "shr", res = FoReco_data$res) #' #' @usage tcsrec(basef, thf_comb, hts_comb, res, avg = "KA", ...) +#' +#' @export tcsrec <- function(basef, thf_comb, hts_comb, res, avg = "KA", ...) { arg_input <- list(...) @@ -101,10 +116,9 @@ tcsrec <- function(basef, thf_comb, hts_comb, res, avg = "KA", ...) { tools <- thf_tools(m) kset <- tools$kset - m <- max(kset) + m <- tools$m kt <- tools$kt - arg_thf <- names(as.list(args(thfrec))) arg_thf <- arg_thf[!(arg_thf %in% c("basef", "keep", "res", "", "comb", "m", "bounds"))] @@ -128,7 +142,7 @@ tcsrec <- function(basef, thf_comb, hts_comb, res, avg = "KA", ...) { ## Step 2: compute time by time cross sectional M matrix arg_hts <- names(as.list(args(htsrec))) - arg_hts <- arg_hts[!(arg_hts %in% c("basef", "keep", "res", "", "comb", "nn", "bounds"))] + arg_hts <- arg_hts[!(arg_hts %in% c("basef", "keep", "res", "", "comb", "nn", "bounds", "nn_type"))] # Create list with lenght p, with time by time temporally reconciled forecasts matrices h <- NCOL(Y1) / kt @@ -138,12 +152,7 @@ tcsrec <- function(basef, thf_comb, hts_comb, res, avg = "KA", ...) { M <- lapply(Y, function(x){ obj <- do.call("htsrec", c(list(basef = t(x), comb = hts_comb), arg_input[which(names(arg_input) %in% arg_hts)])) - out <- obj[["M"]] - if(is.null(out)){ - return(obj[["G"]]) - }else{ - return(out) - } + return(obj[["M"]]) }) } else { # Create list with lenght p, with time by time temporally reconciled residuals matrices @@ -154,21 +163,18 @@ tcsrec <- function(basef, thf_comb, hts_comb, res, avg = "KA", ...) { M <- mapply(function(Y, E){ obj <- do.call("htsrec", c(list(basef = t(Y), comb = hts_comb, res = t(E)), arg_input[which(names(arg_input) %in% arg_hts)])) - out <- obj[["M"]] - if(is.null(out)){ - return(obj[["G"]]) - }else{ - return(out) - } + return(obj[["M"]]) }, Y = Y, E = E, SIMPLIFY = FALSE ) } if (avg == "KA") { - meanM <- apply(simplify2array(lapply(M, as.matrix)), c(1, 2), sum) / length(M) + #meanM <- apply(simplify2array(lapply(M, as.matrix)), c(1, 2), sum) / length(M) + meanM <- Reduce("+", M)/length(M) } else { Mw <- mapply(function(a, A) a * A, A = M, a = split(kset, 1:length(kset)), SIMPLIFY = FALSE) - meanM <- apply(simplify2array(lapply(Mw, as.matrix)), c(1, 2), sum) / kt + #meanM <- apply(simplify2array(lapply(Mw, as.matrix)), c(1, 2), sum) / kt + meanM <- Reduce("+", Mw)/kt } ## Step 3: Cross-Temporal reconciled forecasts with heuristic diff --git a/R/tdrec.R b/R/tdrec.R new file mode 100644 index 0000000..8a7e5b4 --- /dev/null +++ b/R/tdrec.R @@ -0,0 +1,170 @@ +#' Top-down forecast reconciliation for genuine hierarchical/grouped time series +#' +#' @description +#' \loadmathjax +#' Top-down forecast reconciliation for genuine hierarchical/grouped time series, +#' where the forecast of a `Total' (top-level series, expected to be positive) +#' is disaggregated according to a proportional scheme given by a vector +#' of proportions (weights). +#' Besides the fulfillment of any aggregation constraint, +#' the top-down reconciled forecasts should respect two main properties: +#' - the top-level value remains unchanged; +#' - all the bottom time series reconciled forecasts are non-negative. +#' The top-down procedure is extended to deal with both temporal and cross-temporal cases. +#' Since this is a post forecasting function, the weight vector must be given +#' in input by the user, and is not calculated automatically (see Examples). +#' +#' @param topf (\mjseqn{h \times 1}) vector of the top-level base forecast to be +#' disaggregated; \mjseqn{h} is the forecast horizon (for the lowest temporal +#' aggregation order in temporal and cross-temporal cases). +#' @param C (\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. +#' @param m Highest available sampling frequency per seasonal cycle (max. order +#' of temporal aggregation, \mjseqn{m}), or a subset of the \mjseqn{p} factors +#' of \mjseqn{m}. +#' @param weights vector of weights to be used to disaggregate topf: +#' (\mjseqn{n_b \times h}) matrix in the cross-sectional framework; +#' (\mjseqn{m \times h}) matrix in the temporal framework; +#' (\mjseqn{n_b m \times h}) matrix in the cross-temporal framework. +#' +#' @details Fix \mjseqn{h = 1}, then +#' \mjsdeqn{\widetilde{\mathbf{y}} = \mathbf{S}\mathbf{w}\widehat{a}_1} +#' where \mjseqn{\widetilde{\mathbf{y}}} is the vector of reconciled forecasts, +#' \mjseqn{\mathbf{S}} is the summing matrix (whose pattern depends on which type +#' of reconciliation is being performed), \mjseqn{\mathbf{w}} is the vector of weights, +#' and \mjseqn{\widehat{a}_1} is the top-level value to be disaggregated. +#' +#' @return The function returns an (\mjseqn{h \times n}) matrix of +#' cross-sectionally reconciled forecasts, or an (\mjseqn{h(k^\ast + m) \times 1}) +#' vector of top-down temporally reconciled forecasts, or an +#' (\mjseqn{n \times h (k^\ast + m)}) matrix of top-down +#' cross-temporally reconciled forecasts. +#' +#' @references +#' Athanasopoulos, G., Ahmed, R.A., Hyndman, R.J. (2009), Hierarchical +#' forecasts for Australian domestic tourism, \emph{International Journal of +#' Forecasting}, 25, 1, 146–166. +#' +#' @keywords top-down +#' @family reconciliation procedures +#' +#' @examples +#' data(FoReco_data) +#' ### CROSS-SECTIONAL TOP-DOWN RECONCILIATION +#' # Cross sectional aggregation matrix +#' C <- FoReco_data$C +#' # monthly base forecasts +#' id <- which(simplify2array(strsplit(colnames(FoReco_data$base), split = "_"))[1, ] == "k1") +#' mbase <- t(FoReco_data$base[, id]) +#' obs_1 <- FoReco_data$obs$k1 +#' # average historical proportions +#' props <- colMeans(obs_1[1:168,-c(1:3)]/obs_1[1:168,1]) +#' cs_td <- tdrec(topf = mbase[,1], C = C, weights = props) +#' +#' ### TEMPORAL TOP-DOWN RECONCILIATION +#' # top ts base forecasts ([lowest_freq' ... highest_freq']') +#' top_obs12 <- FoReco_data$obs$k12[1:14,1] +#' bts_obs1 <- FoReco_data$obs$k1[1:168,1] +#' # average historical proportions +#' props <- colMeans(matrix(bts_obs1, ncol = 12, byrow = TRUE)/top_obs12) +#' topbase <- FoReco_data$base[1, 1] +#' t_td <- tdrec(topf = topbase, m = 12, weights = props) +#' +#' ### CROSS-TEMPORAL TOP-DOWN RECONCILIATION +#' top_obs <- FoReco_data$obs$k12[1:14,1] +#' bts_obs <- FoReco_data$obs$k1[1:168,-c(1:3)] +#' bts_obs <- lapply(1:5, function(x) matrix(bts_obs[,x], nrow=14, byrow = TRUE)) +#' bts_obs <- do.call(cbind, bts_obs) +#' # average historical proportions +#' props <- colMeans(bts_obs/top_obs) +#' ct_td <- tdrec(topf = topbase, m = 12, C = C, weights = props) +#' +#' @export +#' +#' @import Matrix +tdrec <- function(topf, C, m, weights){ + if(missing(C) & missing(m)){ + stop("Missing C or/and m. Remember:\n - put only C for a cross-sectional top-down reconciliation\n - put only m for a temporal top-down reconciliation\n - put C AND m for a cross-temporal top-down reconciliation.") + }else if(missing(C)){ + obj_thf <- thf_tools(m = m) + m <- obj_thf$m + S <- obj_thf$R + h <- length(topf) + kset <- obj_thf$kset + + Dh <- Dmat(h = h, m = kset, n = 1) + vnames <- paste("k", rep(kset, h * rev(kset)), "h", + do.call("c", as.list(sapply( + rev(kset) * h, + function(x) seq(1:x)))), + sep = "") + }else if(missing(m)){ + S <- hts_tools(C = C)$S + + cnames <- if (is.null(rownames(C)) | is.null(colnames(C))) { + paste("serie", 1:(NROW(C)+NCOL(C)), sep = "") + } else { + c(rownames(C), colnames(C)) + } + rnames <- paste("h", 1:length(topf), sep="") + }else{ + obj_ctf <- ctf_tools(C = C, m = m, sparse = TRUE) + nb <- obj_ctf$hts$nb + na <- obj_ctf$hts$na + kt <- obj_ctf$thf$kt + ks <- obj_ctf$thf$ks + m <- obj_ctf$thf$m + h <- length(topf) + kset <- obj_ctf$thf$kset + S <- obj_ctf$ctf$Stilde + + Dh <- Dmat(h = h, m = kset, n = na+nb) + rnames <- if (is.null(rownames(C)) | is.null(colnames(C))) { + paste("serie", 1:(na+nb), sep = "") + } else { + c(rownames(C), colnames(C)) + } + cnames <- paste("k", rep(kset, h * rev(kset)), "h", + do.call("c", as.list(sapply( + rev(kset) * h, + function(x) seq(1:x)))), + sep = "") + } + + if(is.vector(weights) | NCOL(weights) == 1 | NROW(weights) == 1){ + if(NCOL(weights) == 1 | NROW(weights) == 1){ + weights <- as.vector(weights) + } + + if(NROW(weights) != NCOL(S)){ + stop(paste0("The weights vector must have ", NCOL(S)," elements")) + } + + p <- weights/sum(weights) + recf <- t(S %*% p %*% t(topf)) + }else{ + if(NROW(weights) != NCOL(S)){ + stop(paste0("The weights matrix must have ", NCOL(S)," rows")) + } + + if(any(abs(colSums(weights)-1) 0) } + OUTF <- BASEF %*% t(R) outf <- as.vector(t(Dh) %*% as.vector(t(OUTF))) @@ -264,7 +367,7 @@ thfrec.default <- function(basef, m, comb, res, mse = TRUE, corpcor = FALSE, )) if (keep == "list") { return(list( - recf = outf, S = R, + recf = outf, R = R, M = R %*% cbind(matrix(0, m, ks), diag(m)) )) } else { @@ -334,12 +437,12 @@ thfrec.default <- function(basef, m, comb, res, mse = TRUE, corpcor = FALSE, if (type == "S") { rec_sol <- recoS( basef = BASEF, W = Omega, S = R, sol = sol, nn = nn, keep = keep, - settings = settings, b_pos = b_pos, bounds = bounds + settings = settings, b_pos = b_pos, bounds = bounds, nn_type = nn_type ) } else { rec_sol <- recoM( - basef = BASEF, W = Omega, H = Zt, sol = sol, nn = nn, keep = keep, - settings = settings, b_pos = b_pos, bounds = bounds + basef = BASEF, W = Omega, Ht = Zt, sol = sol, nn = nn, keep = keep, S = R, + settings = settings, b_pos = b_pos, bounds = bounds, nn_type = nn_type ) } @@ -403,7 +506,7 @@ thfrec.list <- function(basef, m, ..., Omega){ stop("You don't need thfrech for h = 1, please use thfrec()", call. = FALSE) } - Dh <- Dmat(h = h, kset = tools$kset, n = 1) + Dh <- Dmat(h = h, m = tools$kset, n = 1) baseh <- matrix(Dh%*%basef, nrow = h, byrow = T) if(h != length(Omega) | !is.list(Omega)){ diff --git a/R/ut2c.R b/R/ut2c.R new file mode 100644 index 0000000..a241332 --- /dev/null +++ b/R/ut2c.R @@ -0,0 +1,169 @@ +#' Cross-sectional 'structural representation' of a general linearly constrained +#' multiple time series +#' +#' @description +#' \loadmathjax +#' Switching from a zero-constraints kernel representation of a linearly constrained +#' multiple time series, \mjseqn{\mathbf{U}'\mathbf{y}=\mathbf{0}}, +#' to a 'structural representation', \mjseqn{\mathbf{y}=\mathbf{S}\mathbf{b}}. +#' (\emph{Experimental version}) +#' +#' @param Ut (\mjseqn{r \times c}) zero constraints cross-sectional +#' (contemporaneous) kernel matrix \mjseqn{(\textbf{U}'\textbf{y} = \mathbf{0})} +#' spanning the null space valid for the target forecasts. +#' @param sparse Option to return sparse object (\emph{default} is \code{TRUE}). +#' @param verbose If \code{TRUE}, print intermediate steps of \code{\link{srref}}. +#' +#' @details +#' Consider the simple example of linearly constrained multiple time series consisting +#' of two hierarchies, each with distinct bottom time series, +#' with a common top-level series (\mjseqn{T}): +#' \mjsdeqn{\begin{array}{ll} +#' 1)\; T = A + B & 4)\; T = X + Y \cr +#' 2)\; A = AA + AB & 5)\; X = XX + XY \cr +#' 3)\; B = BA + BB + BC & 6)\; Y = YX + YY +#' \end{array}.} +#' Given the cross-sectional aggregation matrices of each hierarchy, +#' \mjsdeqn{\mathbf{C}_1 = \left[\begin{array}{ccccc} +#' 1 & 1 & 0 & 0 & 0\cr +#' 0 & 0 & 1 & 1 & 1 +#' \end{array}\right]\quad \mathrm{and} \quad \mathbf{C}_2 = \left[\begin{array}{ccccc} +#' 1 & 1 & 0 & 0\cr +#' 0 & 0 & 1 & 1 +#' \end{array}\right],} +#' the zero constraints cross-sectional kernel matrix \mjseqn{\mathbf{U}'}, +#' which accounts for the constraints of both hierarchies, +#' can be built as follows: +#' \mjsdeqn{\footnotesize\mathbf{U}' = \left[\begin{array}{cccccccccccccc|c} +#' 1 & 0 & 0 &-1 &-1 &-1 &-1 &-1 & 0 & 0 & 0 & 0 & 0 & 0 & T\cr +#' 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &-1 &-1 &-1 &-1 & T\cr +#' 0 & 1 & 0 &-1 &-1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & A\cr +#' 0 & 0 & 1 & 0 & 0 &-1 &-1 &-1 & 0 & 0 & 0 & 0 & 0 & 0 & B\cr +#' 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 &-1 &-1 & 0 & 0 & X\cr +#' 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 &-1 &-1 & Y\cr +#' \hline +#' T & A & B & AA& AB& BA& BB& BC& X & Y & XX& XY& YX& YY& +#' \end{array}\right].} +#' Function \code{\link{ut2c}} returns a matrix +#' \mjseqn{\footnotesize\mathbf{U}'_{rref} = \left[\mathbf{I}_6 \; -\mathbf{C} \right]}: +#' \mjsdeqn{\footnotesize\mathbf{U}'_{rref} = \left[\begin{array}{cccccccccccccc|c} +#' 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &-1 &-1 &-1 &-1 & T\cr +#' 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 &-1 &-1 &-1 &-1 & A\cr +#' 0 & 0 & 1 & 0 & 0 & 0 & 0 &-1 &-1 &-1 & 0 & 0 & 0 & 0 & B\cr +#' 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 &-1 &-1 & 0 & 0 & X\cr +#' 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 &-1 &-1 & Y\cr +#' 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 &-1 &-1 &-1 &-1 & AA\cr +#' \hline +#' T & A & B & X & Y & AA& AB& BA& BB& BC& XX& XY& YX& YY& +#' \end{array}\right],} +#' from which the 'structural sum matrix' +#' \mjsdeqn{\mathbf{S} = \left[\mathbf{C}' \; \mathbf{I}_8 \right]'} +#' may be easily obtained. +#' @family utilities +#' +#' @return A list with +#' \item{\code{Ut_reshape}}{matrix with rows and columns rearranged so that +#' each row has 1 as the first non-null element starting from the left, +#' and that the first r columns have 1 as the first element starting from +#' the bottom, \mjseqn{\mathbf{U}_{reshape}' = \mathbf{U}'[,\mbox{cid}]}.} +#' \item{\code{Ut_rref}}{reduced row echelon form of \mjseqn{\mathbf{U}'_{reshape}} with +#' \link{srref}.} +#' \item{\code{C}}{(\mjseqn{n_a \times n_b}) cross-sectional (contemporaneous) matrix +#' mapping the bottom level series into the higher level ones, +#' \mjseqn{\mathbf{U}_{srref}' = [\mathbf{I} \ -\mathbf{C}]}.} +#' \item{\code{cid}}{(\mjseqn{c \times 1}) vector of the column permutations of the matrix \mjseqn{\mathbf{U}'}.} +#' \item{\code{nb}}{number of bottom time series, \mjseqn{n_b}.} +#' \item{\code{na}}{number of upper time series, \mjseqn{n_a = r}.} +#' \item{\code{n}}{number of time series, \mjseqn{n_a + n_b = n = c}.} +#' +#' @references +#' Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: +#' Optimal Combination Method and Heuristic Alternatives, Department of Statistical +#' Sciences, University of Padua, \href{https://arxiv.org/abs/2006.08570}{arXiv:2006.08570}. +#' +#' @export +#' +#' @examples +#' C1 <- Matrix(c(1,1,0,0,0, +#' 0,0,1,1,1), 2, byrow = TRUE, sparse = TRUE) +#' C2 <- Matrix(c(1,1,0,0, +#' 0,0,1,1), 2, byrow = TRUE, sparse = TRUE) +#' Ut <- rbind(c(1, rep(0, NROW(C1)), rep(-1, NCOL(C1)), rep(0, NROW(C2)), rep(0, NCOL(C2))), +#' c(1, rep(0, NROW(C1)), rep(0, NCOL(C1)), rep(0, NROW(C2)), rep(-1, NCOL(C2))), +#' cbind(0, Diagonal(NROW(C1)), -C1, Matrix(0, NROW(C1), NROW(C2)+NCOL(C2))), +#' cbind(0, Matrix(0, NROW(C2), NROW(C1)+NCOL(C1)), Diagonal(NROW(C2)), -C2)) +#' colnames(Ut) <- c("T", "A", "B", "AA", "AB", "BA", "BB", "BC", +#' "X", "Y", "XX", "XY", "YX", "YY") +#' rownames(Ut) <- c("T", "T", "A", "B", "X", "Y") +#' obj_ut2c <- ut2c(Ut, sparse = FALSE) +#' Ut_struc <- obj_ut2c$Ut_rref +#' S <- rbind(obj_ut2c$C, diag(1, obj_ut2c$nb)) +#' +ut2c <- function(Ut, sparse = TRUE, verbose = FALSE){ + if(!is(Ut, "dgCMatrix")){ + Ut <- as(Ut, "dgCMatrix") + } + + obj <- reshape_mat(mat = Ut) + Ut_reshape <- obj$mat + cid <- obj$cid + + out <- list() + out$Ut_reshape <- Ut_reshape + + if(!isDiagonal(Ut_reshape[,c(1:NROW(Ut_reshape))]) | !all(Ut_reshape@x[1:NROW(Ut_reshape)]==1)){ + Ut_rref <- srref(A = Ut_reshape, sparse = TRUE, verbose = verbose) + + rs <- rowSums(abs(Ut_rref))==0 + if(any(rs)){ + message("rref(Ut) has ", sum(rs), " zeros rows. These rows have been removed") + Ut_rref <- Ut_rref[!rs,,drop=FALSE] + } + + colnames(Ut_rref) <- colnames(Ut_reshape) + rownames(Ut_rref) <- colnames(Ut_reshape)[1:NROW(Ut_rref)] + out$Ut_rref <- Ut_rref + + out$C <- -Ut_rref[,-c(1:NROW(Ut_rref)),drop=FALSE] + colnames(out$C) <- colnames(Ut_reshape)[-c(1:NROW(Ut_rref))] + rownames(out$C) <- colnames(Ut_reshape)[1:NROW(Ut_rref)] + if(!sparse){ + out$Ut_reshape <- as.matrix(out$Ut_reshape) + out$Ut_rref <- as.matrix(out$Ut_rref) + out$C <- as.matrix(out$C) + } + }else{ + out$C <- -Ut_reshape[,-c(1:NROW(Ut_reshape)),drop=FALSE] + colnames(out$C) <- colnames(Ut_reshape)[-c(1:NROW(Ut_reshape))] + rownames(out$C) <- colnames(Ut_reshape)[1:NROW(Ut_reshape)] + + if(!sparse){ + out$Ut_reshape <- as.matrix(out$Ut_reshape) + out$C <- as.matrix(out$C) + } + } + + names(cid) <- colnames(Ut)[cid] + out$cid <- cid + out$nb <- NCOL(out$C) + out$na <- NROW(out$C) + out$n <- NCOL(out$C)+NROW(out$C) + return(out) +} + +reshape_mat <- function(mat){ + if(!is(mat, "dgCMatrix")){ + mat <- as(mat, "dgCMatrix") + } + r <- mat@i+1 + c <- findInterval(seq(mat@x)-1,mat@p[-1])+1 + c_new <- c[which(mat@x==1)][!duplicated(r[which(mat@x==1)])] + id <- c(unique(c_new), c(1:NCOL(mat))[!c(1:NCOL(mat)) %in% unique(c_new)]) + #r_new <- r[which(mat@x==1)][!duplicated(r[which(mat@x==1)])] # not now + out <- list() + out$mat <- mat[, id] + out$cid <- id + return(out) +} + + diff --git a/README.Rmd b/README.Rmd index da95fa4..8c3ab85 100755 --- a/README.Rmd +++ b/README.Rmd @@ -20,28 +20,31 @@ knitr::opts_chunk$set( [![R build status](https://github.com/daniGiro/FoReco/workflows/R-CMD-check/badge.svg)](https://github.com/daniGiro/FoReco/actions) [![CRAN status](https://www.r-pkg.org/badges/version/FoReco)](https://CRAN.R-project.org/package=FoReco) [![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental) -[![devel version](https://img.shields.io/badge/devel%20version-0.1.2-blue.svg)](https://github.com/daniGiro/FoReco) +[![devel version](https://img.shields.io/badge/devel%20version-0.2.0-blue.svg)](https://github.com/daniGiro/FoReco) [![License: GPL-3](https://img.shields.io/badge/license-GPL--3-forestgreen.svg)](https://cran.r-project.org/web/licenses/GPL-3) -The **FoReco** (**Fo**recast **Reco**nciliation) package is designed for point forecast reconciliation, a **post-forecasting** process aimed to improve the quality of the base forecasts for a system of linearly constrained (e.g. hierarchical/grouped) time series. +The **FoReco** (**Fo**recast **Reco**nciliation) package is designed for point forecast reconciliation, a **post-forecasting** process aimed to improve the accuracy of the base forecasts for a system of linearly constrained (e.g. hierarchical/grouped) time series. -It offers classical (bottom-up), optimal and heuristic combination forecast reconciliation procedures by exploiting cross-sectional, temporal, and cross-temporal relationships linking the time series. +It offers 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. The main functions are: * `htsrec()`: cross-sectional (contemporaneous) forecast reconciliation. * `thfrec()`: forecast reconciliation for a single time series through temporal hierarchies. +* `lccrec()`: level conditional forecast reconciliation for genuine hierarchical/grouped time series. +* `tdrec()`: top-down (cross-sectional, temporal, cross-temporal) forecast reconciliation for genuine hierarchical/grouped time series. +* `ctbu()`: bottom-up cross-temporal forecast reconciliation. * `tcsrec()`: heuristic first-temporal-then-cross-sectional cross-temporal forecast reconciliation. * `cstrec()`: heuristic first-cross-sectional-then-temporal cross-temporal forecast reconciliation. * `iterec()`: heuristic iterative cross-temporal forecast reconciliation. * `octrec()`: optimal combination cross-temporal forecast reconciliation. - + ## Installation You can install the **stable** version on [R CRAN](https://cran.r-project.org/) ```r -install.packages('FoReco', dependencies = TRUE) +install.packages("FoReco") ``` You can also install the **development** version from @@ -49,7 +52,7 @@ You can also install the **development** version from ```r # install.packages("devtools") -devtools::install_github("daniGiro/FoReco", build_vignettes = TRUE) +devtools::install_github("daniGiro/FoReco") ``` ## Links * Source code: https://github.com/daniGiro/FoReco diff --git a/README.md b/README.md index 4680a7f..d5aad8f 100755 --- a/README.md +++ b/README.md @@ -1,8 +1,7 @@ -FoReco logo -================================================================================================================== +# FoReco logo @@ -13,19 +12,20 @@ status](https://www.r-pkg.org/badges/version/FoReco)](https://CRAN.R-project.org [![Lifecycle: experimental](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://www.tidyverse.org/lifecycle/#experimental) [![devel -version](https://img.shields.io/badge/devel%20version-0.1.2-blue.svg)](https://github.com/daniGiro/FoReco) +version](https://img.shields.io/badge/devel%20version-0.2.0-blue.svg)](https://github.com/daniGiro/FoReco) [![License: GPL-3](https://img.shields.io/badge/license-GPL--3-forestgreen.svg)](https://cran.r-project.org/web/licenses/GPL-3) The **FoReco** (**Fo**recast **Reco**nciliation) package is designed for point forecast reconciliation, a **post-forecasting** process aimed to -improve the quality of the base forecasts for a system of linearly +improve the accuracy of the base forecasts for a system of linearly constrained (e.g. hierarchical/grouped) time series. -It offers classical (bottom-up), optimal and heuristic combination -forecast reconciliation procedures by exploiting cross-sectional, -temporal, and cross-temporal relationships linking the time series. +It offers 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. The main functions are: @@ -33,6 +33,12 @@ The main functions are: reconciliation. - `thfrec()`: forecast reconciliation for a single time series through temporal hierarchies. +- `lccrec()`: level conditional forecast reconciliation for genuine + hierarchical/grouped time series. +- `tdrec()`: top-down (cross-sectional, temporal, cross-temporal) + forecast reconciliation for genuine hierarchical/grouped time + series. +- `ctbu()`: bottom-up cross-temporal forecast reconciliation. - `tcsrec()`: heuristic first-temporal-then-cross-sectional cross-temporal forecast reconciliation. - `cstrec()`: heuristic first-cross-sectional-then-temporal @@ -42,14 +48,13 @@ The main functions are: - `octrec()`: optimal combination cross-temporal forecast reconciliation. -Installation ------------- +## Installation You can install the **stable** version on [R CRAN](https://cran.r-project.org/) ``` r -install.packages('FoReco', dependencies = TRUE) +install.packages("FoReco") ``` You can also install the **development** version from @@ -57,17 +62,15 @@ You can also install the **development** version from ``` r # install.packages("devtools") -devtools::install_github("daniGiro/FoReco", build_vignettes = TRUE) +devtools::install_github("daniGiro/FoReco") ``` -Links ------ +## Links - Source code: - Site documentation: -Getting help ------------- +## Getting help If you encounter a clear bug, please file a minimal reproducible example on [GitHub](https://github.com/daniGiro/FoReco/issues). diff --git a/_pkgdown.yml b/_pkgdown.yml index 47a9aa4..d6d216d 100755 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -56,12 +56,15 @@ reference: - cstrec - iterec - ctbu + - lccrec + - tdrec - title: Datasets - contents: - FoReco_data - FoReco-hts - FoReco-thief + - title: Utilities - contents: - commat @@ -72,6 +75,9 @@ reference: - hts_tools - thf_tools - ctf_tools + - srref + - ut2c + - oct_bounds - title: Package - contents: - FoReco-package diff --git a/cran-comments.md b/cran-comments.md new file mode 100644 index 0000000..44fef2d --- /dev/null +++ b/cran-comments.md @@ -0,0 +1,11 @@ +## Test environments +* local R installation, R 4.0.5 +* ubuntu 16.04 (on travis-ci), R 4.0.5 +* win-builder (devel) + +## R CMD check results + +0 errors | 0 warnings | 1 note + +* This is a new release of FoReco 0.2. +* Changelog: https://danigiro.github.io/FoReco/news/index.html diff --git a/docs/404.html b/docs/404.html index 6191dca..cd5eedf 100755 --- a/docs/404.html +++ b/docs/404.html @@ -85,7 +85,7 @@ FoReco - 0.1.2 + 0.2.0 @@ -93,21 +93,21 @@