generate() and reconciliation

[robmoss]

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).