`R/thfrec.R`

`thfrec.Rd`

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

- 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) < (k^\ast + m)\)).

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

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: has the hierarchy been respected?

`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`

.

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 = "KAnn"`

(Kourentzes and Athanasopoulos, 2021)`nn_type = "sntz"`

("set-negative-to-zero")`nn_type = "osqp"`

(Stellato et al., 2020)

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

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.

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