htsrec | R Documentation |
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 |
(\mjseqnh \times n) matrix of base forecasts to be reconciled; \mjseqnh is the forecast horizon and \mjseqnn is the total number of time series. |
comb |
Type of the reconciliation. Except for Bottom-up, each option corresponds to a specific (\mjseqnn \times n) covariance matrix:
|
C |
(\mjseqnn_a \times n_b) cross-sectional (contemporaneous) matrix mapping the bottom level series into the higher level ones. |
res |
(\mjseqnN \times n) in-sample residuals matrix needed when |
Ut |
Zero constraints cross-sectional (contemporaneous) kernel matrix
\mjseqn(\mathbfU'\mathbfy = \mathbf0) 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 (\mjseqn\mboxmin(v) > 0 and \mjseqn\mboxmax(v) < n). |
nn |
Logical value: |
nn_type |
"osqp" (default), "KAnn" (only |
settings |
Settings for osqp (object |
bounds |
(\mjseqnn \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:
|
Let \mjseqn\mathbfy be a (\mjseqnn \times 1) vector of target point forecasts which are wished to satisfy the system of linearly independent constraints \mjsdeqn\mathbfU'\mathbfy = \mathbf0_(r \times 1), where \mjseqn\mathbfU' is a (\mjseqnr \times n) matrix, with rank\mjseqn(\mathbfU') = r \leq n, and \mjseqn\mathbf0_(r \times 1) is a (\mjseqnr \times 1) null vector. Let \mjseqn\widehat\mathbfy be a (\mjseqnn \times 1) vector of unbiased point forecasts, not fulfilling the linear constraints (i.e., \mjseqn\mathbfU'\widehat\mathbfy \ne \mathbf0).
We consider a regression-based reconciliation method assuming that \mjseqn\widehat\mathbfy is related to \mjseqn\mathbfy by \mjsdeqn\widehat\mathbfy = \mathbfy + \mathbf\varepsilon, where \mjseqn\mathbf\varepsilon is a (\mjseqnn \times 1) vector of zero mean disturbances, with known p.d. covariance matrix \mjseqn\mathbfW. The reconciled forecasts \mjseqn\widetilde\mathbfy are found by minimizing the generalized least squares (GLS) objective function \mjseqn\left(\widehat\mathbfy - \mathbfy\right)'\mathbfW^-1 \left(\widehat\mathbfy - \mathbfy\right) constrained by \mjseqn\mathbfU'\mathbfy = \mathbf0_(r \times 1):
\mjsdeqn\widetilde\mathbfy = \mboxargmin_\mathbfy \left(\mathbfy - \widehat\mathbfy \right)' \mathbfW^-1 \left(\mathbfy - \widehat\mathbfy \right), \quad \mboxs.t. \mathbfU'\mathbfy = \mathbf0.
The solution is given by
\mjsdeqn\widetilde\mathbfy= \widehat\mathbfy - \mathbfW\mathbfU
\left(\mathbfU’\mathbfWU\right)^-1\mathbfU'\widehat\mathbfy=
\mathbfM\widehat\mathbfy,
where \mjseqn\mathbfM = \mathbfI_n - \mathbfW\mathbfU\left(
\mathbfU’\mathbfWU\right)^-1\mathbfU’
is a (\mjseqnn \times n) projection matrix. This solution is used by
htsrec
when type = "M"
.
Denoting with \mjseqn\mathbfd_\widehat\mathbfy = \mathbf0 - \mathbfU'\widehat\mathbfy the (\mjseqnr \times 1) vector containing the coherency errors of the base forecasts, we can re-state the solution as \mjsdeqn\widetilde\mathbfy= \widehat\mathbfy + \mathbfWU \left( \mathbfU'\mathbfWU\right)^-1\mathbfd_\widehaty, which makes it clear that the reconciliation formula simply adjusts the vector \mjseqn\widehat\mathbfy with a linear combination – according to the smoothing matrix \mjseqn\mathbfL = \mathbfWU \left(\mathbfU’\mathbfWU\right)^-1 – of the coherency errors of the base forecasts.
In addition, if the error term \mjseqn\mathbf\varepsilon is gaussian, the reconciliation error \mjseqn\widetilde\varepsilon = \widetilde\mathbfy - \mathbfy is a zero-mean gaussian vector with covariance matrix \mjsdeqnE\left( \widetilde\mathbfy - \mathbfy\right) \left( \widetilde\mathbfy - \mathbfy\right)' = \mathbfW - \mathbfWU \left(\mathbfU'\mathbfWU\right)^-1\mathbfU' = \mathbfMW.
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
\mjsdeqn\widehat\mathbfy = \mathbfS\mathbf\beta + \mathbf\varepsilon,
where \mjseqn\mathbfS is the structural summation matrix describing
the aggregation relationships operating on \mjseqn\mathbfy, and
\mjseqn\mathbf\beta is a subset of \mjseqn\mathbfy containing the
target forecasts of the bottom level series, such that
\mjseqn\mathbfy = \mathbfS\mathbf\beta. Since the hypotheses on
\mjseqn\mathbf\varepsilon remain unchanged,
\mjsdeqn\widetilde\mathbf\beta = \left(\mathbfS'\mathbfW^-1\mathbfS
\right)^-1\mathbfS'\mathbfW^-1\widehat\mathbfy
is the best linear unbiased estimate of \mjseqn\mathbf\beta, and
the whole reconciled forecasts vector is given by
\mjsdeqn\widetilde\mathbfy = \mathbfS\widetilde\mathbf\beta = \mathbfSG
\widehat\mathbfy,
where \mjseqn\mathbfG = \left(\mathbfS'\mathbfW^-1
\mathbfS\right)^-1\mathbfS'\mathbfW^-1, and \mjseqn\mathbfM=\mathbfSG.
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 (\mjseqn\mathbfa) and
upper (\mjseqn\mathbfb) bounds like \mjseqn\mathbfa \leq
\widetilde\mathbfy \leq \mathbfb such that:
\mjsdeqn \beginarrayc
a_1 \leq \widetildey_1 \leq b_1
...
a_n \leq \widetildey_n \leq b_n
\endarray \Rightarrow
\mboxbounds = [\mathbfa \; \mathbfb] =
\left[\beginarraycc
a_1 & b_1
\vdots & \vdots
a_n & b_n
\endarray\right],
where \mjseqna_i \in [- \infty, + \infty] and \mjseqnb_i \in [- \infty, + \infty].
If \mjseqny_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)
If the parameter keep
is equal to "recf"
, then the function
returns only the (\mjseqnh \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 |
(\mjseqnh \times n) reconciled forecasts matrix, \mjseqn\widetilde\mathbfY. |
W |
Covariance matrix used for forecast reconciliation, \mjseqn\mathbfW. |
nn_check |
Number of negative values (if zero, there are no values below zero). |
rec_check |
Logical value: |
varf (type="direct" ) |
(\mjseqnn \times 1) reconciled forecasts variance vector for \mjseqnh=1, \mjseqn\mboxdiag(\mathbfMW). |
M (type="direct" ) |
Projection matrix, \mjseqn\mathbfM (projection approach). |
G (type="S" and type="direct" ) |
Projection matrix, \mjseqn\mathbfG (structural approach, \mjseqn\mathbfM=\mathbfS\mathbfG). |
S (type="S" and type="direct" ) |
Cross-sectional summing matrix, \mjseqn\mathbfS. |
info (type="osqp" ) |
matrix with information in columns
for each forecast horizon \mjseqnh (rows): run time ( |
Only if comb = "bu"
, the function returns recf
, S
and M
.
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.
Other reconciliation procedures:
cstrec()
,
ctbu()
,
iterec()
,
lccrec()
,
octrec()
,
tcsrec()
,
tdrec()
,
thfrec()
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")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.