Modelos ARIMA

class: center, middle, inverse, title-slide

.title[
# Modelos ARIMA
]
.author[
### Erika R. Badillo <a href="mailto:erika.badilloen@unaula.edu.co" class="email">erika.badilloen@unaula.edu.co</a> 
]
.date[
### Econometría II Programa de Economía Universidad Autónoma Latinoamericana
]

---

body {
text-align: justify;
}

h1{
      margin-top: -1px;
      margin-bottom: -3px;
}

.small-code pre{
  margin-bottom: -10px;
  
}

.medium-code pre{
 margin-bottom: 2px;
 
} 
</style>

# En este tema

- [Introducción](#introduccion)

- [El enfoque Box-Jenkins](#BJ)

- [Modelos ARIMA estacionales](#estacional)

- [Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos](#r1)

- [Ejercicio aplicado en R: desempleo estacional de los Estados Unidos](#r2)

---
# Lecturas
- Maddala, GS. y Lahiri, K. (2009). *Introduction to econometrics*. 4a edición, Willey. Cap 13

- Enders, W. (2014). *Applied econometric time series*. 4th edition, Wiley. Cap 2, sección 8

- Pfaff, B. (2008). *Analysis integrated and cointegrated series with R*. 2th edition, Springer. Part I, sección 1.4

- Hyndman, R.J., y Athanasopoulos, G. (2021). [*Forecasting: principles and practice*](https://otexts.com/fpp3/), 3rd edition, OTexts: Melbourne, Australia.

Páginas webs
- https://finnstats.com/index.php/2021/04/26/timeseries-analysis-in-r/

---
name: introduccion
# Introducción

- Hemos estudiado cómo una serie de tiempo puede ser explicada o bien por su historia (AR(p)) o por choques contemporáneos o pasados (MA(q))

- También estudiamos que estos dos procesos pueden ponerser juntos en un proceso más general ARMA(p,q)

- Ahora estudiaremos los modelos ARIMA o el enfoque de Box-Jenkins para series de tiempo

- Este enfoque consiste en tres etapas:
 1. Identificación
 2. Estimación
 3. Diagnóstico
 
---
name: BJ
# El enfoque Box-Jenkins

- El enfoque de Box-Jenkins (BJ) es una de las metodologías más amplia para el análisis de las series de tiempo

- Los pasos básicos de la metodología BJ son:

 1. diferenciar la serie, de modo que se alcance la estacionariedad
 2. identificar modelos tentativos
 3. estimar los modelos
 4. diagnósticar los modelos y seleccionar el mejor
 5. usar el modelo para pronosticar
 
---
# El enfoque Box-Jenkins

**Primer paso** 

- Determinar si la series es estacionaria

 - correlograma
 - test de raices unitarias
- Si la serie no es estacionaria deferenciarla (cuantas veces sea necesario) para logra su estacionariedad

**Segundo paso** 

- Examinar el correlograma de la serie estacionaria para decidir los ordenes apropiados de los componentes AR y MA

**Tercer paso** 

- Estimación del modelo ARMA

**Cuarto paso** 

- Diagnóstico de los modelos para comprobar la ideonidad del modelo tentativo

**Quinto paso** 

- Realizar el pronóstico con el modelo ARIMA seleccionado

---
name: estacional
# Modelos ARIMA estacionales

Hasta ahora, hemos restringido la atención a datos no estacionales y modelo ARIMA no estacionales. Sin embargo, los modelos ARIMA también son capaces de modelar un amplio rango de dato estacionales. Los datos estacionales tienen la forma:

<img src="ARIMA_files/figure-html/unnamed-chunk-1-1.png" width="50%" style="display: block; margin: auto;" />

Un modelo estacional ARIMA incluye términos estacionales adicionales y tiene la siguiente estructura

`$$ARIMA \underbrace{(p,d,q)}_{\text{No estacional}} \underbrace{(P,D,Q)_m}_{\text{ Estacional}}$$`
donde `$m=$` el periodo estacional (ejemplo, número de observaciones por año)

---
# Modelos ARIMA estacionales
Para definir la estructura estacional de una serie de tiempo se sigue el mismo procedimiento que con la parte nos estacional, es decir viendo el correlograma

Correlogragra de un proceso `$ARIMA(0,0,0)(0,0,1)_{12}$`

- AC: un pico estadísticamente significativo en el rezago 12 pero no significancia de otros picos
- PAC: decae exponencialmente en los rezagos estacionales (es decir, en los rezagos 12, 24, 36,...)

Correlogragra de un proceso `$ARIMA(0,0,0)(1,0,0)_{12}$`

- AC: decae exponencialmente en los rezagos estacionales (es decir, en los rezagos 12, 24, 36,...)
- PAC: un pico estadísticamente significativo en el rezago 12 pero no significancia de otros picos

---
name: r1
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos

``` r
library(urca); library(tidyverse); library(forecast); library(stats); library(lmtest)
data(npext)
data <- npext |> select(year, gnpperca) |> mutate(year = as.numeric(year, format="%Y"))|> 
 filter(year>=1909 & year<=1988)
```

``` r
ggplot(data) +
  geom_line(aes(x=year, y=gnpperca), color="blue") +
  geom_hline(yintercept = 7.728, linetype = "dashed", color = "red") +
  labs(title="PIB real per capita (media = 7.7)", x="Años", y="") +
  theme(axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 12)) + 
  scale_x_continuous(n.breaks = 10)
```

La serie no es estacionaria:

- Existe una tendencia creciente
- Cuando se hace el Dickey-Fuller Aumentado indica que la serie no es estacionaria (al 1%)
- El correlograma muestra una estructura AR(1)

---
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos

Diferenciamos la serie

``` r
data <- data |> 
 mutate(dgnpperca = c(NA,diff(gnpperca)))

ggplot(data) +
  geom_line(aes(x=year, y=dgnpperca), color="blue") +
  geom_hline(yintercept = 0.01681, linetype = "dashed", color = "red") +
  labs(title="Primera diferencia del PIB real per capita (media = 0.017)", x="Años", y="") +
  theme(axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 12)) + 
  scale_x_continuous(n.breaks = 10)
```

La serie ya es estacionaria:

- No tiene tendencia aunque es algo volátil
- Cuando se hace el Dickey-Fuller Aumentado indica que la serie es estacionaria

---
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos

Ahora calculamos el correlograma de la serie estacionaria para determinar las estructuras ARIMA posibles

.small-code[
.pull-left-50[

``` r
ggAcf(data$dgnpperca, main="Función AC dif PIB real per capita", ylab="AC")
```

<img src="ARIMA_files/figure-html/unnamed-chunk-5-1.png" width="100%" style="display: block; margin: auto;" />
]

.pull-right-50[

``` r
ggPacf(data$dgnpperca, main="Función PAC dif PIB real per capita", ylab="PAC")
```

<img src="ARIMA_files/figure-html/unnamed-chunk-6-1.png" width="100%" style="display: block; margin: auto;" />
]
]

Parece ser un proceso ARIMA(1,1,1). Estimemos también un proceso más parsimonioso ARIMA(1,1,0) y utilicemos la función `auto.arima` para definir un tercer proceso y ponerlos a competir

---
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos

``` r
arima111 <- arima(data$gnpperca, order=c(1,1,1))
coeftest(arima111)
```

```

z test of coefficients:

Estimate Std. Error z value Pr(>|z|)
ar1 0.297515   0.232692  1.2786   0.2010
ma1 0.087575   0.233412  0.3752   0.7075
```

``` r
arima110 <- arima(data$gnpperca, order=c(1,1,0))
coeftest(arima110)
```

```

z test of coefficients:

Estimate Std. Error z value  Pr(>|z|)    
ar1  0.37186    0.10359  3.5898 0.0003309 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

``` r
arima_auto <- auto.arima(data$gnpperca, allowdrift = F)
coeftest(arima_auto)
```

```

z test of coefficients:

Estimate Std. Error z value Pr(>|z|) 
ar1 0.824304 0.092250 8.9355 < 2.2e-16 ***
ar2 -0.864219 0.080500 -10.7356 < 2.2e-16 ***
ma1 -0.621339 0.075901 -8.1862 2.697e-16 ***
ma2 0.959728 0.095625 10.0364 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos

Para comparar los modelos se debe analizar los residuales de cada uno y utilizar los criterios de información. Se debe entonces:

- calcular el correlograma
- realizar los tests: `checkresiduals` y Shapiro-Wilk de normalidad

En los tres modelos los residuales se comportan bien. Pasamos a mirar los criterios de información (el mejor modelo es que tenga menor valor)

``` r
AIC(arima111); AIC(arima110); AIC(arima_auto)
```

```
[1] -225.0382
```

```
[1] -226.8975
```

```
[1] -227.649
```

``` r
BIC(arima111); BIC(arima110); BIC(arima_auto)
```

```
[1] -217.9298
```

```
[1] -222.1586
```

```
[1] -215.8017
```

Con el AIC el ARIMA(2,1,2) es el mejor, pero con el BIC el ARIMA(1,1,0) es mejor. Sin embargo, el ARIMA(2,1,2) mostró alta significancia en los coeficientes estimados, por eso se selecciona este último modelo

---
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos

Una vez se ha seleccionado un modelo ARIMA, puede ser usado para predecir valores futuros de la variable de interes

Estas predicciones pueden ser calculadas recursivamente desde el predictor lineal

$$
`\begin{aligned}
Y_T(h) = & \rho_1 \overline{Y}_{T+h-1} + ... + \rho_1 \overline{Y}_{T+h-p} +\\
         & \epsilon_t + \theta_1 \epsilon_{t-T-1} + ... + \epsilon_{t-T-q}  
\end{aligned}`
$$

donde `$\overline{Y}_t = Y_t$` para `$t\leq T$` y `$\overline{Y}_{T+j} = Y_T(j)$` para `$j = 1,...,h-1$`

Este predictor es equivalente a

`$$Y_T(h) = \mu + \psi_h \epsilon_t + \psi_{h+1} \epsilon_{t-1} + \psi_{h+2} \epsilon_{t-2} + ...$$`

Cuando el horizonte de predicción `$h$` es mayor que el orden q del MA, la predicción son determinandas sólo por los términos autorregresivos

---
# Ejercicio aplicado en R: PIB real per cápita de los Estados Unidos
Vamos a predecir 10 años de la tasa de desempleo

``` r
arima212 <- arima(data$gnpperca, order=c(2,1,2))
data.pred <- data |> 
 add_row(year=c(1989:1998)) |> 
 mutate(arima212.pred = c(rep(NA, length(data$gnpperca)-1), data$gnpperca[length(data$gnpperca)], predict(arima212, n.ahead = 10)$pred),
 arima212.se = c(rep(NA, length(data$gnpperca)), predict(arima212, n.ahead = 10)$se),
 ic_ub = arima212.pred + qt(0.025, nobs(arima212)-length(arima212$coefficients), lower.tail = F)*arima212.se,
 ic_lb = arima212.pred - qt(0.025, nobs(arima212)-length(arima212$coefficients), lower.tail = F)*arima212.se)

ggplot(data.pred) + 
  geom_line(aes(x = year, y = gnpperca, color = "Observada"), linewidth = 1.5) + 
  geom_line(aes(x = year, y = arima212.pred, color = "Predicha"), linewidth = 1.5) + 
  geom_ribbon(aes(x = year, y = arima212.pred, ymin = ic_lb, ymax = ic_ub, fill="IC"), alpha = 0.1) +
  theme(legend.text = element_text(size = 14), text = element_text(size=16), legend.spacing.y = unit(-0.4, "cm"), legend.background=element_blank()) + 
  guides(shape = guide_legend(order = 2), col = guide_legend(order = 1)) + 
  scale_color_manual(name = "", values = c("Observada" = "darkblue", "Predicha" = "red")) + 
  scale_fill_manual(values = c(IC = "steelblue"), labels = c(IC = "IC 95%")) + 
  labs(title="Gráfico 1. PIB real per cápita", x = "Años", y ="", fill = "") + 
  scale_x_continuous(n.breaks = 10)
```

---
name: r2
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos

En este ejemplo se decribe la modelación de un ARIMA estacional usando los datos mensuales de desempleo para Estados Unidos. Se toman los datos para el sector de ocio y hostelería desde enero de 2001 a septiembre de 2019

``` r
data("us_employment")

leisure <- us_employment |> filter(Title == "Leisure and Hospitality", year(Month) > 2000) |>
 mutate(Employed = Employed/1000) |> select(Month, Employed)
autoplot(leisure, Employed) + labs(title = "Desempleo US: Sector de ocio", y="# de personas (en millones)")
```

Los datos son claramente no estacionarios, con una fuerte estacionalidad y tendencia no lineal

---
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos
Se toma una diferencia estacional en 12 meses, para eliminar esa estacionalidad

``` r
leisure <- leisure |> mutate(Employed.e12 = difference(leisure$Employed, 12))
ts.plot(leisure$Employed.e12)
```

.small-code[
.pull-left-50[

``` r
ggAcf(leisure$Employed.e12, main="Función AC", ylab="AC")
```

<img src="ARIMA_files/figure-html/unnamed-chunk-12-1.png" width="50%" style="display: block; margin: auto;" />
]

.pull-right-50[

``` r
ggPacf(leisure$Employed.e12, main="Función PAC", ylab="PAC")
```

<img src="ARIMA_files/figure-html/unnamed-chunk-13-1.png" width="50%" style="display: block; margin: auto;" />
]
]

---
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos

Se observa que no hay estacionariedad, así que se toman primeras diferencias y se vuelve a mirar el correlograma

``` r
leisure <- leisure |> mutate(dEmployed.e12 = difference(leisure$Employed.e12, 1))
ts.plot(leisure$dEmployed.e12)
```

.small-code[
.pull-left-50[

``` r
ggAcf(leisure$dEmployed.e12, main="Función AC", ylab="AC")
```

<img src="ARIMA_files/figure-html/unnamed-chunk-15-1.png" width="50%" style="display: block; margin: auto;" />
]

.pull-right-50[

``` r
ggPacf(leisure$dEmployed.e12, main="Función PAC", ylab="PAC")
```

<img src="ARIMA_files/figure-html/unnamed-chunk-16-1.png" width="50%" style="display: block; margin: auto;" />
]
]

---
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos

Potenciales modelos que surgen con base en el anterior correlograma:

- la significancia en el rezago 2 de la AC sugiere un MA(2) en la parte no estacional. La significancia en el rezago 12 en la AC sugiere un MA(1) en la parte estacional. El proceso sería un ARIMA(0,1,2)(0,1,1)
- si se usa la PAC para seleccionar la parte no estacional y la AC para seleccionar la parte estacional, surge un ARIMA(2,1,2)(0,1,1)
- También se incluye la selección atomática

.pull-left-50[

``` r
arima012011 <- arima(leisure$Employed, order=c(0,1,2),
 seasonal = list(order = c(0, 1, 1), period = 12))
arima012011
```

```

Call:
arima(x = leisure$Employed, order = c(0, 1, 2), seasonal = list(order = c(0, 
    1, 1), period = 12))

Coefficients:
         ma1     ma2     sma1
      0.2315  0.2167  -0.5006
s.e.  0.0707  0.0621   0.0814

sigma^2 estimated as 0.001433:  log likelihood = 391.45,  aic = -774.9
```
]

.pull-right-50[

``` r
arima210011 <- arima(leisure$Employed, order=c(2,1,0),
 seasonal = list(order = c(0, 1, 1), period = 12), method="ML")
arima210011
```

```

Call:
arima(x = leisure$Employed, order = c(2, 1, 0), seasonal = list(order = c(0, 
    1, 1), period = 12), method = "ML")

Coefficients:
         ar1     ar2     sma1
      0.2102  0.1941  -0.4969
s.e.  0.0683  0.0679   0.0788

sigma^2 estimated as 0.001425:  log likelihood = 392.09,  aic = -776.19
```
]

``` r
autoarima <- leisure |>
 model(ARIMA(Employed, stepwise = FALSE, approx = FALSE))
report(autoarima)
```

```
Series: Employed 
Model: ARIMA(2,1,0)(1,1,1)[12]

Coefficients:
         ar1     ar2    sar1     sma1
      0.1786  0.1855  0.3295  -0.7507
s.e.  0.0695  0.0679  0.1273   0.0936

sigma^2 estimated as 0.001415:  log likelihood=394.96
AIC=-779.92   AICc=-779.63   BIC=-763.14
```

---
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos

``` r
AIC(arima012011); AIC(arima210011)
```

```
[1] -774.8994
```

```
[1] -776.189
```

``` r
BIC(arima012011); BIC(arima210011)
```

```
[1] -761.4731
```

```
[1] -762.7627
```

``` r
glance(autoarima)
```

```
# A tibble: 1 × 8
 .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
 <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list> 
1 ARIMA(Employed,… 0.00142 395. -780. -780. -763. <cpl> <cpl> 
```

De acuerdo con los criterios de información, el mejor modelo es el de la selección automática

---
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos

Miremos los residuales del modelo

``` r
autoarima |> gg_tsresiduals(lag=36)
```

Un pequeño pero significante pico (en el rezago 11) de 36 es aún consistente con un comportamiento ruido blanco. Para estar seguros, se calcula el test Ljung-Box

``` r
augment(autoarima) |>
  features(.innov, ljung_box, lag=24, dof=4)
```

```
# A tibble: 1 × 3
 .model lb_stat lb_pvalue
 <chr> <dbl> <dbl>
1 ARIMA(Employed, stepwise = FALSE, approx = FALSE) 16.6 0.680
```
El pvalor es muy grande con lo que se confirma que los residuales son similares a un ruido blanco

---
# Ejercicio aplicado en R: desempleo estacional en los Estados Unidos

Haciendo la predicción con el modelo seleccionado

``` r
fabletools::forecast(autoarima, h=36) |>
  autoplot(leisure) +
  labs(title = "Predicción del desempleo en el sector de ocio US",
       y="# de personas (en millones)")
```