Probabilistic forecast reconciliation with bayesRecon

bayesRecon

open-source software for probabilistic forecast reconciliation
same functionalities for both R and Python
contributors:

Probabilistic forecast

At time \(t+h\), given observations
up to time \(t\):

parametric distribution

Probabilistic forecast

At time \(t+h\), given observations
up to time \(t\):

samples

		\(i = 1\)		\(i = 2\)		\(\dots\)		\(i = N\)
\(\hat{Y}_{t+h\|t}\)		2		3		\(\dots\)		1

Hierarchical time series

Constraints:
A = AA + AB, B = BA + BB, Z = A + B

AA, AB, BA, BB: bottom series

Z, A, B: upper series

\[ \mathcal{S} := \left\{\begin{bmatrix}z \\ a \\ \dots \\ bb \end{bmatrix} \;\text{such that }\; \begin{matrix} z = a + b,\\ a = aa+ab,\\ b = ba+bb \end{matrix} \right\} \]

coherence subspace

Probabilistic forecast for hierarchical time series

Base forecasts are incoherent:

Joint forecast distribution

Joint samples

	\(i = 1\)	\(i = 2\)	\(\dots\)	\(i = N\)
\(Z\)	8	6	\(\dots\)	9
\(A \quad\)	3	2	\(\dots\)	4
\(B\)	4	3	\(\dots\)	3
\(AA\)	1	1	\(\dots\)	2
\(AB\)	1	2	\(\dots\)	2
\(BA\)	2	0	\(\dots\)	1
\(BB\)	3	2	\(\dots\)	1

Probabilistic reconciliation

Incoherent base forecasts

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

reconciled forecast distribution

reconciled samples

	\(i = 1\)	\(i = 2\)	\(\dots\)	\(i = N\)
\(Z\)	8	6	\(\dots\)	9
\(A \quad\)	3	3	\(\dots\)	5
\(B\)	5	3	\(\dots\)	4
\(AA\)	1	1	\(\dots\)	3
\(AB\)	2	2	\(\dots\)	2
\(BA\)	2	1	\(\dots\)	1
\(BB\)	3	2	\(\dots\)	3

Probabilistic reconciliation via conditioning

Reconciliation via conditioning

Condition the incoherent base forecast distribution
on the hierarchical constraints

Well-defined for both continuous and discrete distributions¹

If base forecasts are Gaussian²:
reconciliation is analytic and coincides with MinT³

In general:
need to sample from the reconciled distribution

Hierarchy with different types of base forecasts

Case 1: Gaussian

reconc_gaussian(A, base_forecasts_mu, base_forecasts_Sigma)

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 : A \(n \times n\) covariance matrix representing the uncertainties in the base forecasts.

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

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 ¹

reconc_BUIS(A, base_forecasts, in_type, distr, num_samples)

Input:

A: aggregation matrix (\(n_{upp} \times n_{bott}\))
base_forecasts: list of length \(n\)
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

new dataset: we make available both data and base forecasts
data: counts of extreme market events in five economic sectors
daily time series, period: 2005-2018 (3508 trading days)
base forecasts: produced using score-driven model¹
1-step-ahead, negative binomial distribution

Example BUIS: reconciliation of extreme market events

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.03 sec elapsed

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

Daily sales of a single item from a Walmart store

Temporal hierarchy: 1 day - 1 week - 2 weeks - 4 weeks

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)

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

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.11 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.

reconc_MixCond(A, fc_bottom, fc_upper, bottom_in_type, distr, num_samples, return_type, seed)

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_MixCond import reconc_MixCond

np.random.seed(seed)
start = time.time()
mix_cond = reconc_MixCond(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: 17.77 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_TDcond import reconc_TDcond

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

td = reconc_TDcond(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: 22.81 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

This is the initial version for both the packages as there might be changes in the future.

💡 This is an active research project on GitHub — we welcome collaboration and encourage you to report any issues or suggest new methods or algorithms.

R package

Python package

Probabilistic reconciliation via conditioning

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

Reconciled density/pmf can be expressed analytically (up to a normalizing constant):

\[ \tilde{\pi}(\mathbf{b}) \propto \hat{\pi}(\mathbf{A}\mathbf{b},\, \mathbf{b}) \]