ssarima() - State-Space ARIMA

Ivan Svetunkov

2019-07-24

SSARIMA stands for “State-space ARIMA” or “Several Seasonalities ARIMA”. Both names show what happens in the heart of the function: it constructs ARIMA in a state-space form and allows to model several (actually more than several) seasonalities. ssarima() is a function included in smooth package. This vignette covers ssarima() and auto.ssarima() functions. For more details about the underlying model, read (Svetunkov and Boylan 2019).

As usual, we will use data from Mcomp package, so it is advised to install it.

Let’s load the necessary packages:

require(smooth)
require(Mcomp)

The default call constructs ARIMA(0,1,1):

ssarima(M3$N2457$x, h=18, silent=FALSE)
## Time elapsed: 0.02 seconds
## Model estimated: ARIMA(0,1,1)
## Matrix of MA terms:
##         Lag 1
## MA(1) -0.7941
## Initial values were produced using backcasting.
## 
## Loss function type: MSE; Loss function value: 4401089.934
## Error standard deviation: 2097.878
## Sample size: 115
## Number of estimated parameters: 2
## Number of degrees of freedom: 113
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 2089.553 2089.660 2095.042 2095.297

Some more complicated model can be defined using parameter orders the following way:

ssarima(M3$N2457$x, orders=list(ar=c(0,1),i=c(1,0),ma=c(1,1)),lags=c(1,12),h=18)
## Time elapsed: 0.1 seconds
## Model estimated: SARIMA(0,1,1)[1](1,0,1)[12]
## Matrix of AR terms:
##       Lag 12
## AR(1) 0.8298
## Matrix of MA terms:
##         Lag 1  Lag 12
## MA(1) -0.8141 -0.3689
## Initial values were produced using backcasting.
## 
## Loss function type: MSE; Loss function value: 3618199.0582
## Error standard deviation: 1902.156
## Sample size: 115
## Number of estimated parameters: 4
## Number of degrees of freedom: 111
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 2071.027 2071.390 2082.007 2082.869

This would construct us seasonal ARIMA(0,1,1)(1,0,1)\(_{12}\).

We could try selecting orders manually, but this can also be done automatically via auto.ssarima() function:

auto.ssarima(M3$N2457$x, h=18)
## Time elapsed: 1.99 seconds
## Model estimated: SARIMA(0,1,2)[1](0,0,3)[12]
## Matrix of MA terms:
##         Lag 1 Lag 12
## MA(1) -0.5942 0.3920
## MA(2) -0.2769 0.6312
## MA(3)  0.0000 0.0494
## Initial values were produced using backcasting.
## 
## Loss function type: MSE; Loss function value: 3183490.065
## Error standard deviation: 1784.234
## Sample size: 115
## Number of estimated parameters: 6
## Number of degrees of freedom: 109
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 2060.307 2061.085 2076.777 2078.622

Automatic order selection in SSARIMA with optimised initials does not work well and in general is not recommended. This is partially because of the possible high number of parameters in some models and partially because of potential overfitting of first observations when non-zero order of AR is selected. This problem can be seen on example of another time series (which has complicated seasonality):

auto.ssarima(M3$N1683$x, h=18, initial="backcasting")
## Time elapsed: 3.34 seconds
## Model estimated: SARIMA(0,0,3)[1](3,1,3)[12]
## Matrix of AR terms:
##        Lag 12
## AR(1) -0.7001
## AR(2) -0.5211
## AR(3)  0.2211
## Matrix of MA terms:
##        Lag 1  Lag 12
## MA(1) 0.0310 -0.2591
## MA(2) 0.1146  0.2071
## MA(3) 0.3351 -0.6262
## Initial values were produced using backcasting.
## 
## Loss function type: MSE; Loss function value: 79260.4687
## Error standard deviation: 281.5324
## Sample size: 108
## Number of estimated parameters: 10
## Number of degrees of freedom: 98
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1544.784 1547.052 1571.605 1576.915
auto.ssarima(M3$N1683$x, h=18, initial="optimal")
## Time elapsed: 3.99 seconds
## Model estimated: ARIMA(0,0,3) with constant
## Matrix of MA terms:
##        Lag 1
## MA(1) 0.3558
## MA(2) 0.2926
## MA(3) 0.3104
## Constant value is: 3786.8136
## Initial values were optimised.
## 
## Loss function type: MSE; Loss function value: 169345.8886
## Error standard deviation: 411.5166
## Sample size: 108
## Number of estimated parameters: 8
## Number of degrees of freedom: 100
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1622.778 1624.233 1644.235 1647.640

As can be seen from the second graph, ssarima with optimal initial does not select seasonal model and reverts to ARIMA(0,0,3) with constant. In theory this can be due to implemented order selection algorithm, however if we estimate all the model in the pool separately, we will see that this model is optimal for this time series when this type of initials is used.

A power of ssarima() function is that it can estimate SARIMA models with multiple seasonalities. For example, SARIMA(0,1,1)(0,0,1)_6(1,0,1)_12 model can be estimated the following way:

ssarima(M3$N2457$x, orders=list(ar=c(0,0,1),i=c(1,0,0),ma=c(1,1,1)),lags=c(1,6,12),h=18, silent=FALSE)
## Time elapsed: 0.22 seconds
## Model estimated: SARIMA(0,1,1)[1](0,0,1)[6](1,0,1)[12]
## Matrix of AR terms:
##       Lag 12
## AR(1) 0.7505
## Matrix of MA terms:
##         Lag 1   Lag 6  Lag 12
## MA(1) -0.8038 -0.2056 -0.3107
## Initial values were produced using backcasting.
## 
## Loss function type: MSE; Loss function value: 3448766.7318
## Error standard deviation: 1857.085
## Sample size: 115
## Number of estimated parameters: 5
## Number of degrees of freedom: 110
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 2067.512 2068.062 2081.236 2082.542

It probably does not make much sense for this type of data, it would make more sense on high frequency data (for example, taylor series from forecast package). However, keep in mind that multiple seasonal ARIMAs are very slow in estimation and are very capricious. So it is really hard to obtain an appropriate and efficient multiple seasonal ARIMA model.

Now let’s introduce some artificial exogenous variables:

x <- cbind(rnorm(length(M3$N2457$x),50,3),rnorm(length(M3$N2457$x),100,7))

If we save model:

ourModel <- auto.ssarima(M3$N2457$x, h=18, holdout=TRUE, xreg=x, updateX=TRUE)

we can then reuse it:

ssarima(M3$N2457$x, model=ourModel, h=18, holdout=FALSE, xreg=x, updateX=TRUE, interval=TRUE)
## Time elapsed: 0.37 seconds
## Model estimated: SARIMAX(0,0,2)[1](0,0,3)[12] with constant
## Matrix of MA terms:
##        Lag 1 Lag 12
## MA(1) 0.0940 0.0834
## MA(2) 0.0903 0.1032
## MA(3) 0.0000 0.1087
## Constant value is: 3086.6875
## Initial values were provided by user.
## Xreg coefficients were estimated in a crazy style
## 
## Loss function type: MSE; Loss function value: 4413915.0996
## Error standard deviation: 2100.932
## Sample size: 115
## Number of estimated parameters: 1
## Number of provided parameters: 52
## Number of degrees of freedom: 114
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 2087.887 2087.923 2090.632 2090.716 
## 
## 95% parametric prediction interval was constructed

Finally, we can combine several SARIMA models:

ssarima(M3$N2457$x, h=18, holdout=FALSE, interval=TRUE, combine=TRUE)
## Time elapsed: 0.01 seconds
## Model estimated: ARIMA(0,1,1)
## Matrix of MA terms:
##         Lag 1
## MA(1) -0.7941
## Initial values were produced using backcasting.
## 
## Loss function type: MSE; Loss function value: 4401089.934
## Error standard deviation: 2097.878
## Sample size: 115
## Number of estimated parameters: 2
## Number of degrees of freedom: 113
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 2089.553 2089.660 2095.042 2095.297 
## 
## 95% parametric prediction interval was constructed

References

Svetunkov, Ivan, and John E. Boylan. 2019. “State-space ARIMA for supply-chain forecasting.” International Journal of Production Research 0 (0). Taylor & Francis: 1–10. https://doi.org/10.1080/00207543.2019.1600764.