Skip to contents

MCMC for Probabilistic Reconciliation of forecasts via conditioning

Source: R/reconc_MCMC.R
reconc_MCMC.Rd

Uses Markov Chain Monte Carlo algorithm to draw samples from the reconciled forecast distribution, which is obtained via conditioning.

This is a bare-bones implementation of the Metropolis-Hastings algorithm, we suggest the usage of tools to check the convergence. The function only works with Poisson or Negative Binomial base forecasts.

The function reconc_BUIS() is generally faster on most hierarchies.

Usage

reconc_MCMC(
  A,
  base_forecasts,
  distr,
  num_samples = 10000,
  tuning_int = 100,
  init_scale = 1,
  burn_in = 1000,
  seed = NULL
)

Arguments

A

aggregation matrix (n_upper x n_bottom).

base_forecasts

list of the parameters of the base forecast distributions, see details.

distr

a string describing the type of predictive distribution.

num_samples

number of samples to draw using MCMC.

tuning_int

number of iterations between scale updates of the proposal.

init_scale

initial scale of the proposal.

burn_in

number of initial samples to be discarded.

seed

seed for reproducibility.

Value

A list containing the reconciled forecasts. The list has the following named elements:

  • bottom_reconciled_samples: a matrix (n_bottom x num_samples) containing reconciled samples for the bottom time series;

  • upper_reconciled_samples: a matrix (n_upper x num_samples) containing reconciled samples for the upper time series;

  • reconciled_samples: a matrix (n x num_samples) containing the reconciled samples for all time series.

Details

The parameter base_forecast is a list containing n = n_upper + n_bottom elements. Each element is a list containing the estimated:

  • mean and sd for the Gaussian base forecast, see Normal, if distr='gaussian';

  • lambda for the Poisson base forecast, see Poisson, if distr='poisson';

  • size and prob (or mu) for the negative binomial base forecast, see NegBinomial, if distr='nbinom'.

The first n_upper elements of the list are the upper base forecasts, in the order given by the rows of A. The elements from n_upper+1 until the end of the list are the bottom base forecasts, in the order given by the columns of A.

References

Corani, G., Azzimonti, D., Rubattu, N. (2024). Probabilistic reconciliation of count time series. International Journal of Forecasting 40 (2), 457-469. doi:10.1016/j.ijforecast.2023.04.003 .

See also

reconc_BUIS()

Examples


library(bayesRecon)

# Create a minimal hierarchy with 2 bottom and 1 upper variable
rec_mat <- get_reconc_matrices(agg_levels=c(1,2), h=2)
A <- rec_mat$A

#Set the parameters of the Poisson base forecast distributions
lambda1 <- 2
lambda2 <- 4
lambdaY <- 9
lambdas <- c(lambdaY,lambda1,lambda2)

base_forecasts = list()
for (i in 1:length(lambdas)) {
 base_forecasts[[i]] = list(lambda = lambdas[i])
}

#Sample from the reconciled forecast distribution using MCMC
mcmc = reconc_MCMC(A, base_forecasts, distr = "poisson",
                  num_samples = 30000, seed = 42)
samples_mcmc <- mcmc$reconciled_samples

#Compare the reconciled means with those obtained via BUIS
buis = reconc_BUIS(A, base_forecasts, in_type="params",
                   distr="poisson", num_samples=100000, seed=42)
samples_buis <- buis$reconciled_samples

print(rowMeans(samples_mcmc))
#> [1] 7.091733 2.363633 4.728100
print(rowMeans(samples_buis))
#> [1] 7.09171 2.36542 4.72629