predict.tulip: Forecast by drawing sample paths from a fitted object
In timradtke/heuristika: Robust Probabilistic Forecasts to Tinker With

predict.tulip

R Documentation

Forecast by drawing sample paths from a fitted object

Description

Generate a forecast in the form of n sample paths for the forecast horizon h based on the fitted model object.

Usage

## S3 method for class 'tulip'
predict(
  object,
  h = 12,
  n = 10000,
  remove_anomalies = TRUE,
  postprocess_sample = identity,
  ...
)

Arguments

`object`	The fitted model object returned by `tulip()` of class `tulip`.
`h`	The forecast horizon as integer number of periods.
`n`	The integer number of sample paths to draw from the forecast distribution.
`remove_anomalies`	Logical (default `TRUE`); during prediction, some generated samples can have the characteristics of anomalies; as such, they can be identified and interpolated to not adversely affect the state updates when predicting multiple horizos, similar to `remove_anomalies` in `tulip()`. The interpolated values are only used to fit the states. The returned sample paths will still contain the anomalous samples. Which samples are anomalous is determined based on the residuals from the original model fit.
`postprocess_sample`	A function that is applied on a numeric vector of drawn samples for a single step-ahead before the samples are used to update the state of the model, and before outliers are removed (if applicable). By default equal to `identity()`, but could also be something like `function(x) pmax(x, 0)` to enforce a lower bound of 0, or any other transformation of interest. This works best with `family = "bootstrap"` when the applied transformation represents a feature of the data generating process that is non-Gaussian. Note that this can cause arbitrary errors caused by the author of the function provided to `postprocess_sample`.
`...`	arguments passed to or from other methods.

Details

The forecasts of a tulip model are natively represented by sample paths. Instead of providing a single point forecast, or marginal quantiles, possible future outcomes are summarized through samples from the joint distribution across future horizons. These samples can be summarized down to marginal quantiles, but they also contain the dependency between future horizons that is commonly lost by forecast frameworks that solely store marginal quantile forecasts.

Depending on the chosen family, the samples are either based on a distribution that is fitted against the model's residuals during the fitting procedure; or they can be based directly on bootstrapped residuals if the family of the tulip object is overwritten to be bootstrap. In the latter case, the forecast distribution can represent asymmetries that were observed in the input data. You might not want to use bootstrap when the input time series is very short. In such cases, a parametric assumption can help with more reasonable behavior of the error component.

As for all exponential smoothing models, the one-step-ahead sample is fed back into the model to update the state components, before the two-step-ahead forecast sample is drawn given the updated state components. Through this mechanism, the temporal dependencies between samples of the same sample path are maintained. For example, if the initial sample is higher than average, the level and trend components of the model might update to be larger, thus consequently increasing chances that the next sample is again high (compared to the most recently observed training observation and level).

If the model is able to represent the input data well, then sample paths look like reasonable future continuations of the input series. Averaging across sample paths one can derive quantiles of the forecast distribution.

Importantly, because the sample paths maintain the temporal dependencies, one can also aggregate the forecasts across forecast horizons which is not possible when only quantiles are returned by a model. The ability to summarize arbitrary aggregates can be helpful for optimization problems that use the forecast as input.

The returned value is a matrix of dimensions h times n. Each of the n columns represents one sample path and is a random draw from the h-dimensional forecast distribution. See also examples below.

Use the postprocess_sample argument to hook into the sampling process of the predict() method. Provide a function that is applied on the current j-step-ahead vector of samples before it is fed back into the model, before the states are updated given those samples. This allows you to adjust the data generating process–for example, to enforce non-negative count samples– while doing so before the state components (level, trend, seasonality) of the model are updated given the most recent sample. This can prevent the model's level from drifting to zero, for example, because samples never became less than zero, thus affecting subsequent samples to drift even further below zero, and so on. Generated samples stay similar to the input data, allowing the forecast to stay similar to the input data.

Value

An object of class tulip_paths that is a list of:

paths: A matrix of dimensions (h, n), in which each of the n columns represents one sample path over the h-steps-ahead
model: The provided model object

References

Chapter 4.2 of: Syama Sundar Rangapuram, Matthias W. Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, Tim Januschowski (2018). Deep State Space Models for Time Series Forecasting.: https://papers.nips.cc/paper/2018/hash/5cf68969fb67aa6082363a6d4e6468e2-Abstract.html
Chapter 6.1, for example, of: Rob J. Hyndman, Anne B. Koehler, Ralph D. Snyder, and Simone Grose (2002). A State Space Framework for Automatic Forecasting using Exponential Smoothing Methods.: https://doi.org/10.1016/S0169-2070(01)00110-8
Matthias Seeger, Syama Rangapuram, Yuyang Wang, David Salinas, Jan Gasthaus, Tim Januschowski, Valentin Flunkert (2017). Approximate Bayesian Inference in Linear State Space Models for Intermittent Demand Forecasting at Scale.: https://arxiv.org/abs/1709.07638
Matthias Seeger, David Salinas, Valentin Flunkert (2016). Bayesian Intermittent Demand Forecasting for Large Inventories.: https://proceedings.neurips.cc/paper/2016/file/03255088ed63354a54e0e5ed957e9008-Paper.pdf
Alexander Alexandrov, Konstantinos Benidis, Michael Bohlke-Schneider, Valentin Flunkert, Jan Gasthaus, Tim Januschowski, Danielle C. Maddix, Syama Rangapuram, David Salinas, Jasper Schulz, Lorenzo Stella, Ali Caner Türkmen, Yuyang Wang (2019). GluonTS: Probabilistic Time Series Models in Python.: https://arxiv.org/abs/1906.05264

Examples

fitted_model <- tulip(y = tulip::flowers$flowers, m = 12)
forecast <- predict(object = fitted_model, h = 12, n = 10000)

# library(ggplot2)

# visualize the marginal quantile forecast
# autoplot(forecast)

# visualize a handful of sample paths directly
# autoplot(forecast, method = "paths")

# let's take a closer look at the sample paths themselves
m_fc <- forecast$paths

# summarize over draws (columns) to get point forecasts
rowMeans(m_fc)

# marginal quantiles can be derived in much the same way
apply(m_fc, 1, quantile, 0.05)
apply(m_fc, 1, quantile, 0.95)

# we can also summarize the distribution of the sum over horizons 4 to 6
quantile(colSums(m_fc[4:6, ]), c(0.05, 0.5, 0.95))

# sample paths carry autocorrelation, the horizons are not independent of
# another; this is revealed when we drop the autocorrelation by shuffling
# at each margin randomly and try to recompute the same quantiles as above:

m_fc_shuffled <- rbind(
  m_fc[4, sample(x = 1:10000, size = 10000, replace = TRUE), drop = FALSE],
  m_fc[5, sample(x = 1:10000, size = 10000, replace = TRUE), drop = FALSE],
  m_fc[6, sample(x = 1:10000, size = 10000, replace = TRUE), drop = FALSE]
)
quantile(colSums(m_fc_shuffled), c(0.05, 0.5, 0.95))

# one can also look  directly at the correlation
cor(m_fc[4, ], m_fc[6, ])
cor(m_fc_shuffled[1, ], m_fc_shuffled[3, ])

# -------------------

# An example of how the `postprocess_sample` argument can be used to
# enforce non-negative forecasts

set.seed(7759)

y <- rpois(n = 60, lambda = 1)

ls_fit <- tulip(y = y, m = 12, family = "norm", method = "additive")
ls_fit$family <- "bootstrap"

m_fc <- predict(
  object = ls_fit,
  h = 12,
  n = 10000,
  postprocess_sample = function(x) pmax(0, x)
)

# library(ggplot2)
# autoplot(m_fc)

timradtke/heuristika documentation built on April 24, 2023, 1:55 a.m.