t-Rec: reconciliation via conditioning with uncertain covariance via multivariate t-distribution

Reconciles base forecasts in a hierarchy by conditioning on the hierarchical constraints, specified by the aggregation matrix A. The base forecasts are assumed to be jointly Gaussian, conditionally on the covariance matrix of the forecast errors. A Bayesian approach is adopted to account for the uncertainty of the covariance matrix. An Inverse-Wishart prior is specified on the covariance matrix, leading to a multivariate t-distribution for the base forecasts. The reconciliation via conditioning is in closed-form, yielding a multivariate t reconciled distribution.

Usage

reconc_t(
  A,
  base_fc_mean,
  y_train = NULL,
  residuals = NULL,
  ...,
  return_upper = FALSE,
  return_parameters = FALSE
)

Arguments

A: Matrix (n_upp x n_bott) defining the hierarchy (u = Ab).
base_fc_mean: Vector of base forecasts (length n = n_upp + n_bott).
y_train: mts (or matrix) of historical training data (T x n) used for setting prior parameters.
residuals: Matrix (T x n) of base forecast residuals.
...: Additional arguments for advanced usage: see details.
return_upper: Logical; if TRUE, also returns parameters for the upper-level reconciled distribution.
return_parameters: Logical; if TRUE, also returns prior and posterior parameters.

Value

A list containing:

bottom_rec_mean: reconciled bottom-level mean.
bottom_rec_scale_matrix: reconciled bottom-level scale matrix.
bottom_rec_df: reconciled degrees of freedom.

If return_upper is TRUE, also returns:

upper_rec_mean: reconciled upper-level mean.
upper_rec_scale_matrix: reconciled upper-level scale matrix.
upper_rec_df: reconciled upper-level degrees of freedom.

If return_parameters is TRUE, also returns:

prior_nu: prior degrees of freedom.
prior_Psi: prior scale matrix.
posterior_nu: posterior degrees of freedom.
posterior_Psi: posterior scale matrix.

Details

Standard usage.

The standard workflow for this function is to provide the in-sample residuals and the historical training data y_train. The parameters of the Inverse-Wishart prior distribution of the covariance matrix are set as follows:

Prior scale matrix (Psi): set as the covariance of the residuals of naive (or seasonal naive, a criterion is used to choose between the two) forecasts computed on y_train.
Prior degrees of freedom (nu): estimated via Bayesian Leave-One-Out Cross-Validation (LOOCV) to maximize out-of-sample performance.

The posterior distribution of the covariance matrix is still Inverse-Wishart. The parameters of the posterior are computed in closed-form using the sample covariance of the provided residuals.

Advanced Options.

Users can bypass the automated estimation by specifying:

prior: a list with entries 'nu' and 'Psi'. This skips the LOOCV step for \(\nu_0\) and the covariance estimation from y_train. It requires residuals to compute the posterior.
posterior: a list with entries 'nu' and 'Psi'. This skips all internal estimation and updating logic.

Moreover, users can specify:

freq: positive integer, used as frequency of data for the seasonal naive forecast in the specification of the prior scale matrix. By default, if y_train is a multivariate time series, the frequency of the data is used; otherwise, it is set to 1 (no seasonality).
criterion: either 'RSS' (default) or 'seas-test', specifying which criterior is used to choose between the naive and seasonal naive forecasts for the specification of the prior scale matrix. 'RSS' computes the residual sum of squares for both methods and chooses the one with lower RSS, while 'seas-test' uses a statistical test for seasonality (currently implemented using the number of seasonal differences suggested by the forecast package, which must be installed).

Reconciled distribution.

The reconciled distribution is a multivariate t-distribution, specified by a vector of means, a scale matrix, and a number of degrees of freedom. These parameters are computed in closed-form. By default, only the parameters of the reconciled distribution for the bottom-level series are returned. See examples.

References

Carrara, C., Corani, G., Azzimonti, D., & Zambon, L. (2025). Modeling the uncertainty on the covariance matrix for probabilistic forecast reconciliation. arXiv preprint arXiv:2506.19554. https://arxiv.org/abs/2506.19554

Examples


# \donttest{
library(bayesRecon)

if (requireNamespace("forecast", quietly = TRUE)) {

  set.seed(1234)
  n_obs <- 100

  # Simulate 2 bottom series from AR(1) processes
  y1 <- arima.sim(model = list(ar = 0.8), n = n_obs)
  y2 <- arima.sim(model = list(ar = 0.5), n = n_obs)

  y_upper <- y1 + y2   # upper series is the sum of the two bottoms
  A <- matrix(c(1, 1), nrow = 1)  # Aggregation matrix

  # Fit additive ETS models
  fit1 <- forecast::ets(y1, additive.only = TRUE)
  fit2 <- forecast::ets(y2, additive.only = TRUE)
  fit_upper <- forecast::ets(y_upper, additive.only = TRUE)

  # Point forecasts (h = 1)
  fc_upper <- as.numeric(forecast::forecast(fit_upper, h = 1)$mean)
  fc1 <- as.numeric(forecast::forecast(fit1, h = 1)$mean)
  fc2 <- as.numeric(forecast::forecast(fit2, h = 1)$mean)
  base_fc_mean <- c(fc_upper, fc1, fc2)

  # Residuals and training data (n_obs x n matrices, columns in same order as base_fc_mean)
  res <- cbind(residuals(fit_upper), residuals(fit1), residuals(fit2))
  y_train <- cbind(y_upper, y1, y2)

  # --- 1) Generate joint reconciled samples ---
  result <- reconc_t(A, base_fc_mean, y_train = y_train, residuals = res)

  # Sample from the reconciled bottom-level t-distribution
  n_samples <- 2000
  L_chol <- t(chol(result$bottom_rec_scale_matrix))
  z <- matrix(rt(ncol(A) * n_samples, df = result$bottom_rec_df), nrow = ncol(A))
  bottom_samples <- result$bottom_rec_mean + L_chol %*% z  # 2 x n_samples

  # Aggregate bottom samples to get upper samples
  upper_samples <- A %*% bottom_samples              
  joint_samples <- rbind(upper_samples, bottom_samples)
  rownames(joint_samples) <- c("upper", "bottom_1", "bottom_2")

  cat("Reconciled means (from samples):\n")
  print(round(rowMeans(joint_samples), 3))

  cat("Reconciled standard deviations (from samples):\n")
  print(round(apply(joint_samples, 1, sd), 3))

  # --- 2) 95% prediction intervals via t-distribution quantiles ---
  result2 <- reconc_t(A, base_fc_mean, y_train = y_train,
                      residuals = res, return_upper = TRUE)

  alpha <- 0.05
  # Bottom series intervals
  for (i in seq_len(ncol(A))) {
    s_i <- sqrt(result2$bottom_rec_scale_matrix[i, i])
    lo  <- result2$bottom_rec_mean[i] + s_i * qt(alpha / 2,     df = result2$bottom_rec_df)
    hi  <- result2$bottom_rec_mean[i] + s_i * qt(1 - alpha / 2, df = result2$bottom_rec_df)
    cat(sprintf("Bottom %d: 95%% PI = [%.3f, %.3f]\n", i, lo, hi))
  }
  # Upper series interval
  s_u <- sqrt(result2$upper_rec_scale_matrix[1, 1])
  lo  <- result2$upper_rec_mean[1] + s_u * qt(alpha / 2,     df = result2$upper_rec_df)
  hi  <- result2$upper_rec_mean[1] + s_u * qt(1 - alpha / 2, df = result2$upper_rec_df)
  cat(sprintf("Upper:    95%% PI = [%.3f, %.3f]\n", lo, hi))
}
#> Reconciled means (from samples):
#>    upper bottom_1 bottom_2 
#>   -2.370   -0.056   -2.314 
#> Reconciled standard deviations (from samples):
#>    upper bottom_1 bottom_2 
#>    1.574    0.947    1.212 
#> Bottom 1: 95% PI = [-1.917, 1.777]
#> Bottom 2: 95% PI = [-4.640, 0.062]
#> Upper:    95% PI = [-5.377, 0.659]
# }