```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(kableExtra)
set.seed(123)
```

## Topics for today

Review

* White noise  

* Random walks  

Autoregressive (AR) models

Moving average (MA) models

Autoregressive moving average (ARMA) models

Using ACF & PACF for model ID


## White noise (WN)

A time series $\{w_t\}$ is discrete white noise if its values are

1. independent  

2. identically distributed with a mean of zero

The distributional form for $\{w_t\}$ is flexible


## White noise (WN)

```{r, fig.cap="$w_t = 2e_t - 1; e_t \\sim \\text{Bernoulli}(0.5)$"}
par(mfrow = c(1,2), mai = c(1.5,0.9,0.1,0.1), omi = c(0,0,0,0))
tt <- rbinom(100, 1, 0.5) * 2 - 1
plot.ts(tt, ylab = expression(italic(w[t])))
acf(tt)
```


## Gaussian white noise

We often assume so-called _Gaussian white noise_, whereby

$$
w_t \sim \text{N}(0,\sigma^2)
$$

and the following apply as well

&nbsp; &nbsp; autocovariance:&nbsp; $\gamma_k =
    \begin{cases}
      \sigma^2 & \text{if } k = 0 \\
      0 & \text{if } k \geq 1
    \end{cases}$

&nbsp; &nbsp; autocorrelation: &nbsp; $\rho_k =
    \begin{cases}
      1 & \text{if } k = 0 \\
      0 & \text{if } k \geq 1
    \end{cases}$


## Gaussian white noise

```{r ex_gaussian_wn, fig.cap="$w_t \\sim \\text{N}(0,1)$"}
par(mfrow = c(1,2), mai = c(1.5,0.9,0.1,0.1), omi = c(0,0,0,0))
tt <- rnorm(100)
plot.ts(tt, ylab = expression(italic(w[t])))
acf(tt)
```


## Random walk (RW)

A time series $\{x_t\}$ is a random walk if

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

2. $w_t$ is white noise


## Random walk (RW)

```{r ex_rw, fig.cap="$x_t = x_{t-1} + w_t; w_t \\sim \\text{N}(0,1)$"}
par(mfrow = c(1,2), mai = c(1.5,0.9,0.1,0.1), omi = c(0,0,0,0))
tt <- cumsum(rnorm(100))
plot.ts(tt, ylab = expression(italic(x[t])))
acf(tt)
```


## Biased random walk

A _biased random walk_ (or _random walk with drift_) is written as

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

where $u$ is the bias (drift) per time step and $w_t$ is white noise


## Biased random walk

```{r ex_biased_rw, fig.cap="$x_t = x_{t-1} + 1 + w_t; w_t \\sim \\text{N}(0,1)$"}
par(mfrow = c(1,2), mai = c(1.5,0.9,0.1,0.1), omi = c(0,0,0,0))
xx <- ww <- rnorm(100)
uu <- 1
for(t in 2:100) {
  xx[t] <- xx[t-1] + uu + ww[t]
}
plot.ts(xx, ylab = expression(italic(x[t])))
acf(tt)
```


## Differencing a biased random walk

First-differencing a biased random walk yields a constant mean (level) $u$ plus white noise

$$
\begin{align}
  \nabla x_t &= x_{t-1} + u + w_t \\
  x_t - x_{t-1} &= x_{t-1} + u + w_t - x_{t-1} \\
  x_t - x_{t-1} &= u + w_t
\end{align}
$$


## Differencing a biased random walk

```{r ex_diff_biased_rw, fig.cap="$x_t = x_{t-1} + 1 + w_t; w_t \\sim \\text{N}(0,1)$"}
par(mfrow = c(1,2), mai = c(1.5,0.9,0.1,0.1), omi = c(0,0,0,0))
plot.ts(diff(xx), ylab = expression(nabla~italic(x[t])))
acf(tt)
```


## {.flexbox .vcenter}

<font size="10">LINEAR STATIONARY MODELS</font>


## Linear stationary models

We saw last week that linear filters are a useful way of modeling time series

Here we extend those ideas to a general class of models call _autoregressive moving average_ (ARMA) models


## Autoregressive (AR) models

Autoregressive models are widely used in ecology to treat a current state of nature as a function its past state(s)


## Autoregressive (AR) models

An _autoregressive_ model of order _p_, or AR(_p_), is defined as

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

where we assume

1. $w_t$ is white noise

2. $\phi_p \neq 0$ for an order-_p_ process


## Examples of AR(_p_) models

AR(1)

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

<br>
AR(1) with $\phi_1 = 1$ (random walk)

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

<br>
AR(2)

$x_t = -0.2 x_{t-1} + 0.4 x_{t-2} + w_t$


## Examples of AR(_p_) models

```{r}
set.seed(123)
## the 4 AR coefficients
ARp <- c(0.7, 0.2, -0.1, -0.3)
## empty list for storing models
AR_mods <- vector("list", 4L)

par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(p in 1:4) {
  ## assume SD=1, so not specified
  AR_mods[[p]] <- arima.sim(n=50, list(ar=ARp[1:p]))
  plot.ts(AR_mods[[p]], las = 1,
          ylab = expression(italic(x[t])))
  mtext(side = 3, paste0("AR(",p,")"),
        line = 0.5, adj = 0)
}
```


## Stationary AR(_p_) models

Recall that _stationary_ processes have the following properties

1. no systematic change in the mean or variance  
2. no systematic trend  
3. no periodic variations or seasonality

We seek a means for identifying whether our AR(_p_) models are also stationary


## Stationary AR(_p_) models

We can write out an AR(_p_) model using the backshift operator

$$
  x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \dots + \phi_p x_{t-p} + w_t \\
  \Downarrow \\
\begin{align}
  x_t - \phi_1 x_{t-1} - \phi_2 x_{t-2} - \dots - \phi_p x_{t-p} &= w_t \\
  (1 - \phi_1 \mathbf{B} - \phi_2 \mathbf{B}^2 - \dots - \phi_p \mathbf{B}^p) x_t &= w_t \\
  \phi_p (\mathbf{B}) x_t &= w_t \\
\end{align}
$$

## Stationary AR(_p_) models

If we treat $\mathbf{B}$ as a number (or numbers), we can out write the _characteristic equation_ as

$$
\phi_p (\mathbf{B}) x_t = w_t \\
\Downarrow \\
\phi_p (\mathbf{B}) = 0
$$

To be stationary, __all roots__ of the characteristic equation __must exceed 1 in absolute value__


## Stationary AR(_p_) models

For example, consider this AR(1) model from earlier

$$
\begin{align}
  x_t &= 0.5 x_{t-1} + w_t \\
  x_t - 0.5 x_{t-1} &= w_t \\
  (1 - 0.5 \mathbf{B})x_t &= w_t \\
\end{align}
$$


## Stationary AR(_p_) models

For example, consider this AR(1) model from earlier

$$
\begin{align}
  x_t &= 0.5 x_{t-1} + w_t \\
  x_t - 0.5 x_{t-1} &= w_t \\
  (1 - 0.5 \mathbf{B})x_t &= w_t \\
  \Downarrow \\
  1 - 0.5 \mathbf{B} &= 0 \\
  -0.5 \mathbf{B} &= -1 \\
  \mathbf{B} &= 2 \\
\end{align}
$$

This model is indeed stationary because $\mathbf{B} > 1$


## Stationary AR(_p_) models

What about this AR(2) model from earlier?

$$
\begin{align}
  x_t &= -0.2 x_{t-1} + 0.4 x_{t-2} + w_t \\
  x_t + 0.2 x_{t-1} - 0.4 x_{t-2} &= w_t \\
  (1 + 0.2 \mathbf{B} - 0.4 \mathbf{B}^2)x_t &= w_t \\
\end{align}
$$


## Stationary AR(_p_) models

What about this AR(2) model from earlier?

$$
\begin{align}
  x_t &= -0.2 x_{t-1} + 0.4 x_{t-2} + w_t \\
  x_t + 0.2 x_{t-1} - 0.4 x_{t-2} &= w_t \\
  (1 + 0.2 \mathbf{B} - 0.4 \mathbf{B}^2)x_t &= w_t \\
  \Downarrow \\
  1 + 0.2 \mathbf{B} - 0.4 \mathbf{B}^2 &= 0 \\
  \Downarrow \\
  \mathbf{B} \approx -1.35 ~ \text{and}& ~ \mathbf{B} \approx 1.85
\end{align}
$$

This model is _not_ stationary because only one $\mathbf{B} > 1$


## What about random walks?

Consider our random walk model

$$
\begin{align}
  x_t &= x_{t-1} + w_t \\
  x_t - x_{t-1} &= w_t \\
  (1 - 1 \mathbf{B})x_t &= w_t \\
\end{align}
$$


## What about random walks?

Consider our random walk model

$$
\begin{align}
  x_t &= x_{t-1} + w_t \\
  x_t - x_{t-1} &= w_t \\
  (1 - 1 \mathbf{B})x_t &= w_t \\
  \Downarrow \\
  1 - 1 \mathbf{B} &= 0 \\
  -1 \mathbf{B} &= -1 \\
  \mathbf{B} &= 1 \\
\end{align}
$$

Random walks are __not__ stationary because $\mathbf{B} = 1 \ngtr 1$


## Stationary AR(1) models

We can define a space over which all AR(1) models are stationary


## Stationary AR(1) models

For $x_t = \phi x_{t-1} + w_t$, we have

$$
\begin{align}
  1 - \phi \mathbf{B} &= 0 \\
  -\phi \mathbf{B} &= -1 \\
  \mathbf{B} &= \frac{1}{\phi} > 1 \Rightarrow 0 < \phi < 1\\
\end{align}
$$


## Stationary AR(1) models

For $x_t = \phi x_{t-1} + w_t$, we have

$$
\begin{align}
  1 - \phi \mathbf{B} &= 0 \\
  -\phi \mathbf{B} &= -1 \\
  \mathbf{B} &= \frac{1}{\phi} > 1 \Rightarrow 0 < \phi < 1\\
\end{align}
$$

For $x_t = -\phi x_{t-1} + w_t$, we have

$$
\begin{align}
  1 + \phi \mathbf{B} &= 0 \\
  \phi \mathbf{B} &= -1 \\
  \mathbf{B} &= \frac{-1}{\phi} > 1 \Rightarrow -1 < \phi < 0\\
\end{align}
$$


## Stationary AR(1) models

Thus, AR(1) models are stationary if and only if $\lvert \phi \rvert < 1$  


## Coefficients of AR(1) models

```{r ar_comp_pos_neg, fig.height=4}
set.seed(123)
## list description for AR(1) model with small coef
AR_pos <- list(order=c(1,0,0), ar=0.7, sd=0.1)
## list description for AR(1) model with large coef
AR_neg <- list(order=c(1,0,0), ar=-0.7, sd=0.1)
## simulate AR(1)
AR1_pos <- arima.sim(n=500, model=AR_pos)
AR1_neg <- arima.sim(n=500, model=AR_neg)

## get y-limits for common plots
ylm1 <- c(min(AR1_pos[1:50],AR1_neg[1:50]), max(AR1_pos[1:50],AR1_neg[1:50]))

## set the margins & text size
par(mfrow=c(1,2), mai=c(0.8,0.8,0.3,0.2), oma=c(0,0,0,0))
## plot the ts
plot.ts(AR1_pos[1:50], ylim=ylm1, las = 1,
        ylab=expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = 0.7")),
      line = 0.4, adj = 0)
plot.ts(AR1_neg[1:50], ylim=ylm1, las = 1,
        ylab=expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = -0.7")),
      line = 0.4, adj = 0)
```

Same value, but different sign


## Coefficients of AR(1) models

```{r ar_comp_sm_big, fig.height=4}
set.seed(123)
## list description for AR(1) model with small coef
AR_bg <- list(order=c(1,0,0), ar=0.9, sd=0.1)
## list description for AR(1) model with large coef
AR_sm <- list(order=c(1,0,0), ar=0.1, sd=0.1)
## simulate AR(1)
AR1_bg <- arima.sim(n=500, model=AR_bg)
AR1_sm <- arima.sim(n=500, model=AR_sm)

## get y-limits for common plots
ylm2 <- c(min(AR1_bg[1:50],AR1_sm[1:50]), max(AR1_bg[1:50],AR1_sm[1:50]))

## set the margins & text size
par(mfrow = c(1,2), mai = c(0.8,0.8,0.3,0.2), oma = c(0,0,0,0))
## plot the ts
plot.ts(AR1_bg[1:50], ylim = ylm2, las = 1,
        ylab = expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = 0.9")),
      line = 0.4, adj = 0)
plot.ts(AR1_sm[1:50], ylim = ylm2, las = 1,
        ylab = expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = 0.1")),
      line = 0.4, adj = 0)
```

Both positive, but different magnitude


## Autocorrelation function (ACF)

Recall that the _autocorrelation function_ ($\rho_k$) measures the correlation between $\{x_t\}$ and a shifted version of itself $\{x_{t+k}\}$ 


## ACF for AR(1) models

```{r}
## set the margins & text size
par(mfrow=c(2,2), mai=c(0.8,0.8,0.3,0.2), oma=c(0,0,0,0))
## plot the ts
plot.ts(AR1_pos[1:50], ylim=ylm1, las = 1,
        ylab=expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = 0.7")),
      line = 0.4, adj = 0)
acf(AR1_pos, lag.max = 20, las = 1)
plot.ts(AR1_neg[1:50], ylim=ylm1, las = 1,
        ylab=expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = -0.7")),
      line = 0.4, adj = 0)
acf(AR1_neg, lag.max = 20, las = 1)
```

ACF oscillates for model with $-\phi$

## ACF for AR(1) models

```{r}
## set the margins & text size
par(mfrow = c(2,2), mai = c(0.8,0.8,0.3,0.2), oma = c(0,0,0,0))
## plot the ts
plot.ts(AR1_bg[1:50], ylim = ylm2, las = 1,
        ylab = expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = 0.9")),
      line = 0.4, adj = 0)
acf(AR1_bg, lag.max = 20, las = 1)
plot.ts(AR1_sm[1:50], ylim = ylm2, las = 1,
        ylab = expression(italic(x)[italic(t)]),
        main = "")
mtext(side = 3, expression(paste(phi[1]," = 0.1")),
      line = 0.4, adj = 0)
acf(AR1_sm, lag.max = 20, las = 1)
```

For model with large $\phi$, ACF has longer tail

## Partial autocorrelation funcion (PACF)

Recall that the _partial autocorrelation function_ ($\phi_k$) measures the correlation between $\{x_t\}$ and a shifted version of itself $\{x_{t+k}\}$, with the linear dependence of $\{x_{t-1},x_{t-2},\dots,x_{t-k-1}\}$ removed


## ACF & PACF for AR(_p_) models

```{r}
## set 3 AR coefficients
ARp3 <- list(c(0.7, 0.2, -0.1), c(-0.7, 0.2, 0.1))

expr <- list(expression(paste("AR(3) with ", phi[1], " = 0.7, ",
                              phi[2], " = 0.2, ", phi[3], " = -0.1")),
             expression(paste("AR(3) with ", phi[1], " = -0.7, ",
                              phi[2], " = 0.2, ", phi[3], " = 0.1")))

## empty list for storing models
AR3_mods <- vector("list", 2L)

par(mfrow = c(2,3), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(p in 1:2) {
  ## assume SD=1, so not specified
  AR3_mods[[p]] <- arima.sim(n=5000, list(ar=ARp3[[p]]))
  plot.ts(AR3_mods[[p]][1:50], las = 1,
          ylab = expression(italic(x[t])))
  acf(AR3_mods[[p]], lag.max = 20,
      las = 1, main = "")
  mtext(side = 3, expr[[p]],
        line = 0.5, adj = 0.5)
  pacf(AR3_mods[[p]], lag.max = 20,
       las = 1, main = "")
}
```


## PACF for AR(_p_) models

```{r}
## empty list for storing models
pacf_mods <- vector("list", 4L)

par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(p in 1:4) {
  pacf_mods[[p]] <- arima.sim(n=5000, list(ar=ARp[1:p]))
  pacf(pacf_mods[[p]], lag.max = 15,
       las = 1, main = "")
  mtext(side = 3, paste0("AR(",p,")"),
        line = 0.5, adj = 0)
}

```

Do you see the link between the order _p_ and lag _k_?


## Using ACF & PACF for model ID

| Model   | ACF               | PACF                    |
|:-------:|:-----------------:|:-----------------------:|
| AR(_p_) | Tails off slowly  | Cuts off after lag _p_  |


## Moving average (MA) models

Moving average models are most commonly used for forecasting a future state


## Moving average (MA) models

A moving average model of order _q_, or MA(_q_), is defined as

$$
x_t = w_t + \theta_1 w_{t-1} + \theta_2 w_{t-2} + \dots + \theta_q w_{t-q}
$$

where $w_t$ is white noise

Each of the $x_t$ is a sum of the most recent error terms


## Moving average (MA) models

A moving average model of order _q_, or MA(_q_), is defined as

$$
x_t = w_t + \theta_1 w_{t-1} + \theta_2 w_{t-2} + \dots + \theta_q w_{t-q}
$$

where $w_t$ is white noise

Each of the $x_t$ is a sum of the most recent error terms

Thus, _all_ MA processes are stationary because they are finite sums of stationary WN processes


## Examples of MA(_q_) models

```{r}
## the 4 MA coefficients
MAq <- c(0.7, 0.2, -0.1, -0.3)
## empty list for storing models
MA_mods <- vector("list", 4L)
## loop over orders of q
for(q in 1:4) {
  ## assume SD=1, so not specified
  MA_mods[[q]] <- arima.sim(n=500, list(ma=MAq[1:q]))
}

par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(q in 1:4) {
  ## assume SD=1, so not specified
  plot.ts(MA_mods[[q]][1:50], las = 1,
          ylab = expression(italic(x[t])))
  mtext(side = 3, paste0("MA(",q,")"),
        line = 0.5, adj = 0)
}

```


## Examples of MA(_q_) models

```{r}
## compare MA(1) & MA(2) with similar structure
MA1 <- arima.sim(n=50, list(ma=c(0.7)))
MA2 <- arima.sim(n=50, list(ma=c(-1, 0.7)))

par(mfrow = c(1,2), mai = c(1,1,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
plot.ts(MA1, las = 1,
        ylab = expression(italic(x[t])))
mtext(side = 3,
      expression(MA(1):~italic(x[t])==~italic(w[t])+0.7~italic(w[t-1])),
      line = 0.5, adj = 0)
plot.ts(MA2, las = 1,
        ylab = expression(italic(x[t])))
mtext(side = 3,
      expression(MA(2):~italic(x[t])==~italic(w[t])-~italic(w[t-1])+0.7~italic(w[t-2])),
      line = 0.5, adj = 0)

```


## AR(_p_) model as an MA($\infty$) model

It is possible to write an AR(_p_) model as an MA($\infty$) model


## AR(1) model as an MA($\infty$) model

For example, consider an AR(1) model

$$
\begin{align}
x_t &= \phi x_{t-1} + w_t \\
x_t &= \phi (\phi x_{t-2} + w_{t-1}) + w_t \\
x_t &= \phi^2 x_{t-2} + \phi w_{t-1} + w_t \\
x_t &= \phi^3 x_{t-3} + \phi^2 w_{t-2} + \phi w_{t-1} + w_t \\
& \Downarrow \\
x_t &= w_t + \phi w_{t-1}+ \phi^2 w_{t-2} + \dots + \phi^k w_{t-k} + \phi^{k+1} x_{t-k-1}
\end{align}
$$


## AR(1) model as an MA($\infty$) model

If our AR(1) model is stationary, then

$$
\lvert \phi \rvert < 1 ~ \Rightarrow ~ \lim_{k \to \infty} \phi^{k+1} = 0
$$

so

$$
\begin{align}
  x_t &= w_t + \phi w_{t-1}+ \phi^2 w_{t-2} + \dots + \phi^k w_{t-k} + \phi^{k+1} x_{t-k-1} \\
  & \Downarrow \\
  x_t &= w_t + \phi w_{t-1}+ \phi^2 w_{t-2} + \dots + \phi^k w_{t-k}
\end{align}
$$
$$
$$



## Invertible MA(_q_) models

An MA(_q_) process is invertible if it can be written as a _stationary autoregressive process of infinite order without an error term_


## Invertible MA(1) model

For example, consider an MA(1) model

$$
\begin{align}
  x_t &= w_t + \theta w_{t-1} \\
  & \Downarrow \\
  w_t &= x_t - \theta w_{t-1} \\
  w_t &= x_t - \theta (x_{t-1} - \theta w_{t-2}) \\
  w_t &= x_t - \theta x_{t-1} - \theta^2 w_{t-2} \\
  & ~~\vdots \\
  w_t &= x_t - \theta x_{t-1} + \dots + (-\theta)^k x_{t-k}  + (-\theta)^{k+1} w_{t-k-1} \\
\end{align}
$$


## Invertible MA(1) model

If we constrain $\lvert \theta \rvert < 1$, then

$$
\lim_{k \to \infty} (-\theta)^{k+1} w_{t-k-1} = 0
$$

and

$$
\begin{align}
  w_t &= x_t - \theta x_{t-1} + \dots + (-\theta)^k x_{t-k}  + (-\theta)^{k+1} w_{t-k-1} \\
  & \Downarrow \\
  w_t &= x_t - \theta x_{t-1} + \dots + (-\theta)^k x_{t-k} \\
  w_t &= x_t + \sum_{k=1}^\infty(-\theta)^k x_{t-k}
\end{align}
$$


## Invertible MA(_q_) models

Q: Why do we care if an MA(_q_) model is invertible?

A: It helps us identify the model's parameters


## Invertible MA(_q_) models

For example, these MA(1) models are equivalent

$$
x_t = w_t + \frac{1}{5} w_{t-1}, ~\text{with} ~w_t \sim ~\text{N}(0,25)
$$

$$
x_t = w_t + 5 w_{t-1}, ~\text{with} ~w_t \sim ~\text{N}(0,1)
$$


## ACF & PACF for MA(_q_) models

```{r}
## set 3 AR coefficients
MAp3 <- list(c(0.7), c(-0.7, 0.2, 0.1))

expr <- list(expression(paste("MA(1) with ", theta[1], " = 0.7, ")),
             expression(paste("MA(3) with ", theta[1], " = -0.7, ",
                              theta[2], " = 0.2, ", theta[3], " = 0.1")))

## empty list for storing models
MA3_mods <- vector("list", 2L)

par(mfrow = c(2,3), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(p in 1:2) {
  ## assume SD=1, so not specified
  MA3_mods[[p]] <- arima.sim(n=50, list(ma=MAp3[[p]]))
  plot.ts(MA3_mods[[p]][1:50], las = 1,
          ylab = expression(italic(x[t])))
  acf(MA3_mods[[p]], lag.max = 20,
      las = 1, main = "")
  mtext(side = 3, expr[[p]],
        line = 0.5, adj = 0.5)
  pacf(MA3_mods[[p]], lag.max = 20,
       las = 1, main = "")
}
```


## ACF for MA(_q_) models

```{r}
## the 4 MA coefficients
MAq <- c(0.7, 0.2, -0.1, -0.3)

## empty list for storing models
acf_mods <- vector("list", 4L)

par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(q in 1:4) {
  acf_mods[[p]] <- arima.sim(n=5000, list(ma=MAq[1:q]))
  acf(acf_mods[[p]], lag.max = 15,
       las = 1, main = "")
  mtext(side = 3, paste0("MA(",q,")"),
        line = 0.5, adj = 0)
}
```

Do you see the link between the order _q_ and lag _k_?


## Using ACF & PACF for model ID

| Model   | ACF               | PACF                    |
|:-------:|:-----------------:|:-----------------------:|
| AR(_p_) | Tails off slowly  | Cuts off after lag _p_  |
| MA(_q_) | Cuts off after lag _q_  | Tails off slowly  |


## Using ACF & PACF for model ID



## Autoregressive moving average models

An autoregressive moving average, or ARMA(_p_,_q_), model is written as

$$
x_t = \phi_1 x_{t-1} + \dots + \phi_p x_{t-p} + w_t + \theta_1 w_{t-1} + \dots + \theta_q w_{t-q} 
$$


## Autoregressive moving average models

We can write an ARMA(_p_,_q_) model using the backshift operator

$$
\phi_p (\mathbf{B}) x_t=  \theta_q (\mathbf{B}) w_t 
$$

ARMA models are _stationary_ if all roots of $\phi_p (\mathbf{B}) > 1$

ARMA models are _invertible_ if all roots of $\theta_q (\mathbf{B}) > 1$


## Examples of ARMA(_p_,_q_) models

```{r}
arma_mods <- vector("list", 4L)

## ARMA(3,1): phi[1] = 0.7, phi[2] = 0.2, phi[3] = -0.1, theta[1]= 0.5
arma_mods[[1]] <- arima.sim(list(ar=c(0.7, 0.2, -0.1), ma=c(0.5)), n=5000)
## ARMA(2,2): phi[1] = -0.7, phi[2] = 0.2, theta[1] = 0.7, theta[2]= 0.2
arma_mods[[2]] <- arima.sim(list(ar=c(-0.7, 0.2), ma=c(0.7, 0.2)), n=5000)
## ARMA(1,3): phi[1] = 0.7, theta[1] = 0.7, theta[2]= 0.2, theta[3] = 0.5
arma_mods[[3]] <- arima.sim(list(ar=c(0.7), ma=c(0.7, 0.2, 0.5)), n=5000)
## ARMA(2,2): phi[1] = 0.7, phi[2] = 0.2, theta[1] = 0.7, theta[2]= 0.2
arma_mods[[4]] <- arima.sim(list(ar=c(0.7, 0.2), ma=c(0.7, 0.2)), n=5000)

titles <- list(
  expression("ARMA(3,1): "*phi[1]*" = 0.7, "*phi[2]*" = 0.2, "*phi[3]*" = -0.1, "*theta[1]*" = 0.5"),
  expression("ARMA(2,2): "*phi[1]*" = -0.7, "*phi[2]*" = 0.2, "*theta[1]*" = 0.7, "*theta[2]*" = 0.2"),
  expression("ARMA(1,3): "*phi[1]*" = 0.7, "*theta[1]*" = 0.7, "*theta[2]*" = 0.2, "*theta[3]*" = 0.5"),
  expression("ARMA(2,2): "*phi[1]*" = 0.7, "*phi[2]*" = 0.2, "*theta[1]*" = 0.7, "*theta[2]*" = 0.2")
)

par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(i in 1:4) {
  plot.ts(arma_mods[[i]][1:50], las = 1,
          main = "", ylab = expression(italic(x[t])))
  mtext(side = 3, titles[[i]],
        line = 0.5, adj = 0, cex = 0.8)
  
}
```


## ACF for ARMA(_p_,_q_) models

```{r}
par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(i in 1:4) {
  acf(arma_mods[[i]][1:1000], las = 1,
          main = "")
  mtext(side = 3, titles[[i]],
        line = 0.5, adj = 0, cex = 0.8)
  
}
```


## PACF for ARMA(_p_,_q_) models

```{r}
par(mfrow = c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))
## loop over orders of p
for(i in 1:4) {
  pacf(arma_mods[[i]][1:1000], las = 1,
          main = "")
  mtext(side = 3, titles[[i]],
        line = 0.5, adj = 0, cex = 0.8)
  
}
```


## Using ACF & PACF for model ID

| Model   | ACF               | PACF                    |
|:-------:|:-----------------:|:-----------------------:|
| AR(_p_) | Tails off slowly  | Cuts off after lag _p_  |
| MA(_q_) | Cuts off after lag _q_  | Tails off slowly  |
| ARMA(_p_,_q_) | Tails off slowly  | Tails off slowly  |


## {.flexbox .vcenter}

<font size="10">NONSTATIONARY MODELS</font>


## Autoregressive integrated moving average (ARIMA) models

If the data do not appear stationary, differencing can help

This leads to the class of _autoregressive integrated moving average_ (ARIMA) models

ARIMA models are indexed with orders (_p_,_d_,_q_) where _d_ indicates the order of differencing


## ARIMA(_p_,_d_,_q_) models

For $d > 0$, $\{x_t\}$ is an ARIMA(_p_,_d_,_q_) process if $(1-\mathbf{B})^d x_t$ is an ARMA(_p_,_q_) process


## ARIMA(_p_,_d_,_q_) models

For $d > 0$, $\{x_t\}$ is an ARIMA(_p_,_d_,_q_) process if $(1-\mathbf{B})^d x_t$ is an ARMA(_p_,_q_) process

For example, if $\{x_t\}$ is an ARIMA(1,1,0) process then $\nabla \{x_t\}$ is an ARMA(1,0) = AR(1) process


## ARIMA(_p_,_d_,_q_) models

```{r}
xx <- arima.sim(model=list(ar=0.5, sd=0.1), n=100)

yy <- cumsum(xx)

par(mfrow=c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))

plot.ts(yy, las = 1,
        ylab=expression(italic(x[t])))
mtext(side = 3, "ARIMA(1,1,0)", line = 0.5, adj = 0, cex = 0.8)
acf(yy)
```


## ARIMA(_p_,_d_,_q_) models

```{r}
par(mfrow=c(2,2), mai = c(0.7,0.7,0.3,0.1), omi = c(0,0,0,0))

plot.ts(yy, las = 1,
        ylab=expression(italic(x[t])))
mtext(side = 3, "ARIMA(1,1,0)", line = 0.5, adj = 0, cex = 0.8)
acf(yy)

plot.ts(diff(yy), las = 1,
        ylab=expression(paste(symbol("\xd1"), italic(x[t]))))
mtext(side = 3, "ARMA(1,0)", line = 0.5, adj = 0, cex = 0.8)
acf(diff(yy))
```


## Topics for today

Review

* White noise  

* Random walks  

Autoregressive (AR) models

Moving average (MA) models

Autoregressive moving average (ARMA) models

Using ACF & PACF for model ID
