--- title: "t-Rec: Reconciliation of Gaussian forecasts by modelling the uncertainty on the covariance matrix" author: "Dario Azzimonti, Chiara Carrara, Lorenzo Zambon, Giorgio Corani" date: "" lang: "en" output: rmarkdown::html_vignette bibliography: references.bib cite: - '@zambon2024properties' - '@carrara2025' vignette: > %\VignetteIndexEntry{t-Rec: Reconciliation of Gaussian forecasts by modelling the uncertainty on the covariance matrix} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{css, echo=FALSE} p.caption { text-align: left; } caption { text-align: left; } ``` # Introduction This vignette showcases the method *t-Rec*, a Bayesian approach to forecast reconciliation that explicitly accounts for uncertainty in the covariance matrix of the residuals. t-Rec exploits the flexibility of Bayesian models to incorporate parameter uncertainty: it assumes an Inverse-Wishart prior on the covariance matrix. This choice allows the reconciliation to be derived in closed form, resulting in a reconciled predictive distribution that follows a multivariate Student's t. See @carrara2025 for a detailed explanation. ```{r load, message=FALSE, warning=FALSE} # load the packages library(bayesRecon) library(forecast) # base forecasts library(ggplot2) # plots ``` # Loading the data We consider the monthly Swiss overnight tourist stays data divided by canton and aggregated over the whole country. This is a cross-sectional hierarchical structure, consisting of two levels: the national total and the division by canton, for a total of 27 time series (1 upper and 26 bottom). The dataset is available in this package and can be loaded as `bayesRecon::swiss_tourism`. This is a list that contains the `mts` object with the time series (`ts`), the aggregation matrix (`agg_mat`) and the number of upper (`n_upper`) and bottom (`n_bottom`) time series. The raw data is also publicly available on the official website of the Swiss Confederation at this [link](https://www.bfs.admin.ch/asset/en/px-x-1003020000_102). The dataset spans the period from January 2005 to January 2025, comprising 241 monthly observations. The time series exhibit strong seasonality. **Figure 1** shows the top-level (aggregate) time series. ```{r swiss-tourism-plot, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 1**: Swiss tourism: monthly overnight stays in Switzerland.", fig.dim = c(8, 4)} # save all time series Y = swiss_tourism$ts # plot the first (top) time series autoplot(Y[,1], ylab = "Overnight stays in Switzerland",linewidth=0.9)+ scale_y_continuous(labels = function(x) paste0(formatC(x / 1e6, format = "g"), "M")) ``` We load the aggregation matrix and select the length of the training set. The length of the training set can be selected within the range [14, 240] observations: at least 14 observations are needed to initialise the prior from the residuals of the naive forecasts of these monthly series, and at most 240 so that one observation is left for the forecast horizon. ```{r Swiss tourism} # Save aggregation matrix A = swiss_tourism$agg_mat # Number of bottom and upper time series n_b = ncol(A) n_u = nrow(A) n = n_b + n_u # Frequency is monthly: print(frequency(Y)) # Select the length of the training set and the forecast horizon L = 60 h = 1 # the code below assumes h = 1 # Select the training set and the actuals for the forecast horizon train = window(Y, end = time(Y)[L]) actuals = window(Y, start= time(Y)[L + 1], end = time(Y)[L + h]) ``` # Forecasts The base forecasts are generated using the `ets()` function from the `forecast` package, which fits individual exponential smoothing models to each time series. For each time series, we save the base forecast mean (`base_fc`) and the residuals from the fitted model (`res`), which are used to estimate the covariance matrix of the forecast errors. ```{r Base-forecasts and residuals computation} # Compute base forecasts and residuals for each time series base_fc = rep(NA, n) res = matrix(NA, ncol = n, nrow = L) for (k in 1:n){ fit = forecast::ets(train[,k], model = "AZZ") f = forecast(fit, h = h) base_fc[k] = f$mean res[, k] = fit$residuals } ``` # Reconciliation We can now reconcile the forecasts with the function `reconc_t()`, which implements the t-Rec method. This function takes as input the aggregation matrix `A`, the base forecast means `base_fc_mean`, the training data `y_train` and the model residuals `residuals` and returns the parameters of the reconciled forecasts, which are distributed as a multivariate Student's t. The flag `return_upper = TRUE` makes the function return also the parameters of the upper-level reconciled forecasts, while the flag `return_parameters = TRUE` makes it return also the parameters of the posterior distribution of the covariance matrix; they will be used to compare the covariance estimates. ```{r t-reconciliation} t_rec_results = reconc_t(A, base_fc_mean = base_fc, y_train = train, residuals = res, return_parameters = TRUE, return_upper = TRUE) ``` For the selected data window, we save the base forecasts (Base) together with the covariance matrix of their residuals. ```{r Base} # Base forecasts Base_mean = base_fc Base_cov_mat = crossprod(res)/nrow(res) # covariance of the residuals ``` We further compute the reconciled forecasts with the standard Gaussian reconciliation (MinT, @wickramasuriya2019optimal). Note the flag `return_upper = TRUE`, which makes the function return also the parameters of the upper-level reconciled forecasts. ```{r Gaussian} # Gaussian/MinT: compute reconciliation with bayesRecon gauss_results = reconc_gaussian(A, base_fc, residuals = res, return_upper = TRUE) # Reconciled mean for the whole hierarchy: MinT_reconciled_mean = c(gauss_results$upper_rec_mean, gauss_results$bottom_rec_mean) ``` # Comparison of results We now compare the predictive densities for the 1-step ahead forecasts of the upper-level time series (Switzerland) obtained by the three methods above: the base forecasts (Base), Minimum Trace reconciliation (MinT), and t-Rec. ```{r compute uppers, echo=FALSE, eval=TRUE} t_Rec_reconciled_mean <- c(t_rec_results$upper_rec_mean,t_rec_results$bottom_rec_mean) t_Rec_corr_with_upper = A %*% t_rec_results$bottom_rec_scale_matrix t_Rec_scale_par <- rbind(cbind(t_rec_results$upper_rec_scale_matrix, t_Rec_corr_with_upper), cbind(t(t_Rec_corr_with_upper), t_rec_results$bottom_rec_scale_matrix)) # Index of the upper variable to plot i_upper <- 1 # change this if a different index is needed # Extract distribution's parameters for each method # MinT # Build the full reconciled covariance matrix MinT_corr_with_upper = A %*% gauss_results$bottom_rec_cov MinT_cov_mat = rbind(cbind(gauss_results$upper_rec_cov, MinT_corr_with_upper), cbind(t(MinT_corr_with_upper),gauss_results$bottom_rec_cov)) mu_MinT <- as.numeric(MinT_reconciled_mean[i_upper]) sd_MinT <- sqrt(MinT_cov_mat[i_upper, i_upper]) # t-Rec mu_tRec <- as.numeric(t_Rec_reconciled_mean[i_upper]) scale_tRec <- sqrt(t_Rec_scale_par[i_upper, i_upper]) df_tRec <- t_rec_results$upper_rec_df # Base mu_Base <- as.numeric(Base_mean[i_upper]) sd_Base <- sqrt(Base_cov_mat[i_upper, i_upper]) # Create a grid of x values x_vals <- seq(min(mu_MinT, mu_tRec, mu_Base) - 4*max(sd_MinT, scale_tRec, sd_Base), max(mu_MinT, mu_tRec, mu_Base) + 4*max(sd_MinT, scale_tRec, sd_Base), length.out = 1000) # Compute densities dens_MinT <- dnorm(x_vals, mean = mu_MinT, sd = sd_MinT) dens_tRec <- dt((x_vals - mu_tRec) / scale_tRec, df = df_tRec) /scale_tRec dens_Base <- dnorm(x_vals, mean = mu_Base, sd = sd_Base) ``` ```{r comparison-plot, echo=FALSE, eval=TRUE, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 2**: Predictive densities of the upper time series obtained with MinT (purple), t-Rec (green) and Base (blue). The black triangle indicates the actual value.", fig.dim = c(8, 4), warning=FALSE} dens_df <- data.frame( x = rep(x_vals, 3), Method = rep(c("MinT", "t-Rec", "Base"), each = length(x_vals)), Density = c(dens_MinT, dens_tRec, dens_Base) ) method_colors <- c( "MinT" = "#440154", "t-Rec" = "#29AF7F", "Base" = "#3B528B" ) ggplot(dens_df, aes(x = x, y = Density, color = Method)) + geom_area( data = dens_df[dens_df$Method %in% c("MinT", "t-Rec"),], aes(fill = Method), position = "identity", alpha = 0.3, color = NA ) + geom_line(linewidth = 1.5) + scale_color_manual(values = method_colors) + scale_fill_manual(values = method_colors) + geom_point(aes(x = actuals[i_upper], y = 0, shape = "Actual value"), color = "black", size = 3) + scale_shape_manual(values = c("Actual value" = 17)) + scale_x_continuous(labels = function(x) paste0(formatC(x / 1e3, format = "g"), "k")) + guides(fill = "none") + labs( title = "", #"Predictive densities of upper time series", x = "", y = "" ) + theme_minimal(base_size = 12) + theme( legend.position = "bottom", legend.title = element_blank(), plot.title = element_text(size = 18) ) ``` The base forecast and the MinT reconciled forecast are normally distributed, while the t-Rec reconciled forecast follows a Student's t distribution. Both reconciliation methods improve the forecast: their means are closer to the actual value. However, the t-Rec method returns a wider density, which is more likely to contain the actual value. This is because t-Rec accounts for the uncertainty in the covariance matrix of the forecast errors, which leads to a more realistic representation of the forecast uncertainty. # Comparison of the covariance matrix estimates The main strength of t-Rec is that it provides an estimate of the covariance which quantifies its own uncertainty. To illustrate this point, we compare the estimates of the covariance matrix of the base forecast errors used by t-Rec and by the standard Gaussian reconciliation method (MinT), both computed before the reconciliation step. In t-Rec, the covariance matrix of the base forecast errors is estimated in a Bayesian way: we put an Inverse-Wishart prior on the covariance, we assume that the residuals of the base forecasts are Gaussian distributed and we obtain a posterior distribution which is again Inverse-Wishart with known parameters $\nu$ and $\Psi$. The output of the function `reconc_t()` stores those parameters in `t_rec_results$posterior_nu` and `t_rec_results$posterior_Psi`. Note that those are not the parameters of the reconciled forecasts, but the parameters of the posterior distribution of the covariance matrix of the forecast errors. In this section, we denote this estimate as '*IW posterior*' to distinguish it from the reconciled covariance. Since the posterior is an Inverse-Wishart, we can write the marginal distribution of the variances in closed-form as an inverse Gamma distribution with the following parameters: $$ \Sigma_{ii} \sim \text{Inv-Gamma}\left(\frac{\nu - n + 1}{2}, \frac{\Psi_{ii}}{2}\right), $$ where $\nu$ and $\Psi$ denote the posterior Inverse-Wishart parameters and $n$ is the total number of series in the hierarchy. In the second case (MinT), instead, a point estimate of the covariance matrix is obtained by applying the Schäfer Strimmer shrinkage estimator [@schafer2005shrinkage] to the covariance of the residuals; this method is denoted here as '*Schäfer Strimmer*'. For the Swiss tourism forecasts computed above, we focus on the covariance between the upper-level series, denoted as CH, and the bottom-level time series with the largest average number of overnight stays, "Graubünden", denoted as GR. Analogous considerations apply to the other series: the code below is parametrised by the indices `i` and `j`, so that other pairs of series can be inspected by changing them. ```{r compute shrunk matrix} # Full shrinkage covariance matrix with the Schäfer Strimmer shrinkage estimator # The same matrix is computed internally by `reconc_gaussian()` before reconciliation shrink_mat <- bayesRecon::schaferStrimmer_cov(res)$shrink_cov ``` We show here the standard deviation instead of the variance because it is easier to interpret. The density of the standard deviation is obtained from the density of the variance (defined above) through the change-of-variable formula, implemented in the function below: ```{r compute density} # Select which series to plot i = 1 # CH j = 19 # GR # density of the inverse gamma dinvgamma <- function(x, shape, rate) { dgamma(1/x, shape = shape, rate = rate) / x^2 } # density of the standard deviation (square root transform of variance) d_std_dev <- function(x, shape,rate){ dinvgamma(x^2, shape = shape, rate = rate)*2*x } # compute the density of the standard deviation for CH; # mean_ch is the posterior mean of the variance, used only to centre the grid mean_ch <- t_rec_results$posterior_Psi[i, i]/(t_rec_results$posterior_nu - n - 1) x_ch <- sqrt(seq(mean_ch*0.5, mean_ch*1.7, length.out = 1000)) shape_ch <- (t_rec_results$posterior_nu - n + 1) /2 rate_ch <- t_rec_results$posterior_Psi[i, i] / 2 dens_ch <- d_std_dev(x_ch, shape = shape_ch, rate = rate_ch) # compute the density of the standard deviation for GR mean_gr <- t_rec_results$posterior_Psi[j, j]/(t_rec_results$posterior_nu - n - 1) x_gr <- sqrt(seq(mean_gr*0.5, mean_gr*1.7, length.out = 1000)) shape_gr <- (t_rec_results$posterior_nu - n + 1) /2 rate_gr <- t_rec_results$posterior_Psi[j, j] / 2 dens_gr <- d_std_dev(x_gr, shape = shape_gr, rate = rate_gr) ``` **Figure 3** shows the density of the posterior standard deviation of the forecasts for the upper time series (CH) and the bottom time series (GR). In each panel, the vertical dashed line marks the corresponding point estimate obtained with *Schäfer Strimmer*. ```{r density plot, echo=FALSE, eval=TRUE, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 3**: Density of the posterior standard deviation of the forecasts' residuals, estimated with *IW posterior*. The dashed line is the *Schäfer Strimmer* estimate.", fig.dim = c(8, 4), warning=FALSE} df_dens <- data.frame( x = c(x_ch, x_gr), y = c(dens_ch, dens_gr), panel = rep(c("CH", "GR"), each = length(dens_ch)) ) # Schäfer Strimmer point estimates of the standard deviations df_shrink <- data.frame( panel = c("CH", "GR"), x = sqrt(c(shrink_mat[i, i], shrink_mat[j, j])), y = c(max(dens_ch), max(dens_gr)) * 1.03 ) # Plot the density of the standard deviation parameters for CH and GR ggplot(df_dens, aes(x = x, y = y)) + geom_area(fill = "#6B8E23", alpha = 0.4) + geom_line(color = "black") + geom_vline(data = df_shrink, aes(xintercept = x), linetype = "dashed", color = "black") + geom_text(data = df_shrink, aes(x = x, y = y, label = "Schäfer Strimmer"), hjust = -0.05, vjust = 0, size = 4) + # leave room above the curve for the label geom_blank(data = df_shrink, aes(x = x, y = y * 1.12)) + scale_x_continuous(labels = function(x) paste0(formatC(x / 1e3, format = "g"), "k")) + guides(fill = "none") + labs( title = "", #"Posterior standard deviation of the forecasts' residuals", x = "", y = "" ) + facet_wrap(~ panel, nrow = 1, scales = "free") + theme_minimal(base_size = 12) + theme( legend.position = "top", legend.title = element_blank(), strip.text = element_text(size = 16), plot.title = element_text(size = 18) ) ``` The posterior distribution for the covariance and for the correlation values is not available in closed form, but it can be obtained via sampling. Since the posterior distribution is an Inverse-Wishart distribution, we can sample from it with the custom function `rinvwishart()`, defined below, which fixes the seed for reproducibility. ```{r generate IW samples} # generate k samples from an IW(Psi, nu) distribution rinvwishart <- function(k, nu, Psi, seed=42) { p <- nrow(Psi) Sigma <- solve(Psi) set.seed(seed) all_W <- rWishart(k, df = nu, Sigma = Sigma) W <- array(NA, dim = c(p, p, k)) for (i in 1:k) { W[,,i] <- solve(all_W[,,i]) } return(W) } IW_post_samples <- rinvwishart(k = 1000, nu = t_rec_results$posterior_nu, Psi = t_rec_results$posterior_Psi) ``` ```{r compute correlations, echo=FALSE, eval=TRUE} corr_shrink <- cov2cor(shrink_mat) corr_tRec <- cov2cor(t_rec_results$posterior_Psi/(t_rec_results$posterior_nu - n - 1)) ``` **Figure 4** shows the posterior density of the correlation between CH and GR. The value estimated with *Schäfer Strimmer*, plotted as a vertical dashed line, lies in the lower tail of the *IW posterior* and differs from its mode, showing that the two estimates give a different picture of the dependence structure. ```{r plot densities, echo=FALSE, eval=TRUE, dpi=300, out.width = "100%", fig.align='center', fig.cap="**Figure 4**: Density of the posterior correlation between CH and GR obtained with *IW posterior*. The dashed line is the *Schäfer Strimmer* estimate.", fig.dim = c(8, 4), warning=FALSE} # Compute correlation samples for all sample matrices generated above corr_post <- array(apply(IW_post_samples, FUN = function(M) cov2cor(M), MARGIN = c(3)), dim = dim(IW_post_samples)) df_corr <- data.frame(x = corr_post[i, j, ]) ggplot(df_corr, aes(x = x)) + geom_density(fill = "#6B8E23", alpha = 0.4) + geom_vline(xintercept = corr_shrink[i, j], linetype = "dashed", color = "black") + annotate("text", x = corr_shrink[i, j], y = max(density(corr_post[i, j, ])$y) * 0.95, label = "Schäfer Strimmer", hjust = -0.05, vjust = 0, size = 4.5) + labs(x = "", y = "Density", title = "" #expression(paste("Posterior correlation ", rho["CH,GR"])) ) + theme_minimal() + theme( legend.position = "none", plot.title = element_text(size = 18) ) ``` The two estimates therefore differ in both values and interpretation. *IW posterior* returns a density while *Schäfer Strimmer* is a point estimate. Moreover, the standard deviations estimated by *IW posterior* (both mean and mode) are larger than the *Schäfer Strimmer* point estimates (**Figure 3**). The correlation between CH and GR is also estimated differently: `r round(corr_shrink[i, j], 2)` with *Schäfer Strimmer* against a posterior mean of `r round(corr_tRec[i, j], 2)` with *IW posterior* (**Figure 4**). Accounting for this additional uncertainty is what makes the t-Rec predictive distribution wider than the MinT one in **Figure 2**. # References