1

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" )
Village.Idyot
  • 1,359
  • 2
  • 8
  • 1
    You can obtain the parameters from a distribution (in this case the samples from the sample distribution) using the `parameters()` function on the distribution. Try `parameters(fc$StateX)`. – Mitchell O'Hara-Wild Jan 25 '23 at 12:45
  • I tried `parameters()` but I get the error message "Error in parameters(fc$StateX) : could not find function "parameters". Is `parameters()` part of a package? – Village.Idyot Jan 25 '23 at 14:03
  • 1
    Yes, sorry! The `parameters()` function is from the `{distributional}` package (where the distributions are provided). Try `distributional::parameters(fc$StateX)`. – Mitchell O'Hara-Wild Jan 26 '23 at 20:22
  • Ahhh I see you are the author of the `distributional` package. I downloaded the reference manual from CRAN and am reading through it. – Village.Idyot Jan 27 '23 at 06:00

1 Answers1

1

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")
Village.Idyot
  • 1,359
  • 2
  • 8