This function performs optimal (in least squares sense) combination cross-temporal forecast reconciliation (Di Fonzo and Girolimetto 2023a, Girolimetto et al. 2023). The reconciled forecasts are calculated using either a projection approach (Byron, 1978, 1979) or the equivalent structural approach by Hyndman et al. (2011). Non-negative (Di Fonzo and Girolimetto, 2023) and immutable reconciled forecasts can be considered.
Usage
ctrec(base, agg_mat, cons_mat, agg_order, comb = "ols", res = NULL,
tew = "sum", approach = "proj", nn = NULL, settings = NULL,
bounds = NULL, immutable = NULL, ...)
Arguments
- base
A (\(n \times h(k^\ast+m)\)) numeric matrix containing the base forecasts to be reconciled; \(n\) is the total number of variables, \(m\) is the max. order of temporal aggregation, \(k^\ast\) is the sum of (a subset of) (\(p-1\)) factors of \(m\), excluding \(m\), and \(h\) is the forecast horizon for the lowest frequency time series. The row identifies a time series, and the forecasts in each row are ordered from the lowest frequency (most temporally aggregated) to the highest frequency.
- agg_mat
A (\(n_a \times n_b\)) numeric matrix representing the cross-sectional aggregation matrix. It maps the \(n_b\) bottom-level (free) variables into the \(n_a\) upper (constrained) variables.
- cons_mat
A (\(n_a \times n\)) numeric matrix representing the cross-sectional zero constraints. It spans the null space for the reconciled forecasts.
- agg_order
Highest available sampling frequency per seasonal cycle (max. order of temporal aggregation, \(m\)), or a vector representing a subset of \(p\) factors of \(m\).
- comb
A string specifying the reconciliation method. For a complete list, see ctcov.
- res
A (\(n \times N(k^\ast+m)\)) optional numeric matrix containing the in-sample residuals at all the temporal frequencies ordered from the lowest frequency to the highest frequency (columns) for each variable (rows). This matrix is used to compute some covariance matrices.
- tew
A string specifying the type of temporal aggregation. Options include: "
sum
" (simple summation, default), "avg
" (average), "first
" (first value of the period), and "last
" (last value of the period).- approach
A string specifying the approach used to compute the reconciled forecasts. Options include:
- nn
A string specifying the algorithm to compute non-negative reconciled forecasts:
"
osqp
": quadratic programming optimization (osqp solver)."
sntz
": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
- settings
An object of class
osqpSettings
specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020).- bounds
A (\(n(k^\ast + m) \times 2\)) numeric matrix specifying the cross-temporal bounds. The first column represents the lower bound, and the second column represents the upper bound.
- immutable
A matrix with three columns (\(i,k,j\)), such that
- Column 1
Represents the cross-sectional series (\(i = 1, \dots, n\)).
- Column 2
Denotes the temporal aggregation order (\(k = m,\dots,1\)).
- Column 3
Indicates the temporal forecast horizon (\(j = 1,\dots,m/k\)).
For example, when working with a quarterly multivariate time series (\(n = 3\)):
t(c(1, 4, 1))
- Fix the one step ahead annual forecast of the first time series.t(c(2, 1, 2))
- Fix the two step ahead quarterly forecast of the second time series.
- ...
Arguments passed on to
ctcov
mse
If
TRUE
(default) the residuals used to compute the covariance matrix are not mean-corrected.shrink_fun
Shrinkage function of the covariance matrix, shrink_estim (default).
References
Byron, R.P. (1978), The estimation of large social account matrices, Journal of the Royal Statistical Society, Series A, 141, 3, 359-367. doi:10.2307/2344807
Byron, R.P. (1979), Corrigenda: The estimation of large social account matrices, Journal of the Royal Statistical Society, Series A, 142(3), 405. doi:10.2307/2982515
Di Fonzo, T. and Girolimetto, D. (2023a), Cross-temporal forecast reconciliation: Optimal combination method and heuristic alternatives, International Journal of Forecasting, 39, 1, 39-57. doi:10.1016/j.ijforecast.2021.08.004
Di Fonzo, T. and Girolimetto, D. (2023), Spatio-temporal reconciliation of solar forecasts, Solar Energy, 251, 13–29. doi:10.1016/j.solener.2023.01.003
Girolimetto, D., Athanasopoulos, G., Di Fonzo, T. and Hyndman, R.J. (2024), Cross-temporal probabilistic forecast reconciliation: Methodological and practical issues. International Journal of Forecasting, 40, 3, 1134-1151. doi:10.1016/j.ijforecast.2023.10.003
Hyndman, R.J., Ahmed, R.A., Athanasopoulos, G. and Shang, H.L. (2011), Optimal combination forecasts for hierarchical time series, Computational Statistics & Data Analysis, 55, 9, 2579-2589. doi:10.1016/j.csda.2011.03.006
Stellato, B., Banjac, G., Goulart, P., Bemporad, A. and Boyd, S. (2020), OSQP: An Operator Splitting solver for Quadratic Programs, Mathematical Programming Computation, 12, 4, 637-672. doi:10.1007/s12532-020-00179-2
Examples
set.seed(123)
# (3 x 7) base forecasts matrix (simulated), Z = X + Y and m = 4
base <- rbind(rnorm(7, rep(c(20, 10, 5), c(1, 2, 4))),
rnorm(7, rep(c(10, 5, 2.5), c(1, 2, 4))),
rnorm(7, rep(c(10, 5, 2.5), c(1, 2, 4))))
# (3 x 70) in-sample residuals matrix (simulated)
res <- rbind(rnorm(70), rnorm(70), rnorm(70))
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
m <- 4 # from quarterly to annual temporal aggregation
reco <- ctrec(base = base, agg_mat = A, agg_order = m, comb = "wlsv", res = res)
C <- t(c(1, -1, -1)) # Zero constraints matrix for Z - X - Y = 0
reco <- ctrec(base = base, cons_mat = C, agg_order = m, comb = "wlsv", res = res)
# Immutable reconciled forecasts
# Fix all the quarterly forecasts of the second variable.
imm_mat <- expand.grid(i = 2, k = 1, j = 1:4)
immreco <- ctrec(base = base, cons_mat = C, agg_order = m, comb = "wlsv",
res = res, immutable = imm_mat)
# Non negative reconciliation
base[2,7] <- -2*base[2,7] # Making negative one of the quarterly base forecasts for variable X
nnreco <- ctrec(base = base, cons_mat = C, agg_order = m, comb = "wlsv",
res = res, nn = "osqp")
recoinfo(nnreco, verbose = FALSE)$info
#> obj_val run_time iter pri_res status status_polish
#> 1 -635.7645 9.4207e-05 25 1.284524e-15 1 1