oes() - occurrence part of iETS model

The basics

The canonical general iETS model (called iETS_G) can be summarised as: where y_t is the observed values, z_t is the demand size, which is a pure multiplicative ETS model on its own, w(⋅) is the measurement function, r(⋅) is the error function, f(⋅) is the transition function and g(⋅) 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, ϵ_a, t and ϵ_b, t are the mutually independent error terms that follow unknown distribution, 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. μ_a, t and μ_b, t are the conditional expectations for the unobservable variables a_t and b_t. 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: while for the additive error it is: 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 and their expectations. 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): where l_a, t and l_b, t are the levels for each of the shape parameters and α_a and α_b are the smoothing parameters and the error terms 1 + ϵ_a, t and 1 + ϵ_b, t are positive and have means of one. We do not make any other distributional assumptions concerning the error terms. 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: where Y is the vector of all the in-sample observations, θ is the vector of parameters to estimate (initial values and smoothing parameters), T is the number of all observations, T₀ 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 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 = const, μ_b, t = 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 ≤ 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.

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:

The conditional h-steps ahead mean of the demand occurrence probability is calculated as:

The estimate of the probability p is calculated based on the maximisation of the following concentrated log-likelihood function: where T₀ is the number of zero observations and T₁ 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: where T₁ is the number of non-zero observations and 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)
## Vector of initials:
##  level 
## 0.7182 
## 
## Error standard deviation: 1.078
## Sample size: 110
## Number of estimated parameters: 1
## Number of degrees of freedom: 109
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 132.8257 132.8628 135.5262 135.6132

plot(oETSFModel1)

The occurrence part of the model is supported by adam() function. For example, here’s how the iETS(M,M,N)_F can be constructed:

adam(y, "MMN", occurrence="fixed", h=10, holdout=TRUE, silent=FALSE)

## Time elapsed: 0.03 seconds
## Model estimated using adam() function: iETS(MMN)[F]
## With optimal initialisation
## Occurrence model type: Fixed probability
## Distribution assumed in the model: Mixture of Bernoulli and Gamma
## Loss function type: likelihood; Loss function value: 142.8622
## Persistence vector g:
##  alpha   beta 
## 0.0049 0.0000 
## 
## Sample size: 110
## Number of estimated parameters: 6
## Number of degrees of freedom: 104
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 428.5501 429.1270 444.7530 441.4084 
## 
## Forecast errors:
## Asymmetry: 51.237%; sMSE: 124.982%; rRMSE: 1.123; sPIS: -3855.806%; sCE: 447.368%

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):

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)): where â_t is the estimate of μ_a, t.

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. For 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.001 0.001 
## Vector of initials:
##  level  trend 
## 0.5514 1.0107 
## 
## Error standard deviation: 1.1048
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 112.6883 113.0693 123.4903 124.3856

plot(oETSOModel)

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

adam(y, "MMN", occurrence="o", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)

## Time elapsed: 0.07 seconds
## Model estimated using adam() function: iETS(MMN)[O]
## With optimal initialisation
## Occurrence model type: Odds ratio
## Distribution assumed in the model: Mixture of Bernoulli and Gamma
## Loss function type: likelihood; Loss function value: 142.8622
## Persistence vector g:
##  alpha   beta 
## 0.0049 0.0000 
## 
## Sample size: 110
## Number of estimated parameters: 9
## Number of degrees of freedom: 101
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 408.4127 408.9897 432.7171 415.2710 
## 
## Forecast errors:
## Asymmetry: -20.466%; sMSE: 109.145%; rRMSE: 1.05; sPIS: -928.955%; sCE: -101.275%

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

adam(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):

This model resembles the logistic regression, where the probability is obtained from an underlying regression model of x_t′A. 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: where b̂_t is the estimate of μ_b, t.

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 
## 2e-04 2e-04 
## Vector of initials:
##  level  trend 
## 3.9337 0.9469 
## 
## Error standard deviation: 4.8693
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 114.2637 114.6447 125.0656 125.9610

plot(oETSIModel)

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

adam(y, "MMN", occurrence="i", oesmodel="MMN", h=10, holdout=TRUE, silent=FALSE)

## Time elapsed: 0.07 seconds
## Model estimated using adam() function: iETS(MMN)
## With optimal initialisation
## Occurrence model type: Inverse odds ratio
## Distribution assumed in the model: Mixture of Bernoulli and Gamma
## Loss function type: likelihood; Loss function value: 142.8622
## Persistence vector g:
##  alpha   beta 
## 0.0049 0.0000 
## 
## Sample size: 110
## Number of estimated parameters: 9
## Number of degrees of freedom: 101
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 409.9881 410.5650 434.2924 416.8464 
## 
## Forecast errors:
## Asymmetry: -12.872%; sMSE: 108.408%; rRMSE: 1.046; sPIS: -1104.885%; sCE: -68.399%

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

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

iETS_D

This model appears, when a specific restriction is imposed: The pure multiplicative iETS_G(M,N,N) model is then transformed into iETS_D(M,N,N): An option with the additive model in this case has a different, more complicated form:

The estimation of the multiplicative error model is done using the following set of equations: where and κ is a very small number (for example, κ = 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:

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     0 
## Vector of initials:
##  level  trend 
## 0.4194 1.0092 
## 
## Error standard deviation: 0.6991
## Sample size: 110
## Number of estimated parameters: 4
## Number of degrees of freedom: 106
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 112.6148 112.9957 123.4167 124.3120

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:

adam(y, "MMN", occurrence=oETSDModel, h=10, holdout=TRUE, silent=FALSE)

## Time elapsed: 0.04 seconds
## Model estimated using adam() function: iETS(MMN)[D]
## With optimal initialisation
## Occurrence model type: Direct
## Distribution assumed in the model: Mixture of Bernoulli and Gamma
## Loss function type: likelihood; Loss function value: 142.8622
## Persistence vector g:
##  alpha   beta 
## 0.0049 0.0000 
## 
## Sample size: 110
## Number of estimated parameters: 5
## Number of degrees of freedom: 105
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 400.3391 400.9161 413.8415 415.1975 
## 
## Forecast errors:
## Asymmetry: -25.345%; sMSE: 109.589%; rRMSE: 1.052; sPIS: -740.674%; sCE: -133.889%

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: 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 
## 116.7993 117.6148 133.0021 134.9188

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 
## 118.9443 119.3252 129.7462 130.6415

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.12 seconds
## Model estimated using es() function: iETS(MMN)[G]
## With optimal initialisation
## Occurrence model type: General
## Distribution assumed in the model: Mixture of Bernoulli and Normal
## Loss function type: likelihood; Loss function value: 138.2155
## Persistence vector g:
##  alpha   beta 
## 0.0021 0.0015 
## 
## Sample size: 110
## Number of estimated parameters: 13
## Number of degrees of freedom: 97
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 445.4042 445.9812 480.5105 444.2626 
## 
## Forecast errors:
## Asymmetry: 84.077%; sMSE: 196.047%; rRMSE: 1.407; sPIS: -6437.777%; sCE: 1041.916%

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.0821 
## Vector of initials:
##  level 
## 0.4452 
## 
## Error standard deviation: 1.0484
## Sample size: 110
## Number of estimated parameters: 2
## Number of degrees of freedom: 108
## Information criteria: 
##      AIC     AICc      BIC     BICc 
## 114.9443 115.0564 120.3453 120.6088

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.

References