thfrec: Forecast reconciliation through temporal hierarchies...
In FoReco: Forecast Reconciliation
thfrec R Documentation
Forecast reconciliation through temporal hierarchies (temporal reconciliation)
Description
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.
Usage
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)
Arguments
basef
(\mjseqnh(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, \mjseqnm), or a subset of \mjseqnp factors
of \mjseqnm.
comb
Type of the reconciliation. Except for bottom up, all other
options correspond to a different (\mjseqn(k^\ast + m) \times (k^\ast + m))
covariance matrix, \mjseqnk^\ast is the sum of (\mjseqnp-1) factors of
\mjseqnm (excluding \mjseqnm):
-
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 (\mjseqn\mboxmin(v) > 0
and \mjseqn\mboxmax(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
(\mjseqn(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. (\mjseqn(k^\ast + m) \times (k^\ast + m))
matrix or a list of \mjseqnh matrix (one for each forecast horizon);
if comb is different from "omega", Omega is
not used.
Details
\loadmathjax
Let \mjseqnm be the highest available
sampling frequency per seasonal cycle, and denote
\mjseqn\cal K = \left\lbrace k_m, k_p-1, ..., k_2, k_1\right\rbrace
the \mjseqnp factors of \mjseqnm, in descending order, where
\mjseqnk_p=m, and \mjseqnk_1=1. Define \mjseqn\mathbfK
the \mjseqn\left( k^\ast \times m\right) temporal aggregation matrix
converting the high-frequency observations into lower-frequency
(temporally aggregated) ones:
\mjsdeqn \mathbfK = \left[\beginarrayc
\mathbf1_m'
\mathbfI_\fracmk_p-1 \otimes \mathbf1_k_p-1'
\vdots
\mathbfI_\fracmk_2 \otimes \mathbf1_k_2'
\endarray\right].
Denote
\mjseqn\mathbfR = \left[\beginarrayc
\mathbfK
\mathbfI_m
\endarray\right] the \mjseqn\left[(k^\ast+m) \times m \right]
temporal summing matrix, and
\mjseqn\mathbfZ' = \left[ \mathbfI_k^\ast \; -\mathbfK \right]
the zero constraints kernel matrix.
Suppose we have the \mjseqn\left[(k^\ast+m) \times 1\right] vector
\mjseqn\widehat\mathbfy of unbiased base forecasts for the
\mjseqnp temporal aggregates of a single time series \mjseqnY
within a complete time cycle, i.e. at the forecast horizon \mjseqnh=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. \mjseqn\mathbfZ'
\widehat\mathbfy \ne \mathbf0_(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:
\mjsdeqn\widetilde\mathbfy = \widehat\mathbfy -
\mathbf\Omega\mathbfZ\left(\mathbfZ'\mathbf\Omega
\mathbfZ\right)^-1\mathbfZ'\widehat\mathbfy,
where \mjseqn\mathbf\Omega is a \mjseqn\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:
\mjsdeqn\widetilde\mathbfy = \mathbfR\left(\mathbfR'
\mathbf\Omega^-1\mathbfR\right)^-1\mathbfR'
\mathbf\Omega^-1\widehat\mathbfy.
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 (\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_(k^\ast + m) \leq \widetildey_(k^\ast + m) \leq b_(k^\ast + m)
\endarray \Rightarrow
\mboxbounds = [\mathbfa \; \mathbfb] =
\left[\beginarraycc
a_1 & b_1
\vdots & \vdots
a_(k^\ast + m) & b_(k^\ast + m)
\endarray\right],
where \mjseqna_i \in [- \infty, + \infty] and \mjseqnb_i \in [- \infty, + \infty].
If \mjseqny_i is unbounded, the \mjseqni-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 (\mjseqnh(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
(\mjseqnh(k^\ast + m) \times 1) reconciled forecasts vector, \mjseqn\widetilde\mathbfy.
Omega
Covariance matrix used for forecast reconciliation, \mjseqn\mathbf\Omega.
nn_check
Number of negative values (if zero, there are no values below zero).
rec_check
Logical value: rec_check = TRUE when the constraints have been fulfilled.
varf (type="direct")
(\mjseqn(k^\ast + m) \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=\mathbfR\mathbfG).
S (type="S" and type="direct")
Temporal summing matrix, \mjseqn\mathbfR.
info (type="osqp")
matrix with information in columns
for each forecast horizon \mjseqnh (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.
References
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., and Girolimetto, D. (2023), Cross-temporal forecast reconciliation:
Optimal combination method and heuristic alternatives, International Journal
of Forecasting, 39(1), 39-57.
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.
See Also
Other reconciliation procedures:
cstrec(),
ctbu(),
htsrec(),
iterec(),
lccrec(),
octrec(),
tcsrec(),
tdrec()
Examples
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)
FoReco documentation built on May 31, 2023, 5:17 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| thfrec | R Documentation |
Forecast reconciliation through temporal hierarchies (temporal reconciliation)
Description
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.
Usage
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)
Arguments
basef |
(\mjseqnh(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, \mjseqnm), or a subset of \mjseqnp factors of \mjseqnm. |
comb |
Type of the reconciliation. Except for bottom up, all other options correspond to a different (\mjseqn(k^\ast + m) \times (k^\ast + m)) covariance matrix, \mjseqnk^\ast is the sum of (\mjseqnp-1) factors of \mjseqnm (excluding \mjseqnm):
|
res |
vector containing the in-sample residuals at all the temporal
frequencies ordered as |
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) < (k^\ast + m)). |
nn |
Logical value: |
nn_type |
"osqp" (default), "KAnn" (only |
settings |
Settings for osqp (object |
bounds |
(\mjseqn(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:
|
Details
\loadmathjaxLet \mjseqnm be the highest available
sampling frequency per seasonal cycle, and denote
\mjseqn\cal K = \left\lbrace k_m, k_p-1, ..., k_2, k_1\right\rbrace
the \mjseqnp factors of \mjseqnm, in descending order, where
\mjseqnk_p=m, and \mjseqnk_1=1. Define \mjseqn\mathbfK
the \mjseqn\left( k^\ast \times m\right) temporal aggregation matrix
converting the high-frequency observations into lower-frequency
(temporally aggregated) ones:
\mjsdeqn \mathbfK = \left[\beginarrayc
\mathbf1_m'
\mathbfI_\fracmk_p-1 \otimes \mathbf1_k_p-1'
\vdots
\mathbfI_\fracmk_2 \otimes \mathbf1_k_2'
\endarray\right].
Denote
\mjseqn\mathbfR = \left[\beginarrayc
\mathbfK
\mathbfI_m
\endarray\right] the \mjseqn\left[(k^\ast+m) \times m \right]
temporal summing matrix, and
\mjseqn\mathbfZ' = \left[ \mathbfI_k^\ast \; -\mathbfK \right]
the zero constraints kernel matrix.
Suppose we have the \mjseqn\left[(k^\ast+m) \times 1\right] vector
\mjseqn\widehat\mathbfy of unbiased base forecasts for the
\mjseqnp temporal aggregates of a single time series \mjseqnY
within a complete time cycle, i.e. at the forecast horizon \mjseqnh=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. \mjseqn\mathbfZ'
\widehat\mathbfy \ne \mathbf0_(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:
\mjsdeqn\widetilde\mathbfy = \widehat\mathbfy -
\mathbf\Omega\mathbfZ\left(\mathbfZ'\mathbf\Omega
\mathbfZ\right)^-1\mathbfZ'\widehat\mathbfy,
where \mjseqn\mathbf\Omega is a \mjseqn\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:
\mjsdeqn\widetilde\mathbfy = \mathbfR\left(\mathbfR'
\mathbf\Omega^-1\mathbfR\right)^-1\mathbfR'
\mathbf\Omega^-1\widehat\mathbfy.
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 (\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_(k^\ast + m) \leq \widetildey_(k^\ast + m) \leq b_(k^\ast + m)
\endarray \Rightarrow
\mboxbounds = [\mathbfa \; \mathbfb] =
\left[\beginarraycc
a_1 & b_1
\vdots & \vdots
a_(k^\ast + m) & b_(k^\ast + m)
\endarray\right],
where \mjseqna_i \in [- \infty, + \infty] and \mjseqnb_i \in [- \infty, + \infty].
If \mjseqny_i is unbounded, the \mjseqni-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 (\mjseqnh(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 |
(\mjseqnh(k^\ast + m) \times 1) reconciled forecasts vector, \mjseqn\widetilde\mathbfy. |
Omega |
Covariance matrix used for forecast reconciliation, \mjseqn\mathbf\Omega. |
nn_check |
Number of negative values (if zero, there are no values below zero). |
rec_check |
Logical value: |
varf (type="direct") |
(\mjseqn(k^\ast + m) \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=\mathbfR\mathbfG). |
S (type="S" and type="direct") |
Temporal summing matrix, \mjseqn\mathbfR. |
info (type="osqp") |
matrix with information in columns
for each forecast horizon \mjseqnh (rows): run time ( |
Only if comb = "bu", the function returns recf, R and M.
References
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., and Girolimetto, D. (2023), Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives, International Journal of Forecasting, 39(1), 39-57.
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.
See Also
Other reconciliation procedures:
cstrec(),
ctbu(),
htsrec(),
iterec(),
lccrec(),
octrec(),
tcsrec(),
tdrec()
Examples
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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.