|
| 1 | +# bayesRecon: BAyesian reCONciliation of hierarchical forecasts |
| 2 | + |
| 3 | +<!-- |
| 4 | +<img src="./man/figures/logo.png" align="right" /> |
| 5 | +--> |
| 6 | +<!-- badges: start --> |
| 7 | + |
| 8 | +[](https://github.com/IDSIA/bayesRecon/actions/workflows/R-CMD-check.yaml) |
| 9 | +[](https://CRAN.R-project.org/package=bayesRecon) |
| 11 | +[](https://cran.r-project.org/package=bayesRecon) |
| 12 | +[](https://lifecycle.r-lib.org/articles/stages.html#experimental) |
| 14 | +[-yellow.svg)](https://www.gnu.org/licences/lgpl-3.0) |
| 16 | +[](https://coveralls.io/github/IDSIA/bayesRecon) |
| 17 | +<!-- badges: end --> |
| 18 | + |
| 19 | +The package `bayesRecon` implements several methods for probabilistic |
| 20 | +reconciliation of hierarchical time series forecasts. |
| 21 | + |
| 22 | +The main functions are: |
| 23 | + |
| 24 | +- `reconc_gaussian`: reconciliation via conditioning of multivariate |
| 25 | + Gaussian base forecasts; this is done analytically; |
| 26 | +- `reconc_BUIS`: reconciliation via conditioning of any probabilistic |
| 27 | + forecast via importance sampling; this is the recommended option for |
| 28 | + non-Gaussian base forecasts; |
| 29 | +- `reconc_MCMC`: reconciliation via conditioning of discrete |
| 30 | + probabilistic forecasts via Markov Chain Monte Carlo; |
| 31 | +- `reconc_MixCond`: reconciliation via conditioning of mixed |
| 32 | + hierarchies, where the upper forecasts are multivariate Gaussian and |
| 33 | + the bottom forecasts are discrete distributions; |
| 34 | +- `reconc_TDcond`: reconciliation via top-down conditioning of mixed |
| 35 | + hierarchies, where the upper forecasts are multivariate Gaussian and |
| 36 | + the bottom forecasts are discrete distributions. |
| 37 | + |
| 38 | +## Getting started |
| 39 | + |
| 40 | +The starting point for `bayesRecon` functions is the vignette in the "Get Started" |
| 41 | +section of the website. Moreover you can find the documentation in the |
| 42 | +Reference section and additional vignettes in the Articles section. |
| 43 | +Finally a short example is provided below. |
| 44 | + |
| 45 | + |
| 46 | +## Installation |
| 47 | + |
| 48 | +You can install the **stable** version on [R |
| 49 | +CRAN](https://cran.r-project.org/package=bayesRecon) |
| 50 | + |
| 51 | +``` r |
| 52 | +install.packages("bayesRecon", dependencies = TRUE) |
| 53 | +``` |
| 54 | + |
| 55 | +You can also install the **development** version from |
| 56 | +[Github](https://github.com/IDSIA/bayesRecon) |
| 57 | + |
| 58 | +``` r |
| 59 | +# install.packages("devtools") |
| 60 | +devtools::install_github("IDSIA/bayesRecon", build_vignettes = TRUE, dependencies = TRUE) |
| 61 | +``` |
| 62 | + |
| 63 | +## Getting help |
| 64 | + |
| 65 | +If you encounter a clear bug, please file a minimal reproducible example |
| 66 | +on [GitHub](https://github.com/IDSIA/bayesRecon/issues). |
| 67 | + |
| 68 | + |
| 69 | +## Examples |
| 70 | + |
| 71 | +Let us consider the minimal temporal hierarchy in the figure, where the |
| 72 | +bottom variables are the two 6-monthly forecasts and the upper variable |
| 73 | +is the yearly forecast. We denote the variables for the two semesters |
| 74 | +and the year by $S_1, S_2, Y$ respectively. |
| 75 | + |
| 76 | +<img src="./man/figures/minimal_hierarchy.png" alt="Graph of a minimal hierarchy with two bottom and one upper variable." width="50%" style="display: block; margin: auto;" /> |
| 77 | + |
| 78 | +The hierarchy is described by the *aggregation matrix* A, which can be |
| 79 | +obtained using the function `get_reconc_matrices`. |
| 80 | + |
| 81 | +``` r |
| 82 | +library(bayesRecon) |
| 83 | + |
| 84 | +rec_mat <- get_reconc_matrices(agg_levels = c(1, 2), h = 2) |
| 85 | +A <- rec_mat$A |
| 86 | +print(A) |
| 87 | +#> [,1] [,2] |
| 88 | +#> [1,] 1 1 |
| 89 | +``` |
| 90 | + |
| 91 | +### Example 1: Poisson base forecasts |
| 92 | + |
| 93 | +We assume that the base forecasts are Poisson distributed, with |
| 94 | +parameters given by $\lambda_{Y} = 9$, $\lambda_{S_1} = 2$, and |
| 95 | +$\lambda_{S_2} = 4$. |
| 96 | + |
| 97 | +``` r |
| 98 | +lambdaS1 <- 2 |
| 99 | +lambdaS2 <- 4 |
| 100 | +lambdaY <- 9 |
| 101 | +lambdas <- c(lambdaY, lambdaS1, lambdaS2) |
| 102 | +n_tot = length(lambdas) |
| 103 | + |
| 104 | +base_forecasts = list() |
| 105 | +for (i in 1:n_tot) { |
| 106 | + base_forecasts[[i]] = list(lambda = lambdas[i]) |
| 107 | +} |
| 108 | +``` |
| 109 | + |
| 110 | +We recommend using the BUIS algorithm (Zambon et al., 2024) to sample |
| 111 | +from the reconciled distribution. |
| 112 | + |
| 113 | +``` r |
| 114 | +buis <- reconc_BUIS( |
| 115 | + A, |
| 116 | + base_forecasts, |
| 117 | + in_type = "params", |
| 118 | + distr = "poisson", |
| 119 | + num_samples = 100000, |
| 120 | + seed = 42 |
| 121 | +) |
| 122 | + |
| 123 | +samples_buis <- buis$reconciled_samples |
| 124 | +``` |
| 125 | + |
| 126 | +Since there is a positive incoherence in the forecasts |
| 127 | +($\lambda_Y > \lambda_{S_1}+\lambda_{S_2}$), the mean of the bottom |
| 128 | +reconciled forecast increases. We show below this behavior for $S_1$. |
| 129 | + |
| 130 | +``` r |
| 131 | +reconciled_forecast_S1 <- buis$bottom_reconciled_samples[1,] |
| 132 | +range_forecats <- range(reconciled_forecast_S1) |
| 133 | +hist( |
| 134 | + reconciled_forecast_S1, |
| 135 | + breaks = seq(range_forecats[1] - 0.5, range_forecats[2] + 0.5), |
| 136 | + freq = F, |
| 137 | + xlab = "S_1", |
| 138 | + ylab = NULL, |
| 139 | + main = "base vs reconciled" |
| 140 | +) |
| 141 | +points( |
| 142 | + seq(range_forecats[1], range_forecats[2]), |
| 143 | + stats::dpois(seq(range_forecats[1], range_forecats[2]), lambda = |
| 144 | + lambdaS1), |
| 145 | + pch = 16, |
| 146 | + col = 4, |
| 147 | + cex = 2 |
| 148 | +) |
| 149 | +``` |
| 150 | + |
| 151 | +<img src="man/figures/README-unnamed-chunk-6-1.png" alt="PMF of base versus reconciled forecasts." width="80%" style="display: block; margin: auto;" /> |
| 152 | + |
| 153 | +The blue circles represent the probability mass function of a Poisson |
| 154 | +with parameter $\lambda_{S_1}$ plotted on top of the histogram of the |
| 155 | +reconciled bottom forecasts for $S_1$. Note how the histogram is shifted |
| 156 | +to the right. |
| 157 | + |
| 158 | +Moreover, while the base bottom forecast were assumed independent, the |
| 159 | +operation of reconciliation introduced a negative correlation between |
| 160 | +$S_1$ and $S_2$. We can visualize it with the plot below which shows the |
| 161 | +empirical correlations between the reconciled samples of $S_1$ and the |
| 162 | +reconciled samples of $S_2$. |
| 163 | + |
| 164 | +``` r |
| 165 | +AA <- |
| 166 | + xyTable(buis$bottom_reconciled_samples[1, ], |
| 167 | + buis$bottom_reconciled_samples[2, ]) |
| 168 | +plot( |
| 169 | + AA$x , |
| 170 | + AA$y , |
| 171 | + cex = AA$number * 0.001 , |
| 172 | + pch = 16 , |
| 173 | + col = rgb(0, 0, 1, 0.4) , |
| 174 | + xlab = "S_1" , |
| 175 | + ylab = "S_2" , |
| 176 | + xlim = range(buis$bottom_reconciled_samples[1, ]) , |
| 177 | + ylim = range(buis$bottom_reconciled_samples[2, ]) |
| 178 | +) |
| 179 | +``` |
| 180 | + |
| 181 | +<img src="man/figures/README-unnamed-chunk-7-1.png" alt="Bubble plot of correlations between S1 and S2." width="80%" style="display: block; margin: auto;" /> |
| 182 | + |
| 183 | +We also provide a function for sampling using Markov Chain Monte Carlo |
| 184 | +(Corani et al., 2023). |
| 185 | + |
| 186 | +``` r |
| 187 | +mcmc = reconc_MCMC( |
| 188 | + A, |
| 189 | + base_forecasts, |
| 190 | + distr = "poisson", |
| 191 | + num_samples = 30000, |
| 192 | + seed = 42 |
| 193 | +) |
| 194 | + |
| 195 | +samples_mcmc <- mcmc$reconciled_samples |
| 196 | +``` |
| 197 | + |
| 198 | +### Example 2: Gaussian base forecasts |
| 199 | + |
| 200 | +We now assume that the base forecasts are Gaussian distributed, with |
| 201 | +parameters given by |
| 202 | + |
| 203 | +- $\mu_{Y} = 9$, $\mu_{S_1} = 2$, and $\mu_{S_2} = 4$; |
| 204 | +- $\sigma_{Y} = 2$, $\sigma_{S_1} = 2$, and $\sigma_{S_2} = 3$. |
| 205 | + |
| 206 | +``` r |
| 207 | +muS1 <- 2 |
| 208 | +muS2 <- 4 |
| 209 | +muY <- 9 |
| 210 | +mus <- c(muY, muS1, muS2) |
| 211 | + |
| 212 | +sigmaS1 <- 2 |
| 213 | +sigmaS2 <- 2 |
| 214 | +sigmaY <- 3 |
| 215 | +sigmas <- c(sigmaY, sigmaS1, sigmaS2) |
| 216 | + |
| 217 | +base_forecasts = list() |
| 218 | +for (i in 1:n_tot) { |
| 219 | + base_forecasts[[i]] = list(mean = mus[[i]], sd = sigmas[[i]]) |
| 220 | +} |
| 221 | +``` |
| 222 | + |
| 223 | +We use the BUIS algorithm to sample from the reconciled distribution: |
| 224 | + |
| 225 | +``` r |
| 226 | +buis <- reconc_BUIS( |
| 227 | + A, |
| 228 | + base_forecasts, |
| 229 | + in_type = "params", |
| 230 | + distr = "gaussian", |
| 231 | + num_samples = 100000, |
| 232 | + seed = 42 |
| 233 | +) |
| 234 | +samples_buis <- buis$reconciled_samples |
| 235 | +buis_means <- rowMeans(samples_buis) |
| 236 | +``` |
| 237 | + |
| 238 | +If the base forecasts are Gaussian, the reconciled distribution is still |
| 239 | +Gaussian and can be computed in closed form: |
| 240 | + |
| 241 | +``` r |
| 242 | +Sigma <- diag(sigmas ^ 2) #transform into covariance matrix |
| 243 | +analytic_rec <- reconc_gaussian(A, |
| 244 | + base_forecasts.mu = mus, |
| 245 | + base_forecasts.Sigma = Sigma) |
| 246 | +analytic_means_bottom <- analytic_rec$bottom_reconciled_mean |
| 247 | +analytic_means_upper <- A %*% analytic_means_bottom |
| 248 | +analytic_means <- rbind(analytic_means_upper,analytic_means_bottom) |
| 249 | +``` |
| 250 | + |
| 251 | +The base means of $Y$, $S_1$, and $S_2$ are 9, 2, 4. |
| 252 | + |
| 253 | +The reconciled means obtained analytically are 7.41, 2.71, 4.71, while |
| 254 | +the reconciled means obtained via BUIS are 7.41, 2.71, 4.71. |
| 255 | + |
| 256 | +## References |
| 257 | + |
| 258 | +Corani, G., Azzimonti, D., Augusto, J.P.S.C., Zaffalon, M. (2021). |
| 259 | +*Probabilistic Reconciliation of Hierarchical Forecast via Bayes’ Rule*. |
| 260 | +ECML PKDD 2020. Lecture Notes in Computer Science, vol 12459. |
| 261 | +[DOI](https://doi.org/10.1007/978-3-030-67664-3_13) |
| 262 | + |
| 263 | +Corani, G., Azzimonti, D., Rubattu, N. (2024). *Probabilistic |
| 264 | +reconciliation of count time series*. International Journal of |
| 265 | +Forecasting 40 (2), 457-469. |
| 266 | +[DOI](https://doi.org/10.1016/j.ijforecast.2023.04.003) |
| 267 | + |
| 268 | +Zambon, L., Azzimonti, D. & Corani, G. (2024). *Efficient probabilistic |
| 269 | +reconciliation of forecasts for real-valued and count time series*. |
| 270 | +Statistics and Computing 34 (1), 21. |
| 271 | +[DOI](https://doi.org/10.1007/s11222-023-10343-y) |
| 272 | + |
| 273 | +Zambon, L., Agosto, A., Giudici, P., Corani, G. (2024). *Properties of |
| 274 | +the reconciled distributions for Gaussian and count forecasts*. |
| 275 | +International Journal of Forecasting (in press). |
| 276 | +[DOI](https://doi.org/10.1016/j.ijforecast.2023.12.004) |
| 277 | + |
| 278 | +Zambon, L., Azzimonti, D., Rubattu, N., Corani, G. (2024). |
| 279 | +*Probabilistic reconciliation of mixed-type hierarchical time series*. |
| 280 | +The 40th Conference on Uncertainty in Artificial Intelligence, accepted. |
0 commit comments