Cross-sectional (contemporaneous) forecast reconciliation of hierarchical and grouped 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.
htsrec(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FALSE, W,
type = "M", sol = "direct", nn = FALSE, keep = "list",
settings = osqpSettings(verbose = FALSE, eps_abs = 1e-5,
eps_rel = 1e-5, polish_refine_iter = 100, polish=TRUE))(h x n) matrix of base forecasts to be reconciled;
h is the forecast horizon and n is the total number of time series.
Type of the reconciliation. Except for Bottom-up, each
option corresponds to a specified (n x n) covariance matrix:
bu (Bottom-up);
ols (Identity);
struc (Structural variances);
wls (Series variances) - uses res;
shr (Shrunk covariance matrix - MinT-shr) - uses res;
sam (Sample covariance matrix - MinT-sam) - uses res;
w use your personal matrix W in param W.
(na x nb) cross-sectional (contemporaneous) matrix mapping the bottom
level series into the higher level ones.
(N x n) in-sample residuals matrix needed when comb =
{"wls", "shr", "sam"}. The columns must be in
the same order as basef.
Zero constraints cross-sectional (contemporaneous) kernel matrix
\((\textbf{U}'\textbf{Y} = \mathbf{0})\) spanning the null space valid for the reconciled
forecasts. It can be used instead of parameter C, but in this case nb (n = na + nb) is needed. If
the hierarchy admits a structural representation, Ut has dimension (na x n).
Number of bottom time series; if C is present, nb is not used.
Logical value: TRUE (default) calculates the
covariance matrix of the in-sample residuals (when necessary) according to the original
hts and thief formulation: no mean correction, T as denominator.
Logical value: TRUE if corpcor (Sch<U+00E4>fer et
al., 2017) must be used to shrink the sample covariance matrix according to
Sch<U+00E4>fer and Strimmer (2005), otherwise the function uses the same
implementation as package hts.
This option permits to directly enter the covariance matrix:
W must be a p.d. (n x n) matrix;
if comb is different from "w", W is not used.
Approach used to compute the reconciled forecasts: "M" for
the projection approach with matrix M (default), or "S" for the
structural approach with summing matrix S.
Solution technique for the reconciliation problem: either "direct" (default) for the direct
solution or "osqp" for the numerical solution (solving a linearly constrained quadratic
program using solve_osqp).
Logical value: TRUE if non-negative reconciled forecasts are wished.
Return a list object of the reconciled forecasts at all levels.
Settings for osqp (object osqpSettings). The default options
are: verbose = FALSE, eps_abs = 1e-5, eps_rel = 1e-5,
polish_refine_iter = 100 and polish = TRUE. For details, see the
osqp documentation (Stellato et al., 2019).
If the parameter keep is equal to "recf", then the function
returns only the (h x n) reconciled forecasts matrix, otherwise (keep="all")
it returns a list that mainly depends on what type of representation (type)
and methodology (sol) have been used:
recf(h x n) reconciled forecasts matrix.
WCovariance matrix used for forecast reconciliation.
nn_checkNumber of negative values (if zero, there are no values below zero).
rec_checkLogical value: has the hierarchy been respected?
M (type="M" and type="direct")Projection matrix (projection approach)
G (type="S" and type="direct")Projection matrix (structural approach).
S (type="S" and type="direct")Cross-sectional summing matrix, \(\textbf{S}\).
info (type="osqp")matrix with information in columns
for each forecast horizon h (rows): run time (run_time) number of iteration,
norm of primal residual (pri_res), status of osqp's solution (status) and
polish's status (status_polish).
Only if comb = "bu", the function returns recf, S and M.
In case of non-negativity constraints, there are two ways:
sol = "direct" and nn = TRUE: the base forecasts
will be reconciled at first without non-negativity constraints, then, if negative reconciled
values are present, the "osqp" solver is used.
sol = "osqp" and nn = TRUE: the base forecasts will be
reconciled through the "osqp" solver.
Byron, R.P. (1978), The estimation of large social accounts matrices, Journal of the Royal Statistical Society A, 141, 3, 359-367.
Di Fonzo, T., Girolimetto, D. (2020), Cross-Temporal Forecast Reconciliation: Optimal Combination Method and Heuristic Alternatives, Department of Statistical Sciences, University of Padua, arXiv:2006.08570.
Hyndman, R.J., Ahmed, R.A., Athanasopoulos, G., Shang, H.L.(2011), Optimal combination forecasts for hierarchical time series, Computational Statistics & Data Analysis, 55, 9, 2579-2589.
Sch<U+00E4>fer, J.L., Opgen-Rhein, R., Zuber, V., Ahdesmaki, M., Duarte Silva, A.P., Strimmer, K. (2017), Package `corpcor', R package version 1.6.9 (April 1, 2017), https://CRAN.R-project.org/package= corpcor.
Sch<U+00E4>fer, J.L., Strimmer, K. (2005), A Shrinkage Approach to Large-Scale Covariance Matrix Estimation and Implications for Functional Genomics, 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, arXiv:1711.08013.
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), https://CRAN.R-project.org/package=osqp.
van Erven, T., Cugliari, J. (2015), Game-theoretically Optimal Reconciliation of Contemporaneous Hierarchical Time Series Forecasts, in Antoniadis, A., Poggi, J.M., Brossat, X. (eds.), Modeling and Stochastic Learning for Forecasting in High Dimensions, Berlin, Springer, 297-317.
Wickramasuriya, S.L., Athanasopoulos, G., Hyndman, R.J. (2019), Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization, Journal of the American Statistical Association, 114, 526, 804-819.
# NOT RUN {
data(FoReco_data)
# monthly base forecasts
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")
mres <- t(FoReco_data$res[, id])
obj <- htsrec(mbase, C = FoReco_data$C, comb = "shr", res = mres)
# }
Run the code above in your browser using DataLab