```{r ts-setup, include=FALSE, purl=FALSE}
knitr::opts_knit$set(unnamed.chunk.label = "labboxjenkins-")
```

# Fitting ARIMA models in R  {#chap-labboxjenkins-}
\chaptermark{Fitting ARIMA models}

In this chapter, you will practice doing fitting an ARIMA model to catch data using the Box-Jenkins method.  You will prepare simple forecasts using the **forecast** package.  

### Data

Load catch landings from Greek waters [landings.RData](landings.RData) and the Chinook landings  [chinook.RData](chinook.RData) in Washington data.

```{r read_data}
load("landings.RData")
load("chinook.RData")
```

Ensure you have the necessary packages.

```{r load_packages}
library(ggplot2)
library(gridExtra)
library(reshape2)
library(tseries)
library(urca)
library(forecast)
```

## Box-Jenkins method

A. Model form selection

  1. Evaluate stationarity
  2. Selection of the differencing level (d) -- to fix stationarity problems
  3. Selection of the AR level (p)
  4. Selection of the MA level (q)

B. Parameter estimation

C. Model checking

## Stationarity

It is important to test and transform (via differencing) your data to ensure stationarity when fitting an ARMA model using standard algorithms. The standard algorithms for ARIMA models assume stationarity and we will be using those algorithms.  It possible to fit ARMA models without transforming the data.  We will cover that in later chapters. However, that is not commonly done in the literature on forecasting with ARMA models, certainly not in the literature on catch forecasting. 

Keep in mind also that many ARMA models are stationary and you do not want to get in the situation of trying to fit an incompatible process model to your data.  We will see examples of this when we start fitting models to non-stationary data and random walks.

### Look at stationarity in simulated data

We will start by looking at white noise and a stationary AR(1) process from simulated data.  White noise is simply a string of random numbers drawn from a Normal distribution.  `rnorm()` with return random numbers drawn from a Normal distribution.  Use `?rnorm` to understand what the function requires.


```{r white_noise}
TT=100
y = rnorm(TT, mean=0, sd=1) # 100 random numbers
op=par(mfrow=c(1,2))
plot(y, type="l")
acf(y)
par(op)
```

Here we use `ggplot()` to plot 10 white noise time series.

```{r white.noise.ggplot}
dat = data.frame(t=1:TT, y=y)
p1 = ggplot(dat, aes(x=t, y=y)) + geom_line() + 
  ggtitle("1 white noise time series") + xlab("") + ylab("value")
ys = matrix(rnorm(TT*10),TT,10)
ys = data.frame(ys)
ys$id = 1:TT

ys2=melt(ys, id.var="id")
p2 = ggplot(ys2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("10 white noise processes")
grid.arrange(p1, p2, ncol = 1)
```

These are stationary because the variance and mean (level) does not change with time.

An AR(1) process is also stationary.

```{r ar1.plot}
theta=0.8
nsim=10
ar1=arima.sim(TT, model=list(ar=theta))
plot(ar1)
```

We can use ggplot to plot 10 AR(1) time series, but we need to change the data to a data frame.

```{r ar1.ggplot}
dat = data.frame(t=1:TT, y=ar1)
p1 = ggplot(dat, aes(x=t, y=y)) + geom_line() + 
  ggtitle("AR-1") + xlab("") + ylab("value")
ys = matrix(0,TT,nsim)
for(i in 1:nsim) ys[,i]=as.vector(arima.sim(TT, model=list(ar=theta)))
ys = data.frame(ys)
ys$id = 1:TT

ys2=melt(ys, id.var="id")
p2 = ggplot(ys2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("The variance of an AR-1 process is steady")
grid.arrange(p1, p2, ncol = 1)
```

### Stationary around a linear trend

Fluctuating around a linear trend is a very common type of stationarity used in ARMA modeling and forecasting.  This is just a stationary process, like white noise or AR(1), around an linear trend up or down.

```{r wn.w.linear.trend}
intercept = .5
trend=.1
sd=0.5
TT=20
wn=rnorm(TT, sd=sd) #white noise
wni = wn+intercept #white noise witn interept
wnti = wn + trend*(1:TT) + intercept
```

See how the white noise with trend is just the white noise overlaid on a linear trend.

```{r wnt.plot}
op=par(mfrow=c(1,3))
plot(wn, type="l")
plot(trend*1:TT)
plot(wnti, type="l")
par(op)
```

We can make a similar plot with ggplot.

```{r wnt.ggplot}
dat = data.frame(t=1:TT, wn=wn, wni=wni, wnti=wnti)
p1 = ggplot(dat, aes(x=t, y=wn)) + geom_line() + ggtitle("White noise")
p2 = ggplot(dat, aes(x=t, y=wni)) + geom_line() + ggtitle("with non-zero mean")
p3 = ggplot(dat, aes(x=t, y=wnti)) + geom_line() + ggtitle("with linear trend")
grid.arrange(p1, p2, p3, ncol = 3)
```

We can make a similar plot with AR(1) data.  Ignore the warnings about not knowing how to pick the scale.

```{r ar1.trend.plot}
beta1 <- 0.8
ar1 <- arima.sim(TT, model=list(ar=beta1), sd=sd)
ar1i <- ar1 + intercept
ar1ti <- ar1 + trend*(1:TT) + intercept
dat = data.frame(t=1:TT, ar1=ar1, ar1i=ar1i, ar1ti=ar1ti)
p4 = ggplot(dat, aes(x=t, y=ar1)) + geom_line() + ggtitle("AR1")
p5 = ggplot(dat, aes(x=t, y=ar1i)) + geom_line() + ggtitle("with non-zero mean")
p6 = ggplot(dat, aes(x=t, y=ar1ti)) + geom_line() + ggtitle("with linear trend")

grid.arrange(p4, p5, p6, ncol = 3)
```

### Greek landing data

We will look at the anchovy data.  Notice the two `==` in the subset call not one `=`.  We will use the Greek data before 1989 for the lab.

```{r anchovy}
anchovy <- subset(landings, Species=="Anchovy" & Year <= 1989)$log.metric.tons
anchovyts <- ts(anchovy, start=1964)
```

Plot the data.

```{r anchovy.plot}
plot(anchovyts, ylab="log catch")
```

Questions to ask.

* Does it have a trend (goes up or down)? Yes, definitely
* Does it have a non-zero mean?  Yes
* Does it look like it might be stationary around a trend? Maybe

## Augmented Dickey-Fuller test

The null hypothesis for the Dickey-Fuller test is that the data are non-stationary.  We want to REJECT the null hypothesis for this test, so we want a p-value of less that 0.05 (or smaller).

### Test simulated data

Let's start by doing the test on data that we know are stationary, white noise.  We will use `adf.test()` from **tseries**.  Use `?adf.test` to read about this function.  The default values are generally fine.

```{r df.wn}
TT <- 100
wn <- rnorm(TT) # white noise
tseries::adf.test(wn)
```
The null hypothesis is rejected.

Try the test on white noise with a trend and intercept.

```{r df.wn.trend}
intercept <- 1
wnt <- wn + 1:TT + intercept
tseries::adf.test(wnt)
```

The null hypothesis is still rejected.  `adf.test()` uses a model that allows an intercept and trend.

Let's try the test on a random walk (nonstationary).

```{r df.rw}
rw <- cumsum(rnorm(TT))
tseries::adf.test(rw)
```

The null hypothesis is NOT reject as the p-value is greater than 0.05.

### Test the anchovy data

```{r df.anchovy}
tseries::adf.test(anchovyts)
```
The p-value is greater than 0.05.  We cannot reject the null hypothesis.  The null hypothesis is that the data are non-stationary. 

## KPSS test

The null hypothesis for the KPSS test is that the data are stationary.  For this test, we do NOT want to reject the null hypothesis.  In other words, we want the p-value to be greater than 0.05 not less than 0.05.

### Test simulated data

Let's try the KPSS test on white noise with a trend.  The default is a null hypothesis with no trend.  We will change this to `null="Trend"`.

```{r kpss.wnt}
tseries::kpss.test(wnt, null="Trend")
```

The p-value is greater than 0.05.  The null hypothesis of stationarity around a trend is not rejected.

Let's try the KPSS test on white noise with a trend but let's use the default of stationary with no trend.

```{r kpss.wnt.level}
tseries::kpss.test(wnt, null="Level")
```

The p-value is less than 0.05.  The null hypothesis of stationarity around a level is rejected. This is white noise around a trend so it is definitely a stationary process but has a trend.  This illustrates that you need to be thoughtfull when applying stationarity tests.

### Test the anchovy data

Let's try the anchovy data.

```{r kpss.anchovy}
kpss.test(anchovyts, null="Trend")
```

The null is rejected (p-value less than 0.05).  Again stationarity is not supported.

## Optional: Augmented Dickey-Fuller test using `ur.df()`

The `ur.df()` Augmented Dickey-Fuller test in the **urca** package gives us a bit more information on and control over the test.

### Dickey-Fuller test

The Dickey-Fuller test is testing if $\phi=0$ in this model of the data:
$$z_t = \alpha + \beta t + \phi z_{t-1} + e_t$$
which is written as
$$\Delta z_t = z_t-z_{t-1}= \alpha + \beta t + \gamma z_{t-1} + e_t$$
where $z_t$ is your data.  I used $z$ since `ur.df()` denotes the data as `z`.  It is written this way so we can do a linear regression test.

```
lm(z.diff ~ z.lag.1 + 1 + tt)
```
and test if `z.lag.1` is 0.

### Augmented Dickey-Fuller test

The Augmented Dickey-Fuller test allows for higher-order autoregressive processes by including $\Delta z_{t-p}$ in the model.  But our test is still if $\gamma = 0$.
$$\Delta z_t = \alpha + \beta t + \gamma z_{t-1} + \delta_1 \Delta z_{t-1} + \delta_2 \Delta z_{t-2} + \dots$$

### Dickey-Fuller test on white noise

Let's first do the test on data we know is stationary, white noise.  We have to make choose the `type` and `lags`.  If you have no particular reason to not include an intercept and trend, then use `type="trend"`.  This allows both intercept and trend.  Next you need to chose the `lags`.  We will use `lags=0` to do the Dickey-Fuller test.  Note the number of lags you can test will depend on the amount of data that you have.

`lag=0` is fitting this model to the data. You are testing if the effect for `z.lag.1` is 0.  

`z.diff = gamma * z.lag.1 + intercept + trend * tt`

If you see `***` or `**` on the coefficients list for `z.lag.1`, it indicates that the effect of `z.lag.1` is significantly different than 0 and this supports the assumption of stationarity.

The `intercept` and `tt` estimates indicate where there is a non-zero level (intercept) or linear trend (tt).

```{r df.wn2}
wn <- rnorm(TT)
test <- urca::ur.df(wn, type="trend", lags=0)
summary(test)
```

The coefficient part of the summary indicates that `z.lag.1` is different than 0 (so stationary) an no support for intercept or trend.

Notice that the test statistic is LESS than the critical value for `tau3` at 5 percent.  This means the null hypothesis is rejected at $\alpha=0.05$, a standard level for significance testing.

### When you might want to use `ur.df()`

If you remove the trend (and/or level) from your data, the `ur.df()` test allows you to increase the power of the test by removing the trend and/or level from the model.  In addition, `ur.df()` allows you to select the number of lags to use via model selection.  That way you don't use a stationarity test that is too complex to be supported by your data length.

## Dealing with non-stationarity

The anchovy data have failed both tests for the stationarity, the Augmented Dickey-Fuller and the KPSS test.  How to we fix this?  The standard approach is to use differencing.

Let's see how this works with random walk data.  A random walk is non-stationary but the difference is white noise so is stationary:

$$x_t - x_{t-1} = e_t, e_t \sim N(0,\sigma)$$
```{r adf.wn.diff}
adf.test(diff(rw))
kpss.test(diff(rw))
```
If we difference random walk data, the null is rejected for the ADF test and not rejected for the KPSS test.  This is what we want.

Let's try a single difference with the anchovy data.  A single difference means `dat(t)-dat(t-1)`.  We get this using `diff(anchovyts)`.

```{r adf.anchovy.diff}
diff1dat <- diff(anchovyts)
adf.test(diff1dat)
kpss.test(diff1dat)
```

If a first difference were not enough, we would try a second difference which is the difference of a first difference.

```{r second.diff}
diff2dat <- diff(diff1dat)
adf.test(diff2dat)
```
The null hypothesis of a random walk is now rejected so you might think that a 2nd difference is needed for the anchovy data.  However the actual problem is that the default for `adf.test()` includes a trend but we removed the trend with our first difference.  Thus we included an unneeded trend parameter in our test.  Our data are not that long and this affects the result.  

Let's repeat without the trend and we'll see that the null hypothesis is rejected.  The number of lags is set to be what would be used by `adf.test()`.  See `?adf.test`.

```{r}
k <- trunc((length(diff1dat)-1)^(1/3))
test <- urca::ur.df(diff1dat, type="drift", lags=k)
summary(test)
```


### `ndiffs()`

As an alternative to trying many different differences and remembering to include or not include the trend or level, you can use the `ndiffs()` function in the forecast package.  This automates finding the number of differences needed.

```{r ndiff}
forecast::ndiffs(anchovyts, test="kpss")
forecast::ndiffs(anchovyts, test="adf")
```

One difference is required to pass both the ADF and KPSS stationarity tests.

## Summary: stationarity testing

The basic stationarity diagnostics are the following

* Plot your data.  Look for
    - An increasing trend
    - A non-zero level (if no trend)
    - Strange shocks or steps in your data (indicating something dramatic changed like the data collection methodology)
* Apply stationarity tests
    - `adf.test()` p-value should be less than 0.05 (reject null)
    - `kpss.test()` p-value should be greater than 0.05 (do not reject null)
* If stationarity tests are failed, then try differencing to correct
    - Try `ndiffs()` in the forecast package or manually try different differences.
  
  
## Estimating AR and MA parameters

Let's start with fitting to simulated data.

### AR(2) data

Simulate AR(2) data and add a mean level so that the data are not mean 0.

$$x_t = 0.8 x_{t-1} + 0.1 x_{t-2} + e_t$$
$$y_t = x_t + m$$

```{r}
m <- 1
ar2 <- arima.sim(n=100, model=list(ar=c(.8,.1))) + m
```
To see info on `arima.sim()`, type `?arima.sim`.

This is not this model

$$y_t = m + 0.8 y_{t-1} + 0.1 y_{t-2} + e_t$$
It is this model

$$y_t = (m-0.8 m-0.1 m) + 0.8 y_{t-1} + 0.1 y_{t-2} + e_t$$

Fit an ARMA(2) with level to the data.

```{r}
forecast::Arima(ar2, order=c(2,0,0), include.constant=TRUE)
```


Simulate AR(1) data:

$$x_t = 0.8 x_{t-1} + e_t$$

We could also use `arima()`.
```{r}
arima(ar2, order=c(2,0,0), include.mean=TRUE)
```
However we will not be using `arima()` directly because for differenced data, it will not allow us to include and estimated mean level.  Unless we have done transformed our differenced data in a way that ensures it is mean zero, then want to include a mean.

_Try increasing the length of the simulated data (from 100 to 1000 say) and see how that affects your parameter estimates. Run the simulation a few times._

### AR(1) simulated data

```{r}
ar1 = arima.sim(n=100, model=list(ar=c(.8)))+m
forecast::Arima(ar1, order=c(1,0,0), include.constant=TRUE)
```

### ARMA(1,2) simulated data

Simulate ARMA(1,2)
$$x_t = 0.8 x_{t-1} + e_t + 0.8 e_{t-1} + 0.2 e_{t-2}$$

```{r}
arma12 = arima.sim(n=100, model=list(ar=c(0.8), ma=c(0.8, 0.2)))+m
forecast::Arima(arma12, order=c(1,0,2), include.constant=TRUE)
```


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.

### These functions work for data with missing values

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

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

Fit

```{r}
fit <- forecast::Arima(ar2miss, order=c(2,0,0))
fit
```

Note `fitted()` does not return the expected value at time $t$. It is the expected value of $y_t$ given the data up to time $t-1$.

```{r}
plot(ar2miss, type="l")
title("many missing values")
lines(fitted(fit), col="blue")
```
It is easy enough to get the expected value of $y_t$ for all the missing values but we'll learn to do that when we learn the **MARSS** package and can apply the Kalman Smoother in that package.

## Estimating the AR and MA orders

We will use the `auto.arima()` function in **forecast**.  This function will estimate the level of differencing needed to make our data stationary and estimate the AR and MA orders using AICc (or BIC if we choose).

### Example: model selection for AR(2) data

```{r}
forecast::auto.arima(ar2)
```

Works with missing data too though might not estimate very close to the true model form.

```{r}
forecast::auto.arima(ar2miss)
```

## Fitting to 100 simulated data sets

Let's fit to 100 simulated data sets and see how often the true (generating) model form is selected.
```{r}
save.fits = rep(NA,100)
for(i in 1:100){
  a2 = arima.sim(n=100, model=list(ar=c(.8,.1)))
  fit = auto.arima(a2, seasonal=FALSE, max.d=0, max.q=0)
  save.fits[i] = paste0(fit$arma[1], "-", fit$arma[2])
}
table(save.fits)
```

`auto.arima()` uses AICc for selection by default.  You can change that to AIC or BIC using `ic="aic"` or `ic="bic"`.

_Repeat the simulation using AIC and BIC to see how the choice of the information criteria affects the model that is selected._

### Trace=TRUE

We can set `Trace=TRUE` to see what models `auto.arima()` fit.

```{r}
forecast::auto.arima(ar2, trace=TRUE)
```

### stepwise=FALSE

We can set `stepwise=FALSE` to use and exhaustive search.  The model may be different than the result from the non-exhaustive search.

```{r}
forecast::auto.arima(ar2, trace=TRUE, stepwise=FALSE)
```

## Fit to the anchovy data

```{r}
fit <- auto.arima(anchovyts)
fit
```

Note `arima()` writes a MA model like:

$$x_t = e_t + b_1 e_{t-1} + b_2 e_{t-2}$$

while many authors use this 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 (1996) who analyze these same data.  *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.*


## Check residuals

We can do a test of autocorrelation of the residuals with `Box.test()` with `fitdf` adjusted for the number of parameters estimated in the fit.  In our case, MA(1) and drift parameters.

```{r}
res <- resid(fit)
Box.test(res, type="Ljung-Box", lag=12, fitdf=2)
```

`checkresiduals()` in the **forecast** package will automate this test and show some standard diagnostics plots.

```{r}
forecast::checkresiduals(fit)
```

## Forecast

We can create a forecast from our anchovy ARIMA model using `forecast()`.  The shading is the 80% and 95% prediction intervals.

```{r}
fr <- forecast::forecast(fit, h=10)
plot(fr)
```

## Seasonal ARIMA model

The chinook data are monthly and start in January 1990.  To make this into a ts object do

```{r}
load("chinook.RData")
chinookts <- ts(chinook$log.metric.tons, start=c(1990,1), 
                frequency=12)
```
`start` is the year and month and frequency is the number of months in the year.  

Use `?ts` to see more examples of how to set up ts objects.

### Plot seasonal data

```{r}
plot(chinookts)
```

### `auto.arima()` for seasonal ts

`auto.arima()` will recognize that our data has season and fit a seasonal ARIMA model to our data by default.  Let's define the training data up to 1998 and use 1999 as the test data.

```{r}
traindat <- window(chinookts, c(1990,10), c(1998,12))
testdat <- window(chinookts, c(1999,1), c(1999,12))
fit <- forecast::auto.arima(traindat)
fit
```

Use `?window` to understand how subsetting a ts object works.

## Forecast using seasonal model

Forecasting works the same using the `forecast()` function.

```{r}
fr <- forecast::forecast(fit, h=12)
plot(fr)
points(testdat)
```


## Problems

1. Repeat the stationarity tests for another species in the landings data set.  Sardine, Chub.mackerel, and Horse.mackerel have enough data. Here is how to set up the data for another species.
    b. Do a Dickey-Fuller test using `ur.df()` and `adf.test()`. What do you need to set the lags to to achieve that?
    a. Do an Augmented Dickey-Fuller test using `ur.df()`. How did you choose to set the lags?
    b. Do a KPSS test using `kpss.test()`.

```{r get.another.species}
datdf <- subset(landings, Species=="Chub.mackerel" & Year<=1989)
dat <- ts(datdf$log.metric.tons, start=1964)
```

2. Repeat the stationarity tests and differencing tests for anchovy using the full data.  The full data go to 2007.  

```{r read_data_prob2}
load("landings.RData")
datdf <- subset(landings, Species=="Anchovy")
dat <- ts(datdf$log.metric.tons, start=1964)
```
    a. Do the conclusions regarding stationarity and the amount of differencing needed change? 
    b. Repeat with only the data after 1987.  Are the conclusions different for only the recent data?
    c. Fit the using `auto.arima()` to the full data (1964-2007) and recent data (1987 to present).  Do the selected models change?

4. Fit each of the models listed under the `auto.arima()` trace using `Arima()` and show that you can produce the same AICc value that is shown in the trace table.

```{r results='hide'}
forecast::auto.arima(anchovy, trace=TRUE)
```

5. Use `auto.arima()` with `stepwise=FALSE` to find the top 3 models (according to AICc) for the anchovy.
    a. Use `Arima()` to fit the top 3 models to the data.
    a. Create a 5-year forecast for each of the top 3 models using the fits from (a).  
    b. Do the forecasts differ?
    
6. 
    a. Fit a seasonal model to the Chinook data up to 1999 using `auto.arima()`.  
    b. Then create a forecast through 2015.
    c. Plot the forecast with the 2014 and 2015 actual landings added as data points.
