tulip: Fit a robust exponential smoothing model by...

In timradtke/heuristika: Robust Probabilistic Forecasts to Tinker With

Fit a robust exponential smoothing model by maximum-a-posteriori estimation

Description

Given a univariate time series y, tulip fits an exponential smoothing model using maximum-a-posteriori (MAP) estimation. Prior distributions for the smoothing parameters, the error component, and the probability of anomalies can be provided via priors. The set of available parameter combinations is defined by param_grid over which the optimization of the MAP procedure takes place. The error and trend component are additive (if used), whereas the seasonal component can be additive or multiplicative.

Usage

tulip(
  y,
  m,
  method = c("additive", "multiplicative")[1],
  family = c("norm", "student", "cauchy")[1:2],
  param_grid = NULL,
  priors = NULL,
  init_states = NULL,
  seasonality_threshold = 0.5,
  remove_anomalies = TRUE,
  anomaly_budget = 5,
  anomaly_budget_most_recent_k = 1,
  min_obs_anomaly_removal = 12,
  try_fixed_initial_fit = FALSE,
  check_param_grid_unique = TRUE
)

Arguments

y	A time series as numeric vector, may include NAs for some (but not all) of the observations
m	The time series' period length in number of observations (for example, 12 for yearly seasonality with monthly observations, 1 for no seasonality, ...); does not handle multiple seasonality (e.g. weekly and yearly for daily observations)
method	One of additive or multiplicative; determines whether the model uses an additive or multiplicative seasonal component. In the nomenclature of the forecast package, this corresponds to either the "AAA" model or the "AAM" model. The choice of method does not impact the level and trend components. Default is additive. The series y must be strictly positive when method='multiplicative' is chosen. Note that the multiplicative seasonality can be unstable when y is close to 0, especially for y < 10. You might want to switch to additive seasonality in such cases.
family	The distribution used to describe the error component; must be one or multiple of norm, student, or cauchy. The fitting procedure is slower the more are chosen, and fastest for norm. Default is c("norm", "student") as student is a reliable choice when the noise is somewhat fat-tailed, while cauchy can go too far.
param_grid	Matrix defining the grid of parameters to be trialled during grid search optimization; default parameter grid will be used if left NULL. Can be created using initialize_params_random(), initialize_params_naive(), initialize_params_grid(), or manually.
priors	List of priors on the models parameters; default priors will be used if left NULL. Can be created using add_prior_error(), add_prior_level(), add_prior_seasonality(), add_prior_trend(), add_prior_anomaly()
init_states	List of initial states l, b, and s used to start iterative smoothing of the time series for each set of the parameters in the parameter grid. Default initialization via initialize_states() will be used if left NULL. If s is provided but not both l and b, s will be used as initial seasonal component, and initialize_states() will initialize l and b given the user-provided s.
seasonality_threshold	During initialization of states, only use a seasonal component if more than seasonality_threshold share of the overall variation of the time series is explained by the seasonal component. This doesn't quite correspond to the residual variance as the variation reduction is calculated based on robust measures. Reasonable values are somewhere between 0 and 1.
remove_anomalies	Logical; during fitting, anomalies can be identified and interpolated to not adversely affect the fitted states. The interpolated values are only used to fit the states. When it comes to estimating the error's standard deviation sigma and to measuring the likelihood, interpolated values (i.e., fitted values) are compared against the original "anomalies" so that robust distributions like Cauchy and Student are correctly attributed a higher likelihood than the Normal distribution and may be chosen to project the uncertainty due to anomalies into the future. Default is TRUE.
anomaly_budget	Integerish (default 5); the number of anomalies that can be interpolated during fitting of state components. It can be useful to set a somewhat low hard limit on the number of possible interpolations as some parameter grid combinations may be misspecified compared to the "correct" data generating process and would therefore consider every observation an anomaly.
anomaly_budget_most_recent_k	Integerish (default 1); additional budget reserved to remove anomalies from the k most recent observations only. Especially models with larger alpha and beta smoothing parameters adjust the forecast strongly to the most recent observation(s). But even if the smoothing parameters are small, an anomalous most recent observation(s) can have a large impact on the forecast if it is sufficiently large (or small). At the same time, it can be undesired that the forecast adjusts heavily to the one, two, or k most recent observations even if they are very different from the rest of the series. For a discussion of this, see also: Michael Bohlke-Schneider, Shubham Kapoor, Tim Januschowski (2020). Resilient Neural Forecasting Systems. https://www.amazon.science/publications/resilient-neural-forecasting-systems
min_obs_anomaly_removal	Integerish (default 12); the anomaly detection relies on the fitted values' errors' standard deviation. The standard deviation is iteratively updated as the state components are fitted from the first observation to the last observation. This parameter defines after which observation of the time series there are sufficiently many observations available to reliably estimate the error's standard deviation and thus determine whether an observation should be considered an anomaly. The default of 12 is useful for monthly observations and not too low. But it also implies that any anomaly in the first 12 observations will have an impact on the estimated state components.
try_fixed_initial_fit	Logical (default FALSE); should the initial global fit provided via init_states / initialize_states() be evaluated alongside the smoothed model fits? This will evaluate the fit corresponding to init_states$l_global, init_states$b_global, and init_states$s, resulting in fitted values init_states$fitted_global. In contrast, the normal smoothed model fits rely on init_states$l, init_states$b, and init_states$s.
check_param_grid_unique	Logical (default TRUE); should tulip check that the provided param_grid rows are distinct (and remove duplicates otherwise)? Turn off to avoid computation if you know your param_grid has only distinct rows.

Value

An object of class tulip, a list with components:

y_hat

Fitted values, numeric vector of same length as y

y

Input time series

y_cleaned

Copy of the input time series that is used to update the state components during model fitting. The copy is used to achieve update behavior specific to the states when dealing with missing values and anomalies. In contrast to y, missing values and anomalies are replaced to update the states in robust ways. Note that the update behavior for sigma may differ, see also y_na.

n_cleaned

Number of cleaned observations, integer

y_na

Copy of the input time series that is used to update the sigma parameter during model fitting. The copy is used to achieve update behavior specific to sigma when dealing with missing values and anomalies. At the moment, anomalies are not cleaned in y_na, such that sigma will eventually adjust towards them. This behavior might change in the future.

param_grid

Set of smoothing parameters of fitted model

sigma

Fitted scale parameter of likelihood function

l

Fitted level state component

b

Fitted trend state component

s

Fitted seasonality state component

l_init

Initial level state component

b_init

Initial trend state component

s_init

Initial seasonality state component

log_joint

Value of the (log) joint distribution at the chosen parameter values

family

The distribution family of the fitted model

m

The suspected seasonality period length

method

Either additive or multiplicative seasonal component

priors

The list of priors used during the model estimation. This can deviate from the user-provided priors argument when the user did not provide a prior for all parameters. The returned list lists all priors that were effectively used.

comment

A character string; missing value in the normal case, else describing a predefined exception, e.g. in the case of a very short input series

full

List containing all fitted models for the full parameter grid

See Also

predict.tulip(), initialize_states(), initialize_params_grid(), add_prior_level(), add_prior_trend(), add_prior_seasonality(), add_prior_error(), add_prior_anomaly()

Examples

fitted_model <- tulip(y = tulip::flowers$flowers, m = 12)

print(fitted_model$family)
print(fitted_model$param_grid)

plot(tulip::flowers$flowers, type = "l", col = "grey", xlab = NA)
points(tulip::flowers$flowers, pch = 21, bg = "black", col = "white")

# add fitted values
lines(fitted_model$y_hat, col = "blue")

# indicate observations identified as anomalies and consequently interpolated
idx_anomalies <- which(fitted_model$y_cleaned != tulip::flowers$flowers)
points(
  x = idx_anomalies,
  y = fitted_model$y_cleaned[idx_anomalies],
  col = "darkorange", pch = 19
)
points(
  x = idx_anomalies,
  y = fitted_model$y_cleaned[idx_anomalies],
  col = "white", pch = 21
)

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