```{r 09_setup, include=FALSE}
knitr::opts_chunk$set(echo=TRUE, eval=FALSE, error=TRUE)
```

# 9 ARIMA and SARIMA Models {-}

## Learning Goals {-}

- Explain and illustrate how first and second order differences remove linear and quadratic trends.
- Explain and illustrate how seasonal differences remove seasonality.
- Fit ARIMA and SARIMA models to real data and use model selection tools to decide on a final model.
- Explain and illustrate the general theoretical procedure for doing forecasting.

<br>

Slides from today are available [here](https://docs.google.com/presentation/d/1qLB-yqbbb0XXMuoo50eLDwMG5QW40ihi4j6gesDhTFk/edit?usp=sharing).

## Backshift Operators {-}

It can be useful to use the backshift operator, $B$, to simplify the notation of our models, such that 

$$B^k Y_t = Y_{t-k}\text{ for }k = 0,1,2,3,4,...$$


Our $AR(p)$ can be written as 

$$(1 - \phi_1B - \cdots-\phi_pB^p)Y_t = W_t$$

Our $MA(q)$ can be written as 

$$Y_t = (1 + \theta_1B + \cdots+\theta_qB^q)W_t$$

Our Seasonal $AR(P)$ can be written as 

$$(1 - \Phi_1B^S - \cdots-\Phi_PB^{PS})Y_t = W_t$$

Our Seasonal $MA(Q)$ can be written as 

$$Y_t = (1 + \Theta_1B^S + \cdots+\Theta_QB^{QS})W_t$$

All together, an ARMA(p,q) + Seasonal ARMA(P,Q)

$$(1 - \Phi_1B^S - \cdots-\Phi_PB^{PS})(1 - \phi_1B - \cdots-\phi_pB^p)Y_t = (1 + \Theta_1B^S + \cdots+\Theta_QB^{QS})(1 + \theta_1B + \cdots+\theta_qB^q)W_t$$

Another fun fact about the backshift operator:

- A first order difference to remove a linear trend:

$$Y_t - Y_{t-1} = (1-B)Y_t$$

- A second order difference to remove a quadratic trend:

$$(Y_t - Y_{t-1}) - (Y_{t-1} - Y_{t-2}) = (1-B)^2Y_t$$

- A first order seasonal difference to remove seasonality:

$$Y_t - Y_{t-S} = (1-B^S)Y_t$$

- A combination of first order for trend and seasonality:

$$(1-B)(1-B^S)Y_t = (1 - B -  B^{S} +  B^{S+1}) Y_t = (Y_t - Y_{t-1}) - ( Y_{t-S}-Y_{t-S-1})  $$

## Simulate that Model {-}

**Download a template RMarkdown file to start from [here](template_rmds/09-arima-sarima-models.Rmd).**

Below, there is three R code chunks to generate three time series and show the theoretical and estimated ACF.

1. For each, write out the data generation model in random process notation like above (using $\phi$'s and $\theta$'s). 
2. For each, compare and contrast the estimated ACF & PACF and the theoretical ACF & PACF.
3. For each, add R code to fit the data generation model to the data.
4. For each, add R code to forecast into the future 6 time units.


```{r}
library(dplyr)
library(astsa)

y1 <- sarima.sim(ar = c(0.9, -0.1))
plot(y1, type='l')

#estimated ACF and PACF
acf2(y1)

#theoretical ACF and PACF
ARMAacf(ar = c(0.9, -0.1),lag.max = 20) %>% plot(type='h',ylim=c(-1,1))
ARMAacf(ar = c(0.9, -0.1), pacf=TRUE,lag.max = 20) %>% plot(type='h',ylim=c(-1,1))

#fit model to data

#forecast


```



```{r}
y2 <- sarima.sim(ma = c(0.7), d = 1, D = 1, S = 6)
plot(y2, type='l')

#estimated ACF and PACF
y2 %>% diff(lag = 1) %>% diff(lag = 6) %>% acf2()

#theoretical ACF and PACF
ARMAacf(ma = 0.7, lag.max = 20) %>% plot(type='h',ylim=c(-1,1))
ARMAacf(ma = 0.7, pacf=TRUE,lag.max = 20) %>% plot(type='h',ylim=c(-1,1))

#fit model to data

#forecast

```


```{r}
y3 <- sarima.sim(ar = c(0.9), sar = 0.5, S = 12)
plot(y3, type='l')

#estimated ACF and PACF
acf2(y3)

#theoretical ACF and PACF
ARMAacf(ar = c(0.9, rep(0, 10), 0.5, -0.9*0.5)) %>% plot(type='h',ylim=c(-1,1))
ARMAacf(ar = c(0.9, rep(0, 10), 0.5, -0.9*0.5), pacf=TRUE) %>% plot(type='h',ylim=c(-1,1))

#fit model to data

#forecast

```

