size/speed of forecasting sample paths can be prohibitive

[davidtedfordholt]

Not a serious problem, but worth talking about in the documentation. I'm happy to add something to the Details section about uncertainty.samples (which appears as times based on fable::ARIMA). I'd prefer for times to be replaced with something more intuitive, but that's a preference thing and not important.

The most important part is making sure people understand how to deal with using fable.prophet when forecasting a large number of models. A significant issue comes when parallelising a large number of time series/models. The model objects can be large enough that combining the output of futures can fail. If the user is primarily in need of a point forecast (yes, I know!), and the requirements for accuracy allow, the speed increase and size decrease can be pretty huge. For shorter time horizons, it has less of an impact on accuracy, as well.

To super extra prove that I'm aware of the potentially problematic nature of reducing the number of sample paths, I've gone ahead and created an example. Obviously, there are series where it's appropriate and series where it isn't in this dataset.

library(dplyr)
library(ggplot2)
library(fable.prophet)
library(future)
options(future.fork.enable = TRUE)

models <-
    tsibbledata::gafa_stock %>%
    tsibble::group_by_key() %>%
    tsibble::fill_gaps() %>%
    tidyr::fill(Open:Volume, .direction = "down") %>%
    tidyr::pivot_longer(Open:Volume) %>%
    model(data, prophet(value ~ growth() + season("week") + season("year")))

plan("sequential")
system.time(fcast_1000 <- forecast(models, h = 365, times = 1000))
system.time(fcast_100 <- forecast(models, h = 365, times = 100))
system.time(fcast_10 <- forecast(models, h = 365, times = 10))

plan("multisession")
system.time(fcast_1000_ms <- forecast(models, h = 365, times = 1000))
system.time(fcast_100_ms <- forecast(models, h = 365, times = 100))
system.time(fcast_10_ms <- forecast(models, h = 365, times = 10))

plan("multicore")
system.time(fcast_1000_mc <- forecast(models, h = 365, times = 1000))
system.time(fcast_100_mc <- forecast(models, h = 365, times = 100))
system.time(fcast_10_mc <- forecast(models, h = 365, times = 10))

format(object.size(fcast_1000), unit = "auto")
format(object.size(fcast_100), unit = "auto")
format(object.size(fcast_10), unit = "auto")

fcasts <- as_tibble(fcast_1000) %>%
    select(Symbol, name, Date, .mean) %>%
    rename(.mean_1000 = .mean) %>%
    full_join(as_tibble(fcast_100) %>%
                  select(Symbol, name, Date, .mean) %>%
                  rename(.mean_100 = .mean), 
              by = c("Symbol", "name", "Date")) %>%
    full_join(as_tibble(fcast_10) %>%
                  select(Symbol, name, Date, .mean) %>%
                  rename(.mean_10 = .mean), 
              by = c("Symbol", "name", "Date"))

fcasts %>%
    ggplot(aes(x = Date)) +
    geom_point(aes(y = (.mean_1000 - .mean_10)/.mean_1000), colour = "red", size = .5) +
    geom_point(aes(y = (.mean_1000 - .mean_100)/.mean_1000), colour = "blue", size = .5) +
    geom_line(aes(y = 0), colour = "green") +
    scale_y_continuous(labels = scales::percent) +
    ggtitle("Percent Difference from Mean of 1000 samples",
            subtitle = "Mean of 100 samples (blue) and 10 samples (red)") +
    labs(y = "Percent diff") + 
    facet_wrap(name ~ Symbol, scales = "free")

fcasts %>%
    ggplot(aes(x = Date)) +
    geom_line(aes(y = .mean_10), colour = "red", size = .5) +
    geom_line(aes(y = .mean_100), colour = "blue", size = .5) +
    geom_line(aes(y = .mean_1000), colour = "green", size = .5) +
    scale_y_continuous(labels = scales::dollar) +
    ggtitle("Point forecasts with 1000, 100 and 10 samples",
            subtitle = "1000: green, 100: blue, 10: red") +
    labs(y = "Opening Price") + 
    facet_wrap(name ~ Symbol, scales = "free")