generate() and reconciliation

I've been exploring a number of models for a hierarchical time-series and encountered situations where generate() returns samples that are not consistent with the variances returned by forecast(). Here is a hopefully-minimal working example that demonstrates this:

suppressPackageStartupMessages(library(fpp3))

set.seed(12345)

tourism_fit <- tsibble::tourism |>
  aggregate_key(State, Trips = sum(Trips)) |>
  model(ets = ETS(Trips)) |>
  reconcile(bu = bottom_up(ets)) |>
  select(bu)

tourism_fit |>
  forecast(h = 1) |>
  filter(is_aggregated(State)) |>
  pull(Trips) |>
  distributional::variance() |>
  print()

tourism_fit |>
  forecast(h = 1, simulate = TRUE) |>
  filter(is_aggregated(State)) |>
  pull(Trips) |>
  distributional::variance() |>
  print()

tourism_fit |>
  filter(is_aggregated(State)) |>
  generate(h = 1, times = 10000) |>
  pull(.sim) |>
  var() |>
  print()

This outputs:

[1] 480068
[1] 471076.9
[1] 704529.2

For reference, in my original problem generate() returned samples with much smaller variance than the distribution returned by forecast() (roughly 600 vs 3000).

Run the minimal R example using reconcile(bu = bottom_up(ets)), then compare the variances from forecast(), forecast(simulate = TRUE), and generate() for the aggregated tourism series. Trace the generate() and forecast() entry points involved in reconciliation. Done means the generated sample variance is consistent with the forecast distribution for this reproducible case.