ssarima() - State-Space ARIMA

Ivan Svetunkov

2026-09-03

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)

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

ssarima(AirPassengers, h=12, silent=FALSE)
## Time elapsed: 0.01 seconds
## Model estimated using ssarima() function: SSARIMA(0,1,1)
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 700.9087
## ARMA parameters of the model:
##        Lag 1
## MA(1) 0.4064
## 
## Sample size: 144
## Number of estimated parameters: 3
## Number of degrees of freedom: 141
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1407.817 1407.989 1416.727 1417.153

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

ssarima(AirPassengers, orders=list(ar=c(0,1),i=c(1,0),ma=c(1,1)), lags=c(1,12), h=12)
## Time elapsed: 0.04 seconds
## Model estimated using ssarima() function: SSARIMA(0,1,1)[1](1,0,1)[12]
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 558.1803
## ARMA parameters of the model:
##       Lag 12
## AR(1) 1.0888
##         Lag 1  Lag 12
## MA(1) -0.3945 -0.1584
## 
## Sample size: 144
## Number of estimated parameters: 17
## Number of degrees of freedom: 127
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1150.361 1155.218 1200.847 1212.917

This would construct 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(AirPassengers, h=12)
## Time elapsed: 2.02 seconds
## Model estimated using ssarima() function: SSARIMA(0,1,3)[1](0,1,0)[12] with drift
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 558.4325
## Intercept/Drift value: 0.3338
## ARMA parameters of the model:
##         Lag 1
## MA(1) -0.3424
## MA(2)  0.0500
## MA(3) -0.2528
## 
## Sample size: 144
## Number of estimated parameters: 19
## Number of degrees of freedom: 125
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1154.865 1160.994 1211.291 1226.522

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:

auto.ssarima(AirPassengers, h=12, initial="backcasting")
## Time elapsed: 2.38 seconds
## Model estimated using ssarima() function: SSARIMA(0,1,3)[1](0,1,0)[12] with drift
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 558.4325
## Intercept/Drift value: 0.3338
## ARMA parameters of the model:
##         Lag 1
## MA(1) -0.3424
## MA(2)  0.0500
## MA(3) -0.2528
## 
## Sample size: 144
## Number of estimated parameters: 19
## Number of degrees of freedom: 125
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1154.865 1160.994 1211.291 1226.522
auto.ssarima(AirPassengers, h=12, initial="optimal")
## Time elapsed: 12.22 seconds
## Model estimated using ssarima() function: SSARIMA(0,1,3)[1](0,1,0)[12] with drift
## With optimal initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 550.2778
## Intercept/Drift value: 0.2176
## ARMA parameters of the model:
##         Lag 1
## MA(1) -0.3234
## MA(2)  0.0745
## MA(3) -0.1695
## 
## Sample size: 144
## Number of estimated parameters: 18
## Number of degrees of freedom: 126
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1136.556 1142.028 1190.012 1203.610

As can be seen from the example above the model with optimal initials takes more time and we end up with a different model than in the case of backcasting.

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(AirPassengers, orders=list(ar=c(0,0,1),i=c(1,0,0),ma=c(1,1,1)), lags=c(1,6,12), h=12, silent=FALSE)

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. To tackle this issue, I’ve developed an alternative ARIMA model for multiple seasonalities, called msarima().

Now let’s introduce some artificial exogenous variables:

x <- cbind(rnorm(length(AirPassengers),50,3),rnorm(length(AirPassengers),100,7))

If we save model:

ourModel <- auto.ssarima(AirPassengers, h=12, holdout=TRUE, xreg=x)

we can then reuse it:

ssarima(AirPassengers, model=ourModel, h=12, holdout=FALSE, xreg=x, interval=TRUE)
## Time elapsed: 0 seconds
## Model estimated using ssarima() function: SSARIMAX(0,1,3)[1](0,1,0)[12] with drift
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 598.4589
## Intercept/Drift value: 0.5233
## ARMA parameters of the model:
##         Lag 1
## MA(1) -0.2519
## MA(2)  0.0403
## MA(3) -0.2627
## 
## Sample size: 144
## Number of estimated parameters: 3
## Number of degrees of freedom: 141
## Number of provided parameters: 5
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1202.918 1203.089 1211.827 1212.253

Finally, we can combine several SARIMA models:

ssarima(AirPassengers, h=12, holdout=FALSE, interval=TRUE, combine=TRUE)
## Time elapsed: 0.01 seconds
## Model estimated using ssarima() function: SSARIMA(0,1,1)
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 700.9087
## ARMA parameters of the model:
##        Lag 1
## MA(1) 0.4064
## 
## Sample size: 144
## Number of estimated parameters: 3
## Number of degrees of freedom: 141
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1407.817 1407.989 1416.727 1417.153

MSARIMA

While SSARIMA is flexible, it is not fast. In fact, it cannot handle high frequency data well and most probably will take ages to estimate the parameter and produce forecasts. This is because of the transition matrix, which becomes huge in case of multiple seasonalities. The MSARIMA model (Multiple Seasonal ARIMA) is formulated in a different state-space form, which reduces the size of transition matrix, significantly reducing the computational time for cases with high frequency data.

There are auto.msarima() and msarima() function in the package, that do things similar to auto.ssarima() and ssarima(). Here’s just one example of what can be done with it:

msarima(AirPassengers, orders=list(ar=c(0,0,1),i=c(1,0,0),ma=c(1,1,1)),lags=c(1,6,12),h=12, silent=FALSE)
## Time elapsed: 0.07 seconds
## Model estimated using msarima() function: SARIMA(0,1,1)[1](0,0,1)[6](1,0,1)[12]
## With backcasting initialisation
## Distribution assumed in the model: Normal
## Loss function type: likelihood; Loss function value: 555.7288
## ARMA parameters of the model:
##       Lag 12
## AR(1)  1.047
##         Lag 1  Lag 6 Lag 12
## MA(1) -0.3362 0.0416  -0.12
## 
## Sample size: 144
## Number of estimated parameters: 24
## Number of degrees of freedom: 120
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 1159.458 1169.542 1230.733 1255.791

The forecasts of the two models might differ due to the different state space form. The detailed explanation of MSARIMA is given in Chapter 9 of ADAM textbook.

References

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