Fitting and forecasting a seasonal ARIMA model

Today let's understand seasonal ARIMA models, or SARIMA: how to add a seasonal part on top of the ARIMA you already know, and use it to fit and forecast a series with a real yearly pattern.

The running example is the monthly turnover of cafes, restaurants and takeaway food services in Victoria, Australia, measured in millions of Australian dollars, from the tsibbledata package's aus_retail dataset. It runs 441 months, April 1982 to December 2018. Every December, Victorians spend noticeably more eating out than in any other month, and that repeating shape is exactly what an ordinary ARIMA(p, d, q) has no way to reach.

Here is the whole series, plotted month by month.

That climbing line carries a spike every twelve months, some years wider than others, but always there. This lesson builds a model that expects that spike, instead of being surprised by it every December.

The seasonal signature in the ACF, at lags 12, 24 and 36

Building the seasonal part of a model starts with seeing the season clearly first. Build the running example, and split it the way the whole lesson will use it: a training window to fit on, and a held-out test window to check a forecast against later.

RInteractive R
# Build the Victoria cafes, restaurants and takeaway turnover series, and split into train and test library(tsibbledata) library(tsibble) library(dplyr) library(feasts) retail <- tsibbledata::aus_retail |> filter(State == "Victoria", Industry == "Cafes, restaurants and takeaway food services") |> select(Month, Turnover) nrow(retail) #> [1] 441 train <- retail |> filter(Month <= yearmonth("2016 Dec")) test <- retail |> filter(Month > yearmonth("2016 Dec")) nrow(train) #> [1] 417 nrow(test) #> [1] 24

  

train holds 417 months, April 1982 through December 2016. test holds the last 24 months, January 2017 through December 2018, real turnover figures set aside so a fitted model can be checked later against numbers it never saw.

The tool for reading a repeating pattern like this one is the ACF, short for autocorrelation function. Its value at a given lag is the correlation between the series and its own value that many months earlier. A value near 1 at some lag means the series looks almost the same shape that many months apart.

Compute the ACF of the training series up to lag 36, three full years, and read off the values at exactly 12, 24 and 36 months.

RInteractive R
# Compute the ACF of the training series up to lag 36 acf_raw <- train |> ACF(Turnover, lag_max = 36) round(acf_raw$acf[c(12, 24, 36)], 3) #> [1] 0.904 0.817 0.731 round(acf_raw$acf, 3) #> [1] 0.981 0.969 0.962 0.954 0.944 0.932 0.930 0.924 0.916 0.907 0.904 0.904 #> [13] 0.887 0.876 0.870 0.864 0.855 0.844 0.842 0.837 0.830 0.821 0.818 0.817 #> [25] 0.802 0.792 0.786 0.779 0.770 0.759 0.757 0.751 0.743 0.734 0.731 0.731

  

Every one of the 36 values is high, since a series with a strong trend correlates with itself at almost any lag. But look specifically at 12, 24 and 36: 0.904, 0.817 and 0.731, decaying only slowly, and each one landing at an exact multiple of 12. That is the seasonal signature, a series that keeps echoing itself at the period and every whole multiple of it, long after nearby lags would normally have faded.

See the shape of it as a bar chart, with the three seasonal lags marked apart from the rest.

An ordinary ARIMA(p, d, q) has no way to produce a bar this tall out at lag 36. Its AR and MA terms only reach back p or q lags, and nobody sets p or q to 36 just to explain one distant echo. Reaching lag 12, 24 and 36 on purpose, without needing 36 ordinary terms to do it, is exactly what the seasonal part of a SARIMA model is built for.