om() - ADAM-based occurrence model

The om() function is the ADAM-based occurrence model for intermittent demand data. It estimates the probability of demand occurrence \(p_t\) using the full ADAM state-space framework, which means the underlying model for the occurrence dynamics can be an ETS model, an ARIMA model, a regression model, or a combination of all three.

The five link functions (occurrence types) supported by om() are identical to those in oes(): fixed, odds-ratio, inverse-odds-ratio, direct, and general. For the mathematical derivations behind each link function see the oes vignette. This vignette focuses on the added flexibility that om() inherits from ADAM.

The data

We use the same artificial data as in the oes vignette:

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

The series starts sparse (low Poisson rate) and becomes increasingly dense. Any non-zero value in the input is automatically binarised to 1 by om(), so passing a count series is fine.

om(): ETS occurrence model

The simplest call resembles oes(). We fix an ETS(M,N,N) model and use the odds-ratio link:

omETS <- om(y, model="MNN", occurrence="odds-ratio", h=10, holdout=TRUE)
omETS

## Occurrence model
## Time elapsed: 0.02 seconds
## Model estimated using om() function: oETS(MNN)[O]
## With backcasting initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 61.5463
## Persistence vector g:
##  alpha 
## 0.5605 
## 
## Sample size: 110
## Number of estimated parameters: 1
## Number of degrees of freedom: 109
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 125.0925 125.1295 127.7930 127.8800 
## 
## Forecast errors:
## Asymmetry: 100%; sMSE: 0.072%; rRMSE: Inf; sPIS: -147.893%; sCE: 26.89%

plot(omETS, 7)

Model selection across ETS components is supported through the same "Z" / "X" / "Y" wildcards used in adam() and es(). The default model="ZXZ" searches across all sensible error types and seasonal patterns:

omETSAuto <- om(y, occurrence="odds-ratio", h=10, holdout=TRUE)
omETSAuto

## Occurrence model
## Time elapsed: 0.1 seconds
## Model estimated using om() function: oETS(AAN)[O]
## With backcasting initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 58.7418
## Persistence vector g:
##  alpha   beta 
## 0.0320 0.0012 
## 
## Sample size: 110
## Number of estimated parameters: 2
## Number of degrees of freedom: 108
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 121.4836 121.5958 126.8846 127.1481 
## 
## Forecast errors:
## Asymmetry: 100%; sMSE: 0.684%; rRMSE: Inf; sPIS: -472.993%; sCE: 82.384%

The fitted occurrence model can be passed directly to adam() for a full iETS model that jointly describes demand occurrence and demand sizes:

adam(y, "MNN", occurrence=omETS, h=10, holdout=TRUE, silent=FALSE)

## Time elapsed: 0.02 seconds
## Model estimated using adam() function: iETS(MNN)[O]
## With backcasting initialisation
## Occurrence model type: Odds ratio
## Distribution assumed in the model: Mixture of Bernoulli and Gamma
## Loss function type: likelihood; Loss function value: 145.5116
## Persistence vector g:
##  alpha 
## 0.1784 
## 
## Sample size: 110
## Number of estimated parameters: 2
## Number of degrees of freedom: 108
## Number of provided parameters: 1
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 418.1158 418.2279 423.5167 423.7803 
## 
## Forecast errors:
## Asymmetry: 17.581%; sMSE: 47.116%; rRMSE: 0.907; sPIS: 1204.445%; sCE: -8.642%

om(): ARIMA occurrence model

Any ADAM-supported ARIMA structure can drive the occurrence dynamics. Set model="NNN" to suppress ETS components and specify the ARIMA orders directly. Note that automatic ARIMA order selection requires using auto.om() (see below) or passing orders=list(ar=..., i=..., ma=..., select=TRUE):

omARIMA <- om(y, model="NNN",
              orders=list(ar=1, i=0, ma=1),
              occurrence="odds-ratio", h=10, holdout=TRUE)
omARIMA

## Occurrence model
## Time elapsed: 0.02 seconds
## Model estimated using om() function: oARIMA(1,0,1)[O]
## With backcasting initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 60.32
## ARMA parameters of the model:
##       Lag 1
## AR(1)     1
##         Lag 1
## MA(1) -0.6735
## 
## Sample size: 110
## Number of estimated parameters: 2
## Number of degrees of freedom: 108
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 124.6400 124.7522 130.0410 130.3046 
## 
## Forecast errors:
## Asymmetry: 100%; sMSE: 0.243%; rRMSE: Inf; sPIS: -271.342%; sCE: 49.338%

plot(omARIMA, 7)

om(): Regression occurrence model

External regressors are incorporated via cbind(), exactly as in adam(). The first column of the matrix is treated as the response; subsequent columns become regressors. Here we use the Seatbelts dataset from the datasets package: the response is a binary indicator for above-median driver casualties, and the predictor is the seatbelt law indicator:

data(Seatbelts)
highDeaths <- as.integer(Seatbelts[, "DriversKilled"] > median(Seatbelts[, "DriversKilled"]))
ySeat <- ts(cbind(highDeaths, law=Seatbelts[, "law"]), start=c(1969, 1), frequency=12)
omReg <- om(ySeat, model="NNN",
            occurrence="odds-ratio", regressors="use",
            h=12, holdout=TRUE)
omReg

## Occurrence model
## Time elapsed: 0.01 seconds
## Model estimated using om() function: Regression
## With optimal initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 122.0001
## Sample size: 180
## Number of estimated parameters: 2
## Number of degrees of freedom: 178
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 248.0002 248.0680 254.3861 254.5621 
## 
## Forecast errors:
## Asymmetry: 3.655%; sMSE: 24.505%; rRMSE: 0.857; sPIS: 421.551%; sCE: 181.3%

plot(omReg, 7)

Setting regressors="select" would apply automatic variable selection, dropping regressors that do not improve the information criterion.

omg(): the general two-component occurrence model

omg() fits the iETS\(_G\) model directly, estimating separate ETS (or ARIMA, or regression) structures for both the \(a_t\) (numerator) and \(b_t\) (denominator) components of the probability link. Unlike om(..., occurrence="general"), which applies the same model to both sides, omg() lets you specify different model types for each:

omgModel <- omg(y, modelA="MNN", modelB="ANN", h=10, holdout=TRUE)
omgModel

## Time elapsed: 0.03 seconds
## General occurrence model: oETS[G](MNN)(ANN)
## Model A: ETS(MNN)
## Model B: ETS(ANN)
## 
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 59.6297
## 
## Sample size: 110
## Number of estimated parameters: 2
## Number of degrees of freedom: 108
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 123.2595 123.3716 128.6604 128.9240

plot(omgModel, 7)

Calling om(..., occurrence="general") is a shortcut that delegates to omg() with the same model specification on both sides.

auto.om(): automatic occurrence type selection

auto.om() fits om() (or omg()) for each candidate occurrence type and returns the model with the lowest information criterion. By default all five types are tested:

autoOmModel <- auto.om(y, model="MNN", h=10, holdout=TRUE)
autoOmModel

## Occurrence model
## Time elapsed: 0.61 seconds
## Model estimated using auto.om() function: oETS(MNN)[D]
## With backcasting initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 59.2279
## Persistence vector g:
##  alpha 
## 0.1137 
## 
## Sample size: 110
## Number of estimated parameters: 1
## Number of degrees of freedom: 109
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 120.4558 120.4929 123.1563 123.2433 
## 
## Forecast errors:
## Asymmetry: 100%; sMSE: 0.008%; rRMSE: Inf; sPIS: -49.013%; sCE: 8.912%

Passing orders=list(ar=2, i=1, ma=2, select=TRUE) activates automatic ARIMA order selection for each candidate occurrence type, returning the best overall combination of occurrence link and ARIMA structure:

autoOmARIMA <- auto.om(y, model="MNN",
                       orders=list(ar=2, i=1, ma=2, select=TRUE),
                       h=10, holdout=TRUE)
autoOmARIMA

## Occurrence model
## Time elapsed: 0.47 seconds
## Model estimated using auto.om() function: oETS(MNN)[D]
## With backcasting initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 59.2279
## Persistence vector g:
##  alpha 
## 0.1137 
## 
## Sample size: 110
## Number of estimated parameters: 1
## Number of degrees of freedom: 109
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 120.4558 120.4929 123.1563 123.2433 
## 
## Forecast errors:
## Asymmetry: 100%; sMSE: 0.008%; rRMSE: Inf; sPIS: -49.013%; sCE: 8.912%

A subset of occurrence types can be tested by supplying an explicit vector:

autoOmSubset <- auto.om(y, model="MNN",
                        occurrence=c("fixed", "odds-ratio", "inverse-odds-ratio"),
                        h=10, holdout=TRUE)
autoOmSubset

## Occurrence model
## Time elapsed: 0.43 seconds
## Model estimated using auto.om() function: oETS(MNN)[O]
## With backcasting initialisation
## Distribution assumed in the model: Cumulative Logistic
## Loss function type: likelihood; Loss function value: 61.5463
## Persistence vector g:
##  alpha 
## 0.5605 
## 
## Sample size: 110
## Number of estimated parameters: 1
## Number of degrees of freedom: 109
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 125.0925 125.1295 127.7930 127.8800 
## 
## Forecast errors:
## Asymmetry: 100%; sMSE: 0.072%; rRMSE: Inf; sPIS: -147.893%; sCE: 26.89%