The package `bayesRecon`

implements several methods for
probabilistic reconciliation of hierarchical time series forecasts.

The main functions are:

`reconc_gaussian`

: reconciliation via conditioning of multivariate Gaussian base forecasts; this is done analytically;`reconc_BUIS`

: reconciliation via conditioning of any probabilistic forecast via importance sampling; this is the recommended option for non-Gaussian base forecasts;`reconc_MCMC`

: reconciliation via conditioning of discrete probabilistic forecasts via Markov Chain Monte Carlo;`reconc_MixCond`

: reconciliation via conditioning of mixed hierarchies, where the upper forecasts are multivariate Gaussian and the bottom forecasts are discrete distributions;`reconc_TDcond`

: reconciliation via top-down conditioning of mixed hierarchies, where the upper forecasts are multivariate Gaussian and the bottom forecasts are discrete distributions.

:boom: [2024-05-29] Added `reconc_MixCond`

and
`reconc_TDcond`

and the vignette “Reconciliation of M5
hierarchy with mixed-type forecasts”.

:boom: [2023-12-19] Added the vignette “Properties of the reconciled distribution via conditioning”.

:boom: [2023-08-23] Added the vignette “Probabilistic Reconciliation
via Conditioning with bayesRecon”. Added the
`schaferStrimmer_cov`

function.

:boom: [2023-05-26] bayesRecon v0.1.0 is released!

You can install the **stable** version on R CRAN

`install.packages("bayesRecon", dependencies = TRUE)`

You can also install the **development** version from Github

```
# install.packages("devtools")
::install_github("IDSIA/bayesRecon", build_vignettes = TRUE, dependencies = TRUE) devtools
```

Let us consider the minimal temporal hierarchy in the figure, where the bottom variables are the two 6-monthly forecasts and the upper variable is the yearly forecast. We denote the variables for the two semesters and the year by \(S_1, S_2, Y\) respectively.

The hierarchy is described by the *aggregation matrix* A,
which can be obtained using the function
`get_reconc_matrices`

.

```
library(bayesRecon)
<- get_reconc_matrices(agg_levels = c(1, 2), h = 2)
rec_mat <- rec_mat$A
A print(A)
#> [,1] [,2]
#> [1,] 1 1
```

We assume that the base forecasts are Poisson distributed, with parameters given by \(\lambda_{Y} = 9\), \(\lambda_{S_1} = 2\), and \(\lambda_{S_2} = 4\).

```
<- 2
lambdaS1 <- 4
lambdaS2 <- 9
lambdaY <- c(lambdaY, lambdaS1, lambdaS2)
lambdas = length(lambdas)
n_tot
= list()
base_forecasts for (i in 1:n_tot) {
= list(lambda = lambdas[i])
base_forecasts[[i]] }
```

We recommend using the BUIS algorithm (Zambon et al., 2024) to sample from the reconciled distribution.

```
<- reconc_BUIS(
buis
A,
base_forecasts,in_type = "params",
distr = "poisson",
num_samples = 100000,
seed = 42
)
<- buis$reconciled_samples samples_buis
```

Since there is a positive incoherence in the forecasts (\(\lambda_Y > \lambda_{S_1}+\lambda_{S_2}\)), the mean of the bottom reconciled forecast increases. We show below this behavior for \(S_1\).

```
<- buis$bottom_reconciled_samples[1,]
reconciled_forecast_S1 <- range(reconciled_forecast_S1)
range_forecats hist(
reconciled_forecast_S1,breaks = seq(range_forecats[1] - 0.5, range_forecats[2] + 0.5),
freq = F,
xlab = "S_1",
ylab = NULL,
main = "base vs reconciled"
)points(
seq(range_forecats[1], range_forecats[2]),
::dpois(seq(range_forecats[1], range_forecats[2]), lambda =
stats
lambdaS1),pch = 16,
col = 4,
cex = 2
)
```

The blue circles represent the probability mass function of a Poisson with parameter \(\lambda_{S_1}\) plotted on top of the histogram of the reconciled bottom forecasts for \(S_1\). Note how the histogram is shifted to the right.

Moreover, while the base bottom forecast were assumed independent, the operation of reconciliation introduced a negative correlation between \(S_1\) and \(S_2\). We can visualize it with the plot below which shows the empirical correlations between the reconciled samples of \(S_1\) and the reconciled samples of \(S_2\).

```
<-
AA xyTable(buis$bottom_reconciled_samples[1, ],
$bottom_reconciled_samples[2, ])
buisplot(
$x ,
AA$y ,
AAcex = AA$number * 0.001 ,
pch = 16 ,
col = rgb(0, 0, 1, 0.4) ,
xlab = "S_1" ,
ylab = "S_2" ,
xlim = range(buis$bottom_reconciled_samples[1, ]) ,
ylim = range(buis$bottom_reconciled_samples[2, ])
)
```

We also provide a function for sampling using Markov Chain Monte Carlo (Corani et al., 2023).

```
= reconc_MCMC(
mcmc
A,
base_forecasts,distr = "poisson",
num_samples = 30000,
seed = 42
)
<- mcmc$reconciled_samples samples_mcmc
```

We now assume that the base forecasts are Gaussian distributed, with parameters given by

- \(\mu_{Y} = 9\), \(\mu_{S_1} = 2\), and \(\mu_{S_2} = 4\);
- \(\sigma_{Y} = 2\), \(\sigma_{S_1} = 2\), and \(\sigma_{S_2} = 3\).

```
<- 2
muS1 <- 4
muS2 <- 9
muY <- c(muY, muS1, muS2)
mus
<- 2
sigmaS1 <- 2
sigmaS2 <- 3
sigmaY <- c(sigmaY, sigmaS1, sigmaS2)
sigmas
= list()
base_forecasts for (i in 1:n_tot) {
= list(mean = mus[[i]], sd = sigmas[[i]])
base_forecasts[[i]] }
```

We use the BUIS algorithm to sample from the reconciled distribution:

```
<- reconc_BUIS(
buis
A,
base_forecasts,in_type = "params",
distr = "gaussian",
num_samples = 100000,
seed = 42
)<- buis$reconciled_samples
samples_buis <- rowMeans(samples_buis) buis_means
```

If the base forecasts are Gaussian, the reconciled distribution is still Gaussian and can be computed in closed form:

```
<- diag(sigmas ^ 2) #transform into covariance matrix
Sigma <- reconc_gaussian(A,
analytic_rec base_forecasts.mu = mus,
base_forecasts.Sigma = Sigma)
<- analytic_rec$bottom_reconciled_mean
analytic_means_bottom <- A %*% analytic_means_bottom
analytic_means_upper <- rbind(analytic_means_upper,analytic_means_bottom) analytic_means
```

The base means of \(Y\), \(S_1\), and \(S_2\) are 9, 2, 4.

The reconciled means obtained analytically are 7.41, 2.71, 4.71, while the reconciled means obtained via BUIS are 7.41, 2.71, 4.71.

Corani, G., Azzimonti, D., Augusto, J.P.S.C., Zaffalon, M. (2021).
*Probabilistic Reconciliation of Hierarchical Forecast via Bayes’
Rule*. ECML PKDD 2020. Lecture Notes in Computer Science, vol 12459.
DOI

Corani, G., Azzimonti, D., Rubattu, N. (2024). *Probabilistic
reconciliation of count time series*. International Journal of
Forecasting 40 (2), 457-469. DOI

Zambon, L., Azzimonti, D. & Corani, G. (2024). *Efficient
probabilistic reconciliation of forecasts for real-valued and count time
series*. Statistics and Computing 34 (1), 21. DOI

Zambon, L., Agosto, A., Giudici, P., Corani, G. (2024).
*Properties of the reconciled distributions for Gaussian and count
forecasts*. International Journal of Forecasting (in press). DOI

Zambon, L., Azzimonti, D., Rubattu, N., Corani, G. (2024).
*Probabilistic reconciliation of mixed-type hierarchical time
series*. The 40th Conference on Uncertainty in Artificial
Intelligence, accepted.

_{Dario Azzimonti}_{(Maintainer)}dario.azzimonti@gmail.com |
_{Nicolò Rubattu}nicolo.rubattu@idsia.ch |
_{Lorenzo Zambon}lorenzo.zambon@idsia.ch |
_{Giorgio Corani}giorgio.corani@idsia.ch |

If you encounter a clear bug, please file a minimal reproducible example on GitHub.