Overview

In class on Thursday we learned how perform diagnostic tests on an ARIMA model. We also introduced how to model seasonality using an ARIMA model step-by-step.

Box Pierce and Ljung-Box Statistic

The Box Pierce and Ljung-Box statistics allow us to perform residual analysis for the ARIMA model. When an ARIMA model is adequately modeling the data, then the residuals are small. The test stat is then also small. So we prefer a large p-value. The null hypothesis is that the model adequately models the data, and the alternative is that the model inadequately models the data.

The Box Pierce statistic can be used, but in class we covered the Ljung-Box statistic because it is better than the Box Pierce statistic.

Seasonal ARIMA Model Building Steps

I will use the souvenir sales from a gift shop to demonstrate how to build a seasonal ARIMA model.

souvenir <- scan("http://robjhyndman.com/tsdldata/data/fancy.dat")
souvenir <- ts(souvenir, frequency=12, start=c(1987,1))
library(forecast)
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)
head(souvenir)
##          Jan     Feb     Mar     Apr     May     Jun
## 1987 1664.81 2397.53 2840.71 3547.29 3752.96 3714.74

Step 1: Plotting

First we need to visualize the data to check if there is obvious seasonality

plot.ts(souvenir)

A similar pattern each year indicates that we have seasonality. Because the size of the peaks is not constant, we will try a log transform on the data.

logsouvenir<-log(souvenir)
plot.ts(logsouvenir)

The size of the peaks now appear constant, so we will continue to use the log transformed data. We can better check for seasonality by decomposing the data.

plot(decompose(logsouvenir))

The seasonal plot, which is the third plot, indicates that we have a seasonal trend.

Step 2: Differences

The second step is to estimate what appropriate values are for d and D. If there is seasonality and no trend, we will have \(D\neq 0\) and \(d=0\). If there is no seasonality and trend then we will have \(D=0\) and \(d\neq0\). If there is both seasonality and trend, we start with \(D>0\) then check to see if we still require \(d>0\) for any remaining trend.

We will use a lag of 12 (S=12) for differences because our seasonality is yearly and our data is in months. First we will try to take care of seasonality, because we previously identified it.

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

We don’t see that D=1 also improved the trend. We can better determine this with a hypothesis test.

Hypothesis Testing for Trend

We can use the Augmented Dickey-Fuller Test (adf test) to see if there is a trend remaining. A small p-value allows us to reject the null and accept that the plot is stationary.

adf.test(diff1)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff1
## Dickey-Fuller = -2.6959, Lag order = 4, p-value = 0.292
## alternative hypothesis: stationary

We can see that the p-value is .292, which is large. We fail to say that the plot is stationary, so there remains a trend. Now we can move on to examining the d value.

We again use the difference data from above, but lag is set to 1 (S=1) becasue nonseasonal differences are point to point.

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

We now check to see if the trend is removed.

adf.test(diff2)
## Warning in adf.test(diff2): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff2
## Dickey-Fuller = -5.3558, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

The p-value is 0.01, so we can reject the null and accept the alternative that the data is stationary and that there is no remaining trend. We now have our D and d values, which are both 1.

Step 3: Choosing p, q and P, Q

We can choose values for p,q and P,Q by examining the ACF and PACF. We can determine q and Q from ACF because is is the autocorrelation function, which tells us the correlation between two points given a certain lag. PACF is the partial autocorrelation function. It focuses on shorter lag. We can determine p and P with PACF.

Acf(diff2)

For q, we examine spikes for low lages. A spike is important if it crosses the dotted blue line. We have one spike for low lage so we will try q=1. For Q, we check if we have spikes for lags that are a multiple of 12. There’s one that’s close to the blue line, but we will move forward with Q=0.

Pacf(diff2)

For p, we again look at low lag values. We have one spike over the line, and one close. We will choose p=1. For P, we again look at multiples of 12 and will use P=0.

Step 4: Build the Seasonal ARIMA Model

We can now use our p, d, q and P, D, Q values to build an ARIMA model. Our model is ARIMA (1,1,1)x(0,1,0) 12

mod<-arima(logsouvenir, order=c(1,1,1), season=list(order=c(0,1,0), period=12))

Step 5: Examine the Model

The last step is to perform hypothesis tests and use diagostics to see if our model is adequate or can be improved.

Box.test(mod$residuals, type="Ljung")
## 
##  Box-Ljung test
## 
## data:  mod$residuals
## X-squared = 0.10818, df = 1, p-value = 0.7422

For the Ljung Box Test we have a large p-value of .7422, so we don’t have evidence that the model is inadequate. We can check each coefficient by using the test stat in the summary.

summary(mod)
## 
## Call:
## arima(x = logsouvenir, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), 
##     period = 12))
## 
## Coefficients:
##          ar1      ma1
##       0.2020  -0.7403
## s.e.  0.2614   0.2064
## 
## sigma^2 estimated as 0.03695:  log likelihood = 16.09,  aic = -26.17
## 
## Training set error measures:
##                        ME      RMSE       MAE        MPE     MAPE
## Training set -0.007853116 0.1767525 0.1334457 -0.1052367 1.455206
##                   MASE        ACF1
## Training set 0.3309956 -0.03525604

If a test statistic is less than 2, we might want to try removing it. ar1 has a test stat less than two, so we can check the effects of removing it. That is, our model is now ARIMA (0,1,1)x(0,1,0) 12

mod2<-arima(logsouvenir, order = c(0,1,1), season = list( order = c( 0,1,0), period=12))
Box.test(mod2$residuals, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  mod2$residuals
## X-squared = 0.064518, df = 1, p-value = 0.7995

The p-value is again large, so our model is not inadequate. Let’s again check the summary.

summary(mod2)
## 
## Call:
## arima(x = logsouvenir, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0), 
##     period = 12))
## 
## Coefficients:
##           ma1
##       -0.5653
## s.e.   0.1227
## 
## sigma^2 estimated as 0.03722:  log likelihood = 15.89,  aic = -27.79
## 
## Training set error measures:
##                        ME     RMSE      MAE         MPE     MAPE      MASE
## Training set -0.006290552 0.177391 0.132399 -0.08573452 1.443332 0.3283994
##                   ACF1
## Training set 0.0272264

The test statistic for ma1 is close to 2, so we will try increasing q to 2. Our model is now ARIMA (0,1,2) x (0, 1, 0) 12.

mod3<-mod3 <- arima(logsouvenir, order = c(0,1,2), season = list( order = c( 0,1,0), period=12))
mod3
## 
## Call:
## arima(x = logsouvenir, order = c(0, 1, 2), seasonal = list(order = c(0, 1, 0), 
##     period = 12))
## 
## Coefficients:
##           ma1      ma2
##       -0.5637  -0.0350
## s.e.   0.1156   0.1197
## 
## sigma^2 estimated as 0.03716:  log likelihood = 15.93,  aic = -25.87

The ma2 test stat is around 0, so we will move forward with the previous model. For comparison we could check the model AICs, but we instead will compare our models with R’s auto ARIMA model.

Step 6: Use R to Build a Model

Let’s see what the auto ARIMA function in R thinks is the best model.

automod<-auto.arima(logsouvenir)
automod
## Series: logsouvenir 
## ARIMA(2,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##          ar1     ar2     sar1   drift
##       0.3493  0.3602  -0.3278  0.0247
## s.e.  0.1086  0.1159   0.1334  0.0044
## 
## sigma^2 estimated as 0.03182:  log likelihood=23.04
## AIC=-36.09   AICc=-35.18   BIC=-24.71

R returned ARIMA (2,0,0) x (1,1,0) 12. The best model does not include differences for nonseasonality. q=0 is also best. Let’s now check the AIC values between our models.

AIC(automod); AIC(mod); AIC(mod2); AIC(mod3)
## [1] -36.08965
## [1] -26.17089
## [1] -27.78538
## [1] -25.86905

As expected, the automod model has the lowest AIC. We now know how to estimate an ARIMA model, but since building the model is subject to our best judgement, we can see how letting R decide yields the best model.