In this lab, you will fit data with the auto.arima() function in the forecast package.

Set-up

If you have not loaded the setup package with the Greek landings data, do so now. You would have needed this for the TV Regression lab, so you should not need to do this again.

library(devtools)
devtools::install_github("RVerse-Tutorials/RWorkflowsetup")

We will use the data before 1987 for the lab and you will look at the later data in the problems.

load("landings.RData")
landings$log.metric.tons <- log(landings$metric.tons)
landings <- subset(landings, Year <= 1987)
anchovy <- subset(landings, Species=="Anchovy")
sardine <- subset(landings, Species=="Sardine")

Load the necessary packages.

library(forecast)
library(tseries)

Fitting to simulated data

To see info on arima.sim(), type ?arima.sim and ?arima.

Fit to AR(2) data:

\[x_t = 0.8 x_{t-1} + 0.1 x_{t-2} + e_t\]

a1 = arima.sim(n=100, model=list(ar=c(.8,.1)))
auto.arima(a1, seasonal=FALSE, max.d=0)
## Series: a1 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.8681  1.2128
## s.e.  0.0505  0.7470
## 
## sigma^2 estimated as 1.101:  log likelihood=-146.42
## AIC=298.83   AICc=299.08   BIC=306.65

Fit to AR(1) data:

\[x_t = 0.8 x_{t-1} + e_t\]

a1 = arima.sim(n=100, model=list(ar=c(.8)))
auto.arima(a1, seasonal=FALSE, max.d=0)
## Series: a1 
## ARIMA(1,0,0) with zero mean 
## 
## Coefficients:
##          ar1
##       0.7987
## s.e.  0.0615
## 
## sigma^2 estimated as 1.142:  log likelihood=-148.54
## AIC=301.08   AICc=301.2   BIC=306.29

\[x_t = 0.8 x_{t-1} + e_t + 0.8 e_{t-1} + 0.2 e_{t-2}\]

Fit to ARMA(1,2) data:

We will up the number of data points to 1000 because models with a MA component take a lot of data to estimate. Models with MA(>1) are not very practical for fisheries data for that reason.

a1 = arima.sim(n=1000, model=list(ar=c(0.8), ma=c(0.8, 0.2)))
auto.arima(a1, seasonal=FALSE, max.d=0)
## Series: a1 
## ARIMA(1,0,2) with zero mean 
## 
## Coefficients:
##          ar1     ma1     ma2
##       0.8703  0.7218  0.1614
## s.e.  0.0178  0.0354  0.0353
## 
## sigma^2 estimated as 1.027:  log likelihood=-1432.35
## AIC=2872.7   AICc=2872.74   BIC=2892.33

Fitting to 100 simulated data sets

save.fits = rep(NA,100)
for(i in 1:100){
  a1 = arima.sim(n=100, model=list(ar=c(.8,.1)))
  fit = auto.arima(a1, seasonal=FALSE, max.d=0, max.q=0)
  save.fits[i] = paste0(fit$arma[1], "-", fit$arma[2])
}
table(save.fits)
## save.fits
## 1-0 2-0 3-0 
##  76  21   3

These functions work for data with missing values

Create some AR(2) data and then add missing values (NA).

a1 = arima.sim(n=100, model=list(ar=c(.8,.1)))
a1[sample(100,50)]=NA
plot(a1, type="l")
title("many missing values")

As always to understand what the code is doing, run pieces of it and look at the output. Use ?sample (say) to read about functions you don’t know.

Fit to the data with missing values. You need to give auto.arima() a little help by telling it not to include models with differencing. So we set max.d=0. You can remove that and see what happens.

auto.arima(a1, seasonal=FALSE, max.d=0)
## Series: a1 
## ARIMA(1,0,0) with zero mean 
## 
## Coefficients:
##          ar1
##       0.8236
## s.e.  0.0686
## 
## sigma^2 estimated as 1.229:  log likelihood=-84.37
## AIC=172.74   AICc=172.87   BIC=177.93

auto.arima() is nice because it selects the model order for you. But if you know what model you want to fit, you can use Arima():

Arima(a1, order = c(2,0,0))
## Series: a1 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2    mean
##       0.8277  -0.0084  0.4385
## s.e.  0.1655   0.1703  0.6032
## 
## sigma^2 estimated as 1.266:  log likelihood=-84.1
## AIC=176.19   AICc=176.62   BIC=186.57

Fit to the Anchovy landing data

fit <- auto.arima(anchovy$log.metric.tons)
fit
## Series: anchovy$log.metric.tons 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##           ma1   drift
##       -0.5731  0.0641
## s.e.   0.1610  0.0173
## 
## sigma^2 estimated as 0.03583:  log likelihood=6.5
## AIC=-6.99   AICc=-5.73   BIC=-3.58

Compare to Table 8 in Stergiou and Christou. Note arima() writes a MA model like:

\[x_t = e_t + b_1 e_{t-1} + b_2 e_{t-2}\]

while Stergiou and Christou use the following (also common) notation:

\[x_t = e_t - \theta_1 e_{t-1} - \theta_2 e_{t-2}\]

so the MA parameters reported by auto.arima() will be NEGATIVE of that reported in Stergiou and Christou. Note, in Stergiou and Christou, the model is written in backshift notation on page 112. To see the model as the equation above, I translated from backshift to non-backshift notation.

Here are the values for anchovy in Table 8.

ARIMA Model \(\theta_1\) drift (c) R\(^2\) BIC LB
(0,1,1) 0.563 0.064 0.83 1775 5.4

The model selected by auto.arima() is the same. The model selection criteria is similar in spirit but not in detail. Stergiou and Christou used a combination of model fit and BIC (a model selection criteria), which auto.arima() uses AICc by default (you can specifiy BIC).

The parameter estimates are similar.

The BIC value is different, that is expected because there is a constant that is often left off.

The Ljung-Box test is not significant. They don’t give the p-value, but state in the table legend that it is not significant. That is good. That means we don’t reject the null hypothesis of independence.

Here is how to do the test that they did:

res <- resid(fit)
Box.test(res, type="Ljung-Box", lag=12, fitdf=2)
## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 5.3725, df = 10, p-value = 0.8649

Here is how to do the default test with checkresiduals():

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1) with drift
## Q* = 1.4883, df = 3, p-value = 0.685
## 
## Model df: 2.   Total lags used: 5

Problems

  1. Fit to the sardine data.
  • Use auto.arima()
  • Compare to Table 8 in Stergiou and Christou. Do you get the same result?

I am not sure how Stergiou and Christou selected the amount of differencing to use. auto.arima() uses the KPSS test by default. Let’s try the Augmented Dickey-Fuller test instead. Specify this by passing in test="adf".

auto.arima(sardine$log.metric.tons, test="adf")
  • Is the selected model the same?
  1. Now fit to the anchovy data 1989-2007
load("landings.RData")
landings$log.metric.tons <- log(landings$metric.tons)
anchovy <- subset(landings, Species=="Anchovy" & Year>=1989)
  • First plot the data.
  • Fit with auto.arima(). Do you get the same model? You don’t. ARIMA(0,1,1) is very different than ARIMA(0,0,1).
  • Force a fit of ARIMA(0,1,1) with drift using Arima(). This is the best fit model with auto.arima().
Arima(anchovy$log.metric.tons, order=c(0,1,1), include.drift=TRUE)
## Series: anchovy$log.metric.tons 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##           ma1    drift
##       -1.0000  -0.0069
## s.e.   0.1959   0.0070
## 
## sigma^2 estimated as 0.03183:  log likelihood=5.08
## AIC=-4.15   AICc=-2.44   BIC=-1.48
  • Compare the estimates of ma1.