oes() - occurrence part of iETS model

Ivan Svetunkov

2019-07-24

smooth package has a mechanism of treating the data with zero values. This might be useful for cases of intermittent demand (the demand that happens at random). All the univariate functions in the package have a parameter occurrence that allows handling this type of data. The canonical model, used in smooth, is called “iETS” - intermittent exponential smoothing model. This vignette explains how the iETS model and its occurrence part are implemented in the smooth package.

The basics

The canonical general iETS model (called iETS\(_G\)) can be summarised as: \[\begin{equation} \label{eq:iETS} \tag{1} \begin{matrix} y_t = o_t z_t \\ o_t \sim \text{Bernoulli} \left(p_t \right) \\ p_t = f(a_t, b_t) \\ a_t = w_a(v_{a,t-L}) + r_a(v_{a,t-L}) \epsilon_{a,t} \\ v_{a,t} = f_a(v_{a,t-L}) + g_a(v_{a,t-L}) \epsilon_{a,t} \\ (1 + \epsilon_{a,t}) \sim \text{log}\mathcal{N}(0, \sigma_{a}^2) \\ b_t = w_a(v_{b,t-L}) + r_a(v_{b,t-L}) \epsilon_{b,t} \\ v_{b,t} = f_a(v_{b,t-L}) + g_a(v_{b,t-L}) \epsilon_{b,t} \\ (1 + \epsilon_{b,t}) \sim \text{log}\mathcal{N}(0, \sigma_{b}^2) \end{matrix}, \end{equation}\] where \(y_t\) is the observed values, \(z_t\) is the demand size, which is a pure multiplicative ETS model on its own, \(w(\cdot)\) is the measurement function, \(r(\cdot)\) is the error function, \(f(\cdot)\) is the transition function and \(g(\cdot)\) is the persistence function (the subscripts allow separating the functions for different parts of the model). These four functions define how the elements of the vector \(v_{t}\) interact with each other. Furthermore, \(\epsilon_{a,t}\) and \(\epsilon_{b,t}\) are the mutually independent error terms, \(o_t\) is the binary occurrence variable (1 - demand is non-zero, 0 - no demand in the period \(t\)) which is distributed according to Bernoulli with probability \(p_t\) that has a logit-normal distribution (\(o_t \sim \text{Bernoulli} \left(p_t \right)\), \(p_t \sim \text{logit} \mathcal{N}\)). Any ETS model can be used for \(a_t\) and \(b_t\), and the transformation of them into the probability \(p_t\) depends on the type of the error. The general formula for the multiplicative error is: \[\begin{equation} \label{eq:oETS(MZZ)} p_t = \frac{a_t}{a_t+b_t} , \end{equation}\] while for the additive error it is: \[\begin{equation} \label{eq:oETS(AZZ)} p_t = \frac{\exp(a_t)}{\exp(a_t)+\exp(b_t)} . \end{equation}\] This is because both \(a_t\) and \(b_t\) need to be positive, and the additive error models support the real plane. The canonical iETS model assumes that the pure multiplicative model is used for the both \(a_t\) and \(b_t\). This type of model is positively defined for any values of error, trend and seasonality, which is essential for the values of \(a_t\) and \(b_t\). If a combination of additive and multiplicative error models is used, then the additive part is exponentiated prior to the usage of the formulae for the calculation of the probability.

An example of an iETS model is the basic local-level model iETS(M,N,N)\(_G\)(M,N,N)(M,N,N): \[\begin{equation} \label{eq:iETSGExample} \begin{matrix} y_t = o_t z_t \\ z_t = l_{z,t-1} \left(1 + \epsilon_{z,t} \right) \\ l_{z,t} = l_{z,t-1}( 1 + \alpha_{z} \epsilon_{z,t}) \\ (1 + \epsilon_{t}) \sim \text{log}\mathcal{N}(0, \sigma_\epsilon^2) \\ \\ o_t \sim \text{Bernoulli} \left(p_t \right) \\ p_t = \frac{a_t}{a_t+b_t} \\ a_t = l_{a,t-1} \left(1 + \epsilon_{a,t} \right) \\ l_{a,t} = l_{a,t-1}( 1 + \alpha_{a} \epsilon_{a,t}) \\ (1 + \epsilon_{a,t}) \sim \text{log}\mathcal{N}(0, \sigma_{a}^2) \\ b_t = l_{b,t-1} \left(1 + \epsilon_{b,t} \right) \\ l_{b,t} = l_{b,t-1}( 1 + \alpha_{b} \epsilon_{b,t}) \\ (1 + \epsilon_{b,t}) \sim \text{log}\mathcal{N}(0, \sigma_{b}^2) \end{matrix}, \end{equation}\] where \(l_{a,t}\) and \(l_{b,t}\) are the levels for each of the shape parameters and \(\alpha_{a}\) and \(\alpha_{b}\) are the smoothing parameters. More advanced models can be constructing by specifying the ETS models for each part and / or adding explanatory variables.

In the notation of the model iETS(M,N,N)\(_G\)(M,N,N)(M,N,N), the first brackets describe the ETS model for the demand sizes, the underscore letter points out at the specific subtype of model (see below), the second brackets describe the ETS model, underlying the variable \(a_t\) and the last ones stand for the model for the \(b_t\). If only one variable is needed (either \(a_t\) or \(b_t\)), then the redundant brackets are dropped, so that the notation simplifies, for example, to: iETS(M,N,N)\(_O\)(M,N,N). If the same type of model is used for both demand sizes and demand occurrence, then the second brackets can be dropped as well, simplifying the view to: iETS(M,N,N)\(_G\). Furthermore, the notation without any brackets, such as iETS\(_G\) stands for a general class of a specific subtype of iETS model (so any error / trend / seasonality). Also, given that iETS\(_G\) is the most general model of all iETS models, the “\(G\)” can be dropped, when the properties are applicable to all subtypes. Finally, the “oETS” notation is used when the occurrence part of the model is discussed explicitly, skipping the demand sizes.

The concentrated likelihood function for the iETS model is: \[\begin{equation} \label{eq:LogNormalConcentratedLogLikelihood} \tag{2} \ell(\boldsymbol{\theta}, \hat{\sigma}_\epsilon^2 | \textbf{Y}) = - \frac{1}{2} \left( T \log(2 \pi e \hat{\sigma}_\epsilon^2) + T_0 \right) - {\sum_{o_t=1}} \log(z_t) + {\sum_{o_t=1}} \log(\hat{p}_t) + {\sum_{o_t=0}} \log(1-\hat{p}_t) , \end{equation}\] where \(\textbf{Y}\) is the vector of all the in-sample observations, \(\boldsymbol{\theta}\) is the vector of parameters to estimate (initial values and smoothing parameters), \(T\) is the number of all observations, \(T_0\) is the number of zero observations, \(\hat{\sigma_\epsilon}^2 = \frac{1}{T} \sum_{o_t=1} \log^2 \left(1 + \epsilon_t \right)\) is the scale parameter of the one-step-ahead forecast error for the demand sizes and \(\hat{p}_t\) is the estimated probability of a non-zero demand at time \(t\). This likelihood is used for the estimation of all the special cases of the iETS\(_G\) model.

Depending on the restrictions on \(a_t\) and \(b_t\), there can be several iETS models:

  1. iETS\(_F\): \(a_t = \text{const}\), \(b_t = \text{const}\) - the model with the fixed probability of occurrence;
  2. iETS\(_O\): \(b_t = 1\) - the “odds ratio” model, based on the logistic transform of the \(o_t\) variable;
  3. iETS\(_I\): \(a_t = 1\) - the “inverse odds ratio” model, based on the inverse logistic transform of the \(o_t\) variable;
  4. iETS\(_D\): \(a_t + b_t = 1\), \(a_t \leq 1\) - the direct probability model, where the \(p_t\) is calculated directly from the occurrence variable \(o_t\);
  5. iETS\(_G\): No restrictions - the model based on the evolution of both \(a_t\) and \(b_t\).

Depending on the type of the model, there are different mechanisms of the model construction, error calculation, update of the states and the generation of forecasts. The one thing uniting all the subtypes of models is that all the multiplicative error ETS models in smooth package rely on the conditional medians rather than on the conditional means. This is caused by the assumption of log normality of the error term and simplifies some of the calculations, preserving the main idea of ETS models (separation of components of the model).

In this vignette we will use ETS(M,N,N) model as a base for the different parts of the models. Although, this is a simplification, it allows better understanding the basics of the different types of iETS model, without the loss of generality.

We will use an artificial data in order to see how the functions work:

y <- ts(c(rpois(20,0.25),rpois(20,0.5),rpois(20,1),rpois(20,2),rpois(20,3),rpois(20,5)))

All the models, discussed in this vignette, are implemented in the functions oes() and oesg(). The only missing element in all of this at the moment is the model selection mechanism for the demand occurrence part. So neither oes() nor oesg() currently support “ZZZ” ETS models.

iETS\(_F\)

In case of the fixed \(a_t\) and \(b_t\), the iETS\(_G\) model reduces to: \[\begin{equation} \label{eq:ISSETS(MNN)Fixed} \tag{3} \begin{matrix} y_t = o_t z_t \\ o_t \sim \text{Bernoulli}(p) \end{matrix} . \end{equation}\]

The conditional h-steps ahead mean of the demand occurrence probability is calculated as: \[\begin{equation} \label{eq:pt_fixed_expectation} \mu_{o,t+h|t} = \tilde{p}_{t+h|t} = \hat{p} . \end{equation}\]

The estimate of the probability \(p\) is calculated based on the maximisation of the following concentrated log-likelihood function: \[\begin{equation} \label{eq:ISSETS(MNN)FixedLikelihoodSmooth} \ell \left({p} | o_t \right) = T_1 \log {p} + T_0 \log (1-{p}) , \end{equation}\] where \(T_0\) is the number of zero observations and \(T_1\) is the number of non-zero observations in the data. The number of estimated parameters in this case is equal to \(k_z+1\), where \(k_z\) is the number of parameters for the demand sizes part, and 1 is for the estimation of the probability \(p\). Maximising this likelihood deems the analytical solution for the \(p\): \[\begin{equation} \label{eq:ISSETS(MNN)FixedLikelihoodSmoothProbability} \hat{p} = \frac{1}{T} \sum_{t=1}^T o_t , \end{equation}\] where \(T\) is the number of all the available observations.

The occurrence part of the model oETS\(_F\) is constructed using oes() function:

oETSFModel1 <- oes(y, occurrence="fixed", h=10, holdout=TRUE)
oETSFModel1
## Occurrence state space model estimated: Fixed probability
## Underlying ETS model: oETS[F](MNN)
## Smoothing parameters:
## level 
##     0 
## Vector of initials:
##  level 
## 0.6455 
## 
## Error standard deviation: 1.0909
## Sample size: 110
## Number of estimated parameters: 1
## Number of degrees of freedom: 109
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 145.0473 145.0844 147.7478 147.8349
plot(oETSFModel1)

All the smooth forecasting functions support the occurrence part of the model. For example, here’s how the iETS(M,M,N)\(_F\) can be constructed:

es(y, "MMN", occurrence="fixed", h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.1 seconds
## Model estimated: iETS(MMN)
## Occurrence model type: Fixed probability
## Persistence vector g:
##  alpha   beta 
## 0.0028 0.0028 
## Initial values were optimised.
## 
## Loss function type: MSE; Loss function value: 0.1636
## Error standard deviation: 0.4045
## Sample size: 110
## Number of estimated parameters: 6
## Number of degrees of freedom: 104
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 412.7310 413.3079 428.9339 425.5893 
## 
## Forecast errors:
## Bias: 91.3%; sMSE: 112.5%; RelRMSE: 1.7; sPIS: -4995.2%; sCE: -851.8%

iETS\(_O\)

The odds-ratio iETS uses only one model for the occurrence part, for the \(a_t\) variable (setting \(b_t=1\)), which simplifies the iETS\(_G\) model. For example, for the iETS\(_O\)(M,N,N): \[\begin{equation} \label{eq:iETSO} \tag{5} \begin{matrix} y_t = o_t z_t \\ o_t \sim \text{Bernoulli} \left(p_t \right) \\ p_t = \frac{a_t}{a_t+1} \\ a_t = l_{a,t-1} \left(1 + \epsilon_{a,t} \right) \\ l_{a,t} = l_{a,t-1}( 1 + \alpha_{a} \epsilon_{a,t}) \\ (1 + \epsilon_{a,t}) \sim \text{log}\mathcal{N}(0, \sigma_{a}^2) \end{matrix}. \end{equation}\]

In the estimation of the model, the initial level is set to the transformed mean probability of occurrence \(l_{a,0}=\frac{\bar{p}}{1-\bar{p}}\) for multiplicative error model and \(l_{a,0} = \log l_{a,0}\) for the additive one, where \(\bar{p}=\frac{1}{T} \sum_{t=1}^T o_t\), the initial trend is equal to 0 in case of the additive and 1 in case of the multiplicative types. In cases of seasonal models, the regression with dummy variables is fitted, and its parameters are then used for the initials of the seasonal indices after the transformations similar to the level ones.

The construction of the model is done via the following set of equations (example with oETS\(_O\)(M,N,N)): \[\begin{equation} \label{eq:iETSOEstimation} \begin{matrix} \hat{p}_t = \frac{\hat{a}_t}{\hat{a}_t+1} \\ \hat{a}_t = l_{a,t-1} \\ l_{a,t} = l_{a,t-1}( 1 + \alpha_{a} e_{a,t}) \\ 1+e_{a,t} = \frac{u_t}{1-u_t} \\ u_{t} = \frac{1 + o_t - \hat{p}_t}{2} \end{matrix}. \end{equation}\]

Given that the model is estimated using the likelihood (2), it has \(k_z+k_a\) parameters to estimate, where \(k_z\) includes all the initial values, the smoothing parameters and the scale of the error of the demand sizes part of the model, and \(k_a\) includes only initial values and the smoothing parameters of the model for the demand occurrence. In case of iETS\(_O\)(M,N,N) this number is equal to 5.

The occurrence part of the model iETS\(_O\) is constructed using the very same oes() function, but also allows specifying the ETS model to use. Foe example, here’s the ETS(M,M,N) model:

oETSOModel <- oes(y, model="MMN", occurrence="o", h=10, holdout=TRUE)
oETSOModel
## Occurrence state space model estimated: Odds ratio
## Underlying ETS model: oETS[O](MMN)
## Smoothing parameters:
##  level  trend 
## 0.0230 0.0021 
## Vector of initials:
##  level  trend 
## 0.5146 0.9486 
## 
## Error standard deviation: 1.6063
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 114.2799 114.6609 125.0819 125.9772
plot(oETSOModel)

And here’s the full iETS(M,M,N)\(_O\) model:

es(y, "MMN", occurrence="o", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.17 seconds
## Model estimated: iETS(MMN)
## Occurrence model type: Odds ratio
## Persistence vector g:
##  alpha   beta 
## 0.0028 0.0028 
## Initial values were optimised.
## 
## Loss function type: MSE; Loss function value: 0.1636
## Error standard deviation: 0.4045
## Sample size: 110
## Number of estimated parameters: 9
## Number of degrees of freedom: 101
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 381.9636 382.5405 406.2679 388.8219 
## 
## Forecast errors:
## Bias: 21%; sMSE: 43.1%; RelRMSE: 1.053; sPIS: -947.8%; sCE: -70.6%

This should give the same results as running, meaning that we ask explicitly for the es() function to use the earlier estimated model:

es(y, "MMN", occurrence=oETSOModel, h=10, holdout=TRUE, silent=FALSE)

This gives an additional flexibility, because the construction can be done in two steps, with a more refined model for the occurrence part (e.g. including explanatory variables).

iETS\(_I\)

Similarly to the odds-ratio iETS, inverse-odds-ratio model uses only one model for the occurrence part, but for the \(b_t\) variable instead of \(a_t\) (now \(a_t=1\)). Here is an example of iETS\(_I\)(M,N,N): \[\begin{equation} \label{eq:iETSI} \tag{6} \begin{matrix} y_t = o_t z_t \\ o_t \sim \text{Bernoulli} \left(p_t \right) \\ p_t = \frac{1}{1+b_t} \\ b_t = l_{b,t-1} \left(1 + \epsilon_{b,t} \right) \\ l_{b,t} = l_{b,t-1}( 1 + \alpha_{b} \epsilon_{b,t}) \\ (1 + \epsilon_{b,t}) \sim \text{log}\mathcal{N}(0, \sigma_{b}^2) \end{matrix}. \end{equation}\]

In the estimation of the model, the initial level is set to the transformed mean probability of occurrence \(l_{b,0}=\frac{1-\bar{p}}{\bar{p}}\) for multiplicative error model and \(l_{b,0} = \log l_{b,0}\) for the additive one, where \(\bar{p}=\frac{1}{T} \sum_{t=1}^T o_t\), the initial trend is equal to 0 in case of the additive and 1 in case of the multiplicative types. The seasonality is treated similar to the iETS\(_O\) model, but using the inverse-odds transformation.

The construction of the model is done via the set of equations similar to the ones for the iETS\(_O\) model: \[\begin{equation} \label{eq:iETSIEstimation} \begin{matrix} \hat{p}_t = \frac{\hat{b}_t}{\hat{b}_t+1} \\ \hat{b}_t = l_{b,t-1} \\ l_{b,t} = l_{b,t-1}( 1 + \alpha_{b} e_{b,t}) \\ 1+e_{b,t} = \frac{1-u_t}{u_t} \\ u_{t} = \frac{1 + o_t - \hat{p}_t}{2} \end{matrix}. \end{equation}\]

So the model iETS\(_I\) is like a mirror reflection of the model iETS\(_O\). However, it produces different forecasts, because it focuses on the probability of non-occurrence, rather than the probability of occurrence. Interestingly enough, the probability of occurrence \(p_t\) can also be estimated if \(1+b_t\) in the denominator is set to be equal to the demand intervals (between the demand occurrences). The model (6) underlies Croston’s method in this case.

Once again oes() function is used in the construction of the model:

oETSIModel <- oes(y, model="MMN", occurrence="i", h=10, holdout=TRUE)
oETSIModel
## Occurrence state space model estimated: Inverse odds ratio
## Underlying ETS model: oETS[I](MMN)
## Smoothing parameters:
##  level  trend 
## 0.0519 0.0000 
## Vector of initials:
##  level  trend 
## 8.3670 0.9139 
## 
## Error standard deviation: 7.5936
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 118.0928 118.4738 128.8947 129.7901
plot(oETSIModel)

And here’s the full iETS(M,M,N)\(_O\) model:

es(y, "MMN", occurrence="i", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.18 seconds
## Model estimated: iETS(MMN)
## Occurrence model type: Inverse odds ratio
## Persistence vector g:
##  alpha   beta 
## 0.0028 0.0028 
## Initial values were optimised.
## 
## Loss function type: MSE; Loss function value: 0.1636
## Error standard deviation: 0.4045
## Sample size: 110
## Number of estimated parameters: 9
## Number of degrees of freedom: 101
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 385.7765 386.3534 410.0808 392.6348 
## 
## Forecast errors:
## Bias: 33.9%; sMSE: 44.3%; RelRMSE: 1.067; sPIS: -1217.3%; sCE: -115.4%

Once again, an earlier estimated model can be used in the univariate forecasting functions:

es(y, "MMN", occurrence=oETSIModel, h=10, holdout=TRUE, silent=FALSE)

iETS\(_D\)

This model appears, when a specific restriction is imposed: \[\begin{equation} \label{eq:iETSGRestriction} a_t + b_t = 1, a_t \in [0, 1] \end{equation}\] The pure multiplicative iETS model is then transformed into: \[\begin{equation} \label{eq:iETSD} \tag{7} \begin{matrix} y_t = o_t z_t \\ o_t \sim \text{Bernoulli} \left(a_t \right) \\ a_t = \min \left(l_{a,t-1} \left(1 + \epsilon_{a,t} \right) , 1 \right) \\ l_{a,t} = l_{a,t-1}( 1 + \alpha_{a} \epsilon_{a,t}) \\ (1 + \epsilon_{a,t}) \sim \text{log}\mathcal{N}(0, \sigma_{a}^2) \end{matrix}. \end{equation}\] An option with the additive model in this case has a different, more complicated form: \[\begin{equation} \label{eq:iETSDAdditive} \begin{matrix} y_t = o_t z_t \\ o_t \sim \text{Bernoulli} \left(a_t \right) \\ a_t = \max \left( \min \left(l_{a,t-1} + \epsilon_{a,t} , 1 \right), 0 \right) \\ l_{a,t} = l_{a,t-1} + \alpha_{a} \epsilon_{a,t} \\ \epsilon_{a,t} \sim \mathcal{N}(0, \sigma_{a}^2) \end{matrix}. \end{equation}\] The conditional expectation and variance for the demand occurrence part in this case is quite nasty, although it has an analytical solution. In case of the additive model it might cause problems, so the implementation in smooth is not accurate and needs to be fixed. However, in case of multiplicative model, the medians are used, which are much simpler, and, as a result, the point forecasts for several steps ahead correspond to the ones from the ETS models, restricted with the [0, 1] region. So, I would recommend either sticking with the multiplicative error models or using some other iETS subtypes.

The estimation of the multiplicative error model is done using the following set of equations: \[\begin{equation} \label{eq:ISSETS(MNN)_probability_estimate} \begin{matrix} \hat{y}_t = o_t \hat{l}_{z,t-1} \\ \hat{l}_{z,t} = \hat{l}_{z,t-1}( 1 + \alpha e_t) \\ \hat{a}_t = min(\hat{l}_{a,t-1}, 1) \\ \hat{l}_{a,t} = \hat{l}_{a,t-1}( 1 + \alpha_{a} e_{a,t}) \end{matrix}, \end{equation}\] where \[\begin{equation} \label{eq:ISSETS(MNN)_TSB_model_error_approximation} e_{a,t} = \frac{o_t (1 - 2 \kappa) + \kappa - \hat{a}_t}{\hat{a}_t}, \end{equation}\] and \(\kappa\) is a very small number (for example, \(\kappa = 10^{-10}\)), needed only in order to make the model estimable. The estimate of the error term in case of the additive model is much simpler and does not need any specific tricks to work: \[\begin{equation} \label{eq:ISSETS(MNN)_TSB_model_error_approximation2} e_{a,t} = o_t - \hat{a}_t . \end{equation}\]

The initials of the iETS\(_D\) model are calculated directly from the data without any additional transformations

Here’s an example of the application of the model to the same artificial data:

oETSDModel <- oes(y, model="MMN", occurrence="d", h=10, holdout=TRUE)
oETSDModel
## Occurrence state space model estimated: Direct probability
## Underlying ETS model: oETS[D](MMN)
## Smoothing parameters:
##  level  trend 
## 0.0823 0.0000 
## Vector of initials:
## level trend 
## 0.320 1.001 
## 
## Error standard deviation: 1.1514
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 122.7514 123.1324 133.5534 134.4487
plot(oETSDModel)

The usage of the model in case of univariate forecasting functions is the same as in the cases of other occurrence models, discussed above:

es(y, "MMN", occurrence=oETSDModel, h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.1 seconds
## Model estimated: iETS(MMN)
## Occurrence model type: Direct
## Persistence vector g:
##  alpha   beta 
## 0.0028 0.0028 
## Initial values were optimised.
## 
## Loss function type: MSE; Loss function value: 0.1636
## Error standard deviation: 0.4045
## Sample size: 110
## Number of estimated parameters: 5
## Number of provided parameters: 4
## Number of degrees of freedom: 105
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 382.4351 383.0120 395.9375 397.2934 
## 
## Forecast errors:
## Bias: 37.1%; sMSE: 44.2%; RelRMSE: 1.065; sPIS: -1264.1%; sCE: -130.3%

iETS\(_G\)

This model has already been discussed above and was presented in (1). The estimation of iETS(M,N,N)\(_G\) model is done via the following set of equations: \[\begin{equation} \label{eq:ISSETS(MNN)Estimated} \begin{matrix} \hat{y}_t = o_t \hat{z}_t \\ e_t = o_t \frac{y_t - \hat{z}_t}{\hat{z}_t} \\ \hat{z}_t = \hat{l}_{z,t-1} \\ \hat{l}_{z,t} = \hat{l}_{z,t-1}( 1 + \alpha_z e_t) \\ e_{a,t} = \frac{u_t}{1-u_t} -1 \\ \hat{a}_t = \hat{l}_{a,t-1} \\ \hat{l}_{a,t} = \hat{l}_{a,t-1}( 1 + \alpha_{a} e_{a,t}) \\ \hat{b}_t = \hat{l}_{b,t-1} \\ \hat{l}_{b,t} = \hat{l}_{b,t-1}( 1 + \alpha_{b} e_{b,t}) \\ e_{b,t} = \frac{1-u_t}{u_t} -1 \end{matrix} . \end{equation}\] The initialisation of the parameters of the iETS\(_G\) model is done separately for the variables \(a_t\) and \(b_t\), based on the principles, described above for the iETS\(_O\) and iETS\(_I\).

There is a separate function for this model, called oesg(). It has twice more parameters than oes(), because it allows fine tuning of the models for the variables \(a_t\) and \(b_t\). This gives an additional flexibility. For example, here is how we can use ETS(M,N,N) for the \(a_t\) and ETS(A,A,N) for the \(b_t\), resulting in oETS\(_G\)(M,N,N)(A,A,N):

oETSGModel1 <- oesg(y, modelA="MNN", modelB="AAN", h=10, holdout=TRUE)
oETSGModel1
## Occurrence state space model estimated: General
## Underlying ETS model: oETS[G](MNN)(AAN)
## 
## Sample size: 110
## Number of estimated parameters: 6
## Number of degrees of freedom: 104
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 117.6459 118.4614 133.8488 135.7655
plot(oETSGModel1)

The oes() function accepts occurrence="g" and in this case calls for oesg() with the same types of ETS models for both parts:

oETSGModel2 <- oes(y, model="MNN", occurrence="g", h=10, holdout=TRUE)
oETSGModel2
## Occurrence state space model estimated: General
## Underlying ETS model: oETS[G](MNN)(MNN)
## 
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 120.5913 120.9723 131.3932 132.2886
plot(oETSGModel2)

Finally, the more flexible way to construct iETS model would be to do it in two steps: either using oesg() or oes() and then using the es() with the provided model in occurrence variable. But a simpler option is available as well:

es(y, "MMN", occurrence="g", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)
## Time elapsed: 0.24 seconds
## Model estimated: iETS(MMN)
## Occurrence model type: General
## Persistence vector g:
##  alpha   beta 
## 0.0028 0.0028 
## Initial values were optimised.
## 
## Loss function type: MSE; Loss function value: 0.1636
## Error standard deviation: 0.4045
## Sample size: 110
## Number of estimated parameters: 13
## Number of degrees of freedom: 97
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 389.5435 390.1204 424.6498 388.4018 
## 
## Forecast errors:
## Bias: 20.8%; sMSE: 43.1%; RelRMSE: 1.053; sPIS: -945.3%; sCE: -70.2%

iETS\(_A\)

Finally, there is an occurrence type selection mechanism. It tries out all the iETS subtypes of models, discussed above and selects the one that has the lowest information criterion (i.e. AIC). This subtype is called iETS\(_A\) (automatic), although it does not represent any specific model. Here’s an example:

oETSAModel <- oes(y, model="MNN", occurrence="a", h=10, holdout=TRUE)
oETSAModel
## Occurrence state space model estimated: Odds ratio
## Underlying ETS model: oETS[O](MNN)
## Smoothing parameters:
##  level 
## 0.0892 
## Vector of initials:
##  level 
## 0.2235 
## 
## Error standard deviation: 1.4873
## Sample size: 110
## Number of estimated parameters: 2
## Number of degrees of freedom: 108
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 116.5913 116.7035 121.9923 122.2559
plot(oETSAModel)

The main restriction of the iETS models at the moment (smooth v.2.5.0) is that there is no model selection between the ETS models for the occurrence part. This needs to be done manually. Hopefully, this feature will appear in the next release of the package.

The integer-valued iETS

By default, the models assume that the data is continuous, which sounds counter intuitive for the typical intermittent demand forecasting tasks. However, (Svetunkov and Boylan 2017) showed that these models perform quite well in terms of forecasting accuracy for many cases. Still, there is also an option for the rounded up values, which is implemented in the es() function. This is not described in the manual and can be triggered via the rounded=TRUE parameter provided in ellipsis. Here’s an example:

es(rpois(100,0.3), "MNN", occurrence="g", oesmodel="MNN", h=10, holdout=TRUE, silent=FALSE, interval=TRUE, rounded=TRUE)
## Time elapsed: 0.12 seconds
## Model estimated: iETS(MNN)
## Occurrence model type: General
## Persistence vector g:
## alpha 
## 7e-04 
## Initial values were optimised.
## 
## Loss function type: Rounded; Loss function value: 0
## Error standard deviation: 0.074
## Sample size: 90
## Number of estimated parameters: 7
## Number of degrees of freedom: 83
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 120.2603 120.5393 137.7589 120.3876 
## 
## 95% parametric prediction interval was constructed
## 100% of values are in the prediction interval
## Forecast errors:
## Bias: 72.5%; sMSE: 42.2%; RelRMSE: 0.906; sPIS: -2188.6%; sCE: -465.3%

Keep in mind that the model with the rounded up values is estimated differently than it continuous counterpart and produces more adequate results for the highly intermittent data with low level of demand sizes. In all the other cases, the continuous iETS models are recommended. In fact, if you need to produce integer-valued prediction interval, then you can produce the interval from a continuous model and then round them up (see discussion in (Svetunkov and Boylan 2017) for details).

References

Svetunkov, I., and John E. Boylan. 2017. “Multiplicative State-Space Models for Intermittent Time Series.” Department of Management Science, Lancaster University. https://doi.org/10.13140/RG.2.2.35897.06242.