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

## Forecasting with an ARIMA model

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 a ARIMA(2,1,0) model with drift:
$$z_t = \mu + \beta_1 z_{t-1} + \beta_2 z_{t-2} + e_t$$
where $z_t = x_t - x_{t-1}$, the first difference.

`Arima()` would write this model:
$$(z_t-m) = \beta_1 (z_{t-1}-m) + \beta_2 (z_{t-2}-m) + e_t$$
The relationship between $\mu$ and $m$ is $\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:

$$z_t = 0.1065807-0.53850433 z_{t-1} - 0.44732522 z_{t-2} + e_t$$

To get our forecast for 1990, we do this

$$(x_{90}-x_{89}) = 0.106 - 0.538 (x_{89}-x_{88}) - 0.447 (x_{88}-x_{87})$$

Thus

$$x_{90} = x_{89}+0.106 -0.538  (x_{89}-x_{88}) - 0.447 (x_{88}-x_{87})$$

##

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 `forecast::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::Arima(anchovy.miss, order=c(2,1,0), include.constant=TRUE)
fr <- forecast::forecast(fit, h=5)
fr
```

##

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

## Using auto.arima()

We can let forecast to select the ARIMA model:

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

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