ARIMA Testing and Diagnostics

In class, we went over how to perform hypothesis tests and the diagnostics on ARIMA models. In addition, we learned how to model seasonal ARIMAs.

Hypothesis Tests

To continue from last class, we learned how to do hypothesis tests on nonseasonal ARIMA models. Our null and alternative hypotheses are as follows: \[H_0: \theta=0 \text{ or } \phi=0\] \[H_A: \theta\neq0 \text{ or } \phi\neq0\] The test statistics for \(\theta\) and \(\phi\) respectively are: \(\frac{\hat{\theta}}{SE(\hat{\theta})}\) and \(\frac{\hat{\phi}}{SE(\hat{\phi})}\), where they follow a T distribution with n-(number of \(\theta\)s + number of \(\phi\)s) degrees of freedom.

ARIMA Diagnostics

Next, we covered diagnostics on ARIMA Models. To do so, we can do a residual analysis. If the ARIMA model is doing a good job, then the residuals and the test statistic that we calculate should be small. Our hypotheses for the diagnostic test are as follows: \[H_0:\text{ The model adequately models our time series data}\] \[H_A: \text{The model inadequately models our time series data}\] There are two test statistics that we can use to conduct our diagnostic test: 1. Box Pierce Stat 2. Ljung-Box Stat The Box Pierce stat is the original test statistic, whereas the Ljung-Box stat is newer, but also considered better. For this reason, we prefer to use the Ljung-Box test statistic for our diagnostic test. A large p-value while conducting this test implies that the model is adequate, while a small p-value might suggest that the model is inadequate.

Seasonal ARIMA Models

For this type of ARIMA model, we introduce a new component to our model: the seasonal component. This might represent the change in response based on the time of year. The notation takes the form of ARIMA (p,d,q) x (P,D,Q) S , where p,d, and q are the nonseasonal components while P,D,Q, and S are the seasonal components. In this case, P,D, and Q represent the same components as their nonseasonal counterparts p,d,and q respectively. The P is the order of the autoregressive component, the D is the order of integration (seasonal difference), and Q is the order of the moving average component. While d serves to remove the upward/downward trend, D serves to remove the seasonal trend. Finally, S is the time span of the seasonal repeating pattern. For example, if the repeating pattern is every month, then S=12, whereas if it is ever quarter, then S=4.

General Seasonal ARIMA Steps

Next, we want to look at how to put together an ARIMA model with a seasonal component. There are 5 main steps: 1. Plot the data to see if you have seasonality 2. Do differences to identify D and d 3. Examine autocorrelation function (ACF) and partial autocorrelation function (PACF) 4. Estimate the model 5. Perform diagnostics, hypothesis tests, and use AIC to compare

Example

Here is an example of how to walk through the general seasonal ARIMA steps. I used the AirPassengers data set to show how building this model works. First, we need the following packages:

library(forecast)
## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone 'zone/tz/2018c.
## 1.0/zoneinfo/America/Chicago'
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(tseries)
library(timeSeries)
## Loading required package: timeDate
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
library(forecast)
library(xts)

Step 1

The first step is to plot the data to see if there is any seasonality.

data(AirPassengers)
plot(AirPassengers)

From this plot, we can se that there is seasonality, as well as an upward trend. Additionally, it seems as though the seasonality is making more of an impact as time goes on since the vertical spread of the spikes is increasing. We can fix this by changing the mean function with a log transformation.

plot(log(AirPassengers))

Now the vertical spread is much more evened out, which makes it easier to look at the seasonal trend. We can determine what the seasonal time span from the data set

head(AirPassengers)
##      Jan Feb Mar Apr May Jun
## 1949 112 118 132 129 121 135

Evidently, this data is based on a monthly time frame, so S=12.

Step 2

Now we want to do the differences to identify D and d. To do so, we are going to first plot the differences of log(AirPassengers) by using the diff function. Then, we want to conduct a hypothesis test to see whether the differences are significant. This can be done by using the adf.test function on the differences. The null hypothesis states that the differences are not stationary, while the alternative states that the differences are stationary.

diff1 <- diff(log(AirPassengers), lag = 12)
plot.ts(diff1) 

adf.test(diff1)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1
## Dickey-Fuller = -3.1899, Lag order = 5, p-value = 0.09265
## alternative hypothesis: stationary

Our resulting p-value is greater than 0.05 (p=0.09265), meaning that it is not sufficient evidence to reject the null, and therefore we say that the differences are not stationary. This means that we want to take the difference again.

diff2 <- diff(diff1, lag = 12)
plot.ts(diff2) 

adf.test(diff2)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff2
## Dickey-Fuller = -3.5215, Lag order = 4, p-value = 0.04333
## alternative hypothesis: stationary

After looking at the output of the adf.test function, we can conclude that the p-value is less than 0.05 (p=0.04333). This means that we can reject the null in the favor of the alternative, and say that the differences are stationary. So, we want to then use D=2, since this was the second order of differences for seasonality. Also, there doesn’t appear to be a general trend, so we want to use d=0.

Step 3

We need to then examine ACF and PACF in order to choose preliminary values for p,q, P, and Q. The ACF tells you the correlation of two points that are separated by a specific lag, which is the number of units time difference. The PACF is like the ACF, but controls for values of a shorter lag. In order to determine the preliminary values for q and Q, we want to use the autocorrelation function (ACF). For q, we want to look at the spikes from low lags, while with Q, we want to look at the spikes from lags that are mutliples of S (which is 12 in this case). In order to determine the preliminary values for p and P, we want to use the partial autocorrelation function (PACF). For p, we want to look at spikes from low lags, while for P we want to look at the spikes from lags that are multiples of S (which again is 12).

To look at the ACF, we use the Acf function on the second order of seasonal differences that we obtained from above.

Acf(diff2)

The first 5 lags appear to be significant (above the dotted line), and the 6th lag is borderline significant. Since we can’t definitively say what the value of q is from this, we can then say that \(q\geq1\). In terms of multiples of 12, lag 12 is quite significant, but lag 24 is not, so it appears that Q=1.

To look at PACF, we use the Pacf function on diff2.

Pacf(diff2)

The first lag is definitely significant, and the second one is just slightly above the level of significance, so it looks like p=1 or p=2. Looking at multiples of 12, it appears that lag 12 is significant, and lag 24 is slightly significant, meaning P=1 or P=2. We can then say that \(p\geq1\) and \(P\geq1\).

Step 4

Our fourth step is to estimate the model. We are going to do this based on what we found with our ACF and PACF. We are going to start by building the simplest model given our findings. This would mean that our starting model should be ARIMA (1,0,1) x (1,2,1) 12.

mod1<-arima(log(AirPassengers),order=c(1,0,1),season=list(order=c(1,2,1),period=12))

Step 5

Finally, our last step is to run diagnostics, hypothesis tests, and to use AIC to compare. We are going to try and improve our model that we found in the previous step. First, let’s take a look at the summary of our first model

summary(mod1)
## 
## Call:
## arima(x = log(AirPassengers), order = c(1, 0, 1), seasonal = list(order = c(1, 
##     2, 1), period = 12))
## 
## Coefficients:
##          ar1      ma1     sar1     sma1
##       0.9303  -0.4175  -0.4505  -1.0000
## s.e.  0.0448   0.0991   0.0861   0.1038
## 
## sigma^2 estimated as 0.001415:  log likelihood = 203.55,  aic = -397.1
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE
## Training set -0.001134586 0.03439361 0.02451314 -0.01830164 0.4413303
##                   MASE        ACF1
## Training set 0.2706075 -0.02443055

This summary gives us the coefficients for our different components. Let’s take a look at the Box-Ljung test to see if the model is adequate.

Box.test(mod1$residuals, type= "Ljung")
## 
##  Box-Ljung test
## 
## data:  mod1$residuals
## X-squared = 0.08775, df = 1, p-value = 0.7671

The p-value is 0.7671, which means that there is not sufficient evidence to say that the model is inadequate. Now, we want to take a look and see what the values of the test stats are for each coefficient. If its absolute value is \(\leq 2\), then we want to consider removing it from the model.

ar1<-0.9303/.0448
abs(ar1)
## [1] 20.76562
ma1<--0.4175/0.0991
abs(ma1)
## [1] 4.212916
sar1<--0.4505/0.0861
abs(sar1)
## [1] 5.232288
sma1<--1/0.1038
abs(sma1)
## [1] 9.633911

All of them are greater than 2, so there is no need to remove any of them from the model. So, the final model that we built is ARIMA (1,0,1) x (1,2,1) 12. Now, let’s take a look at the ARIMA model that R builds for us with the auto.arima function.

automod<-auto.arima(log(AirPassengers))
automod
## Series: log(AirPassengers) 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.4018  -0.5569
## s.e.   0.0896   0.0731
## 
## sigma^2 estimated as 0.001371:  log likelihood=244.7
## AIC=-483.4   AICc=-483.21   BIC=-474.77

R gives us the model ARIMA (0,1,1) x (0,1,1) 12. We want to use AIC to determine which model would be better to use. The model R gives us is simpler, since the highest order is 1 instead of 2. So, if the AIC of the model that we built is at least 10 less than the one R gives us, then we use our model. However, if not, then we want to use the one R gives us.

AIC(mod1)
## [1] -397.1033
AIC(automod)
## [1] -483.3991

The AIC for automod is much smaller than the one given for the model we built, so we are going to use the one built by R.