htsrec {FoReco} | R Documentation |
Cross-sectional (contemporaneous) forecast reconciliation
Description
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.
Usage
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:
|
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 |
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
|
nb |
Number of bottom time series; if |
mse |
Logical value: |
corpcor |
Logical value: |
type |
Approach used to compute the reconciled forecasts: |
sol |
Solution technique for the reconciliation problem: either
|
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: |
nn_type |
"osqp" (default), "KAnn" (only |
settings |
Settings for osqp (object |
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:
|
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)
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: |
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 ( |
Only if comb = "bu"
, the function returns recf
, S
and M
.
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")