---
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
<div id="refs"></div>

