# Choosing daily seasonal ARIMA model order vs. auto.arima

I have the following output_ts that looks like this:

When I run unit root tests, I receive the following output seeming to indicate that my series is trend stationary:

> adf.test(output_ts)

Augmented Dickey-Fuller Test

data:  output_ts
Dickey-Fuller = -7.3922, Lag order = 14, p-value = 0.01
alternative hypothesis: stationary

Warning message:
In adf.test(output_ts) : p-value smaller than printed p-value
> kpss.test(output_ts, null = c("T"))

KPSS Test for Trend Stationarity

data:  output_ts
KPSS Trend = 0.21269, Truncation lag parameter = 12, p-value = 0.01124

> nsdiffs(output_ts)
[1] 1


This seems to indicate that the suggested order of differencing is d = 1. However, I also know that my data is daily and has a seasonality of S = 7. This is confirmed by spectral decomposition of the time series frequency:

> findfrequency(output_ts)
[1] 7


ACF and PACF plots of output_ts look like this:

Which seems to confirm that this time series has a seasonality of S = 7. Using the guidance given in this tutorial: https://onlinecourses.science.psu.edu/stat510/node/67/

I then applied a seasonal difference of 7 with the following results:

> output_diff7_ts = diff(output_ts, 7)
> plot(output_diff7_ts)
> acf2(output_diff7_ts)


And ACF and PACF plots as follows:

Based on the PACF plot and the spikes at lag 1, it seems like I need to add a non-seasonal AR(1) term. Based on the PACF plot and the spikes at lags 7, 14, 21, etc., it seems like I also need to add a seasonal MA(1) term.

So the final model should look something like: ARIMA (1,0,0)(0,1,1)[7] which I can fit using Arima from Prof. Hyndman's forecast package and check the residuals as follows:

> output_model = Arima(output_ts, order = c(1,0,0), seasonal=list(order=c(0,1,1),period=7))
> checkresiduals(output_model)


Which gives me the following:

    Ljung-Box test

data:  Residuals from ARIMA(1,0,0)(0,1,1)[7]
Q* = 181.25, df = 12, p-value < 2.2e-16

Model df: 2.   Total lags used: 14


Clearly not ideal. Not only is the ACF pretty terrible looking, the Box test scores are not good.

So then I tried fitting output_ts simply using auto.arima:

> output_model = auto.arima(output_ts, ic = "aic")
> summary(output_model)
Series: output_ts
ARIMA(0,0,1)(0,1,1)[7] with drift

Coefficients:
ma1     sma1   drift
0.2831  -0.8924  0.0539
s.e.  0.0162   0.0124  0.0142

sigma^2 estimated as 1558:  log likelihood=-15959.89
AIC=31927.78   AICc=31927.79   BIC=31951.98

Training set error measures:
ME     RMSE      MAE        MPE     MAPE      MASE       ACF1
Training set 0.01322627 39.40545 29.10605 -0.6122858 5.791276 0.7795865 0.03899792


Which seems to indicate that auto.arima is uncovering the same model that we originally selected manually. The residuals look similarly terrible:

> checkresiduals(output_model)

Ljung-Box test

data:  Residuals from ARIMA(0,0,1)(0,1,1)[7] with drift
Q* = 369.05, df = 11, p-value < 2.2e-16

Model df: 3.   Total lags used: 14


Which leads me to hope that someone might be able to help with a few questions:

1. Is the seasonality of the daily data (period = 7) somehow clashing with the auto.arima restrictions on the order of ARIMA model? According to the documentation, the maximum order checked by auto.arima is (5,2,5)x(2,1,2).
2. How would you suggest I go about working towards a better seasonal arima model?

Thanks!

• Post your data to dropbox.com. Specify the beginning date. – Tom Reilly Dec 17 '18 at 20:13
• auto.arima (without x's) definitely clashes with day-of-the-week effects ...which is why one needs to combine both deterministic effects ( like day-of-the-week , week-of-the-month , day-of-the-month, month-of-the-year, week-of-the-year, lead and lag effects around holidays AND arima structure AND of course pulse effects. – IrishStat Dec 17 '18 at 20:29
• File uploaded here: gofile.io/?c=g2pnfm – user50710 Dec 17 '18 at 21:27