thief_vets: Temporal hierarchy reconciliation of vector exponential...

In nicholasjclark/mvforecast: Tools to fit, interrogate and forecast multivariate statistical and machine learning timeseries models

Description

This function fits vector exponential smoothing on the base level and optionally on higher levels of temporal aggregation for a multivariate xts timeseries object

Usage

thief_vets(
  y,
  multi_freq = 12,
  k = 1,
  level = "grouped",
  slope = "none",
  damped = "none",
  seasonal = "common",
  lambda = NULL,
  dependence = "equicorrelation",
  frequency = 52,
  horizon = NULL,
  cores = parallel::detectCores() - 1,
  group = NULL,
  xreg = NULL,
  max_agg = NULL,
  xreg_include = NULL,
  newxreg = NULL,
  save_plots = FALSE,
  fig_path = "",
  discrete = FALSE
)

Arguments

y	xts matrix. The outcome series to be modelled. NAs are currently not supported
multi_freq	integer. Minimum frequency in the temporal hierarchy that will be modelled using a vector exponential smoothing model. Aggregates with frequencies below this threshold will be modelled with univariate models using the forecast function
k	integer specifying the length of the forecast horizon in multiples of frequency. Default is 1, meaning that a final forecast of frequency horizons will be returned
level	character. The dynamics for the level component in the tsvets model. Options are: c("constant", "diagonal", "common", "full", "grouped")
slope	character. The dynamics for the slope component in the tsvets model. Options are: c("none", "constant", "common", "diagonal", "full", "grouped")
damped	character. The dynamics for the dampening component in the tsvets model. Options are: c("none", "common", "diagonal", "full", "grouped")
seasonal	character. The dynamics for the seasonal component in the tsvets model. Options are: c("none", "common", "diagonal", "full", "grouped")
lambda	numeric proportional. The multivariate Box Cox power transformation parameter for all series. Must be between 0 and 1.5 inclusive
dependence	character. The multivariate error dependence structure to impose. Options are: c("diagonal", "full", "equicorrelation", "shrinkage")
frequency	integer. The seasonal frequency in y
horizon	integer. The horizon to forecast. Defaults to frequency
cores	integer. The number of cores to use. This is used to initialize the states of each series using ets_modelspec
group	Optional vector of indices denoting which group the series belongs to (when using the grouped dynamics). Defaults to NULL
xreg	Optional xts matrix of external regressors. Defaults to NULL
max_agg	(optional) integer specifying the maximum number of temporal aggregation levels to use when reconciling, via the structural scaling method. Useful if higher levels of aggregation are unlikely to have 'seen' recent changes in series dynamics and will likely then result in poor forecasts as a result. Default is NULL, meaning that all levels of aggregation are used
xreg_include	optional matrix of dimension ncol(y) by ncol(xreg) populated with either 0, 1 or 2+ (0 = no beta, 1 = individual beta and 2 = grouped beta). It is also possible to have group wise pooling. For instance 2 variables sharing one pooled estimates, and 3 other variables sharing another grouped estimate would have values of (2,2,3,3,3). The index for group wise pooling starts at 2 and should be incremented for each new group added. Defaults to NULL
newxreg	Optional xts matrix of future values for external regressors to be used in forecasting. Defaults to NULL
save_plots	Logical. Plots of fitted and residual values for the unaggregated series will be saved to fig_path if TRUE
fig_path	character. Optional filepath where fitted and residual plots will be saved. Defaults to the current working directory
discrete	logical Is the series in y discrete? If TRUE, use a copula-based method relying on the Probability Integral Transform to map the series to an approximate Gaussian distribution prior to modelling. Forecasts are then back-transformed to the estimated discrete distribution that best fits y. Default is FALSE

Details

Series in y are aggregated at all possible levels up to annual using tsaggregates. estimate.tsvets.spec is used on the unaggregated (original series), and optionally on higher levels of aggregation down to a frequency of multi_freq. At frequencies below multi_freq, the best-fitting univariate model is chosen using automatic ensembles in ensemble_base. Forecasts are reconciled using reconcilethief and are optionally constrained using non-negative optimisation if lambda is provided. Adjustments to the original unaggregated forecast are incorporated and a distribution of 1000 sample paths for each series' forecast are returned

Value

A list containing the reconciled forecast distributions for each series in y. Each element in the list is a horizon x 1000 matrix of forecast predictions

References

Athanasopoulos, G and de Silva, A. (2012),Multivariate Exponential Smoothing for Forecasting Tourist Arrivals, Journal of Travel Research 51(5) 640–652.

de Silva, A., R. Hyndman, and R. D. Snyder. (2010).The Vector Innovations Structural Time Series Framework: A Simple Approach to Multivariate Forecasting, Statistical Modelling (10) 353–74.

Athanasopoulos, G., Hyndman, R., Kourentzes, N., and Petropoulos, F. Forecasting with temporal hierarchies. (2017) European Journal of Operational Research 262(1) 60–74

See Also

vets_modelspec, ets_modelspec, ensemble_base, reconcilethief

Examples

library(mvforecast)
data("ixodes_vets_dat")

# View the returned data
head(ixodes_vets_dat$y_train)
head(ixodes_vets_dat$xreg_train)
ixodes_vets_dat$xreg_include
head(ixodes_vets_dat$future_xreg)

# Fit a vets model with no regressors and common seasonality with the tsvets package
mod1 <- tsvets:::simulate.tsvets.estimate(tsvets:::estimate.tsvets.spec(tsvets:::vets_modelspec(ixodes_vets_dat$y_train,
level = "grouped",
slope = "none",
damped = "none",
seasonal = "common",
lambda = 1,
dependence = "equicorrelation",
frequency = 52,
cores = parallel::detectCores() - 1,
group = ixodes_vets_dat$groups),
solver = "solnp",
control = list(trace = 0)),
nsim = 1000,
h = ixodes_vets_dat$h)

# Calculate a summary of the summed out-of-sample CRPS
calc_crps(simulation = mod1, y_test = ixodes_vets_dat$y_test)

# Explore whether reconciliation of temporal hierarchies improves predictions
mod2 <- thief_vets(y = ixodes_vets_dat$y_train,
multi_freq = 12,
level = "grouped",
slope = "none",
damped = "none",
seasonal = "common",
lambda = 1,
dependence = "equicorrelation",
frequency = 52,
cores = parallel::detectCores() - 1,
group = ixodes_vets_dat$groups,
save_plots = FALSE)

calc_crps(simulation = mod2, y_test = ixodes_vets_dat$y_test)

# Plot one of the forecasts against the true values in the test set
plot_mvforecast(simulation = mod2[[4]])
points(as.vector(ixodes_vets_dat$y_test[,4]))

nicholasjclark/mvforecast documentation built on Dec. 22, 2021, 2:11 a.m.