Level conditional coherent forecast reconciliation for genuine hierarchical/grouped time series

Forecast reconciliation procedure built on and extending the original proposal by Hollyman et al. (2021). Level conditional coherent reconciled forecasts may be computed in cross-sectional, temporal, and cross-temporal frameworks. The reconciled forecasts are conditional to (i.e., constrained by) the base forecasts of a specific upper level of the hierarchy (exogenous constraints). The linear constraints linking the variables may be dealt with endogenously as well (Di Fonzo and Girolimetto, 2022). Combined Conditional Coherent (CCC) forecasts may be calculated as simple averages of LCC and bottom-up reconciled forecasts, with either endogenous or exogenous constraints.

lccrec(basef, m, C, nl, weights, bnaive = NULL, const = "exogenous",
       CCC = TRUE, nn = FALSE, nn_type = "osqp", settings = list(), ...)

Arguments

basef: matrix/vector of base forecasts to be reconciled: (\(h \times n\)) matrix in the cross-sectional framework; (\(h(k^\ast + m) \times 1\)) vector in the temporal framework; (\(n \times h(k^\ast+m)\)) matrix in the cross-temporal framework. \(n\) is the total number of variables, \(m\) is the highest time frequency, \(k^\ast\) is the sum of (a subset of) (\(p-1\)) factors of \(m\), excluding \(m\), and \(h\) is the forecast horizon.
m: Highest available sampling frequency per seasonal cycle (max. order of temporal aggregation, \(m\)), or a subset of the \(p\) factors of \(m\).
C: (\(n_a \times n_b\)) cross-sectional (contemporaneous) matrix mapping the bottom level series into the higher level ones (or a list of matrices forming \(\mathbf{C} = [\mathbf{C}_1' \; \mathbf{C}_2' \; ... \; \mathbf{C}_L']'\), \(1, ..., L\) being the number of the cross-sectional upper levels.
nl: (\(L \times 1\)) vector containing the number of time series in each level of the hierarchy (nl[1] = 1).
weights: covariance matrix or a vector (weights used in the reconciliation: either (\(n_b \times 1\)) for exogenous or (\(n \times 1\)) for endogenous constraints).
bnaive: matrix/vector of naive bts base forecasts (e.g., seasonal averages, as in Hollyman et al., 2021): (\(h \times n_b\)) matrix in the cross-sectional framework; (\(h m \times 1\)) vector in the temporal framework; (\(n_b \times h m\)) matrix in the cross-temporal framework. Ignore it, if only basef are to be used as base forecasts.
const: exogenous (default) or endogenous constraints
CCC: Option to return Combined Conditional Coherent reconciled forecasts (default is TRUE).
nn: Logical value: TRUE if non-negative reconciled forecasts are wished.
nn_type: Non-negative method: "osqp" (default) or "sntz" (set-negative-to-zero, only if CCC = TRUE) with exogenous constraints (const = "exo"); "osqp" (default), "KAnn" (only type == "M") or "sntz" with endogenous constraints (const = "endo").
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).
...: Additional functional arguments passed to htsrec of FoReco.

Value

The function returns the level reconciled forecasts in case of an elementary hierarchy with one level. Otherwise it returns a list with

recf: Level Conditional Coherent (CCC = FALSE) forecasts matrix/vector calculated as simple averages of upper level conditional reconciled forecasts, with either endogenous or exogenous constraints. If CCC = TRUE then it is the Combined Conditional Coherent matrix/vector, as weighted averages of LCC and bottom-up reconciled forecasts.
levrecf: list of level conditional reconciled forecasts (+ BU).
info (type="osqp"): matrix with some useful indicators (columns) for each forecast horizon \(h\) (rows): run time (run_time), number of iteration, norm of primal residual (pri_res), status of osqp's solution (status) and polish's status (status_polish).

Details

Cross-sectional hierarchies

To be as simple as possible, we fix the forecast horizon equal to 1. Let the base forecasts be a vector \(\widehat{\mathbf{y}} = \left[\widehat{\mathbf{a}}' \; \widehat{\mathbf{b}}'\right]'\), where vector \(\widehat{\mathbf{a}}\) consists of the sub-vectors forming each of the upper \(L\) levels of the hierarchy/grouping: \[\widehat{\mathbf{a}} = \left[\begin{array}{c} \widehat{a}_1 \cr \widehat{\mathbf{a}}_2 \cr \vdots \cr \widehat{\mathbf{a}}_L \end{array}\right],\] where \(\widehat{\mathbf{a}}_l\), \(l=1,\ldots, L\), has dimension \((n_l \times 1)\) and \(\sum_{l=1}^{L} n_l = n_a\). Denote \(\mathbf{C}_l\) the \((n_l \times n_b)\) matrix mapping the bts into the level-\(l\) uts, then the aggregation matrix \(\mathbf{C}\) may be written as \[\mathbf{C} = \left[\begin{array}{c} \mathbf{C}_1 \cr\mathbf{C}_2 \cr \vdots \cr \mathbf{C}_L \end{array}\right],\] where the generic matrix \(\mathbf{C}_L\) is (\(n_L \times n_b\)), \(l=1, \ldots, L\). Given a generic level \(l\), the reconciled forecasts coherent with the base forecasts of that level are the solution to a quadratic minimization problem with linear exogenous constraints (const = "exo"): \[\begin{array}{c}\widetilde{\mathbf{b}}_{l} = \arg\min_{\mathbf{b}} \left(\mathbf{b} - \widehat{\mathbf{b}}\right)'\mathbf{W}_b^{-1} \left(\mathbf{b} - \widehat{\mathbf{b}}\right) \quad \mbox{ s.t. } \mathbf{C}_l\mathbf{b} = \widehat{\mathbf{a}}_l, \quad l=1, \ldots, L ,\cr \Downarrow \cr \widetilde{\mathbf{b}}_{l} = \widehat{\mathbf{b}} + \mathbf{W}_b\mathbf{C}_l'\left(\mathbf{C}_l\mathbf{W}_b\mathbf{C}_l' \right)^{-1}\left(\widehat{\mathbf{a}}_l -\mathbf{C}_l\widehat{\mathbf{b}} \right), \quad l=1,\ldots,L,\end{array}\] where \(\mathbf{W}_b\) is a (\(n_b \times n_b\)) p.d. matrix (in Hollyman et al., 2021, \(\mathbf{W}_b\) is diagonal). If endogenous constraints (const = "endo") are considered, denote \(\widehat{\mathbf{y}}_l = \left[\widehat{\mathbf{a}}_l' \; \widehat{\mathbf{b}}'\right]'\) and \(\mathbf{U}_l' = \left[\mathbf{I}_{n_l}' \; \mathbf{C}_l'\right]'\), then \[\begin{array}{c}\widetilde{\mathbf{y}}_{l} = \arg\min_{\mathbf{y}_l} \left(\mathbf{y}_l - \widehat{\mathbf{y}}_l\right)'\mathbf{W}_l^{-1} \left(\mathbf{y}_l - \widehat{\mathbf{y}}_l\right) \quad \mbox{ s.t. } \mathbf{U}_l'\mathbf{y}_l = \mathbf{0}, \quad l=1, \ldots, L ,\cr \Downarrow \cr \widetilde{\mathbf{y}}_{l} = \left(\mathbf{I}_{n_b+n_l} - \mathbf{W}_l\mathbf{U}_l\left(\mathbf{U}_l'\mathbf{W}_l \mathbf{U}_l\right)^{-1}\mathbf{U}_l'\right)\widehat{\mathbf{y}}_{l}, \quad l=1,...,L, \end{array}\] where \(\mathbf{W}_l\) is a (\(n_l + n_b \times n_l + n_b\)) p.d. matrix. Combined Conditional Coherent (CCC) forecasts are obtained by taking the simple average of the \(L\) level conditional, and of the bottom-up reconciled forecasts, respectively (Di Fonzo and Girolimetto, 2022): \[\widetilde{\mathbf{y}}_{CCC}=\frac{1}{L+1}\sum_{l=1}^{L+1} \mathbf{S}\widetilde{\mathbf{b}}_l,\] with \[\widetilde{\mathbf{b}}_{L+1} = \widehat{\mathbf{b}}.\]

Non-negative reconciled forecasts may be obtained by setting nn_type alternatively as:

nn_type = "osqp" to apply non-negative constraints to each level:
nn_type = "sntz" ("set-negative-to-zero") to apply non-negative constraints only to the CCC forecasts:

Temporal hierarchies

The extension to the case of temporal hierarchies is quite simple. Using the same notation as in thfrec(), the following `equivalences' hold between the symbols used for the level conditional cross-sectional reconciliation and the ones used in temporal reconciliation: \[L \equiv p-1, \quad (n_a, n_b, n) \equiv (k^*, m, k^*+m),\] and \[\mathbf{C} \equiv \mathbf{K} , \; \mathbf{C}_1 \equiv \mathbf{K}_1 = \mathbf{1}_m', \; \mathbf{C}_2 \equiv \mathbf{K}_2 = \mathbf{I}_{\frac{m}{k_{p-1}}},\; \ldots, \; \mathbf{C}_L \equiv \mathbf{K}_{p-1} = \mathbf{I}_{\frac{m}{k_{2}}} \otimes \mathbf{1}_{k_{2}}',\; \mathbf{S} \equiv \mathbf{R}.\]

The description of the cross-temporal extension is currently under progress.

References

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., Girolimetto, D. (2022), Forecast combination based forecast reconciliation: insights and extensions, International Journal of Forecasting, in press.

Hollyman, R., Petropoulos, F., Tipping, M.E. (2021), Understanding Forecast Reconciliation, European Journal of Operational Research (in press).

Examples

data(FoReco_data)
### LCC and CCC CROSS-SECTIONAL FORECAST RECONCILIATION
# Cross sectional aggregation matrix
C <- rbind(FoReco_data$C, c(0,0,0,0,1))
# monthly base forecasts
mbase <- FoReco2matrix(FoReco_data$base, m = 12)$k1
mbase <- mbase[, c("T", "A","B","C","AA","AB","BA","BB","C")]
# residuals
mres <- FoReco2matrix(FoReco_data$res, m = 12)$k1
mres <- mres[, c("T", "A","B","C","AA","AB","BA","BB","C")]
# covariance matrix of all the base forecasts, computed using the in-sample residuals
Wres <- cov(mres)
# covariance matrix of the bts base forecasts, computed using the in-sample residuals
Wb <- Wres[c("AA","AB", "BA", "BB", "C"),
           c("AA","AB", "BA", "BB", "C")]
# bts seasonal averages
obs_1 <- FoReco_data$obs$k1
bts_sm <- apply(obs_1, 2, function(x) tapply(x[1:168],rep(1:12, 14), mean))
bts_sm <- bts_sm[,c("AA", "AB", "BA", "BB", "C")]

## EXOGENOUS CONSTRAINTS AND DIAGONAL COVARIANCE MATRIX (Hollyman et al., 2021)
# Forecast reconciliation for an elementary hierarchy:
# 1 top-level series + 5 bottom-level series (Level 2 base forecasts are not considered).
# The input is given by the base forecasts of the top and bottom levels series,
# along with a vector of positive weights for the bts base forecasts
exo_EHd <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")],
                 weights = diag(Wb), bnaive = bts_sm)

# Level conditional reconciled forecasts
# recf/L1: Level 1 reconciled forecasts for the whole hierarchy
# L2: Middle-Out reconciled forecasts hinging on Level 2 conditional reconciled forecasts
# L3: Bottom-Up reconciled forecasts
exo_LCd <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = diag(Wb),
                 CCC = FALSE, bnaive = bts_sm)

# Combined Conditional Coherent (CCC) reconciled forecasts
# recf: CCC reconciled forecasts matrix
# L1: Level 1 conditional reconciled forecasts for the whole hierarchy
# L2: Middle-Out reconciled forecasts hinging on Level 2 conditional reconciled forecasts
# L3: Bottom-Up reconciled forecasts
exo_CCCd <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = diag(Wb))

## EXOGENOUS CONSTRAINTS AND FULL COVARIANCE MATRIX
# Simply substitute weights=diag(Wb) with weights=Wb
exo_EHf <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")], weights = Wb, bnaive = bts_sm)
exo_LCf <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = Wb, CCC = FALSE, bnaive = bts_sm)
exo_CCCf <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = Wb, bnaive = bts_sm)

## ENDOGENOUS CONSTRAINTS AND DIAGONAL COVARIANCE MATRIX
# parameters of function htsrec(), like "type" and "nn_type" may be used

# Forecast reconciliation for an elementary hierarchy with endogenous constraints
W1 <- Wres[c("T","AA","AB", "BA", "BB", "C"),
           c("T","AA","AB", "BA", "BB", "C")]
endo_EHd <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")],
                 weights = diag(W1), const = "endo", CCC = FALSE)

# Level conditional reconciled forecasts with endogenous constraints
endo_LCd <- lccrec(basef = mbase, C = C, nl = c(1,3), weights = diag(Wres),
                  const = "endo", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts with endogenous constraints
endo_CCCd <- lccrec(basef = mbase, C = C, nl = c(1,3),
                    weights = diag(Wres), const = "endo")

## ENDOGENOUS CONSTRAINTS AND FULL COVARIANCE MATRIX
# Simply substitute weights=diag(Wres) with weights=Wres
endo_EHf <- lccrec(basef = mbase[, c("T","AA","AB", "BA", "BB", "C")],
                   weights = W1,
                   const = "endo")
endo_LCf <- lccrec(basef = mbase, C = C, nl = c(1,3),
                   weights = Wres, const = "endo", CCC = FALSE)
endo_CCCf <- lccrec(basef = mbase-40, C = C, nl = c(1,3),
                   weights = Wres, const = "endo")

### LCC and CCC TEMPORAL FORECAST RECONCILIATION
# 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, ]
Om_bt <- cov(matrix(topres[which(simplify2array(strsplit(names(topres),
            "_"))[1,]=="k1")], ncol = 12, byrow = TRUE))
t_exo_LCd <- lccrec(basef = topbase, m = 12, weights = diag(Om_bt), CCC = FALSE)
t_exo_CCCd <- lccrec(basef = topbase, m = 12, weights = diag(Om_bt))

### LCC and CCC CROSS-TEMPORAL FORECAST RECONCILIATION
idr <- which(simplify2array(strsplit(colnames(FoReco_data$res), "_"))[1,]=="k1")
bres <- FoReco_data$res[-c(1:3), idr]
bres <- lapply(1:5, function(x) matrix(bres[x,], nrow=14, byrow = TRUE))
bres <- do.call(cbind, bres)
ctbase <- FoReco_data$base[c("T", "A","B","C","AA","AB","BA","BB","C"),]
ct_exo_LCf <- lccrec(basef = ctbase, m = 12, C = C, nl = c(1,3),
                    weights = diag(cov(bres)), CCC = FALSE)
ct_exo_CCCf <- lccrec(basef = ctbase, m = 12, C = C, nl = c(1,3),
                     weights = diag(cov(bres)), CCC = TRUE)