# Introducing seasonality to an ARIMA model in R

Using the forecast package, I am constructing a dynamic regression / ARIMA model to use X to predict Y (data below). Here is my model:

fit <-auto.arima(y,xreg=x,seasonal=c(0,1,1)[4])
summary(fit)


Question - I'm not sure I am using the seasonal portion of code correctly:

seasonal=c(0,1,1)[4]


I see no difference in the forecasts using the above seasonality code and not including the seasonality code at all.

Question: How do I specify seasonality in this instance? Is there a way the arima code can automatically detect correct seasonality and specify it in code?

Data:

qtrs     y         x
1Q14    41.7    16963482
2Q14    26      8857757
3Q14    22.4    6308120
4Q14    17.6    4872037
1Q15    47.7    20304915
2Q15    30.7    8800350
3Q15    25.6    8065914
4Q15    17.4    8.41E+06
1Q16    59.4    39076150
2Q16    35.6    23348752
3Q16    33.5    26721456
4Q16    23.9    17934787
1Q17    81.7    65821538
2Q17    44.8    35845803
3Q17    42.1    34753415
4Q17    28.5    21449909
1Q18    90      84207978


or also

> dput(y)
c(41.7, 26, 22.4, 17.6, 47.7, 30.7, 25.6, 17.4, 59.4, 35.6, 33.5,
23.9, 81.7, 44.8, 42.1)
dput(x)
c(16963482, 8857757, 6308120, 4872037, 20304915, 8800350, 8065914,
8410000, 39076150, 23348752, 26721456, 17934787, 65821538, 35845803,
34753415)


When I use:

   fit <-auto.arima(y,xreg=x,seasonal=True)


it returns:

  summary(fit)
Series: y
Regression with ARIMA(0,1,0) errors

Coefficients:
x
0e+00
s.e.  3e-04

sigma^2 estimated as 24.53:  log likelihood=-41.75
AIC=87.49   AICc=88.58   BIC=88.77

Training set error measures:
ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
Training set -1.521487 4.611119 3.659424 -7.34176 12.28209 0.1975017 -0.2619158


But when I use the below, it returns the same as above. fit <-auto.arima(y,xreg=x,seasonal=c(0,1,1)[4])

Why is there no seasonal variables in this model when there clearly is seasonality? Why no difference in the two models (seasonality specified and not specified)?

The seasonal parameter expects a simple Boolean input (see ?auto.arima). What you are providing is c(0,1,1)[4], which happens to be a well-formed R expression, namely the fourth entry in the vector c(0,1,1) of length three, or

c(0,1,1)[4]
# NA


And this NA is then coerced to a Boolean value, namely

isTRUE(NA)
# FALSE


So what you did was an elaborate way of specifying a non-seasonal model. Try seasonal=TRUE instead.

(And don't be surprised if the automatic model selection decides to use a non-seasonal model, anyway. If you want to force seasonality, this may be helpful: Seasonality not taken account of in auto.arima().)

EDIT: with the data you posted, the first thing to notice is that x and y are vectors, not time series objects with an associated seasonal frequency. As above, auto.arima does not derive the seasonal period from the seasonal parameter (which is supposed to be, and will be coerced to, a simple Boolean), but from the frequency of the focal time series. So the first thing to do is to identify your quarterly data as such:

y <- ts(c(41.7, 26, 22.4, 17.6, 47.7, 30.7, 25.6, 17.4, 59.4, 35.6, 33.5,
23.9, 81.7, 44.8, 42.1), frequency=4)
x <- ts(c(16963482, 8857757, 6308120, 4872037, 20304915, 8800350, 8065914,
8410000, 39076150, 23348752, 26721456, 17934787, 65821538, 35845803,
34753415),frequency=4)


Next, even if your y is seasonal, this does not necessarily matter. auto.arima() with an xreg parameter fits a regression with ARIMA errors (not ARIMAX), so what matters is the seasonality of the residuals of this regression of y on x. Fortunately, these residuals are also clearly seasonal, as are their first differences (note how we again need to convert the residuals to ts objects):

library(forecast)
seasonplot(ts(residuals(lm(y~x)),frequency=4))
seasonplot(ts(diff(residuals(lm(y~x))),frequency=4))


Interestingly enough, auto.arima() fits a non-seasonal model, as you noticed, even if we don't specify seasonal=FALSE. Then again, if we call auto.arima() on the residuals from the regression y~x, we get a seasonal model, and so do we if we run auto.arima(y,xreg=x,D=1) with enforced seasonality as per this thread - but the seasonal models are different.

> auto.arima(y,xreg=x)
Series: y
Regression with ARIMA(0,1,0) errors

Coefficients:
xreg
0e+00
s.e.  3e-04

sigma^2 estimated as 24.53:  log likelihood=-41.75
AIC=87.49   AICc=88.58   BIC=88.77

> auto.arima(ts(residuals(lm(y~x)),frequency=4))
Series: ts(residuals(lm(y ~ x)), frequency = 4)
ARIMA(0,1,0)(1,0,0)[4]

Coefficients:
sar1
0.9074
s.e.  0.0694

sigma^2 estimated as 12.65:  log likelihood=-40.58
AIC=85.16   AICc=86.25   BIC=86.44

> auto.arima(y,xreg=x,D=1)
Series: y
Regression with ARIMA(0,0,0)(0,1,0)[4] errors

Coefficients:
xreg
0e+00
s.e.  3e-04

sigma^2 estimated as 12.29:  log likelihood=-28.88
AIC=61.76   AICc=63.26   BIC=62.56


Overall, it looks like there is something undocumented going on here - auto.arima() may by default not use seasonal models if there is an xreg present. You may want to ping the author Rob Hyndman about this (best to point him to this thread, he may chime in). Similarly, that the last two calls resulted in different ARIMA models may be due to internals of auto.arima().

It looks like your best bet may be to run the regression outside auto.arima() and then subject the residuals to ARIMA modeling, as I did here. Alternatively, if you can identify whether your residual series is seasonal or not, you can enforce seasonality using D as above.

• I was curious if the added detail above changes your view. Any idea why the same model comes back whether seasonal=True or seasonal=c(0,1,1)[4]? There is clearly seasonality in this series - but both specifications only return a model with a single x coefficient. Thanks!
– ZJAY
Feb 19, 2018 at 14:58
• thanks again, very helpful. I'm reaching out to Rob. Dumb question - Regression with ARIMA(0,0,0)(0,1,0)[4] errors - could you translate this into English? What kind of seasonality was used in this model? Which quarter received a "seasonal pulse." Thanks.
– ZJAY
Feb 20, 2018 at 2:02
• ARIMA(0,0,0)(0,1,0)[4] is actually an extremely simple model. It says that the first seasonal difference (that's the 1 and the [4]), is white noise, $e_t-e_{t-4}=\epsilon_t$ with $\epsilon_t\sim N(0,\sigma^2)$. Note that I'm calling the time series we are looking at $e_t$, because it's the residuals from the regression y~x. Interpret ARIMA models in plain english may be helpful, or the section on seasonal ARIMA in the Hyndman & Athanasopoulos textbook. Feb 20, 2018 at 9:57

Are you sure that's the best seasonal set-up? What I mean is that if you use auto.arima to find the "best" model, you might as use it to find the best seasonal configuration as well? If you leave it out, auto.arima will test automatically for seasonality

• Does auto.arima automatically find the best seasonal configuration? Do I have to specify anything, like fit <-auto.arima(y,xreg=x,seasonal=True) ? Also, how do I see what seasonality auto.arima selects? When I run summary(fit) with no specified seasonality, I only see a x coefficient for the model, no seasonality?
– ZJAY
Feb 19, 2018 at 14:21
• No. You don't have to specify anything. When you have run your model, you enter the name of your model (in your code "fit") in the command line and you will be able to see the model or use fit\$model Feb 19, 2018 at 14:27