tidyverts / fabletools

General fable features useful for extension packages
http://fabletools.tidyverts.org/
89 stars 31 forks source link

generate() and reconciliation #408

Open robmoss opened 3 weeks ago

robmoss commented 3 weeks ago

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

robmoss commented 3 weeks ago

Also, when I include a top-down reconciliation model, forecast(..., simulate = TRUE) returns a fixed value for each horizon:

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),
    td = top_down(ets)
  ) |>
  select(bu, td)

tourism_fit |>
  forecast(h = 1) |>
  filter(is_aggregated(State)) |>
  pull(Trips) |>
  distributional::variance() |>
  print()
# [1] 480068.0 699901.4

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

In this case generate() results in an error (see below), but I don't know if that's me asking it to do something it shouldn't be expected to handle.

Error in `pivot_longer_spec()`:
! Can't combine `bu` <lst_btmup_mdl> and `td` <lst_topdwn_mdl>.
Backtrace:
     ▆
  1. ├─base::print(...)
  2. ├─stats::var(...)
  3. │ └─base::is.data.frame(x)
  4. ├─dplyr::pull(...)
  5. ├─generics::generate(...)
  6. ├─fabletools:::generate.mdl_df(...)
  7. │ ├─tidyr::pivot_longer(...)
  8. │ └─tidyr:::pivot_longer.data.frame(...)
  9. │   └─tidyr::pivot_longer_spec(...)
 10. │     └─vctrs::vec_ptype_common(...)
 11. └─vctrs (local) `<fn>`()
 12.   └─vctrs::vec_default_ptype2(...)
 13.     ├─base::withRestarts(...)
 14.     │ └─base (local) withOneRestart(expr, restarts[[1L]])
 15.     │   └─base (local) doWithOneRestart(return(expr), restart)
 16.     └─vctrs::stop_incompatible_type(...)
 17.       └─vctrs:::stop_incompatible(...)
 18.         └─vctrs:::stop_vctrs(...)
 19.           └─rlang::abort(message, class = c(class, "vctrs_error"), ..., call = call)