-
Notifications
You must be signed in to change notification settings - Fork 1
/
README.Rmd
315 lines (248 loc) · 10.9 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
---
output: github_document
---
<!-- README.md is generated from README.Rmd. Please edit that file -->
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
fig.path = "man/figures/README-",
out.width = "100%"
)
```
# bayesRecon: BAyesian reCONciliation of hierarchical forecasts
<!--
<img src="./man/figures/logo.png" align="right" />
-->
<!-- badges: start -->
[](https://github.com/IDSIA/bayesRecon/actions/workflows/R-CMD-check.yaml)
[](https://CRAN.R-project.org/package=bayesRecon)
[](https://cran.r-project.org/package=bayesRecon)
[](https://lifecycle.r-lib.org/articles/stages.html#experimental)
[-yellow.svg)](https://www.gnu.org/licences/lgpl-3.0)
[](https://coveralls.io/github/IDSIA/bayesRecon)
<!-- badges: end -->
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.
## News
: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!
## Installation
You can install the **stable** version on [R
CRAN](https://cran.r-project.org/package=bayesRecon)
``` r
install.packages("bayesRecon", dependencies = TRUE)
```
You can also install the **development** version from
[Github](https://github.com/IDSIA/bayesRecon)
``` r
# install.packages("devtools")
devtools::install_github("IDSIA/bayesRecon", build_vignettes = TRUE, dependencies = TRUE)
```
## Usage
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.
```{r echo=FALSE, out.width='50%', fig.align='center'}
knitr::include_graphics('./man/figures/minimal_hierarchy.png')
```
The hierarchy is described by the *aggregation matrix* A,
which can be obtained using the function `get_reconc_matrices`.
```{r}
library(bayesRecon)
rec_mat <- get_reconc_matrices(agg_levels = c(1, 2), h = 2)
A <- rec_mat$A
print(A)
```
### Example 1: Poisson base forecasts
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$.
```{r}
lambdaS1 <- 2
lambdaS2 <- 4
lambdaY <- 9
lambdas <- c(lambdaY, lambdaS1, lambdaS2)
n_tot = length(lambdas)
base_forecasts = list()
for (i in 1:n_tot) {
base_forecasts[[i]] = list(lambda = lambdas[i])
}
```
We recommend using the BUIS algorithm (Zambon et al., 2024) to sample from the reconciled distribution.
```{r}
buis <- reconc_BUIS(
A,
base_forecasts,
in_type = "params",
distr = "poisson",
num_samples = 100000,
seed = 42
)
samples_buis <- buis$reconciled_samples
```
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$.
```{r, out.width='80%', fig.align='center'}
reconciled_forecast_S1 <- buis$bottom_reconciled_samples[1,]
range_forecats <- range(reconciled_forecast_S1)
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]),
stats::dpois(seq(range_forecats[1], range_forecats[2]), lambda =
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$.
```{r, out.width='80%', fig.align='center'}
AA <-
xyTable(buis$bottom_reconciled_samples[1, ],
buis$bottom_reconciled_samples[2, ])
plot(
AA$x ,
AA$y ,
cex = 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).
```{r}
mcmc = reconc_MCMC(
A,
base_forecasts,
distr = "poisson",
num_samples = 30000,
seed = 42
)
samples_mcmc <- mcmc$reconciled_samples
```
### Example 2: Gaussian base forecasts
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$.
```{r}
muS1 <- 2
muS2 <- 4
muY <- 9
mus <- c(muY, muS1, muS2)
sigmaS1 <- 2
sigmaS2 <- 2
sigmaY <- 3
sigmas <- c(sigmaY, sigmaS1, sigmaS2)
base_forecasts = list()
for (i in 1:n_tot) {
base_forecasts[[i]] = list(mean = mus[[i]], sd = sigmas[[i]])
}
```
We use the BUIS algorithm to sample from the reconciled distribution:
```{r}
buis <- reconc_BUIS(
A,
base_forecasts,
in_type = "params",
distr = "gaussian",
num_samples = 100000,
seed = 42
)
samples_buis <- buis$reconciled_samples
buis_means <- rowMeans(samples_buis)
```
If the base forecasts are Gaussian, the reconciled distribution is still Gaussian and can be computed in closed form:
```{r}
Sigma <- diag(sigmas ^ 2) #transform into covariance matrix
analytic_rec <- reconc_gaussian(A,
base_forecasts.mu = mus,
base_forecasts.Sigma = Sigma)
analytic_means_bottom <- analytic_rec$bottom_reconciled_mean
analytic_means_upper <- A %*% analytic_means_bottom
analytic_means <- rbind(analytic_means_upper,analytic_means_bottom)
```
The base means of $Y$, $S_1$, and $S_2$ are
`r round(mus,2)`.
The reconciled means obtained analytically are
`r round(analytic_means,2)`,
while the reconciled means obtained via BUIS are
`r round(buis_means,2)`.
## References
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](https://doi.org/10.1007/978-3-030-67664-3_13)
Corani, G., Azzimonti, D., Rubattu, N. (2024).
*Probabilistic reconciliation of count time series*.
International Journal of Forecasting 40 (2), 457-469.
[DOI](https://doi.org/10.1016/j.ijforecast.2023.04.003)
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](https://doi.org/10.1007/s11222-023-10343-y)
Zambon, L., Agosto, A., Giudici, P., Corani, G. (2024).
*Properties of the reconciled distributions for Gaussian and count forecasts*.
International Journal of Forecasting 40 (4), 1438-1448.
[DOI](https://doi.org/10.1016/j.ijforecast.2023.12.004)
Zambon, L., Azzimonti, D., Rubattu, N., Corani, G. (2024).
*Probabilistic reconciliation of mixed-type hierarchical time series*.
Proceedings of the Fortieth Conference on Uncertainty in Artificial Intelligence,
in Proceedings of Machine Learning Research 244:4078-4095. [Available here](https://proceedings.mlr.press/v244/zambon24a.html).
## Contributors
<!-- prettier-ignore-start -->
<!-- markdownlint-disable -->
<table>
<tbody>
<tr>
<td align="center" valign="top" width="14.28%">
<a href="https://sites.google.com/view/darioazzimonti/home">
<img src="https://github.com/dazzimonti.png" width="100px;" alt="Dario Azzimonti" style="border-radius:50%;border:1px solid #646464;"/><br />
<sub><b>Dario Azzimonti</b></sub></a><br />
<sub>(Maintainer)</sub><br />
<a href="mailto:[email protected]?subject=bayesRecon package!">[email protected]</a>
</td>
<td align="center" valign="top" width="14.28%">
<a href="#">
<img src="https://github.com/nicorbtt.png" width="100px;" alt="Nicolò Rubattu" style="border-radius:50%;border:1px solid #646464;"/><br />
<sub><b>Nicolò Rubattu</b></sub></a><br />
<a href="mailto:[email protected]?subject=bayesRecon package!">[email protected]</a>
</td>
<td align="center" valign="top" width="14.28%">
<a href="#">
<img src="https://github.com/LorenzoZambon.png" width="100px;" alt="Lorenzo Zambon" style="border-radius:50%;border:1px solid #646464;"/><br />
<sub><b>Lorenzo Zambon</b></sub></a><br />
<a href="mailto:[email protected]?subject=bayesRecon package!">[email protected]</a>
</td>
<td align="center" valign="top" width="14.28%">
<a href="https://sites.google.com/site/awerbhjkl678214/home">
<img src="https://github.com/gcorani.png" width="100px;" alt="Giorgio Corani" style="border-radius:50%;border:1px solid #646464;"/><br />
<sub><b>Giorgio Corani</b></sub></a><br />
<a href="mailto:[email protected]">[email protected]</a>
</td>
</tr>
</tbody>
</table>
<!-- markdownlint-restore -->
<!-- prettier-ignore-end -->
## Getting help
If you encounter a clear bug, please file a minimal reproducible example
on [GitHub](https://github.com/IDSIA/bayesRecon/issues).