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)(h n) matrix of base forecasts to be reconciled; h is the forecast horizon and n is the total number of time series.
Type of the reconciliation. Except for Bottom-up, each option corresponds to a specific (n 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.
(n_a n_b) cross-sectional (contemporaneous) matrix mapping the bottom level series into the higher level ones.
(N n) in-sample residuals matrix needed when comb =
{"wls", "shr", "sam"}. The columns must be in
the same order as basef.
Zero constraints cross-sectional (contemporaneous) kernel matrix
(U'y = 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
U' [I \ -C]. If the hierarchy
admits a structural representation, U' has dimension
(n_a n).
Number of bottom time series; if C is present, nb
and Ut are not used.
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.
Logical value: TRUE if corpcor (Sch<U+00E4>fer et
al., 2017) must be used to shrink the sample covariance matrix according to
Sch<U+00E4>fer and Strimmer (2005), otherwise the function uses the
same implementation as package hts.
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.
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).
Return a list object of the reconciled forecasts at all levels (if keep = "list") or only the reconciled forecasts matrix (if keep = "recf").
vector index of the fixed base forecast (min(v) > 0 and max(v) < n).
Logical value: TRUE if non-negative reconciled forecasts
are wished.
"osqp" (default), "KAnn" (only type == "M") or "sntz".
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).
(n 2) matrix containing the cross-sectional bounds: the first column is the lower bound, and the second column is the upper bound.
This option permits to directly enter the covariance matrix:
W must be a p.d. (n 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 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 n) reconciled forecasts matrix, Y.
WCovariance matrix used for forecast reconciliation, W.
nn_checkNumber of negative values (if zero, there are no values below zero).
rec_checkLogical value: has the hierarchy been respected?
varf (type="direct")(n 1) reconciled forecasts variance vector for h=1, diag(MW).
M (type="direct")Projection matrix, M (projection approach).
G (type="S" and type="direct")Projection matrix, G (structural approach, M=SG).
S (type="S" and type="direct")Cross-sectional summing matrix, 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 y be a (n 1) vector of target point forecasts which are wished to satisfy the system of linearly independent constraints U'y = 0_(r 1), where U' is a (r n) matrix, with rank(U') = r n, and 0_(r 1) is a (r 1) null vector. Let y be a (n 1) vector of unbiased point forecasts, not fulfilling the linear constraints (i.e., U'y 0).
We consider a regression-based reconciliation method assuming that y is related to y by y = y + , where is a (n 1) vector of zero mean disturbances, with known p.d. covariance matrix W. The reconciled forecasts y are found by minimizing the generalized least squares (GLS) objective function (y - y)'W^-1 (y - y) constrained by U'y = 0_(r 1):
y = argmin_y (y - y )' W^-1 (y - y ), s.t. U'y = 0.
The solution is given by
y= y - WU
(U<U+2019>WU)^-1U'y=
My,
where M = I_n - WU(
U<U+2019>WU)^-1U<U+2019>
is a (n n) projection matrix. This solution is used by
htsrec when type = "M".
Denoting with d_y = 0 - U'y the (r 1) vector containing the coherency errors of the base forecasts, we can re-state the solution as y= y + WU ( U'WU)^-1d_y, which makes it clear that the reconciliation formula simply adjusts the vector y with a linear combination -- according to the smoothing matrix L = WU (U<U+2019>WU)^-1 -- of the coherency errors of the base forecasts.
In addition, if the error term is gaussian, the reconciliation error = y - y is a zero-mean gaussian vector with covariance matrix E( y - y) ( y - y)' = W - WU (U'WU)^-1U' = 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
y = S + ,
where S is the structural summation matrix describing
the aggregation relationships operating on y, and
is a subset of y containing the
target forecasts of the bottom level series, such that
y = S. Since the hypotheses on
remain unchanged,
= (S'W^-1S
)^-1S'W^-1y
is the best linear unbiased estimate of , and
the whole reconciled forecasts vector is given by
y = S = SG
y,
where G = (S'W^-1
S)^-1S'W^-1, and M=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 (a) and
upper (b) bounds like a
y b such that:
arrayc
a_1 y_1 b_1
…
a_n y_n b_n
array
bounds = [a \; b] =
[arraycc
a_1 & b_1
&
a_n & b_n
array],
where a_i [- , + ] and b_i [- , + ].
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., and Girolimetto, D. (2021), Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives, International Journal of Forecasting, in press.
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<U+2013>152.
Sch<U+00E4>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<U+00E4>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()
# NOT RUN {
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")
# }
Run the code above in your browser using DataLab