Package 'FoReco' reference manual

Title:	Forecast Reconciliation
Description:	Classical (bottom-up and top-down), optimal combination and heuristic point (Di Fonzo and Girolimetto, 2023 <doi:10.1016/j.ijforecast.2021.08.004>) and probabilistic (Girolimetto et al. 2023 <doi:10.1016/j.ijforecast.2023.10.003>) forecast reconciliation procedures for linearly constrained time series (e.g., hierarchical or grouped time series) in cross-sectional, temporal, or cross-temporal frameworks.
Authors:	Daniele Girolimetto [aut, cre, fnd] , Tommaso Di Fonzo [aut, fnd]
Maintainer:	Daniele Girolimetto <[email protected]>
License:	GPL-3
Version:	1.0.0.9012
Built:	2025-03-08 05:31:54 UTC
Source:	https://github.com/danigiro/foreco

FoReco: Forecast Reconciliation

Description

Classical (bottom-up and top-down), optimal combination and heuristic point (Di Fonzo and Girolimetto, 2023 doi:10.1016/j.ijforecast.2021.08.004) and probabilistic (Girolimetto et al. 2023 doi:10.1016/j.ijforecast.2023.10.003) forecast reconciliation procedures for linearly constrained time series (e.g., hierarchical or grouped time series) in cross-sectional, temporal, or cross-temporal frameworks.

Author(s)

Maintainer: Daniele Girolimetto [email protected] (ORCID) [funder]

Authors:

Tommaso Di Fonzo (ORCID) [funder]

Non-overlapping temporal aggregation of a time series

Description

Non-overlapping temporal aggregation of a time series according to a specific aggregation order.

Usage

aggts(y, agg_order, tew = "sum", align = "end", rm_na = FALSE)
aggts(y, agg_order, tew = "sum", align = "end", rm_na = FALSE)

Arguments

`y`	Univariate or multivariate time series: a vector/matrix or a `ts` object.
`agg_order`	A numeric vector with the aggregation orders to consider.
`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).
`align`	A string or a vector specifying the alignment of `y`. Options include: "`end`" (end of the series, default), "`start`" (start of the series), an integer (or a vector of integers) indicating the starting period of the temporally aggregated series.
`rm_na`	If `TRUE` the missing values are removed.

Value

A list of vectors or ts objects.

Examples

# Monthly time series (input vector)
y <- ts(rnorm(24), start = 2020, frequency = 12)
# Quarterly time series
x1 <- aggts(y, 3)
# Monthly, quarterly and annual time series
x2 <- aggts(y, c(1, 3, 12))
# All temporally aggregated time series
x3 <- aggts(y)

# Ragged data
y2 <- ts(rnorm(11), start = c(2020, 3), frequency = 4)
# Annual time series: start in 2021
x4 <- aggts(y2, 4, align = 3)
# Semi-annual (start in 2nd semester of 2020) and annual (start in 2021) time series
x5 <- aggts(y2, c(2, 4), align = c(1, 3))

# Monthly time series (input vector)
y <- ts(rnorm(24), start = 2020, frequency = 12)
# Quarterly time series
x1 <- aggts(y, 3)
# Monthly, quarterly and annual time series
x2 <- aggts(y, c(1, 3, 12))
# All temporally aggregated time series
x3 <- aggts(y)

# Ragged data
y2 <- ts(rnorm(11), start = c(2020, 3), frequency = 4)
# Annual time series: start in 2021
x4 <- aggts(y2, 4, align = 3)
# Semi-annual (start in 2nd semester of 2020) and annual (start in 2021) time series
x5 <- aggts(y2, c(2, 4), align = c(1, 3))

Aggregation matrix of a (possibly) unbalanced hierarchy in balanced form

Description

A hierarchy with $L$ upper levels is said to be balanced if each variable at level $l$ has at least one child at level $l+1$ . When this doesn't hold, the hierarchy is unbalanced. This function transforms an aggregation matrix of an unbalanced hierarchy into an aggregation matrix of a balanced one. This function is used to reconcile forecasts with cslcc, which operates exclusively with balanced hierarchies.

Usage

balance_hierarchy(agg_mat, nodes = "auto", sparse = TRUE)
balance_hierarchy(agg_mat, nodes = "auto", sparse = TRUE)

Arguments

`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.
`nodes`	A ( $L \times 1$ ) numeric vector indicating the number of variables in each of the upper $L$ levels of the hierarchy. The default value is the string "`auto`" which calculates the number of variables in each level.
`sparse`	Option to return sparse matrices (default is `TRUE`).

Value

A list containing four elements:

`bam`	The balanced aggregation matrix.
`agg_mat`	The input matrix.
`nodes`	A ( $L \times 1$ ) numeric vector indicating the number of variables in each of the $L$ upper levels of the balanced hierarchy.
`id`	The identification number of each variable in the balanced hierarchy. It may contains duplicated values.

Examples

#    Unbalanced     ->      Balanced
#        T                     T
#    |-------|             |-------|
#    A       |             A       B
#  |---|     |           |---|     |
# AA   AB    B          AA   AB    BA
A <- matrix(c(1, 1, 1,
              1, 1, 0), 2, byrow = TRUE)
obj <- balance_hierarchy(agg_mat = A, nodes = c(1, 1))
obj$bam
#    Unbalanced     ->      Balanced
#        T                     T
#    |-------|             |-------|
#    A       |             A       B
#  |---|     |           |---|     |
# AA   AB    B          AA   AB    BA
A <- matrix(c(1, 1, 1,
              1, 1, 0), 2, byrow = TRUE)
obj <- balance_hierarchy(agg_mat = A, nodes = c(1, 1))
obj$bam

Commutation matrix

Description

This function returns the ( $r c \times r c$ ) commutation matrix $\mathbf{P}$ such that $\mathbf{P} \mbox{vec}(\mathbf{Y}) = \mbox{vec}(\mathbf{Y}'),$ where $\mathbf{Y}$ is a ( $r \times c$ ) matrix (Magnus and Neudecker, 2019).

Usage

commat(r, c)
commat(r, c)

Arguments

`r`	Number of rows of $\mathbf{Y}$ .
`c`	Number of columns of $\mathbf{Y}$ .

Value

A sparse ( $r c \times r c$ ) matrix, $\mathbf{P}$ .

References

Magnus, J.R. and Neudecker, H. (2019), Matrix Differential Calculus with Applications in Statistics and Econometrics, third edition, New York, Wiley, pp. 54-55.

Examples

Y <- matrix(rnorm(30), 5, 6)
P <- commat(5, 6)
P %*% as.vector(Y) == as.vector(t(Y)) # check

Y <- matrix(rnorm(30), 5, 6)
P <- commat(5, 6)
P %*% as.vector(Y) == as.vector(t(Y)) # check

Cross-sectional joint block bootstrap

Description

Joint block bootstrap for generating probabilistic base forecasts that take into account the correlation between different time series (Panagiotelis et al. 2023).

Usage

csboot(model_list, boot_size, block_size, seed = NULL)
csboot(model_list, boot_size, block_size, seed = NULL)

Arguments

`model_list`	A list of all the $n$ base forecasts models. A `simulate()` function for each model has to be available and implemented according to the package forecast, with the following mandatory parameters: object, innov, future, and nsim.
`boot_size`	The number of bootstrap replicates.
`block_size`	Block size of the bootstrap, which is typically equivalent to the forecast horizon.
`seed`	An integer seed.

Value

A list with two elements: the seed used to sample the errors and a 3-d array ( $\text{boot\_size}\times n \times \text{block\_size}$ ).

References

Panagiotelis, A., Gamakumara, P., Athanasopoulos, G. and Hyndman, R.J. (2023), Probabilistic forecast reconciliation: Properties, evaluation and score optimisation, European Journal of Operational Research 306(2), 693–706. doi:10.1016/j.ejor.2022.07.040

Cross-sectional bottom-up reconciliation

Description

This function computes the cross-sectional bottom-up reconciled forecasts (Dunn et al., 1976) for all series by appropriate summation of the bottom base forecasts $\widehat{\mathbf{b}}$ :

$\widetilde{\mathbf{y}} = \mathbf{S}_{cs}\widehat{\mathbf{b}},$

where $\mathbf{S}_{cs}$ is the cross-sectional structural matrix.

Usage

csbu(base, agg_mat, sntz = FALSE)
csbu(base, agg_mat, sntz = FALSE)

Arguments

`base`	A ( $h \times n_b$ ) numeric matrix or multivariate time series (`mts` class) containing bottom base forecasts; $h$ is the forecast horizon, and $n_b$ is the total number of bottom variables.
`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.
`sntz`	If `TRUE`, the negative base forecasts are set to zero before applying bottom-up.

Value

A ( $h \times n$ ) numeric matrix of cross-sectional reconciled forecasts.

References

Dunn, D. M., Williams, W. H. and Dechaine, T. L. (1976), Aggregate versus subaggregate models in local area forecasting, Journal of the American Statistical Association 71(353), 68–71. doi:10.1080/01621459.1976.10481478

Examples

set.seed(123)
# (3 x 2) bottom base forecasts matrix (simulated), Z = X + Y
bts <- matrix(rnorm(6, mean = c(10, 10)), 3, byrow = TRUE)

# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
reco <- csbu(base = bts, agg_mat = A)

# Non negative reconciliation
bts[2,2] <- -bts[2,2] # Making negative one of the base forecasts for variable Y
nnreco <- csbu(base = bts, agg_mat = A, sntz = TRUE)

set.seed(123)
# (3 x 2) bottom base forecasts matrix (simulated), Z = X + Y
bts <- matrix(rnorm(6, mean = c(10, 10)), 3, byrow = TRUE)

# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
reco <- csbu(base = bts, agg_mat = A)

# Non negative reconciliation
bts[2,2] <- -bts[2,2] # Making negative one of the base forecasts for variable Y
nnreco <- csbu(base = bts, agg_mat = A, sntz = TRUE)

Cross-sectional covariance matrix approximation

Description

This function provides an approximation of the cross-sectional base forecasts errors covariance matrix using different reconciliation methods (see Wickramasuriya et al., 2019 and Di Fonzo and Girolimetto, 2023).

Usage

cscov(comb = "ols", n = NULL, agg_mat = NULL, res, mse = TRUE,
      shrink_fun = shrink_estim, ...)
cscov(comb = "ols", n = NULL, agg_mat = NULL, res, mse = TRUE,
      shrink_fun = shrink_estim, ...)

Arguments

`comb`	A string specifying the reconciliation method. Ordinary least squares: "`ols`" (default) - identity error covariance matrix. Weighted least squares: "`str`" - structural variances. "`wls`" - series variances (uses `res`). Generalized least squares (uses `res`): "`shr`" - shrunk covariance (Wickramasuriya et al., 2019). "`oasd`" - oracle shrunk covariance (Ando and Xiao, 2023). "`sam`" - sample covariance.
`n`	Number of variables ( $n = n_a + n_b$ ).
`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.
`res`	An ( $N \times n$ ) optional numeric matrix containing the in-sample residuals. This matrix is used to compute some covariance matrices.
`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).
`...`	Not used.

Value

A ( $n \times n$ ) symmetric positive (semi-)definite matrix.

References

Ando, S., and Xiao, M. (2023), High-dimensional covariance matrix estimation: shrinkage toward a diagonal target. IMF Working Papers, 2023(257), A001. doi:10.5089/9798400260780.001.A001

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. doi:10.1016/j.ijforecast.2021.08.004

Wickramasuriya, S.L., Athanasopoulos, G. and 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. doi:10.1080/01621459.2018.1448825

Examples

# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
# (10 x 3) in-sample residuals matrix (simulated)
res <- t(matrix(rnorm(n = 30), nrow = 3))

cov1 <- cscov("ols", n = 3)          # OLS methods
cov2 <- cscov("str", agg_mat = A)    # STR methods
cov3 <- cscov("wls", res = res)      # WLS methods
cov4 <- cscov("shr", res = res)      # SHR methods
cov5 <- cscov("sam", res = res)      # SAM methods

# Custom covariance matrix
cscov.ols2 <- function(comb, x) diag(x)
cscov(comb = "ols2", x = 3) # == cscov("ols", n = 3)

# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
# (10 x 3) in-sample residuals matrix (simulated)
res <- t(matrix(rnorm(n = 30), nrow = 3))

cov1 <- cscov("ols", n = 3)          # OLS methods
cov2 <- cscov("str", agg_mat = A)    # STR methods
cov3 <- cscov("wls", res = res)      # WLS methods
cov4 <- cscov("shr", res = res)      # SHR methods
cov5 <- cscov("sam", res = res)      # SAM methods

# Custom covariance matrix
cscov.ols2 <- function(comb, x) diag(x)
cscov(comb = "ols2", x = 3) # == cscov("ols", n = 3)

Level conditional coherent reconciliation for genuine hierarchical/grouped time series

Description

This function implements the cross-sectional forecast reconciliation procedure that extends the original proposal by Hollyman et al. (2021). Level conditional coherent reconciled forecasts are conditional on (i.e., constrained by) the base forecasts of a specific upper level in the hierarchy (exogenous constraints). It also allows handling the linear constraints linking the variables endogenously (Di Fonzo and Girolimetto, 2022). The function can calculate Combined Conditional Coherent (CCC) forecasts as simple averages of Level-Conditional Coherent (LCC) and bottom-up reconciled forecasts, with either endogenous or exogenous constraints.

Usage

cslcc(base, agg_mat, nodes = "auto", comb = "ols", res = NULL, CCC = TRUE,
      const = "exogenous", bts = NULL, approach = "proj", nn = NULL,
      settings = NULL, ...)
cslcc(base, agg_mat, nodes = "auto", comb = "ols", res = NULL, CCC = TRUE,
      const = "exogenous", bts = NULL, approach = "proj", nn = NULL,
      settings = NULL, ...)

Arguments

`base`	A ( $h \times n$ ) numeric matrix or multivariate time series (`mts` class) containing the base forecasts to be reconciled; $h$ is the forecast horizon, and $n$ is the total number of time series ( $n = n_a + n_b$ ).
`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.
`nodes`	A ( $L \times 1$ ) numeric vector indicating the number of variables in each of the upper $L$ levels of the hierarchy. The default value is the string "`auto`" which calculates the number of variables in each level.
`comb`	A string specifying the reconciliation method. For a complete list, see cscov.
`res`	An ( $N \times n$ ) optional numeric matrix containing the in-sample residuals. This matrix is used to compute some covariance matrices.
`CCC`	A logical value indicating whether the Combined Conditional Coherent reconciled forecasts reconciliation should include bottom-up forecasts (`TRUE`, default), or not.
`const`	A string specifying the reconciliation constraints: "`exogenous`" (default): Fixes the top level of each sub-hierarchy. "`endogenous`": Coherently revises both the top and bottom levels.
`bts`	A ( $h \times n_b$ ) numeric matrix or multivariate time series (`mts` class) containing bottom base forecasts defined by the user (e.g., seasonal averages, as in Hollyman et al., 2021). This parameter can be omitted if only base forecasts are used (see Di Fonzo and Girolimetto, 2024).
`approach`	A string specifying the approach used to compute the reconciled forecasts. Options include: "`proj`" (default): Projection approach according to Byron (1978, 1979). "`strc`": Structural approach as proposed by Hyndman et al. (2011). "`proj_osqp`": Numerical solution using osqp for projection approach. "`strc_osqp`": Numerical solution using osqp for structural approach.
`nn`	A string specifying the algorithm to compute non-negative forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`bpv`": block principal pivoting algorithm. "`sntz`": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
`settings`	A list of control parameters. `nn = "osqp"` An object of class `osqpSettings` specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020). `nn = "bpv"` It includes: `ptype` for permutation method ("`random`" or "`fixed`", default), `par` for the number of full exchange rules that may be attempted (`10`, default), `tol` for the tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `gtol` for the gradient tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `itmax` for the maximum number of algorithm iterations (`100`, default)
`...`	Arguments passed on to `cscov` `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).

Value

A ( $h \times n$ ) numeric matrix of cross-sectional reconciled forecasts.

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. (2024), Forecast combination-based forecast reconciliation: Insights and extensions, International Journal of Forecasting, 40(2), 490–514. doi:10.1016/j.ijforecast.2022.07.001

Di Fonzo, T. and Girolimetto, D. (2023b) Spatio-temporal reconciliation of solar forecasts. Solar Energy 251, 13–29. doi:10.1016/j.solener.2023.01.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

Hollyman, R., Petropoulos, F. and Tipping, M.E. (2021), Understanding forecast reconciliation. European Journal of Operational Research, 294, 149–160. doi:10.1016/j.ejor.2021.01.017

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)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (2 x 7) base forecasts matrix (simulated)
base <- matrix(rnorm(7*2, mean = c(40, 20, 20, 10, 10, 10, 10)), 2, byrow = TRUE)
# (10 x 7) in-sample residuals matrix (simulated)
res <- matrix(rnorm(n = 7*10), ncol = 7)
# (2 x 7) Naive bottom base forecasts matrix: all forecasts are set equal to 10
naive <- matrix(10, 2, 4)

## EXOGENOUS CONSTRAINTS (Hollyman et al., 2021)
# Level Conditional Coherent (LCC) reconciled forecasts
exo_LC <- cslcc(base = base, agg_mat = A, comb = "wls", bts = naive,
                res = res, nodes = "auto", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts
exo_CCC <- cslcc(base = base, agg_mat = A, comb = "wls", bts = naive,
                 res = res, nodes = "auto", CCC = TRUE)

# Results detailed by level:
# L-1: Level 1 immutable reconciled forecasts for the whole hierarchy
# L-2: Middle-Out reconciled forecasts
# L-3: Bottom-Up reconciled forecasts
info_exo <- recoinfo(exo_CCC, verbose = FALSE)
info_exo$lcc

## ENDOGENOUS CONSTRAINTS (Di Fonzo and Girolimetto, 2024)
# Level Conditional Coherent (LCC) reconciled forecasts
endo_LC <- cslcc(base = base, agg_mat = A, comb = "wls",
                 res = res, nodes = "auto", CCC = FALSE,
                 const = "endogenous")

# Combined Conditional Coherent (CCC) reconciled forecasts
endo_CCC <- cslcc(base = base, agg_mat = A, comb = "wls",
                  res = res, nodes = "auto", CCC = TRUE,
                  const = "endogenous")

# Results detailed by level:
# L-1: Level 1 reconciled forecasts for L1 + L3 (bottom level)
# L-2: Level 2 reconciled forecasts for L2 + L3 (bottom level)
# L-3: Bottom-Up reconciled forecasts
info_endo <- recoinfo(endo_CCC, verbose = FALSE)
info_endo$lcc

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (2 x 7) base forecasts matrix (simulated)
base <- matrix(rnorm(7*2, mean = c(40, 20, 20, 10, 10, 10, 10)), 2, byrow = TRUE)
# (10 x 7) in-sample residuals matrix (simulated)
res <- matrix(rnorm(n = 7*10), ncol = 7)
# (2 x 7) Naive bottom base forecasts matrix: all forecasts are set equal to 10
naive <- matrix(10, 2, 4)

## EXOGENOUS CONSTRAINTS (Hollyman et al., 2021)
# Level Conditional Coherent (LCC) reconciled forecasts
exo_LC <- cslcc(base = base, agg_mat = A, comb = "wls", bts = naive,
                res = res, nodes = "auto", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts
exo_CCC <- cslcc(base = base, agg_mat = A, comb = "wls", bts = naive,
                 res = res, nodes = "auto", CCC = TRUE)

# Results detailed by level:
# L-1: Level 1 immutable reconciled forecasts for the whole hierarchy
# L-2: Middle-Out reconciled forecasts
# L-3: Bottom-Up reconciled forecasts
info_exo <- recoinfo(exo_CCC, verbose = FALSE)
info_exo$lcc

## ENDOGENOUS CONSTRAINTS (Di Fonzo and Girolimetto, 2024)
# Level Conditional Coherent (LCC) reconciled forecasts
endo_LC <- cslcc(base = base, agg_mat = A, comb = "wls",
                 res = res, nodes = "auto", CCC = FALSE,
                 const = "endogenous")

# Combined Conditional Coherent (CCC) reconciled forecasts
endo_CCC <- cslcc(base = base, agg_mat = A, comb = "wls",
                  res = res, nodes = "auto", CCC = TRUE,
                  const = "endogenous")

# Results detailed by level:
# L-1: Level 1 reconciled forecasts for L1 + L3 (bottom level)
# L-2: Level 2 reconciled forecasts for L2 + L3 (bottom level)
# L-3: Bottom-Up reconciled forecasts
info_endo <- recoinfo(endo_CCC, verbose = FALSE)
info_endo$lcc

Cross-sectional middle-out reconciliation

Description

The middle-out forecast reconciliation (Athanasopoulos et al., 2009) combines top-down (cstd) and bottom-up (csbu) for genuine hierarchical/grouped time series. Given the base forecasts of variables at an intermediate level $l$ , it performs

a top-down approach for the levels $<l$ ;
a bottom-up approach for the levels $>l$ .

Usage

csmo(base, agg_mat, id_rows = 1, weights, normalize = TRUE)
csmo(base, agg_mat, id_rows = 1, weights, normalize = TRUE)

Arguments

`base`	A ( $h \times n_l$ ) numeric matrix containing the $l$ -level base forecast; $n_l$ is the number of variables at level $l$ , and $h$ is the forecast horizon.
`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.
`id_rows`	A numeric vector indicating the $l$ -level rows of `agg_mat`.
`weights`	A ( $h \times n_b$ ) numeric matrix containing the proportions for the bottom time series; $h$ is the forecast horizon, and $n_b$ is the total number of bottom variables.
`normalize`	If `TRUE` (default), the `weights` will sum to 1.

Value

A ( $h \times n$ ) numeric matrix of cross-sectional reconciled forecasts.

References

Athanasopoulos, G., Ahmed, R. A. and Hyndman, R.J. (2009) Hierarchical forecasts for Australian domestic tourism. International Journal of Forecasting 25(1), 146–166. doi:10.1016/j.ijforecast.2008.07.004

Examples

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (3 x 2) top base forecasts vector (simulated), forecast horizon = 3
baseL2 <- matrix(rnorm(2*3, 5), 3, 2)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- csmo(base = baseL2, agg_mat = A, id_rows = 2:3, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*3), 3, 4)
recoh <- csmo(base = baseL2, agg_mat = A, id_rows = 2:3, weights = h_weights)

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (3 x 2) top base forecasts vector (simulated), forecast horizon = 3
baseL2 <- matrix(rnorm(2*3, 5), 3, 2)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- csmo(base = baseL2, agg_mat = A, id_rows = 2:3, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*3), 3, 4)
recoh <- csmo(base = baseL2, agg_mat = A, id_rows = 2:3, weights = h_weights)

Projection matrix for optimal combination cross-sectional reconciliation

Description

This function computes the projection or the mapping matrix $\mathbf{M}$ and $\mathbf{G}$ , respectively, such that $\widetilde{\mathbf{y}} = \mathbf{M}\widehat{\mathbf{y}} = \mathbf{S}_{cs}\mathbf{G}\widehat{\mathbf{y}}$ , where $\widetilde{\mathbf{y}}$ is the vector of the reconciled forecasts, $\widehat{\mathbf{y}}$ is the vector of the base forecasts, $\mathbf{S}_{cs}$ is the cross-sectional structural matrix, and $\mathbf{M} = \mathbf{S}_{cs}\mathbf{G}$ . For further information regarding on the structure of these matrices, refer to Girolimetto et al. (2023).

Usage

csprojmat(agg_mat, cons_mat, comb = "ols", res = NULL, mat = "M", ...)
csprojmat(agg_mat, cons_mat, comb = "ols", res = NULL, mat = "M", ...)

Arguments

`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: each row represents a constraint equation, and each column represents a variable. The matrix can be of full rank, meaning the rows are linearly independent, but this is not a strict requirement, as the function allows for redundancy in the constraints.
`comb`	A string specifying the reconciliation method. For a complete list, see cscov.
`res`	An ( $N \times n$ ) optional numeric matrix containing the in-sample residuals. This matrix is used to compute some covariance matrices.
`mat`	A string specifying which matrix to return: "`M`" (default) for $\mathbf{M}$ and "`G`" for $\mathbf{G}$ .
`...`	Arguments passed on to `cscov` `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).

Value

The projection matrix $\mathbf{M}$ (mat = "M") or the mapping matrix $\mathbf{G}$ (mat = "G").

References

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

Examples

# Cross-sectional framework
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
Mcs <- csprojmat(agg_mat = A, comb = "ols")
Gcs <- csprojmat(agg_mat = A, comb = "ols", mat = "G")

# Cross-sectional framework
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
Mcs <- csprojmat(agg_mat = A, comb = "ols")
Gcs <- csprojmat(agg_mat = A, comb = "ols", mat = "G")

Optimal combination cross-sectional reconciliation

Description

This function performs optimal (in least squares sense) combination cross-sectional forecast reconciliation for a linearly constrained (e.g., hierarchical/grouped) multiple time series (Wickramasuriya et al., 2019, Panagiotelis et al., 2022, Girolimetto and Di Fonzo, 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 (including Zhang et al., 2023) reconciled forecasts can be considered.

Usage

csrec(base, agg_mat, cons_mat, comb = "ols", res = NULL, approach = "proj",
      nn = NULL, settings = NULL, bounds = NULL, immutable = NULL, ...)
csrec(base, agg_mat, cons_mat, comb = "ols", res = NULL, approach = "proj",
      nn = NULL, settings = NULL, bounds = NULL, immutable = NULL, ...)

Arguments

`base`	A ( $h \times n$ ) numeric matrix or multivariate time series (`mts` class) containing the base forecasts to be reconciled; $h$ is the forecast horizon, and $n$ is the total number of time series ( $n = n_a + n_b$ ).
`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: each row represents a constraint equation, and each column represents a variable. The matrix can be of full rank, meaning the rows are linearly independent, but this is not a strict requirement, as the function allows for redundancy in the constraints.
`comb`	A string specifying the reconciliation method. For a complete list, see cscov.
`res`	An ( $N \times n$ ) optional numeric matrix containing the in-sample residuals. This matrix is used to compute some covariance matrices.
`approach`	A string specifying the approach used to compute the reconciled forecasts. Options include: "`proj`" (default): Projection approach according to Byron (1978, 1979). "`strc`": Structural approach as proposed by Hyndman et al. (2011). "`proj_osqp`": Numerical solution using osqp for projection approach. "`strc_osqp`": Numerical solution using osqp for structural approach.
`nn`	A string specifying the algorithm to compute non-negative forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`bpv`": block principal pivoting algorithm. "`sntz`": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
`settings`	A list of control parameters. `nn = "osqp"` An object of class `osqpSettings` specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020). `nn = "bpv"` It includes: `ptype` for permutation method ("`random`" or "`fixed`", default), `par` for the number of full exchange rules that may be attempted (`10`, default), `tol` for the tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `gtol` for the gradient tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `itmax` for the maximum number of algorithm iterations (`100`, default)
`bounds`	A matrix (see set_bounds) with 3 columns ( $i,lower,upper$ ), such that Column 1 represents the cross-sectional series ( $i = 1, \dots, n$ ). Columns 2 and 3 indicates the lower and lower bounds, respectively.
`immutable`	A numeric vector containing the column indices of the base forecasts (`base` parameter) that should be fixed.
`...`	Arguments passed on to `cscov` `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).

Value

A ( $h \times n$ ) numeric matrix of cross-sectional reconciled forecasts.

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. (2023), Spatio-temporal reconciliation of solar forecasts, Solar Energy, 251, 13–29. doi:10.1016/j.solener.2023.01.003

Girolimetto, D. and Di Fonzo, T. (2023), Point and probabilistic forecast reconciliation for general linearly constrained multiple time series, Statistical Methods & Applications, 33, 581-607. doi:10.1007/s10260-023-00738-6.

Panagiotelis, A., Athanasopoulos, G., Gamakumara, P. and Hyndman, R.J. (2021), Forecast reconciliation: A geometric view with new insights on bias correction, International Journal of Forecasting, 37, 1, 343–359. doi:10.1016/j.ijforecast.2020.06.004

Zhang, B., Kang, Y., Panagiotelis, A. and Li, F. (2023), Optimal reconciliation with immutable forecasts, European Journal of Operational Research, 308(2), 650–660. doi:10.1016/j.ejor.2022.11.035

Examples

set.seed(123)
# (2 x 3) base forecasts matrix (simulated), Z = X + Y
base <- matrix(rnorm(6, mean = c(20, 10, 10)), 2, byrow = TRUE)
# (10 x 3) in-sample residuals matrix (simulated)
res <- t(matrix(rnorm(n = 30), nrow = 3))

# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
reco <- csrec(base = base, agg_mat = A, comb = "wls", res = res)

# Zero constraints matrix for Z - X - Y = 0
C <- t(c(1, -1, -1))
reco <- csrec(base = base, cons_mat = C, comb = "wls", res = res) # same results

# Non negative reconciliation
base[1,3] <- -base[1,3] # Making negative one of the base forecasts for variable Y
nnreco <- csrec(base = base, agg_mat = A, comb = "wls", res = res, nn = "osqp")
recoinfo(nnreco, verbose = FALSE)$info

set.seed(123)
# (2 x 3) base forecasts matrix (simulated), Z = X + Y
base <- matrix(rnorm(6, mean = c(20, 10, 10)), 2, byrow = TRUE)
# (10 x 3) in-sample residuals matrix (simulated)
res <- t(matrix(rnorm(n = 30), nrow = 3))

# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
reco <- csrec(base = base, agg_mat = A, comb = "wls", res = res)

# Zero constraints matrix for Z - X - Y = 0
C <- t(c(1, -1, -1))
reco <- csrec(base = base, cons_mat = C, comb = "wls", res = res) # same results

# Non negative reconciliation
base[1,3] <- -base[1,3] # Making negative one of the base forecasts for variable Y
nnreco <- csrec(base = base, agg_mat = A, comb = "wls", res = res, nn = "osqp")
recoinfo(nnreco, verbose = FALSE)$info

Cross-sectional top-down reconciliation

Description

Top-down forecast reconciliation for genuine hierarchical/grouped time series (Gross and Sohl, 1990), where the forecast of a ‘Total’ (top-level series, expected to be positive) is disaggregated according to a proportional scheme (weights). Besides fulfilling any aggregation constraint, the top-down reconciled forecasts should respect two main properties:

the top-level value remains unchanged;
all the bottom time series reconciled forecasts are non-negative.

Usage

cstd(base, agg_mat, weights, normalize = TRUE)
cstd(base, agg_mat, weights, normalize = TRUE)

Arguments

`base`	A ( $h \times 1$ ) numeric vector containing the top-level base forecast; $h$ is the forecast horizon.
`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.
`weights`	A ( $h \times n_b$ ) numeric matrix containing the proportions for the bottom time series; $h$ is the forecast horizon, and $n_b$ is the total number of bottom variables.
`normalize`	If `TRUE` (default), the `weights` will sum to 1.

Value

A ( $h \times n$ ) numeric matrix of cross-sectional reconciled forecasts.

References

Gross, C.W. and Sohl, J.E. (1990), Disaggregation methods to expedite product line forecasting. Journal of Forecasting 9(3), 233–254. doi:10.1002/for.3980090304

Examples

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (3 x 1) top base forecasts vector (simulated), forecast horizon = 3
topf <- rnorm(3, 10)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- cstd(base = topf, agg_mat = A, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*3), 3, 4)
recoh <- cstd(base = topf, agg_mat = A, weights = h_weights)

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (3 x 1) top base forecasts vector (simulated), forecast horizon = 3
topf <- rnorm(3, 10)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- cstd(base = topf, agg_mat = A, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*3), 3, 4)
recoh <- cstd(base = topf, agg_mat = A, weights = h_weights)

Cross-sectional reconciliation tools

Description

Some useful tools for the cross-sectional forecast reconciliation of a linearly constrained (e.g., hierarchical/grouped) multiple time series.

Usage

cstools(agg_mat, cons_mat, sparse = TRUE)
cstools(agg_mat, cons_mat, sparse = TRUE)

Arguments

`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: each row represents a constraint equation, and each column represents a variable. The matrix can be of full rank, meaning the rows are linearly independent, but this is not a strict requirement, as the function allows for redundancy in the constraints.
`sparse`	Option to return sparse matrices (default is `TRUE`).

Value

A list with four elements:

`dim`	A vector containing information about the number of series for the complete system (`n`), for upper levels (`na`) and bottom level (`nb`).
`agg_mat`	The cross-sectional aggregation matrix.
`strc_mat`	The cross-sectional structural matrix.
`cons_mat`	The cross-sectional zero constraints matrix.

Examples

# Cross-sectional framework
# One level hierarchy A = [1 1]
A <- matrix(1, 1, 2)
obj <- cstools(agg_mat = A)

# Cross-sectional framework
# One level hierarchy A = [1 1]
A <- matrix(1, 1, 2)
obj <- cstools(agg_mat = A)

Cross-temporal joint block bootstrap

Description

Joint block bootstrap for generating probabilistic base forecasts that take into account the correlation between variables at different temporal aggregation orders (Girolimetto et al. 2023).

Usage

ctboot(model_list, boot_size, agg_order, block_size = 1, seed = NULL)
ctboot(model_list, boot_size, agg_order, block_size = 1, seed = NULL)

Arguments

`model_list`	A list of $n$ elements, one for each cross-sectional series. Each elements is a list with the $(k^\ast+m)$ base forecasts models ordered from the lowest frequency (most temporally aggregated) to the highest frequency. A `simulate()` function for each model has to be available and implemented according to the package forecast, with the following mandatory parameters: object, innov, future, and nsim.
`boot_size`	The number of bootstrap replicates.
`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$ .
`block_size`	Block size of the bootstrap, which is typically equivalent to the forecast horizon for the most temporally aggregated series.
`seed`	An integer seed.

Value

A list with two elements: the seed used to sample the errors and a ( $\text{boot\_size}\times n(k^\ast+m)\text{block\_size}$ ) matrix.

References

Girolimetto, D., Athanasopoulos, G., Di Fonzo, T. and Hyndman, R.J. (2023), 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

Cross-temporal bottom-up reconciliation

Description

Cross-temporal bottom-up reconciled forecasts for all series at any temporal aggregation level are computed by appropriate summation of the high-frequency bottom base forecasts $\widehat{\mathbf{B}^{[1]}}$ :

$\widetilde{\mathbf{X}} = \mathbf{S}_{cs}\widehat{\mathbf{B}^{[1]}}\mathbf{S}'_{te},$

where $\mathbf{S}_{cs}$ and $\mathbf{S}_{te}$ are the cross-sectional and temporal structural matrices, respectively.

Usage

ctbu(base, agg_mat, agg_order, tew = "sum", sntz = FALSE)
ctbu(base, agg_mat, agg_order, tew = "sum", sntz = FALSE)

Arguments

`base`	A ( $n_b \times hm$ ) numeric matrix containing high-frequency bottom base forecasts; $n_b$ is the total number of high-frequency bottom variables, $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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.
`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$ .
`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).
`sntz`	If `TRUE`, the negative base forecasts are set to zero before applying bottom-up.

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

Examples

set.seed(123)
# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
# (2 x 4) high frequency bottom base forecasts matrix (simulated),
# agg_order = 4 (annual-quarterly)
hfbts <- matrix(rnorm(4*2, 2.5), 2, 4)

reco <- ctbu(base = hfbts, agg_mat = A, agg_order = 4)

# Non negative reconciliation
hfbts[1,4] <- -hfbts[1,4] # Making negative one of the quarterly base forecasts for variable X
nnreco <- ctbu(base = hfbts, agg_mat = A, agg_order = 4, sntz = TRUE)

set.seed(123)
# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
# (2 x 4) high frequency bottom base forecasts matrix (simulated),
# agg_order = 4 (annual-quarterly)
hfbts <- matrix(rnorm(4*2, 2.5), 2, 4)

reco <- ctbu(base = hfbts, agg_mat = A, agg_order = 4)

# Non negative reconciliation
hfbts[1,4] <- -hfbts[1,4] # Making negative one of the quarterly base forecasts for variable X
nnreco <- ctbu(base = hfbts, agg_mat = A, agg_order = 4, sntz = TRUE)

Cross-temporal covariance matrix approximation

Description

This function provides an approximation of the cross-temporal base forecasts errors covariance matrix using different reconciliation methods (Di Fonzo and Girolimetto, 2023, Girolimetto et al., 2023).

Usage

ctcov(comb = "ols", n = NULL, agg_mat = NULL, agg_order = NULL, res = NULL,
      tew = "sum", mse = TRUE, shrink_fun = shrink_estim, ...)
ctcov(comb = "ols", n = NULL, agg_mat = NULL, agg_order = NULL, res = NULL,
      tew = "sum", mse = TRUE, shrink_fun = shrink_estim, ...)

Arguments

`comb`	A string specifying the reconciliation method. Ordinary least squares: "`ols`" (default) - identity error covariance. Weighted least squares: "`str`" - structural variances. "`csstr`" - cross-sectional structural variances. "`testr`" - temporal structural variances. "`wlsh`" - hierarchy variances (uses `res`). "`wlsv`" - series variances (uses `res`). Generalized least squares (uses `res`): "`acov`" - series auto-covariance. "`bdshr`"/"`bdsam`" - shrunk/sample block diagonal cross-sectional covariance. "`Sshr`"/"`Ssam`" - series shrunk/sample covariance. "`shr`"/"`sam`" - shrunk/sample covariance. "`hbshr`"/"`hbsam`" - shrunk/sample high frequency bottom time series covariance. "`bshr`"/"`bsam`" - shrunk/sample bottom time series covariance. "`hshr`"/"`hsam`" - shrunk/sample high frequency covariance.
`n`	Cross-sectional number of variables.
`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.
`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$ .
`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).
`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).
`...`	Not used.

Value

A ( $n(k^\ast+m) \times n(k^\ast+m)$ ) symmetric matrix.

References

Examples

set.seed(123)
# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
# (3 x 70) in-sample residuals matrix (simulated),
# agg_order = 4 (annual-quarterly)
res <- rbind(rnorm(70), rnorm(70), rnorm(70))

cov1 <- ctcov("ols", n = 3, agg_order = 4)                     # OLS methods
cov2 <- ctcov("str", agg_mat = A, agg_order = 4)               # STR methods
cov3 <- ctcov("csstr", agg_mat = A, agg_order = 4)             # CSSTR methods
cov4 <- ctcov("testr", n = 3, agg_order = 4)                   # TESTR methods
cov5 <- ctcov("wlsv", agg_order = 4, res = res)                # WLSv methods
cov6 <- ctcov("wlsh", agg_order = 4, res = res)                # WLSh methods
cov7 <- ctcov("shr", agg_order = 4, res = res)                 # SHR methods
cov8 <- ctcov("sam", agg_order = 4, res = res)                 # SAM methods
cov9 <- ctcov("acov", agg_order = 4, res = res)                # ACOV methods
cov10 <- ctcov("Sshr", agg_order = 4, res = res)               # Sshr methods
cov11 <- ctcov("Ssam", agg_order = 4, res = res)               # Ssam methods
cov12 <- ctcov("hshr", agg_order = 4, res = res)               # Hshr methods
cov13 <- ctcov("hsam", agg_order = 4, res = res)               # Hsam methods
cov14 <- ctcov("hbshr", agg_mat = A, agg_order = 4, res = res) # HBshr methods
cov15 <- ctcov("hbsam", agg_mat = A, agg_order = 4, res = res) # HBsam methods
cov16 <- ctcov("bshr", agg_mat = A, agg_order = 4, res = res)  # Bshr methods
cov17 <- ctcov("bsam", agg_mat = A, agg_order = 4, res = res)  # Bsam methods
cov18 <- ctcov("bdshr", agg_order = 4, res = res)              # BDshr methods
cov19 <- ctcov("bdsam", agg_order = 4, res = res)              # BDsam methods

# Custom covariance matrix
ctcov.ols2 <- function(comb, x) diag(x)
cov20 <- ctcov(comb = "ols2", x = 21) # == ctcov("ols", n = 3, agg_order = 4)

set.seed(123)
# Aggregation matrix for Z = X + Y
A <- t(c(1,1))
# (3 x 70) in-sample residuals matrix (simulated),
# agg_order = 4 (annual-quarterly)
res <- rbind(rnorm(70), rnorm(70), rnorm(70))

cov1 <- ctcov("ols", n = 3, agg_order = 4)                     # OLS methods
cov2 <- ctcov("str", agg_mat = A, agg_order = 4)               # STR methods
cov3 <- ctcov("csstr", agg_mat = A, agg_order = 4)             # CSSTR methods
cov4 <- ctcov("testr", n = 3, agg_order = 4)                   # TESTR methods
cov5 <- ctcov("wlsv", agg_order = 4, res = res)                # WLSv methods
cov6 <- ctcov("wlsh", agg_order = 4, res = res)                # WLSh methods
cov7 <- ctcov("shr", agg_order = 4, res = res)                 # SHR methods
cov8 <- ctcov("sam", agg_order = 4, res = res)                 # SAM methods
cov9 <- ctcov("acov", agg_order = 4, res = res)                # ACOV methods
cov10 <- ctcov("Sshr", agg_order = 4, res = res)               # Sshr methods
cov11 <- ctcov("Ssam", agg_order = 4, res = res)               # Ssam methods
cov12 <- ctcov("hshr", agg_order = 4, res = res)               # Hshr methods
cov13 <- ctcov("hsam", agg_order = 4, res = res)               # Hsam methods
cov14 <- ctcov("hbshr", agg_mat = A, agg_order = 4, res = res) # HBshr methods
cov15 <- ctcov("hbsam", agg_mat = A, agg_order = 4, res = res) # HBsam methods
cov16 <- ctcov("bshr", agg_mat = A, agg_order = 4, res = res)  # Bshr methods
cov17 <- ctcov("bsam", agg_mat = A, agg_order = 4, res = res)  # Bsam methods
cov18 <- ctcov("bdshr", agg_order = 4, res = res)              # BDshr methods
cov19 <- ctcov("bdsam", agg_order = 4, res = res)              # BDsam methods

# Custom covariance matrix
ctcov.ols2 <- function(comb, x) diag(x)
cov20 <- ctcov(comb = "ols2", x = 21) # == ctcov("ols", n = 3, agg_order = 4)

Level conditional coherent reconciliation for cross-temporal hierarchies

Description

This function implements a forecast reconciliation procedure inspired by the original proposal by Hollyman et al. (2021) for cross-temporal hierarchies. Level conditional coherent reconciled forecasts are conditional on (i.e., constrained by) the base forecasts of a specific upper level in the hierarchy (exogenous constraints). It also allows handling the linear constraints linking the variables endogenously (Di Fonzo and Girolimetto, 2022). The function can calculate Combined Conditional Coherent (CCC) forecasts as simple averages of Level-Conditional Coherent (LCC) and bottom-up reconciled forecasts, with either endogenous or exogenous constraints.

Usage

ctlcc(base, agg_mat, nodes = "auto", agg_order, comb = "ols", res = NULL,
      CCC = TRUE, const = "exogenous", hfbts = NULL, tew = "sum",
      approach = "proj", nn = NULL, settings = NULL, ...)
ctlcc(base, agg_mat, nodes = "auto", agg_order, comb = "ols", res = NULL,
      CCC = TRUE, const = "exogenous", hfbts = NULL, tew = "sum",
      approach = "proj", nn = NULL, settings = 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.
`nodes`	A ( $L \times 1$ ) numeric vector indicating the number of variables in each of the upper $L$ levels of the hierarchy. The default value is the string "`auto`" which calculates the number of variables in each level.
`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.
`CCC`	A logical value indicating whether the Combined Conditional Coherent reconciled forecasts reconciliation should include bottom-up forecasts (`TRUE`, default), or not.
`const`	A string specifying the reconciliation constraints: "`exogenous`" (default): Fixes the top level of each sub-hierarchy. "`endogenous`": Coherently revises both the top and bottom levels.
`hfbts`	A ( $n \times mh$ ) numeric matrix containing high frequency bottom base forecasts defined by the user. This parameter can be omitted if only base forecasts are used (see Di Fonzo and Girolimetto, 2024).
`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: "`proj`" (default): Projection approach according to Byron (1978, 1979). "`strc`": Structural approach as proposed by Hyndman et al. (2011). "`proj_osqp`": Numerical solution using osqp for projection approach. "`strc_osqp`": Numerical solution using osqp for structural approach.
`nn`	A string specifying the algorithm to compute non-negative forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`bpv`": block principal pivoting algorithm. "`sntz`": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
`settings`	A list of control parameters. `nn = "osqp"` An object of class `osqpSettings` specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020). `nn = "bpv"` It includes: `ptype` for permutation method ("`random`" or "`fixed`", default), `par` for the number of full exchange rules that may be attempted (`10`, default), `tol` for the tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `gtol` for the gradient tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `itmax` for the maximum number of algorithm iterations (`100`, default)
`...`	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).

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

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. (2023b) Spatio-temporal reconciliation of solar forecasts. Solar Energy 251, 13–29. doi:10.1016/j.solener.2023.01.003

Hollyman, R., Petropoulos, F. and Tipping, M.E. (2021), Understanding forecast reconciliation. European Journal of Operational Research, 294, 149–160. doi:10.1016/j.ejor.2021.01.017

Examples

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (7 x 7) base forecasts matrix (simulated), agg_order = 4
base <- rbind(rnorm(7, rep(c(40, 20, 10), c(1, 2, 4))),
              rnorm(7, rep(c(20, 10, 5), c(1, 2, 4))),
              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))),
              rnorm(7, rep(c(10, 5, 2.5), c(1, 2, 4))),
              rnorm(7, rep(c(10, 5, 2.5), c(1, 2, 4))))
# (7 x 70) in-sample residuals matrix (simulated)
res <- matrix(rnorm(70*7), nrow = 7)
# (4 x 4) Naive high frequency bottom base forecasts vector:
# all forecasts are set equal to 2.5
naive <- matrix(2.5, 4, 4)

## EXOGENOUS CONSTRAINTS (Hollyman et al., 2021)
# Level Conditional Coherent (LCC) reconciled forecasts
exo_LC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh", nn = "osqp",
                hfbts = naive, res = res, nodes = "auto", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts
exo_CCC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh",
                hfbts = naive, res = res, nodes = "auto", CCC = TRUE)

# Results detailed by level:
info_exo <- recoinfo(exo_CCC, verbose = FALSE)
# info_exo$lcc

## ENDOGENOUS CONSTRAINTS (Di Fonzo and Girolimetto, 2024)
# Level Conditional Coherent (LCC) reconciled forecasts
endo_LC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh",
                 res = res, nodes = "auto", CCC = FALSE,
                 const = "endogenous")

# Combined Conditional Coherent (CCC) reconciled forecasts
endo_CCC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh",
                  res = res, nodes = "auto", CCC = TRUE,
                  const = "endogenous")

# Results detailed by level:
info_endo <- recoinfo(endo_CCC, verbose = FALSE)
# info_endo$lcc

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (7 x 7) base forecasts matrix (simulated), agg_order = 4
base <- rbind(rnorm(7, rep(c(40, 20, 10), c(1, 2, 4))),
              rnorm(7, rep(c(20, 10, 5), c(1, 2, 4))),
              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))),
              rnorm(7, rep(c(10, 5, 2.5), c(1, 2, 4))),
              rnorm(7, rep(c(10, 5, 2.5), c(1, 2, 4))))
# (7 x 70) in-sample residuals matrix (simulated)
res <- matrix(rnorm(70*7), nrow = 7)
# (4 x 4) Naive high frequency bottom base forecasts vector:
# all forecasts are set equal to 2.5
naive <- matrix(2.5, 4, 4)

## EXOGENOUS CONSTRAINTS (Hollyman et al., 2021)
# Level Conditional Coherent (LCC) reconciled forecasts
exo_LC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh", nn = "osqp",
                hfbts = naive, res = res, nodes = "auto", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts
exo_CCC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh",
                hfbts = naive, res = res, nodes = "auto", CCC = TRUE)

# Results detailed by level:
info_exo <- recoinfo(exo_CCC, verbose = FALSE)
# info_exo$lcc

## ENDOGENOUS CONSTRAINTS (Di Fonzo and Girolimetto, 2024)
# Level Conditional Coherent (LCC) reconciled forecasts
endo_LC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh",
                 res = res, nodes = "auto", CCC = FALSE,
                 const = "endogenous")

# Combined Conditional Coherent (CCC) reconciled forecasts
endo_CCC <- ctlcc(base = base, agg_mat = A, agg_order = 4, comb = "wlsh",
                  res = res, nodes = "auto", CCC = TRUE,
                  const = "endogenous")

# Results detailed by level:
info_endo <- recoinfo(endo_CCC, verbose = FALSE)
# info_endo$lcc

Cross-temporal middle-out reconciliation

Description

The cross-temporal middle-out forecast reconciliation combines top-down (cttd) and bottom-up (ctbu) methods in the cross-temporal framework for genuine hierarchical/grouped time series. Given the base forecasts of an intermediate cross-sectional level $l$ and aggregation order $k$ , it performs

a top-down approach for the aggregation orders $\geq k$ and cross-sectional levels $\geq l$ ;
a bottom-up approach, otherwise.

Usage

ctmo(base, agg_mat, agg_order, id_rows = 1, order = max(agg_order),
     weights, tew = "sum", normalize = TRUE)
ctmo(base, agg_mat, agg_order, id_rows = 1, order = max(agg_order),
     weights, tew = "sum", normalize = TRUE)

Arguments

`base`	A ( $n_l \times hk$ ) numeric matrix containing the $l$ -level base forecasts of temporal aggregation order $k$ ; $n_l$ is the number of variables at level $l$ , $k$ is an aggregation order (a factor of $m$ , and $1<k<m$ ), $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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.
`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$ .
`id_rows`	A numeric vector indicating the $l$ -level rows of `agg_mat`.
`order`	The intermediate fixed aggregation order $k$ .
`weights`	A ( $n_b \times hm$ ) numeric matrix containing the proportions for each high-frequency bottom time series; $n_b$ is the total number of high-frequency bottom variables, $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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).
`normalize`	If `TRUE` (default), the `weights` will sum to 1.

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

Examples

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (2 x 6) base forecasts matrix (simulated), forecast horizon = 3
# and intermediate aggregation order k = 2 (max agg order = 4)
baseL2k2 <- rbind(rnorm(3*2, 5), rnorm(3*2, 5))

# Same weights for different forecast horizons, agg_order = 4
fix_weights <- matrix(runif(4*4), 4, 4)
reco <- ctmo(base = baseL2k2, id_rows = 2:3, agg_mat = A,
             order = 2, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*4*3), 4, 3*4)
recoh <- ctmo(base = baseL2k2, id_rows = 2:3, agg_mat = A,
             order = 2, agg_order = 4, weights = h_weights)

set.seed(123)
# Aggregation matrix for Z = X + Y, X = XX + XY and Y = YX + YY
A <- matrix(c(1,1,1,1,1,1,0,0,0,0,1,1), 3, byrow = TRUE)
# (2 x 6) base forecasts matrix (simulated), forecast horizon = 3
# and intermediate aggregation order k = 2 (max agg order = 4)
baseL2k2 <- rbind(rnorm(3*2, 5), rnorm(3*2, 5))

# Same weights for different forecast horizons, agg_order = 4
fix_weights <- matrix(runif(4*4), 4, 4)
reco <- ctmo(base = baseL2k2, id_rows = 2:3, agg_mat = A,
             order = 2, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*4*3), 4, 3*4)
recoh <- ctmo(base = baseL2k2, id_rows = 2:3, agg_mat = A,
             order = 2, agg_order = 4, weights = h_weights)

Projection matrix for optimal combination cross-temporal reconciliation

Description

This function computes the projection or the mapping matrix $\mathbf{M}$ and $\mathbf{G}$ , respectively, such that $\widetilde{\mathbf{y}} = \mathbf{M}\widehat{\mathbf{y}} = \mathbf{S}_{ct}\mathbf{G}\widehat{\mathbf{y}}$ , where $\widetilde{\mathbf{y}}$ is the vector of the reconciled forecasts, $\widehat{\mathbf{y}}$ is the vector of the base forecasts, $\mathbf{S}_{ct}$ is the cross-temporal structural matrix, and $\mathbf{M} = \mathbf{S}_{ct}\mathbf{G}$ . For further information regarding on the structure of these matrices, refer to Girolimetto et al. (2023).

Usage

ctprojmat(agg_mat, cons_mat, agg_order, comb = "ols", res = NULL,
          mat = "M", tew = "sum", ...)
ctprojmat(agg_mat, cons_mat, agg_order, comb = "ols", res = NULL,
          mat = "M", tew = "sum", ...)

Arguments

`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: each row represents a constraint equation, and each column represents a variable. The matrix can be of full rank, meaning the rows are linearly independent, but this is not a strict requirement, as the function allows for redundancy in the constraints.
`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.
`mat`	A string specifying which matrix to return: "`M`" (default) for $\mathbf{M}$ and "`G`" for $\mathbf{G}$ .
`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).
`...`	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).

Value

The projection matrix $\mathbf{M}$ (mat = "M") or the mapping matrix $\mathbf{G}$ (mat = "G").

References

Examples

# Cross-temporal framework (Z = X + Y, annual-quarterly)
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
Mct <- ctprojmat(agg_mat = A, agg_order = 4, comb = "ols")
Gct <- ctprojmat(agg_mat = A, agg_order = 4, comb = "ols", mat = "G")

# Cross-temporal framework (Z = X + Y, annual-quarterly)
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
Mct <- ctprojmat(agg_mat = A, agg_order = 4, comb = "ols")
Gct <- ctprojmat(agg_mat = A, agg_order = 4, comb = "ols", mat = "G")

Optimal combination cross-temporal reconciliation

Description

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, ...)
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: each row represents a constraint equation, and each column represents a variable. The matrix can be of full rank, meaning the rows are linearly independent, but this is not a strict requirement, as the function allows for redundancy in the constraints.
`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: "`proj`" (default): Projection approach according to Byron (1978, 1979). "`strc`": Structural approach as proposed by Hyndman et al. (2011). "`proj_osqp`": Numerical solution using osqp for projection approach. "`strc_osqp`": Numerical solution using osqp for structural approach.
`nn`	A string specifying the algorithm to compute non-negative forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`bpv`": block principal pivoting algorithm. "`sntz`": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
`settings`	A list of control parameters. `nn = "osqp"` An object of class `osqpSettings` specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020). `nn = "bpv"` It includes: `ptype` for permutation method ("`random`" or "`fixed`", default), `par` for the number of full exchange rules that may be attempted (`10`, default), `tol` for the tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `gtol` for the gradient tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `itmax` for the maximum number of algorithm iterations (`100`, default)
`bounds`	A matrix (see set_bounds) with 5 columns ( $i,k,j,lower,upper$ ), such that Column 1 represents the cross-sectional series ( $i = 1, \dots, n$ ). Column 2 represents the temporal aggregation order ( $k = m,\dots,1$ ). Column 3 represents the temporal forecast horizon ( $j = 1,\dots,m/k$ ). Columns 4 and 5 indicates the lower and lower bounds, respectively.
`immutable`	A matrix with three columns ( $i,k,j$ ), such that Column 1 represents the cross-sectional series ( $i = 1, \dots, n$ ). Column 2 represents the temporal aggregation order ( $k = m,\dots,1$ ). Column 3 represents 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).

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

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

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

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

Cross-temporal top-down reconciliation

Description

Top-down forecast reconciliation for cross-temporal hierarchical/grouped time series, where the forecast of a ‘Total’ (top-level series, expected to be positive) is disaggregated according to a proportional scheme (weights). Besides fulfilling any aggregation constraint, the top-down reconciled forecasts should respect two main properties:

the top-level value remains unchanged;
all the bottom time series reconciled forecasts are non-negative.

Usage

cttd(base, agg_mat, agg_order, weights, tew = "sum", normalize = TRUE)
cttd(base, agg_mat, agg_order, weights, tew = "sum", normalize = TRUE)

Arguments

`base`	A ( $hm \times 1$ ) numeric vector containing top- and $m$ temporal aggregated level base forecasts; $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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.
`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$ .
`weights`	A ( $n_b \times hm$ ) numeric matrix containing the proportions for each high-frequency bottom time series; $n_b$ is the total number of high-frequency bottom variables, $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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).
`normalize`	If `TRUE` (default), the `weights` will sum to 1.

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

Examples

set.seed(123)
# (3 x 1) top base forecasts vector (simulated), forecast horizon = 3
topf <- rnorm(3, 10)
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y

# Same weights for different forecast horizons, agg_order = 4
fix_weights <- matrix(runif(4*2), 2, 4)
reco <- cttd(base = topf, agg_mat = A, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*2*3), 2, 3*4)
recoh <- cttd(base = topf, agg_mat = A, agg_order = 4, weights = h_weights)

set.seed(123)
# (3 x 1) top base forecasts vector (simulated), forecast horizon = 3
topf <- rnorm(3, 10)
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y

# Same weights for different forecast horizons, agg_order = 4
fix_weights <- matrix(runif(4*2), 2, 4)
reco <- cttd(base = topf, agg_mat = A, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- matrix(runif(4*2*3), 2, 3*4)
recoh <- cttd(base = topf, agg_mat = A, agg_order = 4, weights = h_weights)

Cross-temporal reconciliation tools

Description

Some useful tools for the cross-temporal forecast reconciliation of a linearly constrained (e.g., hierarchical/grouped) multiple time series.

Usage

cttools(agg_mat, cons_mat, agg_order, tew = "sum", fh = 1, sparse = TRUE)
cttools(agg_mat, cons_mat, agg_order, tew = "sum", fh = 1, sparse = TRUE)

Arguments

`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: each row represents a constraint equation, and each column represents a variable. The matrix can be of full rank, meaning the rows are linearly independent, but this is not a strict requirement, as the function allows for redundancy in the constraints.
`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$ .
`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).
`fh`	Forecast horizon for the lowest frequency (most temporally aggregated) time series (default is `1`).
`sparse`	Option to return sparse matrices (default is `TRUE`).

Value

A list with four elements:

`dim`	A vector containing information about the number of series for the complete system (`n`), for upper levels (`na`) and bottom level (`nb`), the maximum aggregation order (`m`), the number of factor (`p`), the partial (`ks`) and total sum (`kt`) of factors.
`set`	The vector of the temporal aggregation orders (in decreasing order).
`agg_mat`	The cross-temporal aggregation matrix.
`strc_mat`	The cross-temporal structural matrix.
`cons_mat`	The cross-temporal zero constraints matrix.

Examples

# Cross-temporal framework
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
m <- 4 # from quarterly to annual temporal aggregation
cttools(agg_mat = A, agg_order = m)

# Cross-temporal framework
A <- t(c(1,1)) # Aggregation matrix for Z = X + Y
m <- 4 # from quarterly to annual temporal aggregation
cttools(agg_mat = A, agg_order = m)

Cross-sectional aggregation matrix of a dataframe

Description

This function allows the user to easily build the ( $n_a \times n_b$ ) cross-sectional aggregation matrix starting from a data frame.

Usage

df2aggmat(formula, data, sep = "_", sparse = TRUE, top_label = "Total",
          verbose = TRUE)
df2aggmat(formula, data, sep = "_", sparse = TRUE, top_label = "Total",
          verbose = TRUE)

Arguments

`formula`	Specification of the hierarchical structure: grouped hierarchies are specified using `~ g1 * g2` and nested hierarchies are specified using `~ parent / child`. Mixtures of the two formulations are also possible, like `~ g1 * (grandparent / parent / child)`.
`data`	A dataset in which each column contains the values of the variables in the formula and each row identifies a bottom level time series.
`sep`	Character to separate the names of the aggregated series, (default is "`_`").
`sparse`	Option to return sparse matrices (default is `TRUE`).
`top_label`	Label of the top level variable (default is "`Total`").
`verbose`	If `TRUE` (default), hierarchy informations are printed.

Value

A (na x nb) matrix.

Examples

## Balanced hierarchy
#         T
#    |--------|
#    A        B
#  |---|   |--|--|
# AA   AB  BA BB BC
# Names of the bottom level variables
data_bts <- data.frame(X1 = c("A", "A", "B", "B", "B"),
                       X2 = c("A", "B", "A", "B", "C"),
                       stringsAsFactors = FALSE)
# Cross-sectional aggregation matrix
agg_mat <- df2aggmat(~ X1 / X2, data_bts, sep = "", verbose = FALSE)

## Unbalanced hierarchy
#                 T
#       |---------|---------|
#       A         B         C
#     |---|     |---|     |---|
#    AA   AB   BA   BB   CA   CB
#  |----|         |----|
# AAA  AAB       BBA  BBB
# Names of the bottom level variables
data_bts <- data.frame(X1 = c("A", "A", "A", "B", "B", "B", "C", "C"),
                       X2 = c("A", "A", "B", "A", "B", "B", "A", "B"),
                       X3 = c("A", "B", NA, NA, "A", "B", NA, NA),
                       stringsAsFactors = FALSE)
# Cross-sectional aggregation matrix
agg_mat <- df2aggmat(~ X1 / X2 / X3, data_bts, sep = "", verbose = FALSE)

## Group of two hierarchies
#     T          T         T | A  | B  | C
#  |--|--|  X  |---|  ->  ---+----+----+----
#  A  B  C     M   F       M | AM | BM | CM
#                          F | AF | BF | CF
# Names of the bottom level variables
data_bts <- data.frame(X1 = c("A", "A", "B", "B", "C", "C"),
                       Y1 = c("M", "F", "M", "F", "M", "F"),
                       stringsAsFactors = FALSE)
# Cross-sectional aggregation matrix
agg_mat <- df2aggmat(~ Y1 * X1, data_bts, sep = "", verbose = FALSE)

## Balanced hierarchy
#         T
#    |--------|
#    A        B
#  |---|   |--|--|
# AA   AB  BA BB BC
# Names of the bottom level variables
data_bts <- data.frame(X1 = c("A", "A", "B", "B", "B"),
                       X2 = c("A", "B", "A", "B", "C"),
                       stringsAsFactors = FALSE)
# Cross-sectional aggregation matrix
agg_mat <- df2aggmat(~ X1 / X2, data_bts, sep = "", verbose = FALSE)

## Unbalanced hierarchy
#                 T
#       |---------|---------|
#       A         B         C
#     |---|     |---|     |---|
#    AA   AB   BA   BB   CA   CB
#  |----|         |----|
# AAA  AAB       BBA  BBB
# Names of the bottom level variables
data_bts <- data.frame(X1 = c("A", "A", "A", "B", "B", "B", "C", "C"),
                       X2 = c("A", "A", "B", "A", "B", "B", "A", "B"),
                       X3 = c("A", "B", NA, NA, "A", "B", NA, NA),
                       stringsAsFactors = FALSE)
# Cross-sectional aggregation matrix
agg_mat <- df2aggmat(~ X1 / X2 / X3, data_bts, sep = "", verbose = FALSE)

## Group of two hierarchies
#     T          T         T | A  | B  | C
#  |--|--|  X  |---|  ->  ---+----+----+----
#  A  B  C     M   F       M | AM | BM | CM
#                          F | AF | BF | CF
# Names of the bottom level variables
data_bts <- data.frame(X1 = c("A", "A", "B", "B", "C", "C"),
                       Y1 = c("M", "F", "M", "F", "M", "F"),
                       stringsAsFactors = FALSE)
# Cross-sectional aggregation matrix
agg_mat <- df2aggmat(~ Y1 * X1, data_bts, sep = "", verbose = FALSE)

Reconciled forecasts to matrix/vector

Description

This function splits the temporal vectors and the cross-temporal matrices in a list according to the temporal aggregation order

Usage

FoReco2matrix(x, agg_order, keep_names = FALSE)
FoReco2matrix(x, agg_order, keep_names = FALSE)

Arguments

`x`	An output from any reconciliation function implemented by FoReco.
`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$ .
`keep_names`	If `FALSE` (default), the rownames names of the output matrices are removed.

Value

A list of matrices or vectors distinct by temporal aggregation order.

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))))

reco <- ctrec(base = base, agg_mat = t(c(1,1)), agg_order = 4, comb = "ols")
matrix_list <- FoReco2matrix(reco)
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))))

reco <- ctrec(base = base, agg_mat = t(c(1,1)), agg_order = 4, comb = "ols")
matrix_list <- FoReco2matrix(reco)

Italian Quarterly National Accounts

Description

A subset of the data used by Girolimetto et al. (2023) from the Italian Quarterly National Accounts (output, income and expenditure sides) spanning the period 2000:Q1-2019:Q4.

Usage

# 21 time series of the Italian Quarterly National Accounts
itagdp

# 'agg_mat' and 'cons_mat' for the output side
outside

# 'agg_mat' and 'cons_mat' for the expenditure side
expside

# 'agg_mat' and 'cons_mat' for the income side
incside

# zero constraints matrix encompassing output, expenditure and income sides
gdpconsmat
# 21 time series of the Italian Quarterly National Accounts
itagdp

# 'agg_mat' and 'cons_mat' for the output side
outside

# 'agg_mat' and 'cons_mat' for the expenditure side
expside

# 'agg_mat' and 'cons_mat' for the income side
incside

# zero constraints matrix encompassing output, expenditure and income sides
gdpconsmat

Format

itagdp is a $(80 \times 21)$ ts object, corresponding to 21 time series of the Italian Quarterly National Accounts (2000:Q1-2019:Q4).

outside, income and expenditure are lists with two elements:

agg_mat contains the $(1 \times 2)$ , $(2 \times 4)$ , or $(6 \times 8)$ aggregation matrix according to output, income or expenditure side, respectively.
cons_mat contains the $(1 \times 3)$ , $(2 \times 6)$ , or $(6 \times 14)$ zero constraints matrix according to output, income or expenditure side, respectively.

gdpconsmat is the complete $(9 \times 21)$ zero constraints matrix encompassing output, expenditure and income sides.

Source

https://ec.europa.eu/eurostat/web/national-accounts/

References

Iterative cross-temporal reconciliation

Description

This function performs the iterative procedure described in Di Fonzo and Girolimetto (2023), which produces cross-temporally reconciled forecasts by alternating forecast reconciliation along one single dimension (either cross-sectional or temporal) at each iteration step.

Usage

iterec(base, cslist, telist, res = NULL, itmax = 100, tol = 1e-5,
       type = "tcs", norm = "inf", verbose = TRUE)
iterec(base, cslist, telist, res = NULL, itmax = 100, tol = 1e-5,
       type = "tcs", norm = "inf", verbose = TRUE)

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.
`cslist`	A list of elements for the cross-sectional reconciliation. See csrec for a complete list (excluded `base` and `res`).
`telist`	A list of elements for the temporal reconciliation. See terec for a complete list (excluded `base` and `res`).
`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.
`itmax`	Max number of iteration (`100`, default).
`tol`	Convergence tolerance (`1e-5`, default).
`type`	A string specifying the uni-dimensional reconciliation order: temporal and then cross-sectional ("`tcs`") or cross-sectional and then temporal ("`cst`").
`norm`	Norm used to calculate the temporal and the cross-sectional incoherence: infinity norm ("`inf`", default), one norm ("`one`"), and 2-norm ("`two`").
`verbose`	If `TRUE`, reconciliation information are printed.

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

References

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

rite <- iterec(base = base,
               cslist = list(agg_mat = A, comb = "shr"),
               telist = list(agg_order = m, comb = "wlsv"),
               res = res)

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

rite <- iterec(base = base,
               cslist = list(agg_mat = A, comb = "shr"),
               telist = list(agg_order = m, comb = "wlsv"),
               res = res)

Linear combination (aggregation) matrix for a general linearly constrained multiple time series

Description

This function transforms a general (possibly redundant) zero constraints matrix into a linear combination (aggregation) matrix $\mathbf{A}_{cs}$ . When working with a general linearly constrained multiple ( $n$ -variate) time series, getting a linear combination matrix $\mathbf{A}_{cs}$ is a critical step to obtain a structural-like representation such that

$\mathbf{C}_{cs} = [\mathbf{I} \quad -\mathbf{A}],$

where $\mathbf{C}_{cs}$ is the full rank zero constraints matrix (Girolimetto and Di Fonzo, 2023).

Usage

lcmat(cons_mat, method = "rref", tol = sqrt(.Machine$double.eps),
       verbose = FALSE, sparse = TRUE)
lcmat(cons_mat, method = "rref", tol = sqrt(.Machine$double.eps),
       verbose = FALSE, sparse = TRUE)

Arguments

`cons_mat`	A ( $r \times n$ ) numeric matrix representing the cross-sectional zero constraints.
`method`	Method to use: "`rref`" for the Reduced Row Echelon Form through Gauss-Jordan elimination (default), or "`qr`" for the (pivoting) QR decomposition (Strang, 2019).
`tol`	Tolerance for the "`rref`" or "`qr`" method.
`verbose`	If `TRUE`, intermediate steps are printed (default is `FALSE`).
`sparse`	Option to return a sparse matrix (default is `TRUE`).

Value

A list with two elements: (i) the linear combination (aggregation) matrix (agg_mat) and (ii) the vector of the column permutations (pivot).

References

Strang, G. (2019), Linear algebra and learning from data, Wellesley, Cambridge Press.

Examples

## Two hierarchy sharing the same top-level variable, but not sharing the bottom variables
#        X            X
#    |-------|    |-------|
#    A       B    C       D
#  |---|
# A1   A2
# 1) X = C + D,
# 2) X = A + B,
# 3) A = A1 + A2.
cons_mat <- matrix(c(1,-1,-1,0,0,0,0,
               1,0,0,-1,-1,0,0,
               0,0,0,1,0,-1,-1), nrow = 3, byrow = TRUE)
obj <- lcmat(cons_mat = cons_mat, verbose = TRUE)
agg_mat <- obj$agg_mat # linear combination matrix
pivot <- obj$pivot # Pivot vector

## Two hierarchy sharing the same top-level variable, but not sharing the bottom variables
#        X            X
#    |-------|    |-------|
#    A       B    C       D
#  |---|
# A1   A2
# 1) X = C + D,
# 2) X = A + B,
# 3) A = A1 + A2.
cons_mat <- matrix(c(1,-1,-1,0,0,0,0,
               1,0,0,-1,-1,0,0,
               0,0,0,1,0,-1,-1), nrow = 3, byrow = TRUE)
obj <- lcmat(cons_mat = cons_mat, verbose = TRUE)
agg_mat <- obj$agg_mat # linear combination matrix
pivot <- obj$pivot # Pivot vector

Informations on the reconciliation process

Description

This function extracts reconciliation information from the output of any reconciled function implemented by FoReco.

Usage

recoinfo(x, verbose = TRUE)
recoinfo(x, verbose = TRUE)

Arguments

`x`	An output from any reconciliation function implemented by FoReco.
`verbose`	If `TRUE` (defaults), reconciliation information are printed.

Value

A list containing the following reconciliation process informations:

`rfun`	the reconciliation function.
`cs_n`	the cross-sectional number of variables.
`te_set`	the set of temporal aggregation orders.
`forecast_horizon`	the forecast horizon (in temporal and cross-temporal frameworks, for the most temporally aggregated series).
`framework`	the reconciliation framework (cross-sectional, temporal or cross-temporal).
`info`	non-negative reconciled forecast convergence information.
`lcc`	list of level conditional reconciled forecasts (+ BU) for cslcc, telcc and ctlcc.
`nn`	if `TRUE`, all the forecasts are not negative.
`comb`	the covariance approximation.

One-step and multi-step residuals

Description

These functions can be used to arrange residuals to reconcile temporal or cross-temporal forecasts.

res2matrix takes as input a set of temporal and cross-temporal residuals and re-organizes them into a matrix where the rows correspond to different forecast horizons, capturing the temporal dimension. Meanwhile, the columns are ordered based on the specific arrangement as described in Di Fonzo and Girolimetto (2023).

arrange_hres takes as input a list of multi-step residuals and is designed to organize them in accordance with their time order (Girolimetto et al. 2023). When applied, this function ensures that the sequence of multi-step residuals aligns with the chronological order in which they occurred.

Usage

res2matrix(res, agg_order)

arrange_hres(list_res)
res2matrix(res, agg_order)

arrange_hres(list_res)

Arguments

`res`	A ( $n \times N(k^\ast+m)$ ) numeric matrix (cross-temporal framework) or an ( $N(k^\ast+m) \times 1$ ) numeric vector (temporal framework) representing the in-sample residuals at all the temporal frequencies ordered from the lowest frequency to the highest frequency (columns) for each variable (rows).
`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$ .
`list_res`	A list of $H$ multi-step residuals. Each element in the list can be either a ( $T \times 1$ ) vector (temporal framework) or a ( $T \times n$ ) matrix (cross-temporal framework).

Details

Let $Z_t$ , $t=1,\dots,T$ , be a univariate time series. We can define the multi-step residuals such us

$\widehat{\varepsilon}_{h,t} = Z_{t+h} - \widehat{Z}_{t+h|t} \qquad h \le t \le T-h$

where $\widehat{Z}_{t+h|t}$ is the $h$ -step fitted value, calculated as the $h$ -step ahead forecast condition to the information up to time $t$ . Given the list of errors at different steps

$\left([\widehat{\varepsilon}_{1,1}, \; \dots, \; \widehat{\varepsilon}_{1,T}], \dots, [\widehat{\varepsilon}_{H,1}, \; \dots, \; \widehat{\varepsilon}_{H,T}]\right),$

arrange_hres returns a $T$ -vector with the residuals, organized in the following way:

$[\varepsilon_{1,1} \; \varepsilon_{2,2} \; \dots \; \varepsilon_{H,H} \; \varepsilon_{1,H+1} \; \dots \; \varepsilon_{H,T-H}]'$

A similar organisation can be apply to a multivariate time series.

Value

res2matrix returns a ( $N \times n(k^\ast + m)$ ) matrix, where $n = 1$ for the temporal framework.

arrange_hres returns a ( $N(k^\ast+m) \times 1$ ) vector (temporal framework) or a ( $n \times N(k^\ast+m)$ ) matrix (cross-temporal framework) of multi-step residuals.

References

Set bounds for bounded forecast reconciliation

Description

This function defines the bounds matrix considering cross-sectional, temporal, or cross-temporal frameworks. The output matrix can be used as input for the bounds parameter in functions such as csrec, terec, or ctrec, to perform bounded reconciliations.

Usage

set_bounds(n, k, h, lb = -Inf, ub = Inf, approach = "osqp", bounds = NULL)
set_bounds(n, k, h, lb = -Inf, ub = Inf, approach = "osqp", bounds = NULL)

Arguments

`n`	A ( $b \times 1$ ) vector representing the $i$ th cross-sectional series ( $i = 1, \dots, n$ ), where $b$ is the number of bounds to be set.
`k`	A ( $b \times 1$ ) vector specifying the temporal aggregation orders ( $k = m, \dots, 1$ ).
`h`	A ( $b \times 1$ ) vector representing the forecast horizons ( $j = 1, \dots, m/k$ ).
`lb`, `ub`	A ( $b \times 1$ ) vector of lower and upper bounds.
`approach`	A string specifying the algorithm to compute bounded reconciled forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`sftb`": heuristic "set-forecasts-to-bounds", which adjusts the reconciled forecasts to be within specified bounds without further optimization.
`bounds`	A matrix of previous bounds to be added. If not specified, new bounds will be computed.

Value

A numeric matrix representing the computed bounds, which can be:

Cross-sectional ( $b \times 3$ ) matrix for cross-sectional reconciliation (csrec).
Temporal ( $b \times 4$ ) matrix for temporal reconciliation (terec).
Cross-temporal ( $b \times 5$ ) matrix for cross-temporal reconciliation (ctrec).

Examples

# Example 1
# Two cross-sectional series (i = 2,3),
# with each series required to be between 0 and 1.
n <- c(2, 3)
lb <- c(0, 0)
ub <- c(1,1)
bounds_mat <- set_bounds(n = c(2, 3),
                         lb = rep(0, 2), # or lb = 0
                         ub = rep(1, 2)) # or ub = 1

# Example 2
# All the monthly values are between 0 and 1.
bounds_mat <- set_bounds(k = rep(1, 12),  # or k = 1
                         h = 1:12,
                         lb = rep(0, 12), # or lb = 0
                         ub = rep(1, 12)) # or ub = 1

# Example 3
# For two cross-sectional series (i = 2,3),
# all the monthly values are between 0 and 1.
bounds_mat <- set_bounds(n = rep(c(2, 3), each = 12),
                         k = 1,
                         h = rep(1:12, 2),
                         lb = 0, # or lb = 0
                         ub = 1) # or ub = 1

# Example 1
# Two cross-sectional series (i = 2,3),
# with each series required to be between 0 and 1.
n <- c(2, 3)
lb <- c(0, 0)
ub <- c(1,1)
bounds_mat <- set_bounds(n = c(2, 3),
                         lb = rep(0, 2), # or lb = 0
                         ub = rep(1, 2)) # or ub = 1

# Example 2
# All the monthly values are between 0 and 1.
bounds_mat <- set_bounds(k = rep(1, 12),  # or k = 1
                         h = 1:12,
                         lb = rep(0, 12), # or lb = 0
                         ub = rep(1, 12)) # or ub = 1

# Example 3
# For two cross-sectional series (i = 2,3),
# all the monthly values are between 0 and 1.
bounds_mat <- set_bounds(n = rep(c(2, 3), each = 12),
                         k = 1,
                         h = rep(1:12, 2),
                         lb = 0, # or lb = 0
                         ub = 1) # or ub = 1

Shrinkage of the covariance matrix

Description

Shrinkage of the covariance matrix according to Schäfer and Strimmer (2005).

Usage

shrink_estim(x, mse = TRUE)
shrink_estim(x, mse = TRUE)

Arguments

`x`	A numeric matrix containing the in-sample residuals.
`mse`	If `TRUE` (default), the residuals used to compute the covariance matrix are not mean-corrected.

Value

A shrunk covariance matrix.

References

Schäfer, J.L. and 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

Shrinkage of the covariance matrix using the Oracle approximation

Description

Shrinkage of the covariance matrix according to the Oracle Approximating Shrinkage (OAS) of Chen et al. (2009) and Ando and Xiao (2023).

Usage

shrink_oasd(x, mse = TRUE)
shrink_oasd(x, mse = TRUE)

Arguments

`x`	A numeric matrix containing the in-sample residuals.
`mse`	If `TRUE` (default), the residuals used to compute the covariance matrix are not mean-corrected.

Value

A shrunk covariance matrix.

References

Ando, S., and Xiao, M. (2023), High-dimensional covariance matrix estimation: shrinkage toward a diagonal target. IMF Working Papers, 2023(257), A001. doi:10.5089/9798400260780.001.A001

Chen, Y., Wiesel, A., and Hero, A. O. (2009), Shrinkage estimation of high dimensional covariance matrices, 2009 IEEE international conference on acoustics, speech and signal processing, 2937–2940. IEEE.

Heuristic cross-temporal reconciliation

Description

tcsrec replicates the procedure by Kourentzes and Athanasopoulos (2019): (i) for each time series the forecasts at any temporal aggregation order are reconciled using temporal hierarchies; (ii) time-by-time cross-sectional reconciliation is performed; and (iii) the projection matrices obtained at step (ii) are then averaged and used to cross-sectionally reconcile the forecasts obtained at step (i). In cstrec, the order of application of the two reconciliation steps (temporal first, then cross-sectional), is inverted compared to tcsrec (Di Fonzo and Girolimetto, 2023).

Usage

# First-temporal-then-cross-sectional forecast reconciliation
tcsrec(base, cslist, telist, res = NULL, avg = "KA")

# First-cross-sectional-then-temporal forecast reconciliation
cstrec(base, cslist, telist, res = NULL)
# First-temporal-then-cross-sectional forecast reconciliation
tcsrec(base, cslist, telist, res = NULL, avg = "KA")

# First-cross-sectional-then-temporal forecast reconciliation
cstrec(base, cslist, telist, res = 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.
`cslist`	A list of elements for the cross-sectional reconciliation. See csrec for a complete list (excluded `base` and `res`).
`telist`	A list of elements for the temporal reconciliation. See terec for a complete list (excluded `base` and `res`).
`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.
`avg`	If `avg = "KA"` (default), the final projection matrix $\mathbf{M}$ is the one proposed by Kourentzes and Athanasopoulos (2019), otherwise it is calculated as simple average of all the involved projection matrices at step 2 of the procedure (see Di Fonzo and Girolimetto, 2023).

Value

A ( $n \times h(k^\ast+m)$ ) numeric matrix of cross-temporal reconciled forecasts.

Warning

The two-step heuristic reconciliation allows considering non negativity constraints only in the first step. This means that non-negativity is not guaranteed in the final reconciled values.

References

Kourentzes, N. and Athanasopoulos, G. (2019), Cross-temporal coherent forecasts for Australian tourism, Annals of Tourism Research, 75, 393-409. doi:10.1016/j.annals.2019.02.001

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

rtcs <- tcsrec(base = base,
               cslist = list(agg_mat = A, comb = "shr"),
               telist = list(agg_order = m, comb = "wlsv"),
               res = res)

rcst <- tcsrec(base = base,
               cslist = list(agg_mat = A, comb = "shr"),
               telist = list(agg_order = m, comb = "wlsv"),
               res = res)

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

rtcs <- tcsrec(base = base,
               cslist = list(agg_mat = A, comb = "shr"),
               telist = list(agg_order = m, comb = "wlsv"),
               res = res)

rcst <- tcsrec(base = base,
               cslist = list(agg_mat = A, comb = "shr"),
               telist = list(agg_order = m, comb = "wlsv"),
               res = res)

Temporal joint block bootstrap

Description

Joint block bootstrap for generating probabilistic base forecasts that take into account the correlation between different temporal aggregation orders (Girolimetto et al. 2023).

Usage

teboot(model_list, boot_size, agg_order, block_size = 1, seed = NULL)
teboot(model_list, boot_size, agg_order, block_size = 1, seed = NULL)

Arguments

`model_list`	A list of all the $(k^\ast+m)$ base forecasts models ordered from the lowest frequency (most temporally aggregated) to the highest frequency. A `simulate()` function for each model has to be available and implemented according to the package forecast, with the following mandatory parameters: object, innov, future, and nsim.
`boot_size`	The number of bootstrap replicates.
`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$ .
`block_size`	Block size of the bootstrap, which is typically equivalent to the forecast horizon for the most temporally aggregated series.
`seed`	An integer seed.

Value

A list with two elements: the seed used to sample the errors and a ( $\text{boot\_size}\times (k^\ast+m)\text{block\_size}$ ) matrix.

References

Temporal bottom-up reconciliation

Description

Temporal bottom-up reconciled forecasts at any temporal aggregation level are computed by appropriate aggregation of the high-frequency base forecasts, $\widehat{\mathbf{x}}^{[1]}$ :

$\widetilde{\mathbf{x}} = \mathbf{S}_{te}\widehat{\mathbf{x}}^{[1]},$

where $\mathbf{S}_{te}$ is the temporal structural matrix.

Usage

tebu(base, agg_order, tew = "sum", sntz = FALSE)
tebu(base, agg_order, tew = "sum", sntz = FALSE)

Arguments

`base`	A ( $hm \times 1$ ) numeric vector containing the high-frequency base forecasts; $m$ is the max. temporal aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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$ .
`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).
`sntz`	If `TRUE`, the negative base forecasts are set to zero before applying bottom-up.

Value

A ( $h(k^\ast+m) \times 1$ ) numeric vector of temporal reconciled forecasts.

Examples

set.seed(123)
# (4 x 1) high frequency base forecasts vector (simulated),
# agg_order = 4 (annual-quarterly)
hfts <- rnorm(4, 5)

reco <- tebu(base = hfts, agg_order = 4)

# Non negative reconciliation
hfts[4] <- -hfts[4] # Making negative one of the quarterly base forecasts
nnreco <- tebu(base = hfts, agg_order = 4, sntz = TRUE)

set.seed(123)
# (4 x 1) high frequency base forecasts vector (simulated),
# agg_order = 4 (annual-quarterly)
hfts <- rnorm(4, 5)

reco <- tebu(base = hfts, agg_order = 4)

# Non negative reconciliation
hfts[4] <- -hfts[4] # Making negative one of the quarterly base forecasts
nnreco <- tebu(base = hfts, agg_order = 4, sntz = TRUE)

Temporal covariance matrix approximation

Description

This function provides an approximation of the temporal base forecasts errors covariance matrix using different reconciliation methods (see Di Fonzo and Girolimetto, 2023).

Usage

tecov(comb, agg_order = NULL, res = NULL, tew = "sum",
      mse = TRUE, shrink_fun = shrink_estim, ...)
tecov(comb, agg_order = NULL, res = NULL, tew = "sum",
      mse = TRUE, shrink_fun = shrink_estim, ...)

Arguments

`comb`	A string specifying the reconciliation method. Ordinary least squares: "`ols`" (default) - identity error covariance. Weighted least squares: "`str`" - structural variances. "`wlsh`" - hierarchy variances (uses `res`). "`wlsv`" - series variances (uses `res`). Generalized least squares (uses `res`): "`acov`" - series auto-covariance. "`strar1`" - structural Markov covariance. "`sar1`" - series Markov covariance. "`har1`" - hierarchy Markov covariance. "`shr`"/"`sam`" - shrunk/sample covariance.
`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$ .
`res`	A ( $N(k^\ast+m) \times 1$ ) optional numeric vector containing the in-sample residuals at all the temporal frequencies ordered from the lowest frequency to the highest frequency. This vector is used to compute come 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).
`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)
`...`	Not used.

Value

A ( $(k^\ast+m) \times (k^\ast+m)$ ) symmetric matrix.

References

Examples

# (7 x 70) in-sample residuals matrix (simulated), agg_order = 4
res <- rnorm(70)

cov1 <- tecov("ols", agg_order = 4)                 # OLS methods
cov2 <- tecov("str", agg_order = 4)                 # STRC methods
cov3 <- tecov("wlsv", agg_order = 4, res = res)     # WLSv methods
cov4 <- tecov("wlsh", agg_order = 4, res = res)     # WLSh methods
cov5 <- tecov("acov", agg_order = 4, res = res)     # ACOV methods
cov6 <- tecov("strar1", agg_order = 4, res = res)   # STRAR1 methods
cov7 <- tecov("har1", agg_order = 4, res = res)     # HAR1 methods
cov8 <- tecov("sar1", agg_order = 4, res = res)     # SAR1 methods
cov9 <- tecov("shr", agg_order = 4, res = res)      # SHR methods
cov10 <- tecov("sam", agg_order = 4, res = res)     # SAM methods

# Custom covariance matrix
tecov.ols2 <- function(comb, x) diag(x)
tecov(comb = "ols2", x = 7) # == tecov("ols", agg_order = 4)

# (7 x 70) in-sample residuals matrix (simulated), agg_order = 4
res <- rnorm(70)

cov1 <- tecov("ols", agg_order = 4)                 # OLS methods
cov2 <- tecov("str", agg_order = 4)                 # STRC methods
cov3 <- tecov("wlsv", agg_order = 4, res = res)     # WLSv methods
cov4 <- tecov("wlsh", agg_order = 4, res = res)     # WLSh methods
cov5 <- tecov("acov", agg_order = 4, res = res)     # ACOV methods
cov6 <- tecov("strar1", agg_order = 4, res = res)   # STRAR1 methods
cov7 <- tecov("har1", agg_order = 4, res = res)     # HAR1 methods
cov8 <- tecov("sar1", agg_order = 4, res = res)     # SAR1 methods
cov9 <- tecov("shr", agg_order = 4, res = res)      # SHR methods
cov10 <- tecov("sam", agg_order = 4, res = res)     # SAM methods

# Custom covariance matrix
tecov.ols2 <- function(comb, x) diag(x)
tecov(comb = "ols2", x = 7) # == tecov("ols", agg_order = 4)

Level conditional coherent reconciliation for temporal hierarchies

Description

This function implements a forecast reconciliation procedure inspired by the original proposal by Hollyman et al. (2021) for temporal hierarchies. Level conditional coherent reconciled forecasts are conditional on (i.e., constrained by) the base forecasts of a specific upper level in the hierarchy (exogenous constraints). It also allows handling the linear constraints linking the variables endogenously (Di Fonzo and Girolimetto, 2022). The function can calculate Combined Conditional Coherent (CCC) forecasts as simple averages of Level-Conditional Coherent (LCC) and bottom-up reconciled forecasts, with either endogenous or exogenous constraints.

Usage

telcc(base, agg_order, comb = "ols", res = NULL, CCC = TRUE,
      const = "exogenous", hfts = NULL, tew = "sum",
      approach = "proj", nn = NULL, settings = NULL, ...)
telcc(base, agg_order, comb = "ols", res = NULL, CCC = TRUE,
      const = "exogenous", hfts = NULL, tew = "sum",
      approach = "proj", nn = NULL, settings = NULL, ...)

Arguments

`base`	A ( $h(k^\ast + m) \times 1$ ) numeric vector containing base forecasts to be reconciled ordered from the lowest frequency to the highest frequency; $m$ is the max aggregation order, $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.
`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 tecov.
`res`	A ( $N(k^\ast+m) \times 1$ ) optional numeric vector containing the in-sample residuals at all the temporal frequencies ordered from the lowest frequency to the highest frequency. This vector is used to compute come covariance matrices.
`CCC`	A logical value indicating whether the Combined Conditional Coherent reconciled forecasts reconciliation should include bottom-up forecasts (`TRUE`, default), or not.
`const`	A string specifying the reconciliation constraints: "`exogenous`" (default): Fixes the top level of each sub-hierarchy. "`endogenous`": Coherently revises both the top and bottom levels.
`hfts`	A ( $mh \times 1$ ) numeric vector containing high frequency base forecasts defined by the user. This parameter can be omitted if only base forecasts in `base` are used (see Di Fonzo and Girolimetto, 2024).
`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: "`proj`" (default): Projection approach according to Byron (1978, 1979). "`strc`": Structural approach as proposed by Hyndman et al. (2011). "`proj_osqp`": Numerical solution using osqp for projection approach. "`strc_osqp`": Numerical solution using osqp for structural approach.
`nn`	A string specifying the algorithm to compute non-negative forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`bpv`": block principal pivoting algorithm. "`sntz`": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
`settings`	A list of control parameters. `nn = "osqp"` An object of class `osqpSettings` specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020). `nn = "bpv"` It includes: `ptype` for permutation method ("`random`" or "`fixed`", default), `par` for the number of full exchange rules that may be attempted (`10`, default), `tol` for the tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `gtol` for the gradient tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `itmax` for the maximum number of algorithm iterations (`100`, default)
`...`	Arguments passed on to `tecov` `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)

Value

A ( $h(k^\ast+m) \times 1$ ) numeric vector of temporal reconciled forecasts.

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. (2023b) Spatio-temporal reconciliation of solar forecasts. Solar Energy 251, 13–29. doi:10.1016/j.solener.2023.01.003

Hollyman, R., Petropoulos, F. and Tipping, M.E. (2021), Understanding forecast reconciliation. European Journal of Operational Research, 294, 149–160. doi:10.1016/j.ejor.2021.01.017

Examples

set.seed(123)
# (7 x 1) base forecasts vector (simulated), agg_order = 4
base <- rnorm(7, rep(c(20, 10, 5), c(1, 2, 4)))
# (70 x 1) in-sample residuals vector (simulated)
res <- rnorm(70)
# (4 x 1) Naive high frequency base forecasts vector: all forecasts are set equal to 2.5
naive <- rep(2.5, 4)

## EXOGENOUS CONSTRAINTS
# Level Conditional Coherent (LCC) reconciled forecasts
exo_LC <- telcc(base = base, agg_order = 4, comb = "wlsh", hfts = naive,
                res = res, nodes = "auto", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts
exo_CCC <- telcc(base = base, agg_order = 4, comb = "wlsh", hfts = naive,
                 res = res, nodes = "auto", CCC = TRUE)

# Results detailed by level:
info_exo <- recoinfo(exo_CCC, verbose = FALSE)
# info_exo$lcc

## ENDOGENOUS CONSTRAINTS
# Level Conditional Coherent (LCC) reconciled forecasts
endo_LC <- telcc(base = base, agg_order = 4, comb = "wlsh", res = res,
                 nodes = "auto", CCC = FALSE, const = "endogenous")

# Combined Conditional Coherent (CCC) reconciled forecasts
endo_CCC <- telcc(base = base, agg_order = 4, comb = "wlsh", res = res,
                  nodes = "auto", CCC = TRUE, const = "endogenous")

# Results detailed by level:
info_endo <- recoinfo(endo_CCC, verbose = FALSE)
# info_endo$lcc

set.seed(123)
# (7 x 1) base forecasts vector (simulated), agg_order = 4
base <- rnorm(7, rep(c(20, 10, 5), c(1, 2, 4)))
# (70 x 1) in-sample residuals vector (simulated)
res <- rnorm(70)
# (4 x 1) Naive high frequency base forecasts vector: all forecasts are set equal to 2.5
naive <- rep(2.5, 4)

## EXOGENOUS CONSTRAINTS
# Level Conditional Coherent (LCC) reconciled forecasts
exo_LC <- telcc(base = base, agg_order = 4, comb = "wlsh", hfts = naive,
                res = res, nodes = "auto", CCC = FALSE)

# Combined Conditional Coherent (CCC) reconciled forecasts
exo_CCC <- telcc(base = base, agg_order = 4, comb = "wlsh", hfts = naive,
                 res = res, nodes = "auto", CCC = TRUE)

# Results detailed by level:
info_exo <- recoinfo(exo_CCC, verbose = FALSE)
# info_exo$lcc

## ENDOGENOUS CONSTRAINTS
# Level Conditional Coherent (LCC) reconciled forecasts
endo_LC <- telcc(base = base, agg_order = 4, comb = "wlsh", res = res,
                 nodes = "auto", CCC = FALSE, const = "endogenous")

# Combined Conditional Coherent (CCC) reconciled forecasts
endo_CCC <- telcc(base = base, agg_order = 4, comb = "wlsh", res = res,
                  nodes = "auto", CCC = TRUE, const = "endogenous")

# Results detailed by level:
info_endo <- recoinfo(endo_CCC, verbose = FALSE)
# info_endo$lcc

Temporal middle-out reconciliation

Description

The middle-out forecast reconciliation for temporal hierarchies combines top-down (tetd) and bottom-up (tebu) methods. Given the base forecasts of an intermediate temporal aggregation order $k$ , it performs

a top-down approach for the aggregation orders $<k$ ;
a bottom-up approach for the aggregation orders $>k$ .

Usage

temo(base, agg_order, order = max(agg_order), weights, tew = "sum",
     normalize = TRUE)
temo(base, agg_order, order = max(agg_order), weights, tew = "sum",
     normalize = TRUE)

Arguments

`base`	A ( $hk \times 1$ ) numeric vector containing the temporal aggregated base forecasts of order $k$ ; $k$ is an aggregation order (a factor of $m$ , and $1<k<m$ ), $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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$ .
`order`	The intermediate fixed aggregation order $k$ .
`weights`	A ( $hm \times 1$ ) numeric vector containing the proportions for the high-frequency time series; $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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).
`normalize`	If `TRUE` (default), the `weights` will sum to 1.

Value

A ( $h(k^\ast+m) \times 1$ ) numeric vector of temporal reconciled forecasts.

Examples

set.seed(123)
# (6 x 1) base forecasts vector (simulated), forecast horizon = 3
# and intermediate aggregation order k = 2 (max agg order = 4)
basek2 <- rnorm(3*2, 5)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- temo(base = basek2, order = 2, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- runif(4*3)
recoh <- temo(base = basek2, order = 2, agg_order = 4, weights = h_weights)

set.seed(123)
# (6 x 1) base forecasts vector (simulated), forecast horizon = 3
# and intermediate aggregation order k = 2 (max agg order = 4)
basek2 <- rnorm(3*2, 5)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- temo(base = basek2, order = 2, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- runif(4*3)
recoh <- temo(base = basek2, order = 2, agg_order = 4, weights = h_weights)

Projection matrix for optimal combination temporal reconciliation

Description

This function computes the projection or the mapping matrix $\mathbf{M}$ and $\mathbf{G}$ , respectively, such that $\widetilde{\mathbf{y}} = \mathbf{M}\widehat{\mathbf{y}} = \mathbf{S}_{te}\mathbf{G}\widehat{\mathbf{y}}$ , where $\widetilde{\mathbf{y}}$ is the vector of the reconciled forecasts, $\widehat{\mathbf{y}}$ is the vector of the base forecasts, $\mathbf{S}_{te}$ is the temporal structural matrix, and $\mathbf{M} = \mathbf{S}_{te}\mathbf{G}$ . For further information regarding on the structure of these matrices, refer to Girolimetto et al. (2023).

Usage

teprojmat(agg_order, comb = "ols", res = NULL, mat = "M", tew = "sum", ...)
teprojmat(agg_order, comb = "ols", res = NULL, mat = "M", tew = "sum", ...)

Arguments

`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 tecov.
`res`	A ( $N(k^\ast+m) \times 1$ ) optional numeric vector containing the in-sample residuals at all the temporal frequencies ordered from the lowest frequency to the highest frequency. This vector is used to compute come covariance matrices.
`mat`	A string specifying which matrix to return: "`M`" (default) for $\mathbf{M}$ and "`G`" for $\mathbf{G}$ .
`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).
`...`	Arguments passed on to `tecov` `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)

Value

The projection matrix $\mathbf{M}$ (mat = "M") or the mapping matrix $\mathbf{G}$ (mat = "G").

References

Examples

# Temporal framework (annual-quarterly)
Mte <- teprojmat(agg_order = 4, comb = "ols")
Gte <- teprojmat(agg_order = 4, comb = "ols", mat = "G")

# Temporal framework (annual-quarterly)
Mte <- teprojmat(agg_order = 4, comb = "ols")
Gte <- teprojmat(agg_order = 4, comb = "ols", mat = "G")

Optimal combination temporal reconciliation

Description

This function performs forecast reconciliation for a single time series using temporal hierarchies (Athanasopoulos et al., 2017, Nystrup et al., 2020). The reconciled forecasts can be computed 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

terec(base, agg_order, comb = "ols", res = NULL, tew = "sum",
      approach = "proj", nn = NULL, settings = NULL, bounds = NULL,
      immutable = NULL, ...)
terec(base, agg_order, comb = "ols", res = NULL, tew = "sum",
      approach = "proj", nn = NULL, settings = NULL, bounds = NULL,
      immutable = NULL, ...)

Arguments

`base`	A ( $h(k^\ast + m) \times 1$ ) numeric vector containing base forecasts to be reconciled ordered from the lowest frequency to the highest frequency; $m$ is the max aggregation order, $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.
`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 tecov.
`res`	A ( $N(k^\ast+m) \times 1$ ) optional numeric vector containing the in-sample residuals at all the temporal frequencies ordered from the lowest frequency to the highest frequency. This vector is used to compute come 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: "`proj`" (default): Projection approach according to Byron (1978, 1979). "`strc`": Structural approach as proposed by Hyndman et al. (2011). "`proj_osqp`": Numerical solution using osqp for projection approach. "`strc_osqp`": Numerical solution using osqp for structural approach.
`nn`	A string specifying the algorithm to compute non-negative forecasts: "`osqp`": quadratic programming optimization (osqp solver). "`bpv`": block principal pivoting algorithm. "`sntz`": heuristic "set-negative-to-zero" (Di Fonzo and Girolimetto, 2023).
`settings`	A list of control parameters. `nn = "osqp"` An object of class `osqpSettings` specifying settings for the osqp solver. For details, refer to the osqp documentation (Stellato et al., 2020). `nn = "bpv"` It includes: `ptype` for permutation method ("`random`" or "`fixed`", default), `par` for the number of full exchange rules that may be attempted (`10`, default), `tol` for the tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `gtol` for the gradient tolerance in convergence criteria (`sqrt(.Machine$double.eps)`, default), `itmax` for the maximum number of algorithm iterations (`100`, default)
`bounds`	A matrix (see set_bounds) with 4 columns ( $k,j,lower,upper$ ), such that Column 1 represents the temporal aggregation order ( $k = m,\dots,1$ ). Column 2 represents the temporal forecast horizon ( $j = 1,\dots,m/k$ ). Columns 3 and 4 indicates the lower and lower bounds, respectively.
`immutable`	A matrix with 2 columns ( $k,j$ ), such that Column 1 represents the temporal aggregation order ( $k = m,\dots,1$ ). Column 2 represents the temporal forecast horizon ( $j = 1,\dots,m/k$ ). For example, when working with a quarterly time series: `t(c(4, 1))` - Fix the one step ahead annual forecast. `t(c(1, 2))` - Fix the two step ahead quarterly forecast.
`...`	Arguments passed on to `tecov` `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)

Value

A ( $h(k^\ast+m) \times 1$ ) numeric vector of temporal reconciled forecasts.

References

Athanasopoulos, G., Hyndman, R.J., Kourentzes, N. and Petropoulos, F. (2017), Forecasting with Temporal Hierarchies, European Journal of Operational Research, 262, 1, 60-74. doi:10.1016/j.ejor.2017.02.046

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. (2023), Spatio-temporal reconciliation of solar forecasts, Solar Energy, 251, 13–29. doi:10.1016/j.solener.2023.01.003

Nystrup, P., Lindström, E., Pinson, P. and Madsen, H. (2020), Temporal hierarchies with autocorrelation for load forecasting, European Journal of Operational Research, 280, 1, 876-888. doi:10.1016/j.ejor.2019.07.061

Examples

set.seed(123)
# (7 x 1) base forecasts vector (simulated), m = 4
base <- rnorm(7, rep(c(20, 10, 5), c(1, 2, 4)))
# (70 x 1) in-sample residuals vector (simulated)
res <- rnorm(70)

m <- 4 # from quarterly to annual temporal aggregation
reco <- terec(base = base, agg_order = m, comb = "wlsv", res = res)

# Immutable reconciled forecast
# E.g. fix all the quarterly forecasts
imm_q <- expand.grid(k = 1, j = 1:4)
immreco <- terec(base = base, agg_order = m, comb = "wlsv",
                 res = res, immutable = imm_q)

# Non negative reconciliation
base[7] <- -base[7] # Making negative one of the quarterly base forecasts
nnreco <- terec(base = base, agg_order = m, comb = "wlsv",
                res = res, nn = "osqp")
recoinfo(nnreco, verbose = FALSE)$info

set.seed(123)
# (7 x 1) base forecasts vector (simulated), m = 4
base <- rnorm(7, rep(c(20, 10, 5), c(1, 2, 4)))
# (70 x 1) in-sample residuals vector (simulated)
res <- rnorm(70)

m <- 4 # from quarterly to annual temporal aggregation
reco <- terec(base = base, agg_order = m, comb = "wlsv", res = res)

# Immutable reconciled forecast
# E.g. fix all the quarterly forecasts
imm_q <- expand.grid(k = 1, j = 1:4)
immreco <- terec(base = base, agg_order = m, comb = "wlsv",
                 res = res, immutable = imm_q)

# Non negative reconciliation
base[7] <- -base[7] # Making negative one of the quarterly base forecasts
nnreco <- terec(base = base, agg_order = m, comb = "wlsv",
                res = res, nn = "osqp")
recoinfo(nnreco, verbose = FALSE)$info

Temporal top-down reconciliation

Description

Top-down forecast reconciliation for a univariate time series, where the forecast of the most aggregated temporal level is disaggregated according to a proportional scheme (weights). Besides fulfilling any aggregation constraint, the top-down reconciled forecasts should respect two main properties:

the top-level value remains unchanged;
all the bottom time series reconciled forecasts are non-negative.

Usage

tetd(base, agg_order, weights, tew = "sum", normalize = TRUE)
tetd(base, agg_order, weights, tew = "sum", normalize = TRUE)

Arguments

`base`	A ( $hm \times 1$ ) numeric vector containing the temporal aggregated base forecasts of order $m$ ; $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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$ .
`weights`	A ( $hm \times 1$ ) numeric vector containing the proportions for the high-frequency time series; $m$ is the max aggregation order, and $h$ is the forecast horizon for the lowest frequency time series.
`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).
`normalize`	If `TRUE` (default), the `weights` will sum to 1.

Value

A ( $h(k^\ast+m) \times 1$ ) numeric vector of temporal reconciled forecasts.

Examples

set.seed(123)
# (2 x 1) top base forecasts vector (simulated), forecast horizon = 2
topf <- rnorm(2, 10)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- tetd(base = topf, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- runif(4*2)
recoh <- tetd(base = topf, agg_order = 4, weights = h_weights)

set.seed(123)
# (2 x 1) top base forecasts vector (simulated), forecast horizon = 2
topf <- rnorm(2, 10)
# Same weights for different forecast horizons
fix_weights <- runif(4)
reco <- tetd(base = topf, agg_order = 4, weights = fix_weights)

# Different weights for different forecast horizons
h_weights <- runif(4*2)
recoh <- tetd(base = topf, agg_order = 4, weights = h_weights)

Temporal reconciliation tools

Description

Some useful tools for forecast reconciliation through temporal hierarchies.

Usage

tetools(agg_order, fh = 1, tew = "sum", sparse = TRUE)
tetools(agg_order, fh = 1, tew = "sum", sparse = TRUE)

Arguments

`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$ .
`fh`	Forecast horizon for the lowest frequency (most temporally aggregated) time series (default is `1`).
`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).
`sparse`	Option to return sparse matrices (default is `TRUE`).

Value

A list with five elements:

`dim`	A vector containing information about the maximum aggregation order (`m`), the number of factor (`p`), the partial (`ks`) and total sum (`kt`) of factors.
`set`	The vector of the temporal aggregation orders (in decreasing order).
`agg_mat`	The temporal linear combination or aggregation matrix.
`strc_mat`	The temporal structural matrix.
`cons_mat`	The temporal zero constraints matrix.

Examples

# Temporal framework (quarterly data)
obj <- tetools(agg_order = 4, sparse = FALSE)

# Temporal framework (quarterly data)
obj <- tetools(agg_order = 4, sparse = FALSE)

Aggregation matrix of a balanced hierarchy in (possibly) unbalanced form

Description

A hierarchy with $L$ upper levels is said to be balanced if each variable at level $l$ has at least one child at level $l+1$ . When this doesn't hold, the hierarchy is unbalanced. This function transforms an aggregation matrix of a balanced hierarchy into an aggregation matrix of an unbalanced one, by removing possible duplicated series.

Usage

unbalance_hierarchy(agg_mat, more_info = FALSE, sparse = TRUE)
unbalance_hierarchy(agg_mat, more_info = FALSE, sparse = TRUE)

Arguments

`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.
`more_info`	If `TRUE`, it returns only the aggregation matrix of the unbalanced hierarchy. Default is `FALSE`.
`sparse`	Option to return sparse matrices (default is `TRUE`).

Value

A list containing four elements (more_info = TRUE):

`ubm`	The aggregation matrix of the unbalanced hierarchy.
`agg_mat`	The input matrix.
`idrm`	The identification number of the duplicated variables (row numbers of the aggregation matrix `agg_mat`).
`id`	The identification number of each variable in the balanced hierarchy. It may contains duplicated values.

Examples

#     Balanced     ->     Unbalanced
#        T                    T
#    |-------|            |-------|
#    A       B            A       |
#  |---|     |          |---|     |
# AA   AB    BA        AA   AB    BA
A <- matrix(c(1, 1, 1,
              1, 1, 0,
              0, 0, 1), 3, byrow = TRUE)
obj <- unbalance_hierarchy(agg_mat = A)
obj
#     Balanced     ->     Unbalanced
#        T                    T
#    |-------|            |-------|
#    A       B            A       |
#  |---|     |          |---|     |
# AA   AB    BA        AA   AB    BA
A <- matrix(c(1, 1, 1,
              1, 1, 0,
              0, 0, 1), 3, byrow = TRUE)
obj <- unbalance_hierarchy(agg_mat = A)
obj

Australian Tourism Demand dataset

Description

The Australian Tourism Demand dataset (Wickramasuriya et al. 2019) measures the number of nights Australians spent away from home. It includes 228 monthly observations of Visitor Nights (VNs) from January 1998 to December 2016, and has a cross-sectional grouped structure based on a geographic hierarchy crossed by purpose of travel. The geographic hierarchy comprises 7 states, 27 zones, and 76 regions, for a total of 111 nested geographic divisions. Six of these zones are each formed by a single region, resulting in 105 unique nodes in the hierarchy. The purpose of travel comprises four categories: holiday, visiting friends and relatives, business, and other. To avoid redundancies (Girolimetto et al. 2023), 24 nodes (6 zones are formed by a single region) are not considered, resulting in an unbalanced hierarchy of 525 (304 bottom and 221 upper time series) unique nodes instead of the theoretical 555 with duplicated nodes.

Usage

# 525 time series of the Australian Tourism Demand dataset
vndata

# aggregation matrix
vnaggmat
# 525 time series of the Australian Tourism Demand dataset
vndata

# aggregation matrix
vnaggmat

Format

vndata is a $(228 \times 525)$ ts object, corresponding to 525 time series of the Australian Tourism Demand dataset (1998:01-2016:12).

vnaggmat is the $(221 \times 304)$ aggregation matrix.

Source

https://robjhyndman.com/publications/mint/

Package 'FoReco'

Help Index

FoReco: Forecast Reconciliation

Description

Author(s)

See Also

Non-overlapping temporal aggregation of a time series

Description

Usage

Arguments

Value

See Also

Examples

Aggregation matrix of a (possibly) unbalanced hierarchy in balanced form

Description

Usage

Arguments

Value

See Also

Examples

Commutation matrix

Description

Usage

Arguments

Value

References

See Also

Examples

Cross-sectional joint block bootstrap

Description

Usage

Arguments

Value

References

See Also

Cross-sectional bottom-up reconciliation

Description

Usage

Arguments

Value

References

See Also

Examples

Cross-sectional covariance matrix approximation

Description

Usage

Arguments

Value

References

See Also

Examples

Level conditional coherent reconciliation for genuine hierarchical/grouped time series

Description

Usage

Arguments

Value

References

See Also

Examples

Cross-sectional middle-out reconciliation

Description

Usage

Arguments

Value

References

See Also

Examples

Projection matrix for optimal combination cross-sectional reconciliation

Description

Usage

Arguments

Value

References

See Also

Examples

Optimal combination cross-sectional reconciliation

Description

Usage

Arguments

Value