ves {legion}

Vector Exponential Smoothing in SSOE state space model

Description

Function constructs vector ETS model and returns forecast, fitted values, errors and matrix of states along with other useful variables.

Usage

ves(data, model = "PPP", lags = c(frequency(data)),
  persistence = c("common", "individual", "dependent"),
  transition = c("common", "individual", "dependent"), phi = c("common",
  "individual"), initial = c("individual", "common"),
  initialSeason = c("common", "individual"), loss = c("likelihood",
  "diagonal", "trace"), ic = c("AICc", "AIC", "BIC", "BICc"), h = 10,
  holdout = FALSE, occurrence = c("none", "fixed", "logistic"),
  bounds = c("admissible", "usual", "none"), silent = TRUE, ...)

Arguments

data

The matrix with the data, where series are in columns and observations are in rows.

model

The type of ETS model. Can consist of 3 or 4 chars: ANN, AAN, AAdN, AAA, AAdA, MMdM etc. PPP means that the best pure model will be selected based on the chosen information criteria type. ATTENTION! ONLY PURE ADDITIVE AND PURE MULTIPLICATIVE MODELS ARE AVAILABLE! Pure multiplicative models are done as additive model applied to log(y).

Also model can accept a previously estimated VES model and use all its parameters.

lags

The lags of the model. Needed for seasonal models.

persistence

Persistence matrix G, containing smoothing parameters. Can be:

• "independent" - each series has its own smoothing parameters and no interactions are modelled (all the other values in the matrix are set to zero);

• "dependent" - each series has its own smoothing parameters, but interactions between the series are modelled (the whole matrix is estimated);

• "group" each series has the same smoothing parameters for respective components (the values of smoothing parameters are repeated, all the other values in the matrix are set to zero).

• "seasonal" - each component has its own smoothing parameter, except for the seasonal one, which is common across the time series.

• provided by user as a vector or as a matrix. The value is used by the model.

You can also use the first letter instead of writing the full word.

transition

Transition matrix F. Can be:

• "independent" - each series has its own preset transition matrix and no interactions are modelled (all the other values in the matrix are set to zero);

• "dependent" - each series has its own transition matrix, but interactions between the series are modelled (the whole matrix is estimated). The estimated model behaves similar to VAR in this case;

• "group" each series has the same transition matrix for respective components (the values are repeated, all the other values in the matrix are set to zero).

• provided by user as a vector or as a matrix. The value is used by the model.

You can also use the first letter instead of writing the full word.

phi

In cases of damped trend this parameter defines whether the phi should be estimated separately for each series ("individual") or for the whole set ("common"). If vector or a value is provided here, then it is used by the model.

initial

Can be either character or a vector / matrix of initial states. If it is character, then it can be "individual", individual values of the initial non-seasonal components are used, or "common", meaning that the initials for all the time series are set to be equal to the same value. If vector of states is provided, then it is automatically transformed into a matrix, assuming that these values are provided for the whole group.

initialSeason

Can be either character or a vector / matrix of initial states. Treated the same way as initial. This means that different time series may share the same initial seasonal component.

loss

Type of Loss Function used in optimization. loss can be:

• "likelihood" - which implies the maximisation of likelihood of multivariate normal distribution (or log Normal if the multiplicative model is constructed);

• "diagonal" - similar to "likelihood", but assumes that covariances between the error terms are zero.

• "trace" - the trace of the covariance matrix of errors. The sum of variances is minimised in this case.

• Provided by user as a custom function of actual, fitted and B. Note that internally function transposes the data, so that actual and fitted contain observations in columns and series in rows.

An example of the latter option is: lossFunction <- function(actual,fitted,B){return(mean(abs(actual - fitted)))} followed by loss=lossFunction.

ic

The information criterion used in the model selection procedure.

h

Length of forecasting horizon.

holdout

If TRUE, holdout sample of size h is taken from the end of the data.

occurrence

Defines type of occurrence model used. Can be:

• none, meaning that the data should be considered as non-intermittent;

• fixed, taking into account constant Bernoulli distribution of demand occurrences;

• logistic, based on logistic regression.

In this case, the ETS model inside the occurrence part will correspond to model and probability="dependent". Alternatively, model estimated using oves function can be provided here.

bounds

What type of bounds to use in the model estimation. The first letter can be used instead of the whole word. "admissible" means that the model stability is ensured, while "usual" means that the all the parameters are restricted by the (0, 1) region.

silent

If silent="none", then nothing is silent, everything is printed out and drawn. silent="all" means that nothing is produced or drawn (except for warnings). In case of silent="graph", no graph is produced. If silent="legend", then legend of the graph is skipped. And finally silent="output" means that nothing is printed out in the console, but the graph is produced. silent also accepts TRUE and FALSE. In this case silent=TRUE is equivalent to silent="all", while silent=FALSE is equivalent to silent="none". The parameter also accepts first letter of words ("n", "a", "g", "l", "o").

...

Other non-documented parameters. For example FI=TRUE will make the function also produce Fisher Information matrix, which then can be used to calculated variances of smoothing parameters and initial states of the model. The vector of initial parameter for the optimiser can be provided here as the variable B. The upper bound for the optimiser is provided via ub, while the lower one is lb. Also, the options for nloptr can be passed here:

• maxeval=40*k is the default number of iterations for both optimisers used in the function (k is the number of parameters to estimate).

• algorithm1="NLOPT_LN_BOBYQA" is the algorithm used in the first optimiser, while algorithm2="NLOPT_LN_NELDERMEAD" is the second one.

• xtol_rel1=1e-8 is the relative tolerance in the first optimiser, while xtol_rel2=1e-6 is for the second one. All of this can be amended and passed in ellipsis for finer tuning.

• print_level - the level of output for the optimiser (0 by default). If equal to 41, then the detailed results of the optimisation are returned.

Details

Function estimates vector ETS in a form of the Single Source of Error state space model of the following type:

\mathbf{y}_{t} = (\mathbf{W} \mathbf{v}_{t-l} + \mathbf{x}_t \mathbf{a}_{t-1} + \mathbf{\epsilon}_{t})

\mathbf{v}_{t} = \mathbf{F} \mathbf{v}_{t-l} + \mathbf{G} \mathbf{\epsilon}_{t}

\mathbf{a}_{t} = \mathbf{F_{X}} \mathbf{a}_{t-1} + \mathbf{G_{X}} \mathbf{\epsilon}_{t} / \mathbf{x}_{t}

Where y_{t} is the vector of time series on observation t, \mathbf{v}_{t} is the matrix of states and l is the matrix of lags, \mathbf{x}_t is the vector of exogenous variables. \mathbf{W} is the measurement matrix, \mathbf{F} is the transition matrix and \mathbf{G} is the persistence matrix. Finally, \epsilon_{t} is the vector of error terms.

Conventionally we formulate values as:

\mathbf{y}'_t = (y_{1,t}, y_{2,t}, \dots, y_{m,t})

where m is the number of series in the group.

\mathbf{v}'_t = (v_{1,t}, v_{2,t}, \dots, v_{m,t})

where v_{i,t} is vector of components for i-th time series.

\mathbf{W}' = (w_{1}, \dots , 0; \vdots , \ddots , \vdots; 0 , \vdots , w_{m})

is matrix of measurement vectors.

For the details on the additive model see Hyndman et al. (2008), chapter 17.

In case of multiplicative model, instead of the vector y_t we use its logarithms. As a result the multiplicative model is much easier to work with.

For some more information about the model and its implementation, see the vignette: vignette("ves","legion")

Value

Object of class "legion" is returned. It contains the following list of values:

• model - The name of the fitted model;

• timeElapsed - The time elapsed for the construction of the model;

• states - The matrix of states with components in columns and time in rows;

• persistence - The persistence matrix;

• transition - The transition matrix;

• measurement - The measurement matrix;

• phi - The damping parameter value;

• lagsAll - The vector of the internal lags used in the model;

• B - The vector of all the estimated coefficients;

• initial - The initial values of the non-seasonal components;

• initialSeason - The initial values of the seasonal components;

• nParam - The number of estimated parameters;

• occurrence - The occurrence part of the model estimated with VES;

• data - The matrix with the original data;

• fitted - The matrix of the fitted values;

• holdout - The matrix with the holdout values (if holdout=TRUE in the estimation);

• residuals - The matrix of the residuals of the model;

• Sigma - The covariance matrix of the errors (estimated with the correction for the number of degrees of freedom);

• forecast - The matrix of point forecasts;

• ICs - The values of the information criteria;

• logLik - The log-likelihood function;

• lossValue - The value of the loss function;

• loss - The type of the used loss function;

• lossFunction - The loss function if the custom was used in the process;

• accuracy - the values of the error measures. Currently not available.

• FI - Fisher information if user asked for it using FI=TRUE.

Author(s)

Ivan Svetunkov, [e-mail]

References

• de Silva A., Hyndman R.J. and Snyder, R.D. (2010). The vector innovations structural time series framework: a simple approach to multivariate forecasting. Statistical Modelling, 10 (4), pp.353-374

• Hyndman, R.J., Koehler, A.B., Ord, J.K., and Snyder, R.D. (2008) Forecasting with exponential smoothing: the state space approach, Springer-Verlag.

• Lütkepohl, H. (2005). New Introduction to Multiple Time Series Analysis. New introduction to Multiple Time Series Analysis. Berlin, Heidelberg: Springer Berlin Heidelberg. doi:10.1007/978-3-540-27752-1

• Chen H., Svetunkov I., Boylan J. (2021). A New Taxonomy for Vector Exponential Smoothing and Its Application to Seasonal Time Series.

See Also

vets, es, adam

Examples

Y <- ts(cbind(rnorm(100,100,10),rnorm(100,75,8)),frequency=12)

# The simplest model applied to the data with the default values
ves(Y,model="ANN",h=10,holdout=TRUE)

# Damped trend model with the dependent persistence
ves(Y,model="AAdN",persistence="d",h=10,holdout=TRUE)

# Multiplicative damped trend model with individual phi
ves(Y,model="MMdM",persistence="i",h=10,holdout=TRUE,initialSeason="c")

# Automatic selection between pure models
ves(Y,model="PPP",persistence="i",h=10,holdout=TRUE,initialSeason="c")

# Intermittent demand vector model
Y <- cbind(c(rpois(25,0.1),rpois(25,0.5),rpois(25,1),rpois(25,5)),
           c(rpois(25,0.1),rpois(25,0.5),rpois(25,1),rpois(25,5)))

ves(Y,model="MNN",h=10,holdout=TRUE,occurrence="l")