Optimal combination cross-temporal forecast reconciliation

Source: R/octrec.R

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

octrec(basef, m, C, comb, 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,
       Omega = NULL)

Arguments

basef

(\(n \times h(k^\ast+m)\)) matrix of base forecasts to be reconciled, \(\widehat{\mathbf{Y}}\); \(n\) is the total number of variables, \(m\) is the highest time frequency, \(k^\ast\) is the sum of (a subset of) (\(p-1\)) factors of \(m\), excluding \(m\), and \(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']'.

m

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

C

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

comb

Type of the reconciliation. It corresponds to a specific (\(n(k\ast + m) \times n(k^\ast + m)\)) covariance matrix (Di Fonzo and Girolimetto 2023, Girolimetto et al. 2023), where \(k^\ast\) is the sum of (a subset of) (\(p-1\)) factors of \(m\) (\(m\) is not considered) and \(n\) is the number of variables:

  • ols (Identity);

  • struc (Cross-temporal structural variances);

  • cs_struc (Cross-sectional structural variances and temporal independence);

  • t_struc (Cross-sectional independence and temporal structural variances);

  • wlsh (Hierarchy variances);

  • wlsv (Series variances);

  • bdshr (Shrunk cross-covariance matrix, cross-sectional framework);

  • bdsam (Sample cross-covariance matrix, cross-sectional framework);

  • acov (Series auto-covariance matrix);

  • Sshr (Series shrunk cross-covariance matrix);

  • Ssam (Series cross-covariance matrix);

  • shr (Shrunk cross-covariance matrix);

  • sam (Sample cross-covariance matrix);

  • hbshr (Shrunk high frequency bottom time series cross-covariance matrix);

  • hbsam (Sample high frequency bottom time series cross-covariance matrix);

  • bshr (Shrunk bottom time series cross-covariance matrix);

  • bsam (Sample bottom time series cross-covariance matrix);

  • hshr (Shrunk high frequency cross-covariance matrix);

  • hsam (Sample high frequency cross-covariance matrix);

  • w use your personal matrix W in param W;

  • omega use your personal matrix Omega in param Omega.

res

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

Ut

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 nb (\(n = n_a + n_b\)) is needed if \(\textbf{U}' \neq [\textbf{I} \ -\textbf{C}]\). If the hierarchy admits a structural representation, \(\textbf{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(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.

W, Omega

This option permits to directly enter the covariance matrix:

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

  2. Omega must be a p.d. (\(n(k^\ast + m) \times n(k^\ast + m)\)) matrix or a list of h matrix (one for each forecast horizon);

  3. if comb is different from "w" or "omega", W or Omega is not used.

Value

If the parameter keep is equal to "recf", then the function returns only the (\(n \times h(k^\ast + m)\)) 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

(\(n \times h(k^\ast + m)\)) reconciled forecasts matrix, \(\widetilde{\textbf{Y}}\).

Omega

Covariance matrix used for reconciled forecasts (\(\mbox{vec}(\widehat{\textbf{Y}}')\) representation).

W

Covariance matrix used for reconciled forecasts (\(\mbox{vec}(\widehat{\textbf{Y}})\) representation).

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 (k^\ast + m)\)) reconciled forecasts variance matrix for \(h=1\), \(\mbox{diag}(\mathbf{MW}\)).

M (type="direct")

Projection matrix (projection approach).

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

Projection matrix (structural approach, \(\mathbf{M}=\mathbf{S}\mathbf{G}\)).

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

Cross-temporal summing matrix (\(\widetilde{\textbf{F}}\mbox{vec}(\widehat{\textbf{Y}}')\) representation).

info (type="osqp")

matrix with some useful indicators (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).

Details

Considering contemporaneous and temporal dimensions in the same framework requires to extend and adapt the notations used in htsrec and thfrec. To do that, we define the matrix containing the base forecasts at any considered temporal frequency as \[ \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 \(\cal K\) is a subset of \(p\) factors of \(m\) and, \(\widehat{\textbf{B}}^{[k]}\) and \(\widehat{\textbf{A}}^{[k]}\) are the matrices containing the \(k\)-order temporal aggregates of the bts and uts, of dimension (\(n_b \times h m/k\)) and (\(n_a \times h m/k\)), respectively.

Let us consider the multivariate regression model \[\widehat{\mathbf{Y}} = \mathbf{Y} + \mathbf{E} ,\] where the involved matrices have each dimension \([n \times (k^\ast+m)]\) and contain, respectively, the base (\(\widehat{\mathbf{Y}}\)) and the target forecasts (\(\mathbf{Y}\)), and the coherency errors (\(\mathbf{E}\)) for the \(n\) component variables of the linearly constrained time series of interest. For each variable, \(k^\ast + m\) base forecasts are available, pertaining to all aggregation levels of the temporal hierarchy for a complete cycle of high-frequency observation, \(m\). Consider now two vectorized versions of model, by transforming the matrices either in original form: \[\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: \[\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 \(\mathbf{P}\) the \([n(k^\ast+m) \times n(k^\ast+m)]\) commutation matrix such that \(\mathbf{P}\mbox{vec}(\mathbf{Y}) = \mbox{vec}(\mathbf{Y}')\), \(\mathbf{P}\mbox{vec}(\widehat{\mathbf{Y}}) = \mbox{vec}(\widehat{\mathbf{Y}}')\) and \(\mathbf{P}\mathbf{\varepsilon} = {\bf \eta}\). Let \(\mathbf{W} = \mathrm{E}[\mathbf{\varepsilon\varepsilon}']\) be the covariance matrix of vector \(\mathbf{\varepsilon}\), and \(\mathbf{\Omega} = \mathrm{E}[\mathbf{\eta\eta}']\) the covariance matrix of vector \(\mathbf{\eta}\). Clearly, \(\mathbf{W}\) and \(\mathbf{\Omega}\) are different parameterizations of the same statistical object for which the following relationships hold: \[\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 (type = "M") to a cross-temporal forecast reconciliation problem, we may consider either two vec-forms , e.g. if we follow the first: \[ \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 \(\widehat{\mathbf{y}} = \mbox{vec}(\widehat{\mathbf{Y}}')\) is the row vectorization of the base forecasts matrix \(\widehat{\mathbf{Y}}\) The alternative equivalent solution (type = "S") (following the structural reconciliation approach by Hyndman et al., 2011) is \[\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 \(\widetilde{\mathbf{S}}\) is the cross-temporal summing matrix.

Bounds on the reconciled forecasts

When the reconciliation uses 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}\), where \(\widehat{\mathbf{y}} = \mbox{vec}(\widehat{\mathbf{Y}}')\), such that: \[ \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 \(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.

Girolimetto, D., Athanasopoulos, G., Di Fonzo, T., & Hyndman, R. J. (2023), Cross-temporal Probabilistic Forecast Reconciliation, doi:10.48550/arXiv.2303.17277 .

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.

See also

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

Examples

data(FoReco_data)
obj <- octrec(FoReco_data$base, m = 12, C = FoReco_data$C,
              comb = "bdshr", res = FoReco_data$res)