SARIMA Models

Introduction

So far, we have restricted our attention to non-seasonal data and non-seasonal ARIMA models. However, ARIMA models are also capable of modelling a wide range of seasonal data.

A seasonal ARIMA model is formed by including additional seasonal terms in the ARIMA models we have seen so far. It is written as follows: \[ SARIMA\;\;\;\underbrace{(p, d, q)}_{\text{non-seasonal part}} \times \underbrace{(P, D, Q)_s}_{\text{Seasonal part}} \] where \(s=\) number of observations per year. We use uppercase notation for the seasonal parts of the model, and lowercase notation for the non-seasonal parts of the model.

Seasonal differencing

Almost by definition, it may be necessary to examine differenced data when we have seasonality. Seasonality usually causes the series to be nonstationary because the average values at some particular times within the seasonal span (months, for example) may be different than the average values at other times. For instance, our sales of ice cream and cold drinks will always be higher in the summer months.

Seasonal differencing is defined as a difference between a value and a value with lag that is a multiple of \(s\).

With \(s=12\), which may occur with monthly data, a seasonal difference is \((1-B^{12})y_t = y_t - y_{t-12}\). The differences (from the previous year) may be about the same for each month of the year giving us a stationary series.
With \(s=4\), which may occur with quarterly data, a seasonal difference is \((1-B^4) y_t = y_t - y_{t-4}\)

The SARIMA Notation

\[ SARIMA(p,d,q) \times (P,D,Q)_s \]

\(AR(p)\): Autoregressive component of order \(p\)
\(MA(q)\): Moving average component of order \(q\)
\(I(d)\): Integrated component of order \(d\)
Seasonal \(AR(p)\): Seasonal autoregressive component of order \(P\)
\(MA(Q)\): Seasonal moving average component of order \(Q\)
Seasonal \(I(D)\): ISeasonal integrated component of order \(D\)
\(s\): Seasonal period

The mathematical representation of SARIMA:

\[ (1-\phi_1 B)(1-\Phi_1B^s)(1-B)(1-B^s)y_t = (1+\theta_1 B)(1 + \Theta_1B^s) \epsilon_t \tag{1} \]

\(B\) is the backward shift (lag) operator, , \(BY_t = y_{t-1}\)
\(\phi_1\) is the non-seasonal autoregressive coefficient
\(\Phi_1\) is the seasonal autoregressive coefficient
\(\theta_1\) is the nonseasonal moving average coefficient
\(\Theta_1\) is the seasonal moving average coefficient
\(s\) is the seasonal period
\(\epsilon_t\) is the white noise error at time \(t\)

Lеt’s go through each component of the equation:

Autoregressive (AR) Component: The autoregressive non-seasonal component represented by \((1-\phi_1 B)\) captures the relationship between the current observation and a certain number of lagged observations (previous values in the time series). The \(B\) term represents the backshift operator or the lag operator, which shifts the time series backward by a certain number of time period. The order of the autoregressive component, denoted by (\(p\)), determines the number of past values considered in the model.
Seasonal Autoregressive (SAR) Component: Thе seasonal autoregressive component is represented by \((1 - \Phi_1 B^s)\). This component captures the relationship between the current observation and a certain number of lagged observations at seasonal intervals. Thе \(B^s\) term represents the backshift operator applied to the seasonal lagged observations.
Non-Seasonal Differencing Component: The non-seasonal differencing component is represented by \((1-B)\), where \(d\) is the order of non-seasonal differencing. This component is used to make the time series stationary by differencing it a certain number of times.
Moving Average (MA) Component: The moving average component is represented by \((1 + \theta_1 B )\). This component captures the relationship between the current observation and the residual errors from a moving average model appliеd to lagged observations.
Seasonal Moving Average (SMA) Component: The seasonal moving average component is represented by \((1 + \Theta_1B^s)\). This component captures the relationship between the current observation and the residual errors from a moving average model applied to lagged observations at seasonal intervals.
Error Term: The error term is denoted by \(\epsilon_t\). It represents the random noise or unexplained variation in the time series.

The \(SARIMA(p,d,q) \times (P,D,Q)_s\) model can also be written as follows:

\[ \Phi(B^S) \phi(B)(X_t - \mu) = \Theta(B^S)\theta(B)\epsilon_t \tag{2} \]

The non-seasonal parts are:

AR: \(\phi(B) = 1-\phi_1B- \ldots-\phi_pB^p\)
MA: \(\theta(B) = 1 + \theta_1B + \ldots + \theta_q B^q\)

The seasonal parts are:

Seasonal AR: \(\Phi(B^S)= 1- \Phi_1B^S - \ldots - \Phi_PB^{PS}\)
Seasonal MA: \(\Theta(B^S) = 1 + \Theta_1B^S + \ldots + \Theta_QB^{QS}\)

Note that on the left side of equation (2) the seasonal and non-seasonal AR components multiply each other, and on the right side of equation (2) the seasonal and non-seasonal MA components multiply each other.

Example: \(SARIMA(0,0,1) \times (0, 0, 1)_{12}\)

The model includes a non-seasonal MA(1) term, a seasonal MA(1) term, no differencing, no AR terms and the seasonal period is \(s = 12\).

The non-seasonal MA(1) polynomial is \[ \theta(B) = 1+\theta_1B \].

The seasonal MA(1) polynomial is

\[ \Theta(B^{12}) = 1 + \Theta_1B^{12} \]

The model is \[ (x_t - \mu) = \Theta(B^{12}) \theta(B)\epsilon_t = (1+\Theta_1B^{12})(1+ \theta_1 B)\epsilon_t \]

When we multiply the two polynomials on the right side, we get \[ \begin{aligned} (x_t - \mu) &= (1+\theta_1B + \Theta_1B^{12} + \theta_1\Theta_1B^{13})\epsilon_t \\ &=\epsilon_t + \theta_1\epsilon_{t-1} + \Theta_1\epsilon_{t-12} + \theta_1\Theta_1 \epsilon_{t-13} \end{aligned} \]

Thus, the model has MA terms at lags 1, 12, and 13. This leads many to think that the identifying ACF for the model will have non-zero autocorrelations only at lags 1, 12, and 13.

Example: \(SARIMA(1,0,0) \times (1, 0, 0)_{12}\)

The model includes a non-seasonal AR(1) term, a seasonal AR(1) term, no differencing, no MA terms and the seasonal period is \(s = 12\).

The non-seasonal AR(1) polynomial is \[ \phi(B) = 1 + \phi_1B \] The seasonal AR(1) polynomial is

\[ \Phi(B) = 1 + \Phi_1B^{12} \]

Hence, the model is \[ (1 + \phi_1B)(1 + \Phi_1B^{12})(y_t - \mu) = \epsilon_t \]

If we let \(z_t = y_t - \mu\), (for simplicity), multiplying the two AR components and pushing all but \(z_t\) to the right side we get \[ z_t = \phi_1z_{t-1} + \Phi_1z_{t-12} - \phi_1 \Phi_1z_{t-13} + \epsilon_t \]

This is an AR model with predictors at lags 1, 12, and 13. An example PACF plot for this model is

Some guidelines in identifying a seasonal model

Step 1: Do a time series plot of the data.

Examine it for features such as trend and seasonality. You’ll know that you’ve gathered seasonal data (months, quarters, etc.,) so look at the pattern across those time units (months, etc.) to see if there is indeed a seasonal pattern.

Step 2: Do any necessary differencing.

If there is seasonality and no trend, then take a difference of lag \(s\). For instance, take a 12th difference for monthly data with seasonality. Seasonality will appear in the ACF by tapering slowly at multiples of \(s\). Below is an example of a time series that requires a seasonal difference

Note that the time series plot shows a clear seasonal pattern that repeats over 12 time points. Seasonal differences are supported in the ACF/PACF of the original data because the first seasonal lag in the ACF is close to 1 and decays slowly over multiples of \(s=12\).

Once seasonal differences are taken, the ACF/PACF plots of twelfth differences support a seasonal AR(1) pattern. See figure below.

If there is linear trend and no obvious seasonality, then take a first difference. If there is a curved trend, consider a transformation of the data before differencing.

If there is both trend and seasonality, apply a seasonal difference to the data and then re-evaluate the trend. If a trend remains, then take first differences. For instance, if the series is called \(y\), the commands in R would be:

diff12=diff(y, 12)
plot(diff12)
acf2(diff12)
diff1and12 = diff(diff12, 1)

If there is neither obvious trend nor seasonality, don’t take any differences.

Step 3: Examine the ACF and PACF of the differenced data (if differencing is necessary).

We’re using this information to determine possible models. This can be tricky going involving some (educated) guessing. Some basic guidance:

Non-seasonal terms: Examine the early lags (1, 2, 3, …) to judge non-seasonal terms. Spikes in the ACF (at low lags) with a tapering PACF indicate non-seasonal MA terms. Spikes in the PACF (at low lags) with a tapering ACF indicate possible non-seasonal AR terms.
Seasonal terms: Examine the patterns across lags that are multiples of S. For example, for monthly data, look at lags 12, 24, 36, and so on (probably won’t need to look at much more than the first two or three seasonal multiples). Judge the ACF and PACF at the seasonal lags in the same way you do for the earlier lags.

Step 4: Estimate the model(s) that might be reasonable on the basis of step 3.

Don’t forget to include any differencing that you did before looking at the ACF and PACF. In the software, specify the original series as the data and then indicate the desired differencing when specifying parameters in the arima command that you’re using.

Step 5: Examine the residuals (with ACF, Box-Pierce, and any other means) to see if the model seems good.

Compare AIC or BIC values to choose among several models.

If things don’t look good here, it’s back to Step 3 (or maybe even Step 2)

An example: European quarterly retail trade

autoplot(euretail) + ylab("Retail index") + xlab("Year")

The data are clearly non-stationary, with some seasonality, so we will first take a seasonal difference. The seasonally differenced data are shown below.

euretail %>% diff(lag=4) %>% ggtsdisplay()

These also appear to be non-stationary, so we take an additional first difference and these are shown below.

euretail %>% diff(lag=4) %>% diff() %>% ggtsdisplay()

Our aim now is to find an appropriate ARIMA model based on the ACF and PACF plots. The significant spike at lag 1 in the ACF suggests a non-seasonal MA(1) component, and the significant spike at lag 4 in the ACF suggests a seasonal MA(1) component. Consequently, we begin with an \(ARIMA(0,1,1)\times(0,1,1)_4\) model, indicating a first and seasonal difference, and non-seasonal and seasonal MA(1) components.

fit1 <- euretail %>%
  Arima(order=c(0,1,1), seasonal=c(0,1,1))
summary(fit1)

## Series: . 
## ARIMA(0,1,1)(0,1,1)[4] 
## 
## Coefficients:
##          ma1     sma1
##       0.2903  -0.6913
## s.e.  0.1118   0.1193
## 
## sigma^2 = 0.188:  log likelihood = -34.64
## AIC=75.28   AICc=75.72   BIC=81.51
## 
## Training set error measures:
##                      ME     RMSE       MAE         MPE      MAPE      MASE
## Training set -0.0514931 0.409237 0.2843862 -0.04949815 0.2945738 0.2313337
##                    ACF1
## Training set 0.06136677

The residuals for the fitted model are shown in the plot below.

ggtsdisplay(residuals(fit1))

checkresiduals(fit1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[4]
## Q* = 10.654, df = 6, p-value = 0.09968
## 
## Model df: 2.   Total lags used: 8

By analogous logic applied to the PACF, we could also have started with an \(ARIMA(1,1,0) \times (1,1,0)_4\) model.

fit2 <- euretail %>%
  Arima(order=c(1,1,0), seasonal=c(1,1,0))
summary(fit2)

## Series: . 
## ARIMA(1,1,0)(1,1,0)[4] 
## 
## Coefficients:
##          ar1     sar1
##       0.4224  -0.5297
## s.e.  0.1214   0.1095
## 
## sigma^2 = 0.1962:  log likelihood = -35.29
## AIC=76.57   AICc=77.01   BIC=82.8
## 
## Training set error measures:
##                       ME      RMSE       MAE          MPE      MAPE      MASE
## Training set -0.01214714 0.4179911 0.3092987 -0.009643066 0.3205153 0.2515987
##                     ACF1
## Training set -0.08624114

ggtsdisplay(residuals(fit2))

checkresiduals(fit2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)(1,1,0)[4]
## Q* = 10.827, df = 6, p-value = 0.09387
## 
## Model df: 2.   Total lags used: 8

Next we fit \(ARIMA(0,1,2) \times (0,1,1)_4\).

fit3 <- euretail %>%
  Arima(order=c(0,1,2), seasonal=c(0,1,1))
summary(fit3)

## Series: . 
## ARIMA(0,1,2)(0,1,1)[4] 
## 
## Coefficients:
##          ma1     ma2     sma1
##       0.2303  0.2502  -0.6991
## s.e.  0.1484  0.1188   0.1284
## 
## sigma^2 = 0.1789:  log likelihood = -32.76
## AIC=73.53   AICc=74.27   BIC=81.84
## 
## Training set error measures:
##                      ME      RMSE       MAE         MPE      MAPE      MASE
## Training set -0.0448743 0.3956301 0.2816785 -0.04296557 0.2916949 0.2291311
##                    ACF1
## Training set 0.06982161

ggtsdisplay(residuals(fit3))

checkresiduals(fit3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,1)[4]
## Q* = 7.946, df = 5, p-value = 0.1592
## 
## Model df: 3.   Total lags used: 8

Finally, we fit a \(ARIMA(0,1,3) \times (0,1,1)_4\) model.

fit4 <- euretail  |> 
  Arima(order=c(0,1,3), seasonal=c(0,1,1))
summary(fit4)

## Series: euretail 
## ARIMA(0,1,3)(0,1,1)[4] 
## 
## Coefficients:
##          ma1     ma2     ma3     sma1
##       0.2630  0.3694  0.4200  -0.6636
## s.e.  0.1237  0.1255  0.1294   0.1545
## 
## sigma^2 = 0.156:  log likelihood = -28.63
## AIC=67.26   AICc=68.39   BIC=77.65
## 
## Training set error measures:
##                       ME      RMSE       MAE         MPE      MAPE      MASE
## Training set -0.02965298 0.3661147 0.2787802 -0.02795377 0.2885545 0.2267735
##                     ACF1
## Training set 0.006455781

ggtsdisplay(residuals(fit4))

checkresiduals(fit4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,3)(0,1,1)[4]
## Q* = 0.51128, df = 4, p-value = 0.9724
## 
## Model df: 4.   Total lags used: 8

Thus, we now have a seasonal ARIMA model that passes the required checks and is ready for forecasting. Forecasts from the model for the next three years are shown below.

fit4 |> 
  forecast(h=12) |> 
  autoplot()

The forecasts follow the recent trend in the data, because of the double differencing. The large and rapidly increasing prediction intervals show that the retail trade index could start increasing or decreasing at any time, while the point forecasts trend downwards, the prediction intervals allow for the data to trend upwards during the forecast period.

We could have used auto.arima() function to do most of this work for us. It would have given the same result.

auto.arima(euretail)

## Series: euretail 
## ARIMA(0,1,3)(0,1,1)[4] 
## 
## Coefficients:
##          ma1     ma2     ma3     sma1
##       0.2630  0.3694  0.4200  -0.6636
## s.e.  0.1237  0.1255  0.1294   0.1545
## 
## sigma^2 = 0.156:  log likelihood = -28.63
## AIC=67.26   AICc=68.39   BIC=77.65

The auto.arima() function uses nsdiffs() to determine the number of seasonal differences to use (\(D\)), and ndiffs() to determine the number of ordinary differences to use (\(d\)). The selection of the other model parameters \(p\), \(q\), \(P\), and \(Q\) are all determined by minimizing the AICc, as with non-seasonal ARIMA models.