diff --git a/.Rbuildignore b/.Rbuildignore
index 36d15b7..993be48 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -9,3 +9,4 @@
^codecov\.yml$
^doc$
^Meta$
+^data-raw$
diff --git a/R/bayesRecon-package.R b/R/bayesRecon-package.R
index 34a8905..a62ba47 100644
--- a/R/bayesRecon-package.R
+++ b/R/bayesRecon-package.R
@@ -21,7 +21,7 @@
#'
#' Corani, G., Azzimonti, D., Rubattu, N. (2023). *Probabilistic reconciliation of count time series*. \doi{10.1016/j.ijforecast.2023.04.003}.
#'
-#' Zambon, L., Azzimonti, D. & Corani, G. (2022). *Efficient probabilistic reconciliation of forecasts for real-valued and count time series*. \doi{10.48550/arXiv.2210.02286}.
+#' Zambon, L., Azzimonti, D. & Corani, G. (2024). *Efficient probabilistic reconciliation of forecasts for real-valued and count time series*. \doi{10.1007/s11222-023-10343-y}.
#'
#' Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). *Properties of the reconciled distributions for Gaussian and count forecasts*. \doi{10.48550/arXiv.2303.15135}.
#'
diff --git a/R/extr_mkt_events.R b/R/extr_mkt_events.R
new file mode 100644
index 0000000..71ef9bf
--- /dev/null
+++ b/R/extr_mkt_events.R
@@ -0,0 +1,26 @@
+#' Extreme market events dataset
+#'
+#' Count time series of extreme market events in five economic sectors.
+#' The data refer to the trading days between 2004/12/31 and 2018/12/19 (3508 trading days in total).
+#'
+#' The counts are computed by considering 29 companies included in the Euro Stoxx
+#' 50 index and observing if the value of the CDS spread on a given day exceeds
+#' the 90-th percentile of its distribution in the last trading year.
+#' The companies are divided in the following sectors: Financial (FIN), Information
+#' and Communication Technology (ICT), Manufacturing (MFG), Energy (ENG), and Trade (TRD).
+#'
+#' There are 6 time series:
+#' - 5 bottom time series, corresponding to the daily counts for each sector
+#' - 1 upper time series, which is the sum of all the bottom (ALL)
+#'
+#' @format
+#' A multivariate time series of class \link[stats]{ts}.
+#'
+#' @source
+#' Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). *Properties of the reconciled distributions for Gaussian and count forecasts*. \doi{10.48550/arXiv.2303.15135}.
+#'
+#' @references
+#' Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). *Properties of the reconciled distributions for Gaussian and count forecasts*. \doi{10.48550/arXiv.2303.15135}.
+#'
+#' Agosto, A. (2022). *Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion*. \doi{10.2139/ssrn.4119895}
+"extr_mkt_events"
diff --git a/R/extr_mkt_events_basefc.R b/R/extr_mkt_events_basefc.R
new file mode 100644
index 0000000..7cbfd7c
--- /dev/null
+++ b/R/extr_mkt_events_basefc.R
@@ -0,0 +1,27 @@
+#' Base forecasts for the extreme market events dataset
+#'
+#' Base forecasts for the `extr_mkt_events` dataset, computed using the model by
+#' Agosto, A. (2022). *Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion*. \doi{10.2139/ssrn.4119895}.
+#'
+#' The predictive distribution for the bottom time series is a multivariate negative
+#' binomial with a static vector of dispersion parameters and a time-varying vector
+#' of location parameters following a score-driven dynamics.
+#' The base forecasts for the upper time series are computed using a univariate version of this model.
+#' They are in-sample forecasts: for each training instant, they are computed for
+#' time t+1 by conditioning on the counts observed up to time t.
+#'
+#' @format
+#' A list `extr_mkt_events_basefc` containing
+#' \describe{
+#' \item{`extr_mkt_events_basefc$mu`}{data frame of the base forecast means, for each day}
+#' \item{`extr_mkt_events_basefc$size`}{data frame of the static base forecast size parameters}
+#' }
+#'
+#' @source
+#' Agosto, A. (2022). *Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion*. \doi{10.2139/ssrn.4119895}
+#'
+#' @references
+#' Agosto, A. (2022). *Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion*. \doi{10.2139/ssrn.4119895}
+#'
+#' Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). *Properties of the reconciled distributions for Gaussian and count forecasts*. \doi{10.48550/arXiv.2303.15135}.
+"extr_mkt_events_basefc"
diff --git a/R/reconc.R b/R/reconc.R
index 3ba2eb6..9659c95 100644
--- a/R/reconc.R
+++ b/R/reconc.R
@@ -240,7 +240,7 @@
#'print(rowMeans(samples_buis))
#'
#' @references
-#' Zambon, L., Azzimonti, D. & Corani, G. (2022). *Efficient probabilistic reconciliation of forecasts for real-valued and count time series*. \doi{10.48550/arXiv.2210.02286}.
+#' Zambon, L., Azzimonti, D. & Corani, G. (2024). *Efficient probabilistic reconciliation of forecasts for real-valued and count time series*. \doi{10.1007/s11222-023-10343-y}.
#'
#'
#' @seealso
diff --git a/README.Rmd b/README.Rmd
index c64206f..d76b0c1 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -236,8 +236,8 @@ Corani, G., Azzimonti, D., Augusto, J.P.S.C., Zaffalon, M. (2021). *Probabilisti
Corani, G., Azzimonti, D., Rubattu, N. (2023). *Probabilistic reconciliation of count
time series*. [DOI](https://doi.org/10.1016/j.ijforecast.2023.04.003)
-Zambon, L., Azzimonti, D. & Corani, G. (2022). *Efficient probabilistic reconciliation of forecasts
-for real-valued and count time series*. [DOI](https://doi.org/10.48550/arXiv.2210.02286)
+Zambon, L., Azzimonti, D. & Corani, G. (2024). *Efficient probabilistic reconciliation of forecasts
+for real-valued and count time series*. [DOI](https://doi.org/10.1007/s11222-023-10343-y)
Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). *Properties of the reconciled distributions for Gaussian and count forecasts*. [DOI](https://doi.org/10.48550/arXiv.2303.15135)
diff --git a/README.md b/README.md
index ed0d8c9..10e2d7d 100644
--- a/README.md
+++ b/README.md
@@ -256,9 +256,9 @@ Corani, G., Azzimonti, D., Rubattu, N. (2023). *Probabilistic
reconciliation of count time series*.
[DOI](https://doi.org/10.1016/j.ijforecast.2023.04.003)
-Zambon, L., Azzimonti, D. & Corani, G. (2022). *Efficient probabilistic
+Zambon, L., Azzimonti, D. & Corani, G. (2024). *Efficient probabilistic
reconciliation of forecasts for real-valued and count time series*.
-[DOI](https://doi.org/10.48550/arXiv.2210.02286)
+[DOI](https://doi.org/10.1007/s11222-023-10343-y)
Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). *Properties of
the reconciled distributions for Gaussian and count forecasts*.
diff --git a/data-raw/NB_score_driven_forecasts.xlsx b/data-raw/NB_score_driven_forecasts.xlsx
new file mode 100644
index 0000000..1aff083
Binary files /dev/null and b/data-raw/NB_score_driven_forecasts.xlsx differ
diff --git a/data-raw/extr_mkt_events.R b/data-raw/extr_mkt_events.R
new file mode 100644
index 0000000..494052f
--- /dev/null
+++ b/data-raw/extr_mkt_events.R
@@ -0,0 +1,9 @@
+library(readxl)
+
+dir_path <- dirname(rstudioapi::getSourceEditorContext()$path)
+setwd(dir_path)
+
+data <- read_excel("NB_score_driven_forecasts.xlsx", sheet = "Data")
+extr_mkt_events <- ts(data)
+
+usethis::use_data(extr_mkt_events, overwrite = TRUE)
diff --git a/data-raw/extr_mkt_events_basefc.R b/data-raw/extr_mkt_events_basefc.R
new file mode 100644
index 0000000..793e828
--- /dev/null
+++ b/data-raw/extr_mkt_events_basefc.R
@@ -0,0 +1,12 @@
+library(readxl)
+
+dir_path <- dirname(rstudioapi::getSourceEditorContext()$path)
+setwd(dir_path)
+
+predictions <- read_excel("NB_score_driven_forecasts.xlsx", sheet = "Predictions")
+alpha <- read_excel("NB_score_driven_forecasts.xlsx", sheet = "alpha", n_max = 1)
+
+extr_mkt_events_basefc <- list(mu = data.frame(predictions),
+ size = data.frame(1 / alpha))
+
+usethis::use_data(extr_mkt_events_basefc, overwrite = TRUE)
diff --git a/data/extr_mkt_events.rda b/data/extr_mkt_events.rda
new file mode 100644
index 0000000..1fe33d1
Binary files /dev/null and b/data/extr_mkt_events.rda differ
diff --git a/data/extr_mkt_events_basefc.rda b/data/extr_mkt_events_basefc.rda
new file mode 100644
index 0000000..f57f2ab
Binary files /dev/null and b/data/extr_mkt_events_basefc.rda differ
diff --git a/man/bayesRecon-package.Rd b/man/bayesRecon-package.Rd
index 8d08e3e..f2b1a79 100644
--- a/man/bayesRecon-package.Rd
+++ b/man/bayesRecon-package.Rd
@@ -39,7 +39,7 @@ Corani, G., Azzimonti, D., Augusto, J.P.S.C., Zaffalon, M. (2021). \emph{Probabi
Corani, G., Azzimonti, D., Rubattu, N. (2023). \emph{Probabilistic reconciliation of count time series}. \doi{10.1016/j.ijforecast.2023.04.003}.
-Zambon, L., Azzimonti, D. & Corani, G. (2022). \emph{Efficient probabilistic reconciliation of forecasts for real-valued and count time series}. \doi{10.48550/arXiv.2210.02286}.
+Zambon, L., Azzimonti, D. & Corani, G. (2024). \emph{Efficient probabilistic reconciliation of forecasts for real-valued and count time series}. \doi{10.1007/s11222-023-10343-y}.
Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). \emph{Properties of the reconciled distributions for Gaussian and count forecasts}. \doi{10.48550/arXiv.2303.15135}.
}
diff --git a/man/extr_mkt_events.Rd b/man/extr_mkt_events.Rd
new file mode 100644
index 0000000..14c88f1
--- /dev/null
+++ b/man/extr_mkt_events.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/extr_mkt_events.R
+\docType{data}
+\name{extr_mkt_events}
+\alias{extr_mkt_events}
+\title{Extreme market events dataset}
+\format{
+A multivariate time series of class \link[stats]{ts}.
+}
+\source{
+Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). \emph{Properties of the reconciled distributions for Gaussian and count forecasts}. \doi{10.48550/arXiv.2303.15135}.
+}
+\usage{
+extr_mkt_events
+}
+\description{
+Count time series of extreme market events in five economic sectors.
+The data refer to the trading days between 2004/12/31 and 2018/12/19 (3508 trading days in total).
+}
+\details{
+The counts are computed by considering 29 companies included in the Euro Stoxx
+50 index and observing if the value of the CDS spread on a given day exceeds
+the 90-th percentile of its distribution in the last trading year.
+The companies are divided in the following sectors: Financial (FIN), Information
+and Communication Technology (ICT), Manufacturing (MFG), Energy (ENG), and Trade (TRD).
+
+There are 6 time series:
+\itemize{
+\item 5 bottom time series, corresponding to the daily counts for each sector
+\item 1 upper time series, which is the sum of all the bottom (ALL)
+}
+}
+\references{
+Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). \emph{Properties of the reconciled distributions for Gaussian and count forecasts}. \doi{10.48550/arXiv.2303.15135}.
+
+Agosto, A. (2022). \emph{Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion}. \doi{10.2139/ssrn.4119895}
+}
+\keyword{datasets}
diff --git a/man/extr_mkt_events_basefc.Rd b/man/extr_mkt_events_basefc.Rd
new file mode 100644
index 0000000..7de3b76
--- /dev/null
+++ b/man/extr_mkt_events_basefc.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/extr_mkt_events_basefc.R
+\docType{data}
+\name{extr_mkt_events_basefc}
+\alias{extr_mkt_events_basefc}
+\title{Base forecasts for the extreme market events dataset}
+\format{
+A list \code{extr_mkt_events_basefc} containing
+\describe{
+\item{\code{extr_mkt_events_basefc$mu}}{data frame of the base forecast means, for each day}
+\item{\code{extr_mkt_events_basefc$size}}{data frame of the static base forecast size parameters}
+}
+}
+\source{
+Agosto, A. (2022). \emph{Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion}. \doi{10.2139/ssrn.4119895}
+}
+\usage{
+extr_mkt_events_basefc
+}
+\description{
+Base forecasts for the \code{extr_mkt_events} dataset, computed using the model by
+Agosto, A. (2022). \emph{Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion}. \doi{10.2139/ssrn.4119895}.
+}
+\details{
+The predictive distribution for the bottom time series is a multivariate negative
+binomial with a static vector of dispersion parameters and a time-varying vector
+of location parameters following a score-driven dynamics.
+The base forecasts for the upper time series are computed using a univariate version of this model.
+They are in-sample forecasts: for each training instant, they are computed for
+time t+1 by conditioning on the counts observed up to time t.
+}
+\references{
+Agosto, A. (2022). \emph{Multivariate Score-Driven Models for Count Time Series to Assess Financial Contagion}. \doi{10.2139/ssrn.4119895}
+
+Zambon, L., Agosto, A., Giudici, P., Corani, G. (2023). \emph{Properties of the reconciled distributions for Gaussian and count forecasts}. \doi{10.48550/arXiv.2303.15135}.
+}
+\keyword{datasets}
diff --git a/man/reconc_BUIS.Rd b/man/reconc_BUIS.Rd
index 5539108..909303f 100644
--- a/man/reconc_BUIS.Rd
+++ b/man/reconc_BUIS.Rd
@@ -148,7 +148,7 @@ print(rowMeans(samples_buis))
}
\references{
-Zambon, L., Azzimonti, D. & Corani, G. (2022). \emph{Efficient probabilistic reconciliation of forecasts for real-valued and count time series}. \doi{10.48550/arXiv.2210.02286}.
+Zambon, L., Azzimonti, D. & Corani, G. (2024). \emph{Efficient probabilistic reconciliation of forecasts for real-valued and count time series}. \doi{10.1007/s11222-023-10343-y}.
}
\seealso{
\code{\link[=reconc_gaussian]{reconc_gaussian()}}
diff --git a/vignettes/bayesRecon.Rmd b/vignettes/bayesRecon.Rmd
index d8454c1..25ff595 100644
--- a/vignettes/bayesRecon.Rmd
+++ b/vignettes/bayesRecon.Rmd
@@ -1,450 +1,450 @@
----
-title: "Probabilistic Reconciliation via Conditioning with `bayesRecon`"
-author: "Nicolò Rubattu, Giorgio Corani, Dario Azzimonti, Lorenzo Zambon"
-date: "2023-11-26"
-lang: "en"
-output: html_vignette
-bibliography: references.bib
-cite:
-- '@zambon2022efficient'
-- '@zambon2023properties'
-- '@corani2023probabilistic'
-- '@corani2021probabilistic'
-vignette: >
- %\VignetteIndexEntry{Probabilistic Reconciliation via Conditioning with `bayesRecon`}
- %\VignetteEngine{knitr::rmarkdown}
- %\VignetteEncoding{UTF-8}
----
-
-```{r setup, include = FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- eval=TRUE ### !!!! set to FALSE here to render only the text !!!!
-)
-set.seed(42)
-```
-
-```{r klippy, echo=FALSE, include=TRUE, eval=FALSE}
-klippy::klippy(position = c('top', 'right'), tooltip_message = 'Copy', tooltip_success = 'Done', color="black")
-```
-
-# Introduction
-
-This vignette shows how to perform *probabilistic reconciliation* with
-the `bayesRecon` package. We provide three examples:
-
-1. *Temporal hierarchy for a count time series*: we build a temporal hierarchy over a count time series, produce the base forecasts using `glarma` and reconcile them via Bottom-Up Importance Sampling (BUIS).
-
-2. *Temporal hierarchy for a smooth time series*: we build a temporal hierarchy over a smooth time series, compute the base forecasts using `ets` and we reconcile them in closed form using Gaussian reconciliation. The covariance matrix is diagonal.
-
-3. *Hierarchical of smooth time series*: this is an example of a cross-sectional hierarchy. We generate the base forecasts using `ets` and we reconcile them via Gaussian reconciliation.
-The covariance matrix is full and estimated via shrinkage.
-
-
-
-# Installation
-
-The package, available at this
-[CRAN page](https://cran.r-project.org/package=bayesRecon), can be installed and loaded with the usual commands:
-```{r install, eval=FALSE}
-install.packages('bayesRecon', dependencies = TRUE)
-```
-
-Load the package:
-
-```{r load}
-library(bayesRecon)
-```
-
-# Temporal hierarchy over a count time series
-We select a monthly time series of counts from the *carparts* dataset, available from
-the expsmooth package [@expsmooth_pkg].
-The data set contains time series of sales of cars part from Jan. 1998 to Mar. 2002.
-For this example we select time series #2655, which we make available as `bayesRecon::carparts_example`.
-
-This time series has a skewed distribution of values.
-```{r carpart-plot, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 1**: Carpart - monthly car part sales.", fig.dim = c(6, 3)}
-layout(mat = matrix(c(1, 2), nrow = 1, ncol = 2), widths = c(2, 1))
-plot(carparts_example, xlab = "Time", ylab = "Car part sales", main = NULL)
-hist(carparts_example, xlab = "Car part sales", main = NULL)
-```
-
-We divide the time series into train and test; the test set contains the last 12 months.
-```{r train-test}
-train <- window(carparts_example, end = c(2001, 3))
-test <- window(carparts_example, start = c(2001, 4))
-```
-
-
-
-
-
-
-We build the temporal hierarchy using the `temporal aggregation` function.
-We specify the aggregation levels using the
-`agg_levels` argument; in this case they are
-*2-Monthly*, *Quarterly*, *4-Monthly*, *Biannual*, and *Annual*.
-
-``` {r temp-agg}
-train.agg <- bayesRecon::temporal_aggregation(train, agg_levels = c(2, 3, 4, 6, 12))
-levels <- c("Annual", "Biannual", "4-Monthly", "Quarterly", "2-Monthly", "Monthly")
-names(train.agg) <- levels
-```
-
-The function returns a list of aggregated time series, ordered from the most aggregated (top of the hierarchy) to the most disagreggated (bottom of the hierarchy). We plot them below.
-``` {r temp-agg-plot, dpi=300, fig.show="hold", out.width="100%", out.heigth="100%", fig.align='center', fig.cap="**Figure 2**: The aggregated time series of the temporal hierarchy.", fig.dim=c(6,3.5)}
-par(mfrow = c(2, 3), mai = c(0.6, 0.6, 0.5, 0.5))
-for (l in levels) {
- plot(train.agg[[l]], xlab = "Time", ylab = "Car part sales", main = l)
-}
-```
-
-We compute the *base forecasts* using the package [`glarma`](https://cran.r-project.org/package=glarma),
-a package specific for forecasting count time series.
-We forecast 12 steps ahead at the monthly level, 4 steps ahead at the quarterly level, etc. by iterating over the levels of the hierarchy,
-At each level, we fit a `glarma` model
-with Poisson predictive distribution and we forecast;
-each forecast is constituted by 20000 integer samples.
-
-
-Eventually we collect the samples of the 28 predictive distributions (one at the *Annual* level, two at the *Biannual* level, etc.) in a list.
-The code below takes about 30 seconds on a standard computer.
-
-``` {r hier-fore, cache=TRUE}
-# install.packages("glarma", dependencies = TRUE)
-#library(glarma)
-
-fc.samples <- list()
-D <- 20000
-fc.count <- 1
-
-# iterating over the temporal aggregation levels
-for (l in seq_along(train.agg)) {
- f.level <- frequency(train.agg[[l]])
- print(paste("Forecasting at ", levels[l], "...", sep = ""))
- # fit an independent model for each aggregation level
- model <- glarma::glarma(
- train.agg[[l]],
- phiLags = if (f.level == 1) 1 else 1:(min(6, f.level - 1)),
- thetaLags = if (f.level == 1) NULL else f.level,
- X = cbind(intercept = rep(1, length(train.agg[[l]]))),
- offset = cbind(intercept = rep(0, length(train.agg[[l]]))),
- type = "Poi"
- )
- # forecast 1 year ahead
- h <- f.level
- tmp <- matrix(data = NA, nrow = h, ncol = D)
- for (s in (1:D)) {
- # each call to 'forecast.glarma' returns a simulation path
- tmp[, s] <- glarma::forecast(
- model,
- n.ahead = h,
- newdata = cbind(intercept = rep(1, h)),
- newoffset = rep(0, h)
- )$Y
- }
- # collect the forecasted samples
- for (i in 1:h) {
- fc.samples[[fc.count]] <- tmp[i, ]
- fc.count <- fc.count + 1
- }
-}
-
-```
-
-
-Reconciliation requires the summing matrix $\mathbf{S}$, which we obtain using the function `get_reconc_matrices`.
-It requires:
-
-* the aggregation factors of the hierarchy, which in this example are $\{2, 3, 4, 6, 12\}$;
-* the length of the forecasting horizon at the bottom level, which is 12 in this example.
-
-``` {r S}
-recon.matrices <- bayesRecon::get_reconc_matrices(agg_levels = c(2, 3, 4, 6, 12), h = 12)
-# Summing matrix
-S <- recon.matrices$S
-# A matrix
-A <- recon.matrices$A
-```
-
-To reconcile using Bottom-Up Important Sampling (BUIS) we
-we use the function `reconc_BUIS`, passing to it
-the $\mathbf{S}$ matrix, the *base forecasts*, the type of the base forecasts (`in_type`="samples") and whether the samples are discrete or integer (`distr` = "discrete").
-
-
-``` {r reconc}
-recon.res <- bayesRecon::reconc_BUIS(
- S,
- base_forecasts = fc.samples,
- in_type = "samples",
- distr = "discrete",
- seed = 42
-)
-```
-
-Here we obtain samples from the reconciled forecast distribution.
-```{r res}
-reconciled_samples <- recon.res$reconciled_samples
-dim(reconciled_samples)
-
-```
-
-We now compute the Mean Absolute Error (MAE) and the Continuous Ranked Probability Score (CRPS) for the bottom (i.e., *monthly*) time series. For computing CRPS, we use the package [`scoringRules`](https://cran.r-project.org/package=scoringRules).
-
-```{r metrics}
-# install.packages("scoringRules", dependencies = TRUE)
-library(scoringRules)
-
-ae.fc <- list()
-ae.reconc <- list()
-crps.fc <- list()
-crps.reconc <- list()
-for (h in 1:length(test)) {
- y.hat_ <- median(fc.samples[[nrow(A) + h]])
- y.reconc_ <- median(recon.res$bottom_reconciled_samples[, h])
- # Compute Absolute Errors
- ae.fc[[h]] <- abs(test[h] - y.hat_)
- ae.reconc[[h]] <- abs(test[h] - y.reconc_)
- # Compute Continuous Ranked Probability Score (CRPS)
- crps.fc[[h]] <-
- scoringRules::crps_sample(y = test[h], dat = fc.samples[[nrow(A) + h]])
- crps.reconc[[h]] <-
- scoringRules::crps_sample(y = test[h], dat = recon.res$bottom_reconciled_samples[, h])
-}
-
-mae.fc <- mean(unlist(ae.fc))
-mae.reconc <- mean(unlist(ae.reconc))
-crps.fc <- mean(unlist(crps.fc))
-crps.reconc <- mean(unlist(crps.reconc))
-metrics <- data.frame(
- row.names = c("MAE", "CRPS"),
- base.forecasts = c(mae.fc, crps.fc),
- reconciled.forecasts = c(mae.reconc, crps.reconc)
-)
-metrics
-
-```
-
-
-# Temporal hierarchy over a smooth time series
-
-In this second example, we select a smooth monthly time series (N1485) from the M3 forecasting competition [@makridakis2000m3]. The data set is available in
-the Mcomp package [@mcomp_pkg] and we make available this time series as `bayesRecon::m3_example`.
-
-
-```{r m3-plot, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 3**: M3 - N1485 time series.", fig.dim = c(6, 3)}
-plot(M3_example$train, xlab = "Time", ylab = "y", main = "N1485")
-```
-
-We build the temporal hierarchy using the `temporal_aggregation` function.
-
-Without specifying `agg_levels`, the function generates by default all the feasible aggregation: 2-Monthly, Quarterly, 4-Monthly, Biannual, and Annual.
-
-```{r m3-agg}
-train.agg <- bayesRecon::temporal_aggregation(M3_example$train)
-levels <- c("Annual", "Biannual", "4-Monthly", "Quarterly", "2-Monthly", "Monthly")
-names(train.agg) <- levels
-```
-
-
-We generate the base forecasts using `ets` from the [forecast](https://cran.r-project.org/package=forecast) package [@pkg_forecast].
-At every level we predict 18 months ahead.
-All the predictive distributions are Gaussian.
-```{r m3-fore}
-# install.packages("forecast", dependencies = TRUE)
-library(forecast)
-
-H <- length(M3_example$test)
-H
-
-fc <- list()
-level.idx <- 1
-fc.idx <- 1
-for (level in train.agg) {
- level.name <- names(train.agg)[level.idx]
- # fit an ETS model for each temporal level
- model <- ets(level)
- # generate forecasts for each level within 18 months
- h <- floor(H / (12 / frequency(level)))
- print(paste("Forecasting at ", level.name, ", h=", h, "...", sep = ""))
- level.fc <- forecast(model, h = h)
- # save mean and sd of the gaussian predictive distribution
- for (i in 1:h) {
- fc[[fc.idx]] <- c(level.fc$mean[[i]],
- (level.fc$upper[, "95%"][[i]] - level.fc$mean[[i]]) / qnorm(0.975))
- fc.idx <- fc.idx + 1
- }
- level.idx <- level.idx + 1
-}
-```
-
-Using the function `get_reconc_matrices`, we get matrix $\mathbf{S}$.
-
-```{r m3-rmat, dpi=300, out.width = '70%', fig.align='center', fig.cap="**Figure 4**: M3 - The aggregating matrix A (red=1, yellow=0).", fig.dim = c(8, 8)}
-rmat <- get_reconc_matrices(agg_levels = c(2, 3, 4, 6, 12), h = 18)
-
-par(mai = c(1,1,0.5,0.5))
-image(1:ncol(rmat$A), 1:nrow(rmat$A),
- t(apply(t(rmat$A),1,rev)),
- xaxt='n', yaxt='n', ylab = "", xlab = levels[6])
-axis(1, at=1:ncol(rmat$A), label=1:ncol(rmat$A), las=1)
-axis(2, at=c(23,22,19,15,9), label=levels[1:5], las=2)
-```
-
-The function `reconc_gaussian` implements the exact Gaussian reconciliation.
-We also run `reconc_BUIS`, to check the consistency
-between the two approaches.
-
-
-```{r m3-reco}
-recon.gauss <- bayesRecon::reconc_gaussian(
- S = rmat$S,
- base_forecasts.mu = sapply(fc, "[[", 1),
- base_forecasts.Sigma = diag(sapply(fc, "[[", 2)) ^ 2
-)
-
-reconc.buis <- bayesRecon::reconc_BUIS(
- S = rmat$S,
- base_forecasts = fc,
- in_type = "params",
- distr = "gaussian",
- num_samples = 20000,
- seed = 42
-)
-
-# check that the algorithms return consistent results
-round(rbind(
- c(rmat$S %*% recon.gauss$bottom_reconciled_mean),
- rowMeans(reconc.buis$reconciled_samples)
-))
-```
-
-We now compare *base forecasts* and *reconciled forecasts*:
-
-```{r m3-plotfore, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 5**: M3 - Base and reconciled forecasts. The black line shows the actual data (dashed in the test). The orange line is the mean of the base forecasts, the blu line is the reconciled mean. We also show the 95% prediction intervals.", fig.dim = c(6, 4)}
-yhat.mu <- tail(sapply(fc, "[[", 1), 18)
-yhat.sigma <- tail(sapply(fc, "[[", 2), 18)
-yhat.hi95 <- qnorm(0.975, mean = yhat.mu, sd = yhat.sigma)
-yhat.lo95 <- qnorm(0.025, mean = yhat.mu, sd = yhat.sigma)
-yreconc.mu <- rowMeans(reconc.buis$bottom_reconciled_samples)
-yreconc.hi95 <- apply(reconc.buis$bottom_reconciled_samples, 1,
- function(x) quantile(x, 0.975))
-yreconc.lo95 <- apply(reconc.buis$bottom_reconciled_samples, 1,
- function(x) quantile(x, 0.025))
-
-ylim_ <- c(min(M3_example$train, M3_example$test, yhat.lo95, yreconc.lo95) - 1,
- max(M3_example$train, M3_example$test, yhat.hi95, yreconc.hi95) + 1)
-
-plot(M3_example$train, xlim = c(1990, 1995.6), ylim = ylim_,
- ylab = "y", main = "N1485 Forecasts")
-lines(M3_example$test, lty = "dashed")
-lines(ts(yhat.mu, start = start(M3_example$test), frequency = 12),
- col = "coral", lwd = 2)
-lines(ts(yreconc.mu, start = start(M3_example$test), frequency = 12),
- col = "blue2", lwd = 2)
-xtest <- time(M3_example$test)
-polygon(c(xtest, rev(xtest)), c(yhat.mu, rev(yhat.hi95)),
- col = "#FF7F5066", border = "#FF7F5066")
-polygon(c(xtest, rev(xtest)), c(yhat.mu, rev(yhat.lo95)),
- col = "#FF7F5066", border = "#FF7F5066")
-polygon(c(xtest, rev(xtest)), c(yreconc.mu, rev(yreconc.hi95)),
- col = "#0000EE4D", border = "#0000EE4D")
-polygon(c(xtest, rev(xtest)), c(yreconc.mu, rev(yreconc.lo95)),
- col = "#0000EE4D", border = "#0000EE4D")
-```
-
-
-
-# Gaussian reconciliation of a cross-sectional hierarchy
-
-In this example, we consider the hierarchical time series *infantgts*, which is available
-from the `hts` package [@hts_pkg]
-We make it available also in our package as `bayesRecon::infantMortality`.
-
-It contains counts of infant mortality (deaths) in Australia, disaggregated by state and sex (male and female).
-
-
-
-We forecast one year ahead using `auto.arima` from the [`forecast` ](https://cran.r-project.org/package=forecast) package.
-We collect the residuals, which we will later use to compute the covariance matrix.
-
-```{r infants-forecasts}
-# install.packages("forecast", dependencies = TRUE)
-library(forecast)
-
-fc <- list()
-residuals <- matrix(NA,
- nrow = length(infantMortality$total),
- ncol = length(infantMortality))
-fc.idx <- 1
-for (s in infantMortality) {
- s.name <- names(infantMortality)[fc.idx]
- print(paste("Forecasting at ", s.name, "...", sep = ""))
- # fit an auto.arima model and forecast with h=1
- model <- auto.arima(s)
- s.fc <- forecast(model, h = 1)
- # save mean and sd of the gaussian predictive distribution
- fc[[s.name]] <- c(s.fc$mean,
- (s.fc$upper[, "95%"][[1]] - s.fc$mean) / qnorm(0.975))
- residuals[, fc.idx] <- s.fc$residuals
- fc.idx <- fc.idx + 1
-}
-```
-
-Now we build the $\mathbf{S}$ matrix.
-
-```{r infants-s, dpi=300, out.width = '70%', fig.align='center', fig.cap="**Figure 6**: Infants mortality - The aggregating matrix A (red=1, yellow=0).", fig.dim = c(8, 8)}
-# we have 16 bottom time series, and 11 upper time series
-A <- matrix(data = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
- 1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,
- 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,
- 1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,
- 0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1), byrow=TRUE, ncol = 16)
-S <- rbind(A, diag(16))
-
-# plot of A
-par(mai = c(1.5,1,0.5,0.5))
-image(1:ncol(A), 1:nrow(A),
- t(apply(t(A),1,rev)),
- xaxt='n', yaxt='n', ann=FALSE)
-axis(1, at=1:ncol(A), label=names(infantMortality)[12:27], las=2)
-axis(2, at=c(1:11), label=rev(names(infantMortality)[1:11]), las=2)
-```
-
-We use `bayesRecon::schaferStrimmer_cov` to estimate the covariance matrix of the residuals with shrinkage [@schafer2005shrinkage].
-
-```{r infants reconc}
-# means
-mu <- sapply(fc, "[[", 1)
-# Shrinkage covariance
-shrink.res <- bayesRecon::schaferStrimmer_cov(residuals)
-print(paste("The estimated shrinkage intensity is", round(shrink.res$lambda_star, 3)))
-Sigma <- shrink.res$shrink_cov
-```
-
-We now perform Gaussian reconciliation:
-
-```{r infants-recon}
-recon.gauss <- bayesRecon::reconc_gaussian(S,
- base_forecasts.mu = mu,
- base_forecasts.Sigma = Sigma)
-
-bottom_mu_reconc <- recon.gauss$bottom_reconciled_mean
-bottom_Sigma_reconc <- recon.gauss$bottom_reconciled_covariance
-
-# Obtain reconciled mu and Sigma for the upper variable
-upper_mu_reconc <- A %*% bottom_mu_reconc
-upper_Sigma_reconc <- A %*% bottom_Sigma_reconc %*% t(A)
-
-upper_mu_reconc
-```
-
-# References
-
+---
+title: "Probabilistic Reconciliation via Conditioning with `bayesRecon`"
+author: "Nicolò Rubattu, Giorgio Corani, Dario Azzimonti, Lorenzo Zambon"
+date: "2023-11-26"
+lang: "en"
+output: html_vignette
+bibliography: references.bib
+cite:
+- '@zambon2022efficient'
+- '@zambon2023properties'
+- '@corani2023probabilistic'
+- '@corani2021probabilistic'
+vignette: >
+ %\VignetteIndexEntry{Probabilistic Reconciliation via Conditioning with `bayesRecon`}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>",
+ eval=TRUE ### !!!! set to FALSE here to render only the text !!!!
+)
+set.seed(42)
+```
+
+```{r klippy, echo=FALSE, include=TRUE, eval=FALSE}
+klippy::klippy(position = c('top', 'right'), tooltip_message = 'Copy', tooltip_success = 'Done', color="black")
+```
+
+# Introduction
+
+This vignette shows how to perform *probabilistic reconciliation* with
+the `bayesRecon` package. We provide three examples:
+
+1. *Temporal hierarchy for a count time series*: we build a temporal hierarchy over a count time series, produce the base forecasts using `glarma` and reconcile them via Bottom-Up Importance Sampling (BUIS).
+
+2. *Temporal hierarchy for a smooth time series*: we build a temporal hierarchy over a smooth time series, compute the base forecasts using `ets` and we reconcile them in closed form using Gaussian reconciliation. The covariance matrix is diagonal.
+
+3. *Hierarchical of smooth time series*: this is an example of a cross-sectional hierarchy. We generate the base forecasts using `ets` and we reconcile them via Gaussian reconciliation.
+The covariance matrix is full and estimated via shrinkage.
+
+
+
+# Installation
+
+The package, available at this
+[CRAN page](https://cran.r-project.org/package=bayesRecon), can be installed and loaded with the usual commands:
+```{r install, eval=FALSE}
+install.packages('bayesRecon', dependencies = TRUE)
+```
+
+Load the package:
+
+```{r load}
+library(bayesRecon)
+```
+
+# Temporal hierarchy over a count time series
+We select a monthly time series of counts from the *carparts* dataset, available from
+the expsmooth package [@expsmooth_pkg].
+The data set contains time series of sales of cars part from Jan. 1998 to Mar. 2002.
+For this example we select time series #2655, which we make available as `bayesRecon::carparts_example`.
+
+This time series has a skewed distribution of values.
+```{r carpart-plot, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 1**: Carpart - monthly car part sales.", fig.dim = c(6, 3)}
+layout(mat = matrix(c(1, 2), nrow = 1, ncol = 2), widths = c(2, 1))
+plot(carparts_example, xlab = "Time", ylab = "Car part sales", main = NULL)
+hist(carparts_example, xlab = "Car part sales", main = NULL)
+```
+
+We divide the time series into train and test; the test set contains the last 12 months.
+```{r train-test}
+train <- window(carparts_example, end = c(2001, 3))
+test <- window(carparts_example, start = c(2001, 4))
+```
+
+
+
+
+
+
+We build the temporal hierarchy using the `temporal aggregation` function.
+We specify the aggregation levels using the
+`agg_levels` argument; in this case they are
+*2-Monthly*, *Quarterly*, *4-Monthly*, *Biannual*, and *Annual*.
+
+``` {r temp-agg}
+train.agg <- bayesRecon::temporal_aggregation(train, agg_levels = c(2, 3, 4, 6, 12))
+levels <- c("Annual", "Biannual", "4-Monthly", "Quarterly", "2-Monthly", "Monthly")
+names(train.agg) <- levels
+```
+
+The function returns a list of aggregated time series, ordered from the most aggregated (top of the hierarchy) to the most disagreggated (bottom of the hierarchy). We plot them below.
+``` {r temp-agg-plot, dpi=300, fig.show="hold", out.width="100%", out.heigth="100%", fig.align='center', fig.cap="**Figure 2**: The aggregated time series of the temporal hierarchy.", fig.dim=c(6,3.5)}
+par(mfrow = c(2, 3), mai = c(0.6, 0.6, 0.5, 0.5))
+for (l in levels) {
+ plot(train.agg[[l]], xlab = "Time", ylab = "Car part sales", main = l)
+}
+```
+
+We compute the *base forecasts* using the package [`glarma`](https://cran.r-project.org/package=glarma),
+a package specific for forecasting count time series.
+We forecast 12 steps ahead at the monthly level, 4 steps ahead at the quarterly level, etc. by iterating over the levels of the hierarchy,
+At each level, we fit a `glarma` model
+with Poisson predictive distribution and we forecast;
+each forecast is constituted by 20000 integer samples.
+
+
+Eventually we collect the samples of the 28 predictive distributions (one at the *Annual* level, two at the *Biannual* level, etc.) in a list.
+The code below takes about 30 seconds on a standard computer.
+
+``` {r hier-fore, cache=TRUE}
+# install.packages("glarma", dependencies = TRUE)
+#library(glarma)
+
+fc.samples <- list()
+D <- 20000
+fc.count <- 1
+
+# iterating over the temporal aggregation levels
+for (l in seq_along(train.agg)) {
+ f.level <- frequency(train.agg[[l]])
+ print(paste("Forecasting at ", levels[l], "...", sep = ""))
+ # fit an independent model for each aggregation level
+ model <- glarma::glarma(
+ train.agg[[l]],
+ phiLags = if (f.level == 1) 1 else 1:(min(6, f.level - 1)),
+ thetaLags = if (f.level == 1) NULL else f.level,
+ X = cbind(intercept = rep(1, length(train.agg[[l]]))),
+ offset = cbind(intercept = rep(0, length(train.agg[[l]]))),
+ type = "Poi"
+ )
+ # forecast 1 year ahead
+ h <- f.level
+ tmp <- matrix(data = NA, nrow = h, ncol = D)
+ for (s in (1:D)) {
+ # each call to 'forecast.glarma' returns a simulation path
+ tmp[, s] <- glarma::forecast(
+ model,
+ n.ahead = h,
+ newdata = cbind(intercept = rep(1, h)),
+ newoffset = rep(0, h)
+ )$Y
+ }
+ # collect the forecasted samples
+ for (i in 1:h) {
+ fc.samples[[fc.count]] <- tmp[i, ]
+ fc.count <- fc.count + 1
+ }
+}
+
+```
+
+
+Reconciliation requires the summing matrix $\mathbf{S}$, which we obtain using the function `get_reconc_matrices`.
+It requires:
+
+* the aggregation factors of the hierarchy, which in this example are $\{2, 3, 4, 6, 12\}$;
+* the length of the forecasting horizon at the bottom level, which is 12 in this example.
+
+``` {r S}
+recon.matrices <- bayesRecon::get_reconc_matrices(agg_levels = c(2, 3, 4, 6, 12), h = 12)
+# Summing matrix
+S <- recon.matrices$S
+# A matrix
+A <- recon.matrices$A
+```
+
+To reconcile using Bottom-Up Important Sampling (BUIS) we
+we use the function `reconc_BUIS`, passing to it
+the $\mathbf{S}$ matrix, the *base forecasts*, the type of the base forecasts (`in_type`="samples") and whether the samples are discrete or integer (`distr` = "discrete").
+
+
+``` {r reconc}
+recon.res <- bayesRecon::reconc_BUIS(
+ S,
+ base_forecasts = fc.samples,
+ in_type = "samples",
+ distr = "discrete",
+ seed = 42
+)
+```
+
+Here we obtain samples from the reconciled forecast distribution.
+```{r res}
+reconciled_samples <- recon.res$reconciled_samples
+dim(reconciled_samples)
+
+```
+
+We now compute the Mean Absolute Error (MAE) and the Continuous Ranked Probability Score (CRPS) for the bottom (i.e., *monthly*) time series. For computing CRPS, we use the package [`scoringRules`](https://cran.r-project.org/package=scoringRules).
+
+```{r metrics}
+# install.packages("scoringRules", dependencies = TRUE)
+library(scoringRules)
+
+ae.fc <- list()
+ae.reconc <- list()
+crps.fc <- list()
+crps.reconc <- list()
+for (h in 1:length(test)) {
+ y.hat_ <- median(fc.samples[[nrow(A) + h]])
+ y.reconc_ <- median(recon.res$bottom_reconciled_samples[, h])
+ # Compute Absolute Errors
+ ae.fc[[h]] <- abs(test[h] - y.hat_)
+ ae.reconc[[h]] <- abs(test[h] - y.reconc_)
+ # Compute Continuous Ranked Probability Score (CRPS)
+ crps.fc[[h]] <-
+ scoringRules::crps_sample(y = test[h], dat = fc.samples[[nrow(A) + h]])
+ crps.reconc[[h]] <-
+ scoringRules::crps_sample(y = test[h], dat = recon.res$bottom_reconciled_samples[, h])
+}
+
+mae.fc <- mean(unlist(ae.fc))
+mae.reconc <- mean(unlist(ae.reconc))
+crps.fc <- mean(unlist(crps.fc))
+crps.reconc <- mean(unlist(crps.reconc))
+metrics <- data.frame(
+ row.names = c("MAE", "CRPS"),
+ base.forecasts = c(mae.fc, crps.fc),
+ reconciled.forecasts = c(mae.reconc, crps.reconc)
+)
+metrics
+
+```
+
+
+# Temporal hierarchy over a smooth time series
+
+In this second example, we select a smooth monthly time series (N1485) from the M3 forecasting competition [@makridakis2000m3]. The data set is available in
+the Mcomp package [@mcomp_pkg] and we make available this time series as `bayesRecon::m3_example`.
+
+
+```{r m3-plot, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 3**: M3 - N1485 time series.", fig.dim = c(6, 3)}
+plot(M3_example$train, xlab = "Time", ylab = "y", main = "N1485")
+```
+
+We build the temporal hierarchy using the `temporal_aggregation` function.
+
+Without specifying `agg_levels`, the function generates by default all the feasible aggregation: 2-Monthly, Quarterly, 4-Monthly, Biannual, and Annual.
+
+```{r m3-agg}
+train.agg <- bayesRecon::temporal_aggregation(M3_example$train)
+levels <- c("Annual", "Biannual", "4-Monthly", "Quarterly", "2-Monthly", "Monthly")
+names(train.agg) <- levels
+```
+
+
+We generate the base forecasts using `ets` from the [forecast](https://cran.r-project.org/package=forecast) package [@pkg_forecast].
+At every level we predict 18 months ahead.
+All the predictive distributions are Gaussian.
+```{r m3-fore}
+# install.packages("forecast", dependencies = TRUE)
+library(forecast)
+
+H <- length(M3_example$test)
+H
+
+fc <- list()
+level.idx <- 1
+fc.idx <- 1
+for (level in train.agg) {
+ level.name <- names(train.agg)[level.idx]
+ # fit an ETS model for each temporal level
+ model <- ets(level)
+ # generate forecasts for each level within 18 months
+ h <- floor(H / (12 / frequency(level)))
+ print(paste("Forecasting at ", level.name, ", h=", h, "...", sep = ""))
+ level.fc <- forecast(model, h = h)
+ # save mean and sd of the gaussian predictive distribution
+ for (i in 1:h) {
+ fc[[fc.idx]] <- c(level.fc$mean[[i]],
+ (level.fc$upper[, "95%"][[i]] - level.fc$mean[[i]]) / qnorm(0.975))
+ fc.idx <- fc.idx + 1
+ }
+ level.idx <- level.idx + 1
+}
+```
+
+Using the function `get_reconc_matrices`, we get matrix $\mathbf{S}$.
+
+```{r m3-rmat, dpi=300, out.width = '70%', fig.align='center', fig.cap="**Figure 4**: M3 - The aggregating matrix A (red=1, yellow=0).", fig.dim = c(8, 8)}
+rmat <- get_reconc_matrices(agg_levels = c(2, 3, 4, 6, 12), h = 18)
+
+par(mai = c(1,1,0.5,0.5))
+image(1:ncol(rmat$A), 1:nrow(rmat$A),
+ t(apply(t(rmat$A),1,rev)),
+ xaxt='n', yaxt='n', ylab = "", xlab = levels[6])
+axis(1, at=1:ncol(rmat$A), label=1:ncol(rmat$A), las=1)
+axis(2, at=c(23,22,19,15,9), label=levels[1:5], las=2)
+```
+
+The function `reconc_gaussian` implements the exact Gaussian reconciliation.
+We also run `reconc_BUIS`, to check the consistency
+between the two approaches.
+
+
+```{r m3-reco}
+recon.gauss <- bayesRecon::reconc_gaussian(
+ S = rmat$S,
+ base_forecasts.mu = sapply(fc, "[[", 1),
+ base_forecasts.Sigma = diag(sapply(fc, "[[", 2)) ^ 2
+)
+
+reconc.buis <- bayesRecon::reconc_BUIS(
+ S = rmat$S,
+ base_forecasts = fc,
+ in_type = "params",
+ distr = "gaussian",
+ num_samples = 20000,
+ seed = 42
+)
+
+# check that the algorithms return consistent results
+round(rbind(
+ c(rmat$S %*% recon.gauss$bottom_reconciled_mean),
+ rowMeans(reconc.buis$reconciled_samples)
+))
+```
+
+We now compare *base forecasts* and *reconciled forecasts*:
+
+```{r m3-plotfore, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 5**: M3 - Base and reconciled forecasts. The black line shows the actual data (dashed in the test). The orange line is the mean of the base forecasts, the blu line is the reconciled mean. We also show the 95% prediction intervals.", fig.dim = c(6, 4)}
+yhat.mu <- tail(sapply(fc, "[[", 1), 18)
+yhat.sigma <- tail(sapply(fc, "[[", 2), 18)
+yhat.hi95 <- qnorm(0.975, mean = yhat.mu, sd = yhat.sigma)
+yhat.lo95 <- qnorm(0.025, mean = yhat.mu, sd = yhat.sigma)
+yreconc.mu <- rowMeans(reconc.buis$bottom_reconciled_samples)
+yreconc.hi95 <- apply(reconc.buis$bottom_reconciled_samples, 1,
+ function(x) quantile(x, 0.975))
+yreconc.lo95 <- apply(reconc.buis$bottom_reconciled_samples, 1,
+ function(x) quantile(x, 0.025))
+
+ylim_ <- c(min(M3_example$train, M3_example$test, yhat.lo95, yreconc.lo95) - 1,
+ max(M3_example$train, M3_example$test, yhat.hi95, yreconc.hi95) + 1)
+
+plot(M3_example$train, xlim = c(1990, 1995.6), ylim = ylim_,
+ ylab = "y", main = "N1485 Forecasts")
+lines(M3_example$test, lty = "dashed")
+lines(ts(yhat.mu, start = start(M3_example$test), frequency = 12),
+ col = "coral", lwd = 2)
+lines(ts(yreconc.mu, start = start(M3_example$test), frequency = 12),
+ col = "blue2", lwd = 2)
+xtest <- time(M3_example$test)
+polygon(c(xtest, rev(xtest)), c(yhat.mu, rev(yhat.hi95)),
+ col = "#FF7F5066", border = "#FF7F5066")
+polygon(c(xtest, rev(xtest)), c(yhat.mu, rev(yhat.lo95)),
+ col = "#FF7F5066", border = "#FF7F5066")
+polygon(c(xtest, rev(xtest)), c(yreconc.mu, rev(yreconc.hi95)),
+ col = "#0000EE4D", border = "#0000EE4D")
+polygon(c(xtest, rev(xtest)), c(yreconc.mu, rev(yreconc.lo95)),
+ col = "#0000EE4D", border = "#0000EE4D")
+```
+
+
+
+# Gaussian reconciliation of a cross-sectional hierarchy
+
+In this example, we consider the hierarchical time series *infantgts*, which is available
+from the `hts` package [@hts_pkg]
+We make it available also in our package as `bayesRecon::infantMortality`.
+
+It contains counts of infant mortality (deaths) in Australia, disaggregated by state and sex (male and female).
+
+
+
+We forecast one year ahead using `auto.arima` from the [`forecast` ](https://cran.r-project.org/package=forecast) package.
+We collect the residuals, which we will later use to compute the covariance matrix.
+
+```{r infants-forecasts}
+# install.packages("forecast", dependencies = TRUE)
+library(forecast)
+
+fc <- list()
+residuals <- matrix(NA,
+ nrow = length(infantMortality$total),
+ ncol = length(infantMortality))
+fc.idx <- 1
+for (s in infantMortality) {
+ s.name <- names(infantMortality)[fc.idx]
+ print(paste("Forecasting at ", s.name, "...", sep = ""))
+ # fit an auto.arima model and forecast with h=1
+ model <- auto.arima(s)
+ s.fc <- forecast(model, h = 1)
+ # save mean and sd of the gaussian predictive distribution
+ fc[[s.name]] <- c(s.fc$mean,
+ (s.fc$upper[, "95%"][[1]] - s.fc$mean) / qnorm(0.975))
+ residuals[, fc.idx] <- s.fc$residuals
+ fc.idx <- fc.idx + 1
+}
+```
+
+Now we build the $\mathbf{S}$ matrix.
+
+```{r infants-s, dpi=300, out.width = '70%', fig.align='center', fig.cap="**Figure 6**: Infants mortality - The aggregating matrix A (red=1, yellow=0).", fig.dim = c(8, 8)}
+# we have 16 bottom time series, and 11 upper time series
+A <- matrix(data = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
+ 1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,
+ 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,
+ 1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,
+ 0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1), byrow=TRUE, ncol = 16)
+S <- rbind(A, diag(16))
+
+# plot of A
+par(mai = c(1.5,1,0.5,0.5))
+image(1:ncol(A), 1:nrow(A),
+ t(apply(t(A),1,rev)),
+ xaxt='n', yaxt='n', ann=FALSE)
+axis(1, at=1:ncol(A), label=names(infantMortality)[12:27], las=2)
+axis(2, at=c(1:11), label=rev(names(infantMortality)[1:11]), las=2)
+```
+
+We use `bayesRecon::schaferStrimmer_cov` to estimate the covariance matrix of the residuals with shrinkage [@schafer2005shrinkage].
+
+```{r infants reconc}
+# means
+mu <- sapply(fc, "[[", 1)
+# Shrinkage covariance
+shrink.res <- bayesRecon::schaferStrimmer_cov(residuals)
+print(paste("The estimated shrinkage intensity is", round(shrink.res$lambda_star, 3)))
+Sigma <- shrink.res$shrink_cov
+```
+
+We now perform Gaussian reconciliation:
+
+```{r infants-recon}
+recon.gauss <- bayesRecon::reconc_gaussian(S,
+ base_forecasts.mu = mu,
+ base_forecasts.Sigma = Sigma)
+
+bottom_mu_reconc <- recon.gauss$bottom_reconciled_mean
+bottom_Sigma_reconc <- recon.gauss$bottom_reconciled_covariance
+
+# Obtain reconciled mu and Sigma for the upper variable
+upper_mu_reconc <- A %*% bottom_mu_reconc
+upper_Sigma_reconc <- A %*% bottom_Sigma_reconc %*% t(A)
+
+upper_mu_reconc
+```
+
+# References
+
\ No newline at end of file
diff --git a/vignettes/img/finTS_hier.jpg b/vignettes/img/finTS_hier.jpg
new file mode 100644
index 0000000..a221502
Binary files /dev/null and b/vignettes/img/finTS_hier.jpg differ
diff --git a/vignettes/reconciliation_properties.Rmd b/vignettes/reconciliation_properties.Rmd
new file mode 100644
index 0000000..3dc2025
--- /dev/null
+++ b/vignettes/reconciliation_properties.Rmd
@@ -0,0 +1,277 @@
+---
+title: "Properties of the reconciled distribution via conditioning"
+author: "Lorenzo Zambon, Giorgio Corani"
+# date: "2023-12-05"
+lang: "en"
+output: rmarkdown::html_vignette
+bibliography: references.bib
+vignette: >
+ %\VignetteIndexEntry{Properties of the reconciled distribution via conditioning}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>"
+)
+```
+
+
+# Introduction
+
+This vignette reproduces the results of the paper *Properties of the reconciled distributions for Gaussian and count forecasts* [@zambon2023properties],
+accepted by the International Journal of Forecasting for publication.
+We replicate here the main experiment of the paper, where we reconcile base forecasts constituted by negative binomial distributions.
+
+We use the R package `BayesRecon`.
+
+```{r setup}
+library(bayesRecon)
+```
+
+# Data and base forecasts
+
+We release a new data set, containing time series of counts of extreme market events
+in five economic sectors in the period 2005-2018 (3508 trading days).
+The counts are computed by considering 29 companies included in the Euro Stoxx 50 index
+and observing if the value of the CDS spread on a given day exceeds the 90-th percentile of its distribution in the last trading year.
+The companies are divided in the following sectors: Financial (FIN), Information and Communication Technology (ICT),
+Manufacturing (MFG), Energy (ENG), and Trade (TRD).
+
+The hierarchy is composed by 5 bottom time series, the daily number of extreme market
+events in each sector, and 1 upper time series (the sum of the different sectors).
+Data are stored in `extr_mkt_events`.
+
+```{r out.width = '50%', echo = FALSE}
+knitr::include_graphics("img/finTS_hier.jpg")
+```
+
+We also make available as `extr_mkt_events_basefc` the base forecasts to be reconciled.
+As discussed in the paper, they are produced
+using the model by [@agosto2022multivariate];
+the predictive distributions
+are negative binomial.
+
+```{r}
+# Hierarchy composed by 6 time series: 5 bottom and 1 upper
+n_b <- 5
+n_u <- 1
+n <- n_b + n_u
+
+A <- matrix(1, ncol = n_b, nrow = n_u) # aggregating matrix
+S <- rbind(A, diag(n_b)) # summing matrix
+
+# Actual values:
+actuals <- data.frame(extr_mkt_events) # convert to data frame
+# Base forecasts:
+base_fc <- extr_mkt_events_basefc
+
+N <- nrow(actuals) # number of days (3508)
+
+# # If you want to run only N reconciliations (instead of 3508):
+# N <- 200
+# actuals <- actuals[1:N,]
+# base_fc$mu <- base_fc$mu[1:N,]
+```
+
+# Reconciliation via conditioning
+
+We reconcile the base forecasts via conditioning, using importance sampling.
+We use the `reconc_BUIS` function, which implements the BUIS algorithm [@zambon2022efficient];
+since there is only one upper time series in the hierarchy, the BUIS algorithm is
+equivalent to importance sampling.
+We perform 3508 reconciliations, one for each day, drawing each time 10,000 samples from the reconciled distribution.
+We use 10,000 samples instead of 100,000 (as in the paper) to speed up the computation.
+
+For each day, we save the empirical mean, median, and quantiles of the reconciled distribution.
+
+```{r}
+# We need to save the mean and median of the reconciled distribution
+# in order to compute the squared error and the absolute error:
+rec_means <- matrix(NA, ncol = n, nrow = N)
+rec_medians <- matrix(NA, ncol = n, nrow = N)
+
+# We need to save the lower and upper quantiles of the reconciled distribution
+# in order to compute the interval score:
+rec_L <- matrix(NA, ncol = n, nrow = N)
+rec_U <- matrix(NA, ncol = n, nrow = N)
+int_cov = 0.9 # use 90% interval
+q1 <- (1 - int_cov) / 2
+q2 <- (1 + int_cov) / 2
+
+# Set the number of samples to draw from the reconciled distribution:
+N_samples <- 1e4
+
+start <- Sys.time()
+for (j in 1:N) {
+ # Base forecasts:
+ base_fc_j <- c()
+ for (i in 1:n) {
+ base_fc_j[[i]] <- c(base_fc$mu[[j,i]], base_fc$size[[i]])
+ }
+
+ # Reconcile via importance sampling:
+ buis <- reconc_BUIS(S, base_fc_j, "params", "nbinom",
+ num_samples = N_samples, seed = 42)
+ samples_y <- buis$reconciled_samples
+
+ # Save mean, median, and lower and upper quantiles:
+ rec_means[j,] <- rowMeans(samples_y)
+ rec_medians[j,] <- apply(samples_y, 1, median)
+ rec_L[j,] <- apply(samples_y, 1, quantile, q1)
+ rec_U[j,] <- apply(samples_y, 1, quantile, q2)
+}
+stop <- Sys.time()
+cat("Computational time for ", N, " reconciliations: ",
+ round(difftime(stop, start, units = "secs"), 2), "s")
+```
+
+We compute the median and the quantiles of the negative binomial base forecasts.
+
+```{r}
+base_means <- base_fc$mu
+base_medians <- matrix(NA, ncol = n, nrow = N)
+base_L <- matrix(NA, ncol = n, nrow = N)
+base_U <- matrix(NA, ncol = n, nrow = N)
+
+for (i in 1:n) {
+ base_medians[,i] <- sapply(base_means[,i],
+ function(mu) qnbinom(p=0.5, size=base_fc$size[[i]], mu=mu))
+ base_L[,i] <- sapply(base_means[,i],
+ function(mu) qnbinom(p=q1, size=base_fc$size[[i]], mu=mu))
+ base_U[,i] <- sapply(base_means[,i],
+ function(mu) qnbinom(p=q2, size=base_fc$size[[i]], mu=mu))
+}
+```
+
+For each day and for each time series, we compute the absolute error, the squared error,
+and the interval score for the base and reconciled forecasts.
+
+```{r}
+# Compute the squared errors
+SE_base <- (base_means - actuals)^2
+SE_rec <- (rec_means - actuals)^2
+
+# Compute the absolute errors
+AE_base <- abs(base_medians - actuals)
+AE_rec <- abs(rec_medians - actuals)
+
+# Define the function for the interval score
+int_score <- function(l, u, actual, int_cov = 0.9) {
+ is <- (u - l) +
+ 2 / (1-int_cov) * (actual - u) * (actual>u) +
+ 2 / (1-int_cov) * (l - actual) * (l>actual)
+ return(is)
+}
+
+# Compute the interval scores
+IS_base <- mapply(int_score, base_L, base_U, data.matrix(actuals))
+IS_base <- matrix(IS_base, nrow = N)
+IS_rec <- mapply(int_score, rec_L, rec_U, data.matrix(actuals))
+IS_rec <- matrix(IS_rec, nrow = N)
+```
+
+We compute and show the skill scores, which measure the improvement of the
+reconciled forecasts over the base forecasts.
+The skill score is symmetric and bounded between -2 and 2.
+
+```{r}
+SS_AE <- (AE_base - AE_rec) / (AE_base + AE_rec) * 2
+SS_SE <- (SE_base - SE_rec) / (SE_base + SE_rec) * 2
+SS_IS <- (IS_base - IS_rec) / (IS_base + IS_rec) * 2
+SS_AE[is.na(SS_AE)] <- 0
+SS_SE[is.na(SS_SE)] <- 0
+SS_IS[is.na(SS_IS)] <- 0
+
+mean_skill_scores <- c(round(colMeans(SS_IS), 2),
+ round(colMeans(SS_SE), 2),
+ round(colMeans(SS_AE), 2))
+mean_skill_scores <- data.frame(t(matrix(mean_skill_scores, nrow = n)))
+colnames(mean_skill_scores) <- names(actuals)
+rownames(mean_skill_scores) <- c("Interval score", "Squared error", "Absolute error")
+knitr::kable(mean_skill_scores)
+```
+
+The table closely matches Table 4 of the paper.
+In order to exactly reproduce the paper table, it is necessary increase the number of samples drawn from the reconciled distribution to 100,000.
+
+# Reconciled mean and variance
+
+We now show the effects of the reconciliation on the mean and variance of the
+forecast distribution. For more details, we refer to Section 3.2 of the paper.
+
+We observe two different behaviors for the reconciled upper mean:
+it can be between the base and the bottom-up mean (*combination* effect)
+or it can be lower than both (*concordant-shift* effect).
+We show them for two different days.
+
+```{r}
+# Now we draw a larger number of samples:
+N_samples <- 1e5
+
+### Example of concordant-shift effect
+j <- 124
+base_fc_j <- c()
+for (i in 1:n) base_fc_j[[i]] <- c(extr_mkt_events_basefc$mu[[j,i]],
+ extr_mkt_events_basefc$size[[i]])
+# Reconcile
+buis <- reconc_BUIS(S, base_fc_j, "params", "nbinom",
+ num_samples = N_samples, seed = 42)
+samples_y <- buis$reconciled_samples
+
+# The reconciled mean is lower than both the base and bottom-up means:
+means <- c(round(extr_mkt_events_basefc$mu[[j,1]], 2),
+ round(sum(extr_mkt_events_basefc$mu[j,2:n]), 2),
+ round(mean(samples_y[1,]), 2))
+col_names <- c("Base upper mean", "Bottom-up upper mean", "Reconciled upper mean")
+knitr::kable(matrix(means, nrow=1), col.names = col_names)
+
+### Example of combination effect
+j <- 1700
+base_fc_j <- c()
+for (i in 1:n) base_fc_j[[i]] <- c(extr_mkt_events_basefc$mu[[j,i]],
+ extr_mkt_events_basefc$size[[i]])
+# Reconcile
+buis <- reconc_BUIS(S, base_fc_j, "params", "nbinom",
+ num_samples = N_samples, seed = 42)
+samples_y <- buis$reconciled_samples
+
+# The reconciled mean is between the base and the bottom-up mean:
+means <- c(round(extr_mkt_events_basefc$mu[[j,1]], 2),
+ round(sum(extr_mkt_events_basefc$mu[j,2:n]), 2),
+ round(mean(samples_y[1,]), 2))
+col_names <- c("Base upper mean", "Bottom-up upper mean", "Reconciled upper mean")
+knitr::kable(matrix(means, nrow=1), col.names = col_names)
+```
+
+Finally, we show an example in which the variance of the bottom time series increases
+after reconciliation.
+This is a major difference with the Gaussian reconciliation, for which the reconciled
+variance is always smaller than the base variance.
+
+```{r}
+j <- 2308
+base_fc_j <- c()
+for (i in 1:n) base_fc_j[[i]] <- c(extr_mkt_events_basefc$mu[[j,i]],
+ extr_mkt_events_basefc$size[[i]])
+# Reconcile
+buis <- reconc_BUIS(S, base_fc_j, "params", "nbinom", num_samples = N_samples, seed = 42)
+samples_y <- buis$reconciled_samples
+
+# Compute the variance of the base and reconciled bottom forecasts
+base_bottom_var <- mapply(function(mu, size) var(rnbinom(n = 1e5, size = size, mu = mu)),
+ extr_mkt_events_basefc$mu[j,2:n],
+ extr_mkt_events_basefc$size[2:n])
+rec_bottom_var <- apply(samples_y[2:n,], MARGIN = 1, var)
+
+# The reconciled variance is greater than the base variance:
+bottom_var <- rbind(base_bottom_var, rec_bottom_var)
+rownames(bottom_var) <- c("var base", "var reconc")
+knitr::kable(round(bottom_var, 2))
+```
+
+# References
+
+
diff --git a/vignettes/references.bib b/vignettes/references.bib
index e654256..3afd588 100644
--- a/vignettes/references.bib
+++ b/vignettes/references.bib
@@ -141,7 +141,8 @@ @article{zambon2022efficient
number={1},
pages={21},
year={2024},
- publisher={Springer}
+ publisher={Springer},
+ doi={10.1007/s11222-023-10343-y}
}
@@ -151,9 +152,14 @@ @article{zambon2023properties
title={Properties of the reconciled distributions for Gaussian and count forecasts},
author={Zambon, Lorenzo and Agosto, Arianna and Giudici, Paolo and Corani, Giorgio},
journal={arXiv preprint arXiv:2303.15135},
- year={2023}
+ year={2023},
+ eprint={2303.15135},
+ archivePrefix={arXiv},
+ primaryClass={stat.AP},
+ doi={10.48550/arXiv.2303.15135}
}
+
@article{corani2023probabilistic,
title={Probabilistic reconciliation of count time series},
author={Corani, Giorgio and Azzimonti, Dario and Rubattu, Nicol{\`o}},
@@ -171,3 +177,11 @@ @book{hyndman2008forecasting
year={2008},
publisher={Springer Science \& Business Media}
}
+
+@article{agosto2022multivariate,
+ title={Multivariate Score-Driven Models for Count Time Series To Assess Financial Contagion},
+ author={Agosto, Arianna},
+ doi = {10.2139/ssrn.4119895},
+ year={2022}
+}
+