rforecastingfable-r

How to extract simulation data from Fable package forecast?


In running the R code posted at the bottom, I derive a forecast for the next 10 periods at the 80% and 95% confidence levels, using forecast() function from the fable package and running 1000 simulation sample paths, as illustrated here:

enter image description here

The resulting fable object looks like this, in the R Studio console: enter image description here

I'd like to access the simulation paths from the above Fable object so I can plot a distribution of forecasts for example at period 20, as conceptually shown in the example in the below. Any ideas how to do this? enter image description here

Code:

library(feasts)
library(fable)
library(fabletools)
library(ggplot2)
library(tsibble)

tmp <- data.frame(
  Month = c(1,2,3,4,5,6,7,8,9,10),
  StateX = c(1527,1297,933,832,701,488,424,353,302,280)
  ) %>%
  as_tsibble(index = Month)

fit <-  tmp %>% model(NAIVE(StateX))

fc <- fit %>% forecast(h = 10, bootstrap = TRUE, times = 1000)

autoplot(fc, tmp) +
  labs(title="Transitions to Dead State X", y="Units" )

Solution

  • Please see the post at this link which more completely resolves this question: https://community.rstudio.com/t/how-to-derive-a-forecast-density-for-a-specified-point-in-time-with-the-r-fable-package-forecast-function/158329. Package ggdist provides the necessary functionality. Adapting to the OP code example, see solution code below:

    library(dplyr)
    library(fable)
    library(ggdist)
    library(ggplot2)
    library(tsibble)
    
    tmp <- data.frame(
      Month = c(1,2,3,4,5,6,7,8,9,10),
      StateX = c(1527,1297,933,832,701,488,424,353,302,280)
      ) %>%
      as_tsibble(index = Month)
    
    # Fit model
    fit <- tmp |>
      model(naive = NAIVE(StateX))
    
    # Generate forecasts
    fc <- bind_rows(
      forecast(fit, h = 10) |> mutate(.model = "Theoretical"),
      forecast(fit, h = 10, bootstrap = TRUE, times = 1000) |> mutate(.model = "Bootstrapped")
    )
    
    # Density + histogram plot (alpha below sets opacity)
    ggplot() +
      stat_dist_slab(aes(dist = StateX, y = 0, fill = "Theoretical"),
                     data = fc |> filter(.model == "Theoretical"),
                     colour = NA, alpha = 0.3
      ) +
      geom_histogram(aes(x = x, y = after_stat(ndensity), fill = "Bootstrapped"),
                     data = tibble(x = fc |>
                                     filter(.model == "Bootstrapped") |>
                                     pull(StateX) |>
                                     unlist()),
                     colour = NA, bins = 50, alpha = 0.3
      ) +
      scale_color_brewer(palette = "Set1") +
      scale_fill_brewer(palette = "Set1") +
      labs(x = "Units reaching dead state X (at h=10)", y = "Forecast density") +
      theme(legend.position = "none")