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:
The resulting fable object looks like this, in the R Studio console:
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?
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" )
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")