```{r setup, include=FALSE}
knitr::opts_chunk$set(R.options=list(width=60))
set.seed(123)
```


```{r load_data, message=FALSE, warning=FALSE, echo=FALSE}
load("landings.RData")
landings$log.metric.tons = log(landings$metric.tons)
landings = subset(landings, Year <= 1989)
anchovy = subset(landings, Species=="Anchovy")$log.metric.tons
anchovyts = ts(anchovy, start=c(1964,1), frequency = 1)
library(ggplot2)
library(gridExtra)
library(reshape2)
library(tseries)
library(urca)
```


## 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

## Good news

Much of the Box-Jenkins method will be automated with the **forecast** package functions, which we will use in the lab.

## Stationarity

Stationarity means 'not changing in time' in the context of time-series models.  Typically we test the trend and variance, however more generally all statistical properties of a time-series is time-constant if the time series is 'stationary'.

## Example

Many ARMA models exhibit stationarity.  White noise is one type:
$$x_t = e_t, e_t \sim N(0,\sigma)$$

```{r fig.stationarity, fig.height = 3.5, fig.width = 7, fig.align = "center", echo=FALSE}
require(gridExtra)
require(reshape2)

TT=100
y = rnorm(TT)
dat = data.frame(t=1:TT, y=y)
p1 = ggplot(dat, aes(x=t, y=y)) + geom_line() + 
  ggtitle("White Noise") + 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("The variance of a white noise process is steady")
grid.arrange(p1, p2, ncol = 1)
```

## Example

An AR-1 process with $-1<b<1$
$$x_t = \phi x_{t-1} + e_t$$
is also stationary.

```{r fig.stationarity2, fig.height = 3.5, fig.width = 7, fig.align = "center", echo=FALSE}
require(ggplot2)
require(reshape2)
theta=0.8
nsim=10
ar1=as.vector(arima.sim(TT, model=list(ar=theta)))
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)
```

## Stationarity around non-zero mean

We can also have stationarity around a non-zero level or around a linear trend.

```{r fig.stationarity3, fig.height = 4.5, fig.width = 7, fig.align = "center", echo=FALSE}
require(ggplot2)
require(gridExtra)
intercept = .5
trend=.1
dat = data.frame(t=1:TT, wn=rnorm(TT))
dat$wni = dat$wn+intercept
dat$wnti = dat$wn + trend*(1:TT) + intercept
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")

ar1 = rep(0,TT)
err = rnorm(TT)
for(i in 2:TT) ar1[i]=theta*ar1[i-1]+err[i]
dat = data.frame(t=1:TT, ar1=ar1)
dat$ar1i = dat$ar1+intercept
dat$ar1ti = dat$ar1 + trend*(1:TT) + intercept
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(p1, p4, p2, p5, p3, p6, ncol = 2)
```

## Mathematically it looks like this

AR-1

1. Non-zero mean adds $\mu$: $x_t = \mu + \phi x_{t-1} + e_t$
2. Linear trend adds $at$: $x_t = \mu + at + \phi x_{t-1} + e_t$

White noise ($b=0$)

1. Non-zero mean: $x_t = \mu + e_t$
2. Linear trend: $x_t = \mu + at + e_t$


## Non-stationarity

One of the most common forms of non-stationarity that is tested for is that the process is a random walk $x_t = x_{t-1} + e_t$.  A test for an underlying random walk is called a 'unit root' test.


```{r fig.nonstationarity, fig.height = 4, fig.width = 8, fig.align = "center", echo=FALSE}
require(ggplot2)
require(reshape2)

rw = rep(0,TT)
for(i in 2:TT) rw[i]=rw[i-1]+err[i]
dat = data.frame(t=1:TT, rw=rw)
p1 = ggplot(dat, aes(x=t, y=rw)) + geom_line() + 
  ggtitle("Random Walk") + xlab("") + ylab("value")
rws = apply(matrix(rnorm(TT*nsim),TT,nsim),2,cumsum)
rws = data.frame(rws)
rws$id = 1:TT

rws2=melt(rws, id.var="id")
p2 = ggplot(rws2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("The variance of a random walk process grows in time")
grid.arrange(p1, p2, ncol = 1)
```

## Random walk with $\mu$ and $at$ added

Similar to the way we added an intecept and linear trend to the stationarity process equations, we can do the same to the random walk equation.

1. Non-zero mean or intercept: $x_t = \mu + x_{t-1} + e_t$

2. Linear trend: $x_t = \mu + at + x_{t-1} + e_t$

## Random walk with $\mu$ and $at$ added

The effects are fundamentally different however.  The addition of $\mu$ leads to a upward mean linear trend while the addition of $at$ leads to exponential growth (or decline).

```{r fig.stationarity4, fig.height = 4, fig.width = 7, fig.align = "center", echo=FALSE}
require(ggplot2)
require(gridExtra)
dat = data.frame(t=1:TT, y=cumsum(rnorm(TT)))
dat$yi = cumsum(rnorm(TT,intercept,1))
dat$yti = cumsum(rnorm(TT,intercept+trend*1:TT,1))
p1 = ggplot(dat, aes(x=t, y=y)) + geom_line() + ggtitle("Random Walk")
p2 = ggplot(dat, aes(x=t, y=yi)) + geom_line() + ggtitle("with non-zero mean added")
p3 = ggplot(dat, aes(x=t, y=yti)) + geom_line() + ggtitle("with linear trend added")

grid.arrange(p1, p2, p3, ncol = 1)
```

## Testing for stationarity

Why is evaluating stationarity important? 

- Many AR models have a flat level or trend and time-constant variance.  If your data do not have those properties, you are fitting a model that is fundamentally inconsistent with your data.
- Many standard algorithms for fitting ARIMA models assume stationarity.  Note, you can fit ARIMA models without making this assumption, but you need to use the appropriate algorithm.

## Testing for stationarity

We will discuss three common approaches to evaluating stationarity:

- Visual test
- (Augmented) Dickey-Fuller test
- KPSS test


## Visual test

The visual test is simply looking at a plot of the data versus time.  Look for

- Change in the level over time.  Is the time series increasing or decreasing?  Does it appear to cycle?
- Change in the variance over time.  Do deviations away from the mean change over time, increase or decrease?


## Example anchovy and sardine catch in Greek waters

```{r fig.vis, fig.height = 4, fig.width = 8, fig.align = "center", echo=FALSE}
require(ggplot2)
dat = subset(landings, Species %in% c("Anchovy", "Sardine") & 
               Year <= 1989)
dat$log.metric.tons = log(dat$metric.tons)
ggplot(dat, aes(x=Year, y=log.metric.tons)) +
  geom_line()+facet_wrap(~Species)
```

## Dickey-Fuller test

The Dickey=Fuller test (and Augmented Dickey-Fuller test) look for evidence that the time series has a unit root.  

The **null hypothesis** is that the time series has a unit root, that is, it has a random walk component.  

The **alternative hypothesis** is some variation of stationarity.  The test has three main verisons. 

## Dickey-Fuller nulls and alternatives

It is hard to see but in the panels on the left, the variance around the trend is increasing and on the right, it is not.

```{r fig.df, fig.height = 4.5, fig.width = 8, fig.align = "center", echo=FALSE}
require(ggplot2)
require(gridExtra)
#####
ys = matrix(0,TT,nsim)
for(i in 2:TT) ys[i,]=ys[i-1,]+rnorm(nsim)
rws = data.frame(ys)
rws$id = 1:TT
library(reshape2)
rws2=melt(rws, id.var="id")
p1 = ggplot(rws2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("Null Non-stationary", subtitle="Random walk")
ys = matrix(0,TT,nsim)
for(i in 2:TT) ys[i,]=theta*ys[i-1,]+rnorm(nsim)
ar1s = data.frame(ys)
ar1s$id = 1:TT
library(reshape2)
ar1s2=melt(ar1s, id.var="id")
p2 = ggplot(ar1s2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("Alternate Stationary", subtitle="AR1")
#####
ys = matrix(intercept,TT,nsim)
for(i in 2:TT) ys[i,]=intercept+ys[i-1,]+rnorm(nsim)
rws = data.frame(ys)
rws$id = 1:TT
library(reshape2)
rws2=melt(rws, id.var="id")
p3 = ggplot(rws2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("", subtitle="Random walk + constant (drift)")
ys = matrix(intercept/(1-theta),TT,nsim)
for(i in 2:TT) ys[i,]=intercept+theta*ys[i-1,]+rnorm(nsim)
ar1s = data.frame(ys)
ar1s$id = 1:TT
library(reshape2)
ar1s2=melt(ar1s, id.var="id")
p4 = ggplot(ar1s2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("Alternate", subtitle="AR1 + constant (non-zero level)")
#####
ys = matrix(intercept+trend*1,TT,nsim)
for(i in 2:TT) ys[i,]=intercept+trend*i+ys[i-1,]+rnorm(nsim)
rws = data.frame(ys)
rws$id = 1:TT
library(reshape2)
rws2=melt(rws, id.var="id")
p5 = ggplot(rws2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("", subtitle="Random walk + linear trend (a*t)")
ys = matrix((intercept+trend*1)/(1-theta),TT,nsim)
for(i in 2:TT) ys[i,]=intercept+trend*i+theta*ys[i-1,]+rnorm(nsim)
ar1s = data.frame(ys)
ar1s$id = 1:TT
library(reshape2)
ar1s2=melt(ar1s, id.var="id")
p6 = ggplot(ar1s2, aes(x=id,y=value,group=variable)) +
  geom_line() + xlab("") + ylab("value") +
  ggtitle("Alternate", subtitle="AR1 + linear trend")
#####
grid.arrange(p1, p2, p3, p4, p5, p6, ncol = 2)
```

<!--
## Dickey-Fuller null and alternatives

In math, here are the null and alternative hypotheses as equations.  In each, we are testing if $(\phi-1)=0$. 

None: Null $x_t = x_{t-1} + e_t$; Alternative $x_t = \phi x_{t-1} +e_t$

Drift (constant): Null $x_t = \mu + x_{t-1} + e_t$; Alternative $x_t = \mu + \phi x_{t-1} + e_t$

Trend: Null $x_t = \mu + at + x_{t-1} + e_t$; Alternative $x_t = \mu + at + \phi x_{t-1} + e_t$

Note typically in the alternative hypothesis we restrict $-1<\phi<1$ which is stationary.  If $|\phi|>1$, the alternative process is 'explosive'.

-->

## Dickey-Fuller test using `tseries::adf.test()`

`adf.test()` in the tseries package will apply the Augmented Dickey-Fuller **with a constant and trend** and report the p-value.  We want to reject the Dickey=Fuller null hypothesis of non-stationarity.  We will set `k=0` to apply the Dickey-Fuller test which tests for AR(1) stationarity.  The Augmented Dickey-Fuller tests for more general lag-p stationarity.

```
adf.test(x, alternative = c("stationary", "explosive"),
         k = trunc((length(x)-1)^(1/3)))
```

## Example: Dickey-Fuller tests on the anchovy time series

Here is how to apply this test to the anchovy data. The null hypothesis is not rejected.  That is not what we want.

```{r adf.test.anchovy1}
adf.test(anchovyts, k=0)
```


## Dickey-Fuller test with `urca::ur.df`

The `urca` R package can also be used to apply the Dickey-Fuller tests.  Use `lags=0` for Dickey-Fuller which tests for AR(1) stationarity.  We will set `type="trend"` to deal with the trend seen in the anchovy data.  Note, `adf.test()` uses this type by default.

```
ur.df(y, type = c("none", "drift", "trend"), lags = 0)
```

## Dickey-Fuller test with 'ur.df'

```{r dickey.fuller, message=FALSE, warning=FALSE}
test = urca::ur.df(anchovyts, type="trend", lags=0)
test
```

## Dickey-Fuller test with 'ur.df'

The test statistic and the critical value at $\alpha = 0.05$ are
```{r teststat}
attr(test, "teststat")
```
```{r cval}
attr(test,"cval")
```
The statistic is larger than the critical value and thus the null hypothesis of non-stationarity is not rejected. That's not what we want.

##

The $\tau_3$ is the one we want. This is the stationarity parameter.

$x_t = \phi x_{t-1} + \mu + a t + e_t$

$x_t - x_{t-1} = \tau_3 x_{t-1} + \phi_2 + \phi_3 t + e_t$

<!--

## Augmented Dickey-Fuller test

The Dickey-Fuller test assumes that the stationary process is AR(1) (autoregressive lag-1).  The Augmented Dickey-Fuller test allows a general stationary process.  The idea of the test however is the same.

We can apply the Augmented Dickey-Fuller test with the `ur.df()` function or the `adf.test()` function in the `tseries` package.

```
adf.test(x, alternative = c("stationary", "explosive"),
         k = trunc((length(x)-1)^(1/3)))
```

The alternative is either stationary like $x_t = \delta x_{t-1} + \eta_t$ with $\delta<1$ or 'explosive' with $\delta>1$.  `k` is the number of lags which determines the number of time lags allowed in the autoregression.  `k` is generally determined by the length of your time series.

## Example: Augmented Dickey-Fuller tests with adf.test()

With the `tseries` package, we apply the Augmented Dickey-Fuller test with `adf.test()`.  This function uses the test where the alternative model is stationary around a linear trend: $x_t = \mu + at + \delta x_{t-1} + e_t$.

```{r dickey.fuller2, message=FALSE, warning=FALSE}
require(tseries)
adf.test(anchovy)
```

In both cases, we do not reject the null hypothesis that the data have a random walk.  Thus there is not support for these time series being stationary.


## Example: Augmented Dickey-Fuller tests with ur.df()

With the `urca` package, we apply the Augmented Dickey-Fuller test with `ur.df()`.  The defaults for `ur.df()` are different than for `adf.test()`.  

`ur.df()` allows you to specify which of the 3 alternative hypotheses you want: none (stationary around 0), drift (stationary around a non-zero intercept), trend (stationary around a linear trend).

Another difference is that by default, `ur.df()` uses a fixed lag of 1 while by default `adf.test()` selects the lag based on the length of the time series.

## Example: Augmented Dickey-Fuller tests with ur.df()

We will specify "trend" to make the test similar to `adf.test()`.  We will set the lags like `adf.test()` does also.

```{r dickey.fuller.ur.df, message=FALSE, warning=FALSE}
require(urca)
k = trunc((length(anchovy)-1)^(1/3))
test = ur.df(anchovy, type="trend", lags=k)
test
```

The test statistic values are the same, but we need to look up the critical values with `summary(test)`.

-->

## KPSS test

The Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test has as the null hypothesis that the time series is stationary around a level trend (or a linear trend). The alternative hypothesis for the KPSS test is a random walk.

The stationarity assumption is general; it does not assume a specific type of stationarity such as white noise.

If both KPSS and Dickey-Fuller tests support non-stationarity, then the stationarity assumption is not supported.

## Example: KPSS tests

```{r kpss.test, message=FALSE, warning=FALSE}
tseries::kpss.test(anchovyts, null="Trend")
```

Here `null="Trend"` was included to account for the increasing trend in the data.  The null hypothesis of stationarity is rejected.  Thus both the KPSS and Dickey-Fuller tests support the hypothesis that the anchovy time series is non-stationary.  That's not what we want.


## Differencing the data to make the mean stationary

Differencing means to create a new time series  $z_t = x_t - x_{t-1}$.  First order differencing means you do this once (so $z_t$) and second order differencing means you do this twice (so $z_t - z_{t-1}$).

The `diff()` function takes the first difference:

```{r diff1}
x <- diff(c(1,2,4,7,11))
x
```

The second difference is the first difference of the first difference.

```{r diff2}
diff(x)
```

## Anchovy catch first differenced

Here is a plot of the anchovy data and its first difference.

```{r diff1.plot, fig.align="center", fig.height = 4, fig.width = 8, echo=FALSE}
par(mfrow=c(1,2))
plot(anchovyts,xlab="")
title("Anchovy")
plot(diff(anchovyts),xlab="")
title("Anchovy 1st diff")
```


## Stationarity of the first differences

Let's test the anchovy data with one difference using the KPSS test.

```{r diff.anchovy.test, message=FALSE, warning=FALSE}
diff.anchovy = diff(anchovyts)
kpss.test(diff.anchovy)
```

The null hypothesis of stationairity is not rejected. That is good.

## Stationarity of the first differences

Let's test the first difference of the anchovy data using the Augmented Dickey-Fuller test.  We do the default test and allow it to chose the number of lags.

```{r test.dickey.fuller.diff1}
adf.test(diff.anchovy)
```

## 

The null hypothesis of non-stationarity is not rejected.  That is not what we want.  However, we differenced which removed the trend thus we are testing against a more general model than we need.  Let's test with an alternative hypothesis that has a non-zero mean and no trend.  We can do this with `ur.df()` and `type='drift'`.

```{r test.dickey.fuller.diff2}
test <- ur.df(diff.anchovy, type="drift")
```

## 

The test statistic and the critical values are
```{r teststat.diff}
attr(test, "teststat")
```

```{r cval.diff}
attr(test,"cval")
```

The null hypothesis of NON-stationairity IS rejected. That is good.


## `forecast::ndiffs()` function

As an alternative to trying many different differences, 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")
```

## Summary

Test stationarity before you fit a ARMA model.

Visual test: Do the data fluctuate around a level or do they have a trend or look like a random walk?

_Yes or maybe?_  -> Apply a "unit root" test. ADF or KPSS

_No or fails the unit root test?_ -> Apply differencing and re-test.

_Still not passing?_ -> Try a second difference or you may need to transform the data (if say it has an exponential trend).

_Still not passing?_ -> ARMA model might not be the best choice. Or you may need to use an adhoc detrend.

**These steps are automated by the forecast package**

## 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

## ACF and PACF

On Tuesday, you learned how to use ACF and PACF to visually infer the AR and MA lags for a ARMA model.

```{r pacf-acf-anchovy, echo=FALSE, fig.cap="ACF and PACF of anchovy time series."}
par(mfrow=c(1,2))
acf(diff(anchovy))
pacf(diff(anchovy))
```

## Formal model selection

This weighs how well the model fits against how many parameters your model has. Basic idea is to fit (many) models and use AIC, AICc or BIC to select.

The `auto.arima()` function in the forecast package in R allows you to easily do this and will also select the level of differencing (via ADF or KPSS tests).

```{r auto.arima, eval=FALSE}
forecast::auto.arima(anchovy)
```

Type `?forecast::auto.arima` to see a full description of the function.

## Model selection with `auto.arima()`

```{r auto.arima2}
forecast::auto.arima(anchovy)
```

The output indicates that the 'best' model is a MA(1) with first difference.  "with drift" means that the mean of the anchovy first differences (the data for the model) is not zero. 

## Trace = TRUE

By default, step-wise selection is used for the model search. You can see what models that `auto.arima()` tried using `trace=TRUE`.  The models are selected on AICc by default and the AICc value is shown next to the model.

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

## Selected model

First difference of the data is MA(1) with drift

$$x_t - x_{t-1} = \mu + w_t + \theta_1 w_{t-1}$$

where $w_t$ is white noise.

## Example: Fit to simulated AR(2) data

```{r fitting.example.1}
set.seed(100)
a1 = arima.sim(n=100, model=list(ar=c(.8,.1)))
forecast::auto.arima(a1, seasonal=FALSE, max.d=0)
```

The 'best-fit' model is AR(1) not AR(2). 

## How often is the 'true' model is chosen

Let's run 100 simulations of a AR(2) process and record the best fits.

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

##

Overwhelmingly the correct type of model (AR) is selected, but usually a simpler model of AR(2) is chosen over AR(2).

Table heading is AR order - MA order.

```{r fit.1000.table}
table(save.fits)
```

## stepwise=FALSE

By default, step-wise selection is used and an approximation is used for the models tried in the model selection step.  For a final model selection, you should turn these off.

```{r}
forecast::auto.arima(anchovy, stepwise=FALSE, 
                     approximation=FALSE)
```

## Summary: model selection and fitting

- Once you have dealt with stationarity, you need to determine the order of the model: the AR part and the MA part.

- Although you could simply use `auto.arima()`, it is best to run `acf()` and `pacf()` on your data to understand it better.  Definitely you want to plot your data and visually look for stationarity issues.

- Also evaluate if there are reasons to assume a particular structure.  

  - Are you using an established model form, from say another paper?
  - Are you fitting to a process that is fundamentally AR only or AR + MA?

## 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**

## Check the residuals

Residuals = difference between the expected (fitted) value of $x_t$ and the data

There is no observation error in an ARMA model. The expected value is the $x_t$ expected from data up to $t-1$.

$x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + w_t$

$\hat{x}_t = \phi_1 x_{t-1} + \phi_2 x_{t-2}$


## `residuals()` function in R

The `residuals()` function will return the residuals for fitted models.

```{r}
fit <- forecast::auto.arima(anchovy)
residuals(fit)
```

## `fitted()` function in R

The `fitted()` function will return the expected values. Remember that for a ARMA model, these are the expected values conditioned on the data up to time $t-1$.

```{r}
fitted(fit)
```

The residuals are `data` minus `fitted`.

## Standard residuals tests

* Plot the residuals. They should look roughly like white noise.
* Look at the ACF of the residuals.  They should be uncorrelated.
* Look at the histogram.  They should be normally distributed (if that is your error assumption).

## Residuals check with forecast package

```{r results="hide"}
forecast::checkresiduals(fit)
```

## Test for autocorrelation

The standard test for autocorrelated time-series residuals is the Ljung-Box test.  The null hypothesis for this test is **no autocorrelation**.  We do not want to reject the null. 

```{r}
forecast::checkresiduals(fit, plot=FALSE)
```

$p>0.05$ would be interpreted as not enough statistical evidence to reject the null hypothesis.

## Forecasting (brief)

The basic idea of forecasting with an ARIMA model to estimate the parameters and forecast forward.

For example, let's say we want to forecast with this ARIMA(2,1,0) model:
$$y_t = \mu + \beta_1 y_{t-1} + \beta_2 y_{t-2} + e_t$$
or as `Arima()` would write this model:
$$(y_t-m) = \beta_1 (y_{t-1}-m) + \beta_2 (y_{t-2}-m) + e_t$$
where $y_t$ is the first difference of our anchovy data. $\mu = m(1 - \beta_1 - \beta_2)$.

##

Let's estimate the $\beta$'s for this model from the anchovy data.

```{r}
fit <- forecast::Arima(anchovyts, order=c(2,1,0), include.constant=TRUE)
coef(fit)
```

```{r}
mu <- coef(fit)[3]*(1-coef(fit)[1]-coef(fit)[2])
mu
```

##

So we will forecast with this model:

$$y_t = 0.1065807-0.53850433 y_{t-1} - 0.44732522 y_{t-2} + e_t$$

To get our forecast for 1990, we do this

$$(x_{1990}-x_{1989}) = 0.106 - 0.538 (x_{1989}-x_{1988}) - 0.447 (x_{1988}-x_{1987})$$

Thus

$$x_{1990} = x_{1989}+0.106 -0.538  (x_{1989}-x_{1988}) - 0.447 (x_{1988}-x_{1987})$$

##

Here is R code to do that:

```{r}
anchovyts[26]+mu+coef(fit)[1]*(anchovyts[26]-anchovyts[25])+
  coef(fit)[2]*(anchovyts[25]-anchovyts[24])
```

## Forecasting with `forecast()`

`forecast(fit, h=h)` automates the forecast calculations for us and computes the upper and lower prediction intervals.  Prediction intervals include uncertainty in parameter estimates plus the process error uncertainty.

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

## Plotting our forecasts

```{r}
plot(fr, xlab="Year")
```

## Missing values

Missing values are allowed for `Arima()` and `arima()`. We can produce forecasts with the same code.

```{r fig.show=FALSE}
anchovy.miss <- anchovyts
anchovy.miss[10:11] <- NA
anchovy.miss[20:21] <- NA
fit <- forecast::auto.arima(anchovy.miss)
fr <- forecast::forecast(fit, h=5)
fr
```

##

```{r}
plot(fr)
```


## Seasonality

### Load the chinook salmon data set

```{r}
load("chinook.RData")
head(chinook)
```

##

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

```{r}
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)
```

## Seasonal ARIMA model

Seasonally differenced data:

$$z_t = x_t - x_{t+s} - m$$

Basic structure of a seasonal AR model

$z_t$ = AR(p) + AR(season) + AR(p+season)


Example AR(1) non-seasonal part + AR(1) seasonal part

$$z_t = \phi_1 z_{t-1} + \Phi_1 z_{t-12} - \phi_1\Phi_1 z_{t-13}$$

## Notation

ARIMA (p,d,q)(ps,ds,qs)S


ARIMA (1,0,0)(1,1,0)[12]

## `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.  We will 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
```

## Forecast using seasonal model

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

## Missing values

Missing values are ok when fitting a seasonal ARIMA model

```{r echo=FALSE}
plot(fr)
```

## Summary

Basic steps for identifying a seasonal model.  **forecast** automates most of this.

* Check that you have specified your season correctly in your ts object.

* Plot your data.  Look for trend, seasonality and random walks.

## Summary

* Use differencing to remove season and trend.
    * Season and no trend.  Take a difference of lag = season
    * No seasonality but a trend.  Try a first difference
    * Both. Do both types of differences
    * Neither. No differencing.
    * Random walk. First difference.
    * Parametric looking curve. Tranform.

## Summary

* Examine the ACF and PACF of the differenced data.
    * Look for patterns (spikes) at seasonal lags

* Estimate likely models and compare with model selection criteria (or cross-validation). Use `TRACE=TRUE`

* Do residual checks. 

### Intercepts and drift in `Arima()`

If $d=0$, 

* $d=0$: Arima(x, order=c(1,0,0), include.drift=FALSE, include.mean=TRUE)

$m$ is estimated and called `intercept`.

$$(x_t-m) = \phi_1 (x_{t-1}-m) + w_t$$

* $d=0$: Arima(x, order=c(1,0,0), include.drift=TRUE, include.mean=FALSE)

$\mu$ is estimated and called `drift`.

$$x_t = \mu + \phi_1 x_{t-1} + w_t$$
* $d=0$: Arima(x, order=c(1,0,0), include.drift=TRUE, include.mean=TRUE)

$\mu$ and $m$ are estimated and called `drift` and `intercept`.

$$(x_t-m) = \mu + \phi_1 (x_{t-1}-m) + w_t$$

If $d=1$, then `include.mean` is ignored in `Arima()` and `include.drift` estimates an intercept like `include.mean` in the $d=0$ case, but it is called `drift` in the output. $y_t = x_t-x_{t-1}$.

* $d=1$: Arima(x, order=c(1,1,0), include.drift=TRUE)

$m$ is estimated and called `drift`. 

$$(y_t-m) = \phi_1 (y_{t-1}-m) + w_t$$
* $d=1$: Arima(x, order=c(1,1,0), include.drift=FALSE)

$$y_t = \phi_1 y_{t-1} + w_t$$

* $d=1$: Arima(x, order=c(0,1,0), include.drift=TRUE)

This is a random walk with drift.

$$(y_t-m) = w_t$$
which is

$$x_t = m + x_{t-1} + w_t$$
If $d\ge2$, then both `include.mean` and `include.drift` are ignored. $z_t = y_t-y_{t-1}$.

* $d=2$: Arima(x, order=c(1,2,0))

$$z_t = \phi_1 z_{t-1} + w_t$$

### Intercepts in `arima()`

If $d=0$, 

* $d=0$: arima(x, order=c(1,0,0), include.mean=TRUE)

$m$ is estimated and called `intercept`.

$$(x_t-m) = \phi_1 (x_{t-1}-m) + w_t$$


If $d=1$, then `include.mean` is ignored and no intercept can be estimated. 

* $d=1$: arima(x, order=c(1,1,0), include.mean=TRUE)

$$y_t = \phi_1 y_{t-1}+ w_t$$
* $d=1$: arima(x, order=c(0,1,0))

Because an intercept cannot be estimated, this means that a random walk with drift cannot be estimated by `arima()`.

$$y_t = w_t$$
only can be estimated which is random walk without drift.

$$x_t = x_{t-1} + w_t$$





