Cross-sectional (contemporaneous) forecast reconciliation

Source: R/htsrec.R

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.

htsrec(basef, comb, C, res, Ut, nb, mse = TRUE, corpcor = FALSE,
       type = "M", sol = "direct", keep = "list",  v = NULL, nn = FALSE,
       nn_type = "osqp", settings = list(), bounds = NULL, W = NULL)

Arguments

basef

(\(h \times n\)) matrix of base forecasts to be reconciled; \(h\) is the forecast horizon and \(n\) is the total number of time series.

comb

Type of the reconciliation. Except for Bottom-up, each option corresponds to a specific (\(n \times 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.

C

(\(n_a \times n_b\)) cross-sectional (contemporaneous) matrix mapping the bottom level series into the higher level ones.

res

(\(N \times n\)) in-sample residuals matrix needed when comb = {"wls", "shr", "sam"}. The columns must be in the same order as basef.

Ut

Zero constraints cross-sectional (contemporaneous) kernel matrix \((\mathbf{U}'\mathbf{y} = \mathbf{0})\) spanning the null space valid for the reconciled forecasts. It can be used instead of parameter C, but nb (\(n = n_a + n_b\)) is needed if \(\mathbf{U}' \neq [\mathbf{I} \ -\mathbf{C}]\). If the hierarchy admits a structural representation, \(\mathbf{U}'\) has dimension (\(n_a \times n\)).

nb

Number of bottom time series; if C is present, nb and Ut are not used.

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 summing matrix S.

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

(\(n \times 2\)) matrix containing the cross-sectional bounds: the first column is the lower bound, and the second column is the upper bound.

W

This option permits to directly enter the covariance matrix:

  1. W must be a p.d. (\(n \times n\)) matrix or a list of \(h\) matrix (one for each forecast horizon);

  2. if comb is different from "w", W is not used.

Value

If the parameter keep is equal to "recf", then the function returns only the (\(h \times n\)) reconciled forecasts matrix, 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 \times n\)) reconciled forecasts matrix, \(\widetilde{\mathbf{Y}}\).

W

Covariance matrix used for forecast reconciliation, \(\mathbf{W}\).

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

(\(n \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{S}\mathbf{G}\)).

S (type="S" and type="direct")

Cross-sectional summing matrix, \(\mathbf{S}\).

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

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

Details

Let \(\mathbf{y}\) be a (\(n \times 1\)) vector of target point forecasts which are wished to satisfy the system of linearly independent constraints \[\mathbf{U}'\mathbf{y} = \mathbf{0}_{(r \times 1)},\] where \(\mathbf{U}'\) is a (\(r \times n\)) matrix, with rank\((\mathbf{U}') = r \leq n\), and \(\mathbf{0}_{(r \times 1)}\) is a (\(r \times 1\)) null vector. Let \(\widehat{\mathbf{y}}\) be a (\(n \times 1\)) vector of unbiased point forecasts, not fulfilling the linear constraints (i.e., \(\mathbf{U}'\widehat{\mathbf{y}} \ne \mathbf{0}\)).

We consider a regression-based reconciliation method assuming that \(\widehat{\mathbf{y}}\) is related to \(\mathbf{y}\) by \[\widehat{\mathbf{y}} = \mathbf{y} + \mathbf{\varepsilon},\] where \(\mathbf{\varepsilon}\) is a (\(n \times 1\)) vector of zero mean disturbances, with known p.d. covariance matrix \(\mathbf{W}\). The reconciled forecasts \(\widetilde{\mathbf{y}}\) are found by minimizing the generalized least squares (GLS) objective function \(\left(\widehat{\mathbf{y}} - \mathbf{y}\right)'\mathbf{W}^{-1} \left(\widehat{\mathbf{y}} - \mathbf{y}\right)\) constrained by \(\mathbf{U}'\mathbf{y} = \mathbf{0}_{(r \times 1)}\):

\[\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 \[\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 \(\mathbf{M} = \mathbf{I}_n - \mathbf{W}\mathbf{U}\left( \mathbf{U}’\mathbf{WU}\right)^{-1}\mathbf{U}’\) is a (\(n \times n\)) projection matrix. This solution is used by htsrec when type = "M".

Denoting with \(\mathbf{d}_{\widehat{\mathbf{y}}} = \mathbf{0} - \mathbf{U}'\widehat{\mathbf{y}}\) the (\(r \times 1\)) vector containing the coherency errors of the base forecasts, we can re-state the solution as \[\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 \(\widehat{\mathbf{y}}\) with a linear combination -- according to the smoothing matrix \(\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 \(\mathbf{\varepsilon}\) is gaussian, the reconciliation error \(\widetilde{\varepsilon} = \widetilde{\mathbf{y}} - \mathbf{y}\) is a zero-mean gaussian vector with covariance matrix \[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 \[\widehat{\mathbf{y}} = \mathbf{S}\mathbf{\beta} + \mathbf{\varepsilon},\] where \(\mathbf{S}\) is the structural summation matrix describing the aggregation relationships operating on \(\mathbf{y}\), and \(\mathbf{\beta}\) is a subset of \(\mathbf{y}\) containing the target forecasts of the bottom level series, such that \(\mathbf{y} = \mathbf{S}\mathbf{\beta}\). Since the hypotheses on \(\mathbf{\varepsilon}\) remain unchanged, \[\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 \(\mathbf{\beta}\), and the whole reconciled forecasts vector is given by \[\widetilde{\mathbf{y}} = \mathbf{S}\widetilde{\mathbf{\beta}} = \mathbf{SG} \widehat{\mathbf{y}},\] where \(\mathbf{G} = \left(\mathbf{S}'\mathbf{W}^{-1} \mathbf{S}\right)^{-1}\mathbf{S}'\mathbf{W}^{-1}\), and \(\mathbf{M}=\mathbf{SG}\). This solution is used by htsrec when type = "S".

Bounds on the reconciled forecasts

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_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 \(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

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.

Di Fonzo, T., Marini, M. (2011), Simultaneous and two-step reconciliation of systems of time series: methodological and practical issues, 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, Computational Statistics & Data Analysis, 55, 9, 2579-2589.

Kourentzes, N., Athanasopoulos, G. (2021), Elucidate structure in intermittent demand series, European Journal of Operational Research, 288, 1, pp. 141–152.

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.

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.

See also

Other reconciliation procedures: cstrec(), ctbu(), iterec(), lccrec(), octrec(), tcsrec(), tdrec(), thfrec()

Examples

data(FoReco_data)
# monthly base forecasts
mbase <- FoReco2matrix(FoReco_data$base, m = 12)$k1
# monthly residuals
mres <- FoReco2matrix(FoReco_data$res, m = 12)$k1
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")