Probabilistic forecast reconciliation with bayesRecon

bayesRecon

Probabilistic forecast

Probabilistic forecast

Hierarchical time series

Constraints: \(U_2 = B_1 + B_2\), \(U_3 = B_3 + B_4\), \(U_1 = U_2 + U_3\)

Let \(\mathbf{Y} = \begin{bmatrix}\mathbf{U}^T \quad \mathbf{B}^T\end{bmatrix}^T = \begin{bmatrix}U_1 \\ U_2 \\ \dots \\ B_4 \end{bmatrix}\) then \(\mathbf{U} = A \mathbf{B} \quad \text{ where } A= \begin{bmatrix}1 & 1 & 1 & 1\\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1\end{bmatrix}\)

Hierarchical time series

Probabilistic forecast for hierarchical time series

Probabilistic forecast reconciliation

\[ \Big\Downarrow %\;\text{ reconciliation} \]

Probabilistic reconciliation via conditioning

Well-defined for both continuous and discrete distributions¹

Probabilistic reconciliation via conditioning

Gaussian base forecasts¹: reconciliation is analytic and coincides with MinT²

In general: need to sample from the reconciled distribution

Case 1: Gaussian

Input:

• A: aggregation matrix (\(n_{upp} \times n_{bott}\))

• base_forecasts_mu: list of length \(n\) representing the mean of the base forecast distribution

• base_forecasts_Sigma : \(n \times n\) covariance matrix for the base forecast distribution

Output:

• reconc_gaussian['bottom_reconciled_mean'] : mean of the bottom reconciled distribution

• reconc_gaussian['bottom_reconciled_covariance'] : covariance of the bottom reconciled distribution

Python Example: Infant mortality dataset

Dataset considered: yearly counts of infant mortality in Australia, disaggregated by state and sex [‘hts’ package (Hyndman et al. 2021)]

Gaussian reconciliation

python

print(sys.version)

3.13.5 | packaged by conda-forge | (main, Jun 16 2025, 08:25:50) [Clang 18.1.8 ]

python

recon_gauss = reconc_gaussian(A=A, base_forecasts_mu=mu, base_forecasts_Sigma=Sigma)

# Extract reconciled means and covariances
bottom_mu_reconc = recon_gauss['bottom_reconciled_mean']
bottom_Sigma_reconc = recon_gauss['bottom_reconciled_covariance']

upper_mu_reconc = A @ bottom_mu_reconc  
upper_Sigma_reconc = A @ bottom_Sigma_reconc @ A.T

Gaussian reconciliation

python

recon_gauss = reconc_gaussian(A=A, base_forecasts_mu=mu, base_forecasts_Sigma=Sigma)

# Extract reconciled means and covariances
bottom_mu_reconc = recon_gauss['bottom_reconciled_mean']
bottom_Sigma_reconc = recon_gauss['bottom_reconciled_covariance']

upper_mu_reconc = A @ bottom_mu_reconc  
upper_Sigma_reconc = A @ bottom_Sigma_reconc @ A.T

Case 2: non-Gaussian

Bottom-up importance sampling (BUIS) algorithm ¹

Input:

• A: aggregation matrix (\(n_{upp} \times n_{bott}\))

• base_forecasts: list of length \(n\) representing the base forecast distribution

• in_type: either “params” or “samples”

• distr: specifies the parametric distribution (e.g. Poisson)
  or the type of samples (continuous or discrete)

Output:

• samples from the reconciled distribution

Example BUIS: reconciliation of extreme market events

• Counts of extreme market events in 5 economic sectors, daily time series,
  period: 2005-2018 (3508 trading days)

Example BUIS: reconciliation of extreme market events

• Counts of extreme market events in 5 economic sectors, daily time series,
  period: 2005-2018 (3508 trading days)

Example BUIS: reconciliation of extreme market events

r

j = 123  # pick one of the trading days
A <- matrix(1, ncol = 5, nrow = 1)  # aggregation matrix (5 bottom, 1 upper)

# base forecasts: list of length n with the parameters of the negbin fc distr
base_fc_j <- c()   
for (i in 1:6) {
  base_fc_j[[i]] <- list(size = extr_mkt_events_basefc$size[[i]], 
                         mu   = extr_mkt_events_basefc$mu[[j,i]])
}

# Reconcile via importance sampling:
tic()
buis <- reconc_BUIS(A, base_fc_j, "params", "nbinom", num_samples = 1e4, seed = 42)
toc()

0.018 sec elapsed

r

samples_y <- buis$reconciled_samples   # matrix: n x N_samples
rownames(samples_y) <- c("ALL", "FIN", "ICT", "MFG", "ENG", "TRD")
head(t(samples_y), n=3)

     ALL FIN ICT MFG ENG TRD
[1,]   5   3   0   1   0   1
[2,]   7   2   1   3   0   1
[3,]   3   2   1   0   0   0

Example BUIS: temporal reconciliation

r

aggr_ts <- bayesRecon::temporal_aggregation(item_ts, agg_levels = c(7, 14, 28))

Example BUIS: temporal reconciliation

Compute base forecasts using ‘smooth’ package (Svetunkov, 2023)

r

n_samples <- 1e4
base_fc_samples <- matrix(nrow = 0, ncol = n_samples)
  
model <- smooth::adam(aggr_ts[[1]], lags=c(1, 365%/%28), distribution=c("dgamma"))
fc <- forecast(model, h = 1, interval="simulated", scenarios=TRUE, nsim = n_samples)
base_fc_samples <- rbind(base_fc_samples, round((fc$scenarios)))

model <- smooth::adam(aggr_ts[[2]], lags=c(1, 365%/%14), distribution = c("dgamma"))
fc <- forecast(model, h = 2, interval="simulated", scenarios=TRUE, nsim = n_samples)
base_fc_samples <- rbind(base_fc_samples, round(fc$scenarios))

model <- smooth::adam(aggr_ts[[3]], lags=c(1, 365%/%7), distribution = c("dgamma"), occurrence="auto")
fc <- forecast(model, h = 4, interval="simulated", scenarios=TRUE, nsim = n_samples)
base_fc_samples <- rbind(base_fc_samples, round(fc$scenarios))

model <- smooth::adam(aggr_ts[[4]], lags=c(1,7,365), distribution=c("dgamma"), occurrence="auto")
fc <- forecast(model, h = 28, interval="simulated", scenarios=TRUE, nsim = n_samples)
base_fc_samples <- rbind(base_fc_samples, round(fc$scenarios))

dim(base_fc_samples)

[1]    35 10000

Example BUIS: temporal reconciliation

r

A <- bayesRecon::get_reconc_matrices(agg_levels = c(7, 14, 28), h = 28)$A
base_fc = apply(base_fc_samples, 1, as.vector, simplify = F)  # list of vectors

tic()
buis <- reconc_BUIS(A, base_forecasts = base_fc, in_type = "samples", distr = "discrete")
toc()

0.046 sec elapsed

Case 3: Mixed: discrete bottom & continuous upper

• Both the reconciling functions reconc_MixCond and reconc_TDcond ¹ have the same input and output structure.

Input:

• A: aggregation matrix (\(n_{upp} \times n_{bott}\))

• fc_bottom: base bottom forecasts

• fc_upper : base upper forecasts, represented as parameters of a multivariate Gaussian

• bottom_in_type : either “params”, “pmf”, or “samples”

• distr : specifies the parametric distribution (only used for “params”)

• return_type : specifies the return type - either “pmf”, “samples”, or “all”.

Python Example: Mixed-case for the M5 dataset

• Dataset considered: M5 [Makridakis et. al 2022] containing daily time series for sales of stores (non-negative).
• Considered store: CA_1; forecast horizon: h=1; \(n_b=3049\) and \(n_u =11\).

python

# Preparing the base forecasts

Sigma_u_np = np.array(Sigma_u['Sigma_u'])
fc_upper_4rec = {'mu': mu_u, 'Sigma': Sigma_u_np}  # Dictionary for upper forecasts
fc_bottom_4rec = {k: np.array(fc['pmf']) for k, fc in base_fc_bottom.items()}

Mixed-conditioning

• Base forecasts are stored in the package as M5_CA1_basefc which were originally obtained in the R counterpart using ADAM [Svetunkov et. al 2023].

python

from bayesreconpy.reconc_mix_cond import reconc_mix_cond

np.random.seed(seed)
start = time.time()
mix_cond = reconc_mix_cond(A, fc_bottom_4rec, fc_upper_4rec, bottom_in_type="pmf",
                           num_samples=N_samples_IS, return_type="pmf", seed=seed)

stop = time.time()
rec_fc['Mixed_cond'] = {
  'bottom': mix_cond['bottom_reconciled']['pmf'],  # Bottom-level reconciled PMFs
  'upper': mix_cond['upper_reconciled']['pmf'],    # Upper-level reconciled PMFs
  'ESS': mix_cond['ESS']                           # Effective Sample Size (ESS)
}

MixCond_time = round(stop - start, 2)
print(f"Computational time for Mix-cond reconciliation: {MixCond_time} seconds")

Computational time for Mix-cond reconciliation: 9.5 seconds

• The return_type="pmf" returns the pmf of the reconciled distribution with non-negative integer support in a memory-efficient way.

Top-down conditioning

python

from bayesreconpy.reconc_td_cond import reconc_td_cond

N_samples_TD = int(1e4)
start = time.time()

td = reconc_td_cond(A, fc_bottom_4rec, fc_upper_4rec,
                    bottom_in_type="pmf", num_samples=N_samples_TD,
                    return_type="pmf", seed=seed)

python

  
stop = time.time()
rec_fc['TD_cond'] = {
  'bottom': td['bottom_reconciled']['pmf'],
  'upper': td['upper_reconciled']['pmf']
}

TDCond_time = round(stop - start, 2)
print(f"Computational time for TD-cond reconciliation: {TDCond_time} seconds")

Computational time for TD-cond reconciliation: 10.28 seconds

bayesRecon / bayesreconpy

• fully probabilistic: produces reconciled distribution

• versatile: deals with many different cases: discrete, continuous, mixed

• efficient: very fast, able to scale to large hierarchies
  

Note

• Both packages are at version 0.x, there might be changes in the future.

• A fable wrapper for the R package might come soon!

bayesRecon / bayesreconpy

• fully probabilistic: produces reconciled distribution

• versatile: deals with many different cases: discrete, continuous, mixed

• efficient: very fast, able to scale to large hierarchies

💡 Those are active research projects on GitHub — please report any issues or suggestions.

Probabilistic reconciliation via conditioning

Geometric intuition: slice the joint forecast distribution along the coherent subspace \(\mathcal{S}\)