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

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

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

`W`

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

`comb`

is different from "`w`

",`W`

is not used.

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

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

.

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

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

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

(Stellato et al., 2020)

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.

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.

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