Forecast reconciliation through temporal hierarchies (temporal reconciliation)

Forecast reconciliation of one time series through temporal hierarchies (Athanasopoulos et al., 2017). The reconciled forecasts are calculated either through a projection approach (Byron, 1978), or the equivalent structural approach by Hyndman et al. (2011). Moreover, the classic bottom-up approach is available.

thfrec(basef, m, comb, res, mse = TRUE, corpcor = FALSE,
       type = "M", sol = "direct", keep = "list", v = NULL, nn = FALSE,
        nn_type = "osqp", settings = list(), bounds = NULL, Omega = NULL)

Arguments

basef

(\(h(k^\ast + m) \times 1\)) vector of base forecasts to be reconciled, containing the forecasts at all the needed temporal frequencies ordered as [lowest_freq' ... highest_freq']'.

m

Highest available sampling frequency per seasonal cycle (max. order of temporal aggregation, \(m\)), or a subset of \(p\) factors of \(m\).

comb

Type of the reconciliation. Except for bottom up, all other options correspond to a different (\((k^\ast + m) \times (k^\ast + m)\)) covariance matrix, \(k^\ast\) is the sum of (\(p-1\)) factors of \(m\) (excluding \(m\)):

bu (Bottom-up);
ols (Identity);
struc (Structural variances);
wlsv (Series variances);
wlsh (Hierarchy variances);
acov (Auto-covariance matrix);
strar1 (Structural Markov);
sar1 (Series Markov);
har1 (Hierarchy Markov);
shr (Shrunk cross-covariance matrix);
sam (Sample cross-covariance matrix);
omega use your personal matrix Omega in param Omega.

res

vector containing the in-sample residuals at all the temporal frequencies ordered as basef, i.e. [lowest_freq' ... highest_freq']', needed to estimate the covariance matrix when comb = {"wlsv", "wlsh", "acov", "strar1", "sar1", "har1", "shr", "sam"}.

mse

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.

corpcor

Logical value: TRUE if corpcor (Schäfer et al., 2017) must be used to shrink the sample covariance matrix according to Schäfer and Strimmer (2005), otherwise the function uses the same implementation as package hts.

type

Approach used to compute the reconciled forecasts: "M" for the projection approach with matrix M (default), or "S" for the structural approach with temporal summing matrix R.

sol

Solution technique for the reconciliation problem: either "direct" (default) for the closed-form matrix solution, or "osqp" for the numerical solution (solving a linearly constrained quadratic program using solve_osqp).

keep

Return a list object of the reconciled forecasts at all levels (if keep = "list") or only the reconciled forecasts matrix (if keep = "recf").

v

vector index of the fixed base forecast (\(\mbox{min}(v) > 0\) and \(\mbox{max}(v)

nn

Logical value: TRUE if non-negative reconciled forecasts are wished.

nn_type

"osqp" (default), "KAnn" (only type == "M") or "sntz".

settings

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).

bounds

(\((k^\ast + m) \times 2\)) matrix with temporal bounds: the first column is the lower bound, and the second column is the upper bound.

Omega

This option permits to directly enter the covariance matrix:

Omega must be a p.d. (\((k^\ast + m) \times (k^\ast + m)\)) matrix or a list of \(h\) matrix (one for each forecast horizon);
if comb is different from "omega", Omega is not used.

Value

If the parameter keep is equal to "recf", then the function returns only the (\(h(k^\ast + m) \times 1\)) reconciled forecasts vector, otherwise (keep="all") it returns a list that mainly depends on what type of representation (type) and solution technique (sol) have been used:

recf: (\(h(k^\ast + m) \times 1\)) reconciled forecasts vector, \(\widetilde{\mathbf{y}}\).
Omega: Covariance matrix used for forecast reconciliation, \(\mathbf{\Omega}\).
nn_check: Number of negative values (if zero, there are no values below zero).
rec_check: Logical value: rec_check = TRUE when the constraints have been fulfilled.
varf (type="direct"): (\((k^\ast + m) \times 1\)) reconciled forecasts variance vector for \(h=1\), \(\mbox{diag}(\mathbf{MW}\)).
M (type="direct"): Projection matrix, \(\mathbf{M}\) (projection approach).
G (type="S" and type="direct"): Projection matrix, \(\mathbf{G}\) (structural approach, \(\mathbf{M}=\mathbf{R}\mathbf{G}\)).
S (type="S" and type="direct"): Temporal summing matrix, \(\mathbf{R}\).
info (type="osqp"): matrix with information in columns for each forecast horizon \(h\) (rows): run time (run_time), number of iteration (iter), norm of primal residual (pri_res), status of osqp's solution (status) and polish's status (status_polish). It will also be returned with nn = TRUE if a solver (see nn_type) will be use.

Only if comb = "bu", the function returns recf, R and M.

Details

Let \(m\) be the highest available sampling frequency per seasonal cycle, and denote \({\cal K} = \left\lbrace k_m, k_{p-1}, \ldots, k_{2}, k_1\right\rbrace\) the \(p\) factors of \(m\), in descending order, where \(k_p=m\), and \(k_1=1\). Define \(\mathbf{K}\) the \(\left( k^\ast \times m\right)\) temporal aggregation matrix converting the high-frequency observations into lower-frequency (temporally aggregated) ones: \[ \mathbf{K} = \left[\begin{array}{c} \mathbf{1}_m' \cr \mathbf{I}_{\frac{m}{k_{p-1}}} \otimes \mathbf{1}_{k_{p-1}}' \cr \vdots \cr \mathbf{I}_{\frac{m}{k_{2}}} \otimes \mathbf{1}_{k_{2}}' \cr \end{array}\right].\] Denote \(\mathbf{R} = \left[\begin{array}{c} \mathbf{K} \cr \mathbf{I}_m \end{array}\right]\) the \(\left[(k^\ast+m) \times m \right]\) temporal summing matrix, and \(\mathbf{Z}' = \left[ \mathbf{I}_{k^\ast} \; -\mathbf{K} \right]\) the zero constraints kernel matrix.

Suppose we have the \(\left[(k^\ast+m) \times 1\right]\) vector \(\widehat{\mathbf{y}}\) of unbiased base forecasts for the \(p\) temporal aggregates of a single time series \(Y\) within a complete time cycle, i.e. at the forecast horizon \(h=1\) for the lowest (most aggregated) time frequency. If the base forecasts have been independently obtained, generally they do not fulfill the temporal aggregation constraints, i.e. \(\mathbf{Z}' \widehat{\mathbf{y}} \ne \mathbf{0}_{(k^\ast \times 1)}\). By adapting the general point forecast reconciliation according to the projection approach (type = "M"), the vector of temporally reconciled forecasts is given by: \[\widetilde{\mathbf{y}} = \widehat{\mathbf{y}} - \mathbf{\Omega}\mathbf{Z}\left(\mathbf{Z}'\mathbf{\Omega} \mathbf{Z}\right)^{-1}\mathbf{Z}'\widehat{\mathbf{y}},\] where \(\mathbf{\Omega}\) is a \(\left[(k^\ast+m) \times (k^\ast+m)\right]\) p.d. matrix, assumed known. The alternative equivalent solution (type = "S") following the structural reconciliation approach by Athanasopoulos et al. (2017) is given by: \[\widetilde{\mathbf{y}} = \mathbf{R}\left(\mathbf{R}' \mathbf{\Omega}^{-1}\mathbf{R}\right)^{-1}\mathbf{R}' \mathbf{\Omega}^{-1}\widehat{\mathbf{y}}.\]

Bounds on the reconciled forecasts

When the reconciliation makes use of the optimization package osqp, the user may impose bounds on the reconciled forecasts. The parameter bounds permits to consider lower (\(\mathbf{a}\)) and upper (\(\mathbf{b}\)) bounds like \(\mathbf{a} \leq \widetilde{\mathbf{y}} \leq \mathbf{b}\) such that: \[ \begin{array}{c} a_1 \leq \widetilde{y}_1 \leq b_1 \cr \dots \cr a_{(k^\ast + m)} \leq \widetilde{y}_{(k^\ast + m)} \leq b_{(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_{(k^\ast + m)} & b_{(k^\ast + m)} \cr \end{array}\right],\] where \(a_i \in [- \infty, + \infty]\) and \(b_i \in [- \infty, + \infty]\). If \(y_i\) is unbounded, the \(i\)-th row of bounds would be equal to c(-Inf, +Inf). Notice that if the bounds parameter is used, sol = "osqp" must be used. This is not true in the case of non-negativity constraints:

sol = "direct": first the base forecasts are reconciled without non-negativity constraints, then, if negative reconciled values are present, the "osqp" solver is used;
sol = "osqp": the base forecasts are reconciled using the "osqp" solver.

In this case it is not necessary to build a matrix containing the bounds, and it is sufficient to set nn = "TRUE".

Non-negative reconciled forecasts may be obtained by setting nn_type alternatively as:

nn_type = "sntz" ("set-negative-to-zero")
nn_type = "osqp" (Stellato et al., 2020)

References

Athanasopoulos, G., Hyndman, R.J., Kourentzes, N., Petropoulos, F. (2017), Forecasting with Temporal Hierarchies, European Journal of Operational Research, 262, 1, 60-74.

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., and Girolimetto, D. (2023), Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives, International Journal of Forecasting, 39(1), 39-57.

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.

Nystrup, P., Lindström, E., Pinson, P., Madsen, H. (2020), Temporal hierarchies with autocorrelation for load forecasting, European Journal of Operational Research, 280, 1, 876-888.

Schä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ä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. (2020). OSQP: An Operator Splitting Solver for Quadratic Programs, 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), https://CRAN.R-project.org/package=osqp.

Examples

data(FoReco_data)
# 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, ]
obj <- thfrec(topbase, m = 12, comb = "acov", res = topres)