## 9.9 Seasonal ARIMA models

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:

ARIMA | \(\underbrace{(p, d, q)}\) | \(\underbrace{(P, D, Q)_{m}}\) |

\({\uparrow}\) | \({\uparrow}\) | |

Non-seasonal part | Seasonal part | |

of the model | of the model |

where \(m =\) the seasonal period (e.g., 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.

The seasonal part of the model consists of terms that are similar to the non-seasonal components of the model, but involve backshifts of the seasonal period. For example, an ARIMA(1,1,1)(1,1,1)\(_{4}\) model (without a constant) is for quarterly data (\(m=4\)), and can be written as \[ (1 - \phi_{1}B)~(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} = (1 + \theta_{1}B)~ (1 + \Theta_{1}B^{4})\varepsilon_{t}. \]

The additional seasonal terms are simply multiplied by the non-seasonal terms.

### ACF/PACF

The seasonal part of an AR or MA model will be seen in the seasonal lags of the PACF and ACF. For example, an ARIMA(0,0,0)(0,0,1)\(_{12}\) model will show:

- a spike at lag 12 in the ACF but no other significant spikes;
- exponential decay in the seasonal lags of the PACF (i.e., at lags 12, 24, 36, …).

Similarly, an ARIMA(0,0,0)(1,0,0)\(_{12}\) model will show:

- exponential decay in the seasonal lags of the ACF;
- a single significant spike at lag 12 in the PACF.

In considering the appropriate seasonal orders for a seasonal ARIMA model, restrict attention to the seasonal lags.

The modelling procedure is almost the same as for non-seasonal data, except that we need to select seasonal AR and MA terms as well as the non-seasonal components of the model. The process is best illustrated via examples.

### Example: Monthly US leisure and hospitality employment

We will describe seasonal ARIMA modelling using monthly US employment data for leisure and hospitality jobs from January 2000 to September 2019, shown in Figure 9.18.

```
<- us_employment %>%
leisure filter(Title == "Leisure and Hospitality",
year(Month) > 2000) %>%
mutate(Employed = Employed/1000) %>%
select(Month, Employed)
autoplot(leisure, Employed) +
labs(title = "US employment: leisure and hospitality",
y="Number of people (millions)")
```

The data are clearly non-stationary, with strong seasonality and a nonlinear trend, so we will first take a seasonal difference. The seasonally differenced data are shown in Figure 9.19.

```
%>%
leisure gg_tsdisplay(difference(Employed, 12),
plot_type='partial', lag=36) +
labs(title="Seasonally differenced", y="")
```

These are also clearly non-stationary, so we take a further first difference in Figure 9.20.

```
%>%
leisure gg_tsdisplay(difference(Employed, 12) %>% difference(),
plot_type='partial', lag=36) +
labs(title = "Double differenced", y="")
```

Our aim now is to find an appropriate ARIMA model based on the ACF and PACF shown in Figure 9.20. The significant spike at lag 2 in the ACF suggests a non-seasonal MA(2) component. The significant spike at lag 12 in the ACF suggests a seasonal MA(1) component. Consequently, we begin with an ARIMA(0,1,2)(0,1,1)\(_{12}\) model, indicating a first difference, a seasonal difference, and non-seasonal MA(2) and seasonal MA(1) component. If we had started with the PACF, we may have selected an ARIMA(2,1,0)(0,1,1)\(_{12}\) model — using the PACF to select the non-seasonal part of the model and the ACF to select the seasonal part of the model. We will also include an automatically selected model. By setting `stepwise=FALSE`

and `approximation=FALSE`

, we are making R work extra hard to find a good model. This takes much longer, but with only one series to model, the extra time taken is not a problem.

```
<- leisure %>%
fit model(
arima012011 = ARIMA(Employed ~ pdq(0,1,2) + PDQ(0,1,1)),
arima210011 = ARIMA(Employed ~ pdq(2,1,0) + PDQ(0,1,1)),
auto = ARIMA(Employed, stepwise = FALSE, approx = FALSE)
)%>% pivot_longer(everything(), names_to = "Model name",
fit values_to = "Orders")
#> # A mable: 3 x 2
#> # Key: Model name [3]
#> `Model name` Orders
#> <chr> <model>
#> 1 arima012011 <ARIMA(0,1,2)(0,1,1)[12]>
#> 2 arima210011 <ARIMA(2,1,0)(0,1,1)[12]>
#> 3 auto <ARIMA(2,1,0)(1,1,1)[12]>
glance(fit) %>% arrange(AICc) %>% select(.model:BIC)
#> # A tibble: 3 x 6
#> .model sigma2 log_lik AIC AICc BIC
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 auto 0.00142 395. -780. -780. -763.
#> 2 arima210011 0.00145 392. -776. -776. -763.
#> 3 arima012011 0.00146 391. -775. -775. -761.
```

The `ARIMA()`

function uses `unitroot_nsdiffs()`

to determine \(D\) (the number of seasonal differences to use), and `unitroot_ndiffs()`

to determine \(d\) (the number of ordinary differences to use), when these are not specified. 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.

The three fitted models have similar AICc values, with the automatically selected model being a little better. Our second “guess” of ARIMA(2,1,0)(0,1,1)\(_{12}\) turned out to be very close to the automatically selected model of ARIMA\((2,1,0)(1,1,1)_{12}\).

The residuals for the best model are shown in Figure 9.21.

`%>% select(auto) %>% gg_tsresiduals(lag=36) fit `

One small but significant spike (at lag 11) out of 36 is still consistent with white noise. To be sure, we use a Ljung-Box test, which has a large p-value, confirming that the residuals are similar to white noise. Note that the alternative models also pass this test.

```
augment(fit) %>% features(.innov, ljung_box, lag=24, dof=4)
#> # A tibble: 3 x 3
#> .model lb_stat lb_pvalue
#> <chr> <dbl> <dbl>
#> 1 arima012011 22.4 0.320
#> 2 arima210011 18.9 0.527
#> 3 auto 16.6 0.680
```

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 in Figure 9.22. The forecasts have captured the seasonal pattern very well, and the increasing trend extends the recent pattern. The trend in the forecasts is induced by the double differencing.

```
forecast(fit, h=36) %>%
filter(.model=='auto') %>%
autoplot(leisure) +
labs(title = "US employment: leisure and hospitality",
y="Number of people (millions)")
```

### Example: Corticosteroid drug sales in Australia

For our second example, we will try to forecast monthly corticosteroid drug sales in Australia. These are known as H02 drugs under the Anatomical Therapeutic Chemical classification scheme.

```
<- PBS %>%
h02 filter(ATC2 == "H02") %>%
summarise(Cost = sum(Cost)/1e6)
%>%
h02 mutate(log(Cost)) %>%
pivot_longer(-Month) %>%
ggplot(aes(x = Month, y = value)) +
geom_line() +
facet_grid(name ~ ., scales = "free_y") +
labs(y="", title="Corticosteroid drug scripts (H02)")
```

Data from July 1991 to June 2008 are plotted in Figure 9.23. There is a small increase in the variance with the level, so we take logarithms to stabilise the variance.

The data are strongly seasonal and obviously non-stationary, so seasonal differencing will be used. The seasonally differenced data are shown in Figure 9.24. It is not clear at this point whether we should do another difference or not. We decide not to, but the choice is not obvious.

The last few observations appear to be different (more variable) from the earlier data. This may be due to the fact that data are sometimes revised when earlier sales are reported late.

```
%>% gg_tsdisplay(difference(log(Cost), 12),
h02 plot_type='partial', lag_max = 24)
```

In the plots of the seasonally differenced data, there are spikes in the PACF at lags 12 and 24, but nothing at seasonal lags in the ACF. This may be suggestive of a seasonal AR(2) term. In the non-seasonal lags, there are three significant spikes in the PACF, suggesting a possible AR(3) term. The pattern in the ACF is not indicative of any simple model.

Consequently, this initial analysis suggests that a possible model for these data is an ARIMA(3,0,0)(2,1,0)\(_{12}\). We fit this model, along with some variations on it, and compute the AICc values shown in Table 9.2.

Model | AICc |
---|---|

ARIMA(3,0,1)(0,1,2)\(_{12}\) | -485.5 |

ARIMA(3,0,1)(1,1,1)\(_{12}\) | -484.2 |

ARIMA(3,0,1)(0,1,1)\(_{12}\) | -483.7 |

ARIMA(3,0,1)(2,1,0)\(_{12}\) | -476.3 |

ARIMA(3,0,0)(2,1,0)\(_{12}\) | -475.1 |

ARIMA(3,0,2)(2,1,0)\(_{12}\) | -474.9 |

ARIMA(3,0,1)(1,1,0)\(_{12}\) | -463.4 |

Of these models, the best is the ARIMA(3,0,1)(0,1,2)\(_{12}\) model (i.e., it has the smallest AICc value).

```
<- h02 %>%
fit model(ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2)))
report(fit)
#> Series: Cost
#> Model: ARIMA(3,0,1)(0,1,2)[12]
#> Transformation: log(Cost)
#>
#> Coefficients:
#> ar1 ar2 ar3 ma1 sma1 sma2
#> -0.1603 0.5481 0.5678 0.3827 -0.5222 -0.1768
#> s.e. 0.1636 0.0878 0.0942 0.1895 0.0861 0.0872
#>
#> sigma^2 estimated as 0.004278: log likelihood=250
#> AIC=-486.1 AICc=-485.5 BIC=-463.3
```

The innovation residuals from this model are shown in Figure 9.25.

`%>% gg_tsresiduals(lag_max=36) fit `

```
augment(fit) %>%
features(.innov, ljung_box, lag = 36, dof = 6)
#> # A tibble: 1 x 3
#> .model lb_stat lb_pvalue
#> <chr> <dbl> <dbl>
#> 1 ARIMA(log(Cost) ~ 0 + pdq(3, 0, 1) + PDQ(0, 1, 2)) 50.7 0.0104
```

There are a few significant spikes in the ACF, and the model fails the Ljung-Box test. The model can still be used for forecasting, but the prediction intervals may not be accurate due to the correlated residuals.

Next we will try using the automatic ARIMA algorithm. Running `ARIMA()`

with all arguments left at their default values led to an ARIMA(2,1,0)(0,1,1)\(_{12}\) model. Running `ARIMA()`

with `stepwise=FALSE`

and `approximation=FALSE`

gives an ARIMA(2,1,3)(0,1,1)\(_{12}\) model. However, both models still fail the Ljung-Box test for 36 lags. Sometimes it is just not possible to find a model that passes all of the tests.

### Test set evaluation:

We will compare some of the models fitted so far using a test set consisting of the last two years of data. Thus, we fit the models using data from July 1991 to June 2006, and forecast the script sales for July 2006 – June 2008. The results are summarised in Table 9.3.

.model | RMSE |
---|---|

ARIMA(3,0,1)(1,1,1)\(_{12}\) | 0.0619 |

ARIMA(3,0,1)(0,1,2)\(_{12}\) | 0.0621 |

ARIMA(2,1,1)(0,1,1)\(_{12}\) | 0.0622 |

ARIMA(2,1,2)(0,1,1)\(_{12}\) | 0.0623 |

ARIMA(2,1,4)(0,1,1)\(_{12}\) | 0.0627 |

ARIMA(2,1,3)(0,1,1)\(_{12}\) | 0.0628 |

ARIMA(3,0,1)(0,1,1)\(_{12}\) | 0.0630 |

ARIMA(3,0,2)(0,1,1)\(_{12}\) | 0.0630 |

ARIMA(2,1,0)(0,1,1)\(_{12}\) | 0.0630 |

ARIMA(3,0,1)(0,1,3)\(_{12}\) | 0.0630 |

ARIMA(3,0,3)(0,1,1)\(_{12}\) | 0.0631 |

ARIMA(3,0,2)(2,1,0)\(_{12}\) | 0.0651 |

ARIMA(3,0,1)(2,1,0)\(_{12}\) | 0.0653 |

ARIMA(2,1,0)(1,1,0)\(_{12}\) | 0.0666 |

ARIMA(3,0,1)(1,1,0)\(_{12}\) | 0.0666 |

ARIMA(3,0,0)(2,1,0)\(_{12}\) | 0.0668 |

ARIMA(2,1,0)(0,1,0)\(_{12}\) | 0.0740 |

The models chosen manually are close to the best model over this test set based on the RMSE values, while those models chosen automatically with `ARIMA()`

are not far behind.

When models are compared using AICc values, it is important that all models have the same orders of differencing. However, when comparing models using a test set, it does not matter how the forecasts were produced — the comparisons are always valid. Consequently, in the table above, we can include some models with only seasonal differencing and some models with both first and seasonal differencing, while in the earlier table containing AICc values, we only compared models with seasonal differencing but no first differencing.

None of the models considered here pass all of the residual tests. In practice, we would normally use the best model we could find, even if it did not pass all of the tests.

Forecasts from the ARIMA(3,0,1)(0,1,2)\(_{12}\) model (which has the second lowest RMSE value on the test set, and the best AICc value amongst models with only seasonal differencing) are shown in Figure 9.26.

```
%>%
h02 model(ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2))) %>%
forecast() %>%
autoplot(h02) +
labs(y=" $AU (millions)",
title="Corticosteroid drug scripts (H02) sales")
```