Capitulo_17.Rmd

# Temas adicionales en la regresión de series temporales {#TARST}

```{r, echo = F}
options(knitr.duplicate.label = "allow")
```

```{r, 882, child="_setup.Rmd"}
```

```{r, 883, eval=knitr::opts_knit$get("rmarkdown.pandoc.to") == "html", results='asis', echo=FALSE}
cat('<hr style="background-color:#03193b;height:2px">')
```

Este capítulo analiza los siguientes temas avanzados en la regresión de series de tiempo y demuestra cómo se pueden aplicar las técnicas básicas usando **R**:

+ *Autorregresiones vectoriales* (VAR). Se enfoca en usar VAR para pronosticar. Otra rama de la literatura se ocupa de los llamados *VAR estructurales*, que están más allá del alcance de este capítulo.
+ Previsiones multiperiodo. Esto incluye una discusión de pronósticos iterados y directos (multivariados).
+ La prueba *DF-GLS*, una modificación de la prueba *ADF* que tiene más potencia que esta última cuando la serie tiene componentes deterministas y está cerca de ser no estacionaria.
+ Análisis de cointegración con aplicación a tipos de interés a corto y largo plazo. Se demuestra cómo estimar un modelo de corrección de errores vectoriales.
+ Modelos de *Heterocedasticidad condicional autorregresiva* (ARCH). Se muestra cómo un modelo *ARCH* generalizado simple (*GARCH*) puede ser útil para cuantificar el riesgo asociado con la inversión en el mercado de valores en términos de estimación y pronóstico de la volatilidad de los rendimientos de los activos.

Para reproducir los ejemplos de código, se necesita instalar los paquetes **R** que se enumeran a continuación. Asegúrese que el siguiente fragmento de código se ejecute sin errores.

+ **AER** [@R-AER]
+ **dynlm** [@R-dynlm]
+ **fGarch** [@R-fGarch]
+ **quantmod** [@R-quantmod]
+ **readxl** [@R-readxl]
+ **scales** [@R-scales]
+ **vars** [@R-vars]

```{r, 884, warning=FALSE, message=FALSE}
library(AER)
library(readxl)
library(dynlm)
library(vars)
library(quantmod)
library(scales)
library(fGarch)
```

## Autorregresiones vectoriales

Un modelo de *Vector autorregresivo* (VAR) es útil cuando se está interesado en predecir múltiples variables de series de tiempo usando un solo modelo. En esencia, el modelo VAR es una extensión del modelo autorregresivo univariante que se ha tratado en los capítulos \@ref(IRSTP) y \@ref(EECD). El concepto clave 16.1 resume los fundamentos del VAR.

```{r, 885, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC16.1">
<h3 class = "right"> Concepto clave <br> 16.1</h3>          
<h3 class = "left">Autoregresiones vectoriales</h3>
<p>

El modelo de autorregresión vectorial (VAR) extiende la idea de autorregresión univariante a regresiones de series de tiempo $k$, donde los valores rezagados de *todas* $k$ series aparecen como regresores. Dicho de otra manera, en un modelo VAR se hace una regresión de un *vector* de variables de series de tiempo en vectores rezagados de estas variables. En cuanto a los modelos AR ($p$), el orden de retraso se denota por $p$ por lo que el modelo VAR ($p$) de dos variables $X_t$ y $Y_t$ ($k = 2$) viene dado por las ecuaciones:

\\begin{align*}
  Y_t =& \\, \\beta_{10} + \\beta_{11} Y_{t-1} + \\dots + \\beta_{1p} Y_{t-p} + \\gamma_{11} X_{t-1} + \\dots + \\gamma_{1p} X_{t-p} + u_{1t}, \\\\
  X_t =& \\, \\beta_{20} + \\beta_{21} Y_{t-1} + \\dots + \\beta_{2p} Y_{t-p} + \\gamma_{21} X_{t-1} + \\dots + \\gamma_{2p} X_{t-p} + u_{2t}.
\\end{align*}

Los $\\beta$s y $\\gamma$s se pueden estimar usando MCO en cada ecuación. Los supuestos para los VAR son los supuestos de series de tiempo presentados en el Concepto clave 14.6 aplicados a cada una de las ecuaciones.

Es sencillo estimar modelos VAR en <tt>R</tt>. Un enfoque factible es simplemente usar <tt>lm()</tt> para la estimación de las ecuaciones individuales. Además, el paquete <tt>R</tt> <tt>vars</tt> proporciona herramientas estándar para estimación, pruebas de diagnóstico y predicción utilizando este tipo de modelos.

</p>
</div>
')
```

```{r, 886, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Autoregresiones vectoriales]{16.1}

El modelo de autorregresión vectorial (VAR) extiende la idea de autorregresión univariante a regresiones de series de tiempo $k$, donde los valores rezagados de *todas* $k$ series aparecen como regresores. Dicho de otra manera, en un modelo VAR se hace una regresión de un \\textit{vector} de variables de series de tiempo en vectores rezagados de estas variables. En cuanto a los modelos AR ($p$), el orden de retraso se denota por $p$ por lo que el modelo VAR ($p$) de dos variables $X_t$ y $Y_t$ ($k = 2$) viene dado por las ecuaciones:

\\begin{align*}
  Y_t =& \\, \\beta_{10} + \\beta_{11} Y_{t-1} + \\dots + \\beta_{1p} Y_{t-p} + \\gamma_{11} X_{t-1} + \\dots + \\gamma_{1p} X_{t-p} + u_{1t}, \\\\
  X_t =& \\, \\beta_{20} + \\beta_{21} Y_{t-1} + \\dots + \\beta_{2p} Y_{t-p} + \\gamma_{21} X_{t-1} + \\dots + \\gamma_{2p} X_{t-p} + u_{2t}.
\\end{align*}

Los $\\beta$s y $\\gamma$s se pueden estimar usando MCO en cada ecuación. Los supuestos para los VAR son los supuestos de series de tiempo presentados en el Concepto clave 14.6 aplicados a cada una de las ecuaciones.\\newline

Es sencillo estimar modelos VAR en \\texttt{R}. Un enfoque factible es simplemente usar \\texttt{lm()} para la estimación de las ecuaciones individuales. Además, el paquete \\texttt{R} \\texttt{vars} proporciona herramientas estándar para estimación, pruebas de diagnóstico y predicción utilizando este tipo de modelos.

\\end{keyconcepts}
')
```

Cuando se cumplen los supuestos del Concepto clave 16.1, los estimadores MCO de los coeficientes VAR son consistentes y conjuntamente normales en muestras grandes, de modo que se pueden utilizar los métodos inferenciales habituales, como los intervalos de confianza y los estadísticos $t$.

La estructura de los VAR también permite probar conjuntamente las restricciones en múltiples ecuaciones. Por ejemplo, puede ser de interés probar si los coeficientes de todos los regresores del rezago $p$ son cero. Esto corresponde a probar el nulo de que el orden de retraso $p-1$ es correcto. La normalidad conjunta de muestras grandes de las estimaciones de coeficientes es conveniente porque implica que simplemente se puede usar una prueba $F$ para este problema de prueba. La fórmula explícita para un estadístico de prueba de este tipo es bastante complicada, pero afortunadamente estos cálculos se realizan fácilmente utilizando las funciones **R** con las que se trabaja en este capítulo. Otra forma de determinar las longitudes de retraso óptimas son los criterios de información como $BIC$ que se ha introducido para las regresiones de series de tiempo univariadas en el Capítulo \@ref(SLRUCI). Al igual que en el caso de una sola ecuación, para un modelo de ecuaciones múltiples se elige la especificación que tiene el menor $BIC(p)$, donde

\begin{align*}
  BIC(p) =& \, \log\left[\text{det}(\widehat{\Sigma}_u)\right] + k(kp+1) \frac{\log(T)}{T}.
\end{align*}

donde $\widehat{\Sigma}_u$ denota la estimación de la matriz de covarianza $k \times k$ de los errores VAR y $\text{det}(\cdot)$ denota el determinante.

En cuanto a los modelos de rezagos distribuidos univariados, se debe pensar detenidamente en las variables que se incluirán en un VAR, ya que agregar variables no relacionadas reduce la precisión del pronóstico al aumentar el error de estimación. Esto es particularmente importante porque el número de parámetros a estimar crece cuadráticamente al número de variables modeladas por el VAR. En la aplicación siguiente se verá que la teoría económica y la evidencia empírica son útiles para esta decisión.

#### Un modelo VAR de la tasa de crecimiento del PIB y el margen temporal {-}

Ahora se muestra cómo estimar un modelo VAR de la tasa de crecimiento del PIB, $GDPGR$, y el diferencial de plazo, $TSpread$. Como sigue la discusión sobre la no estacionariedad del crecimiento del PIB en el Capítulo \@ref(NEIT) (recuerde la posible ruptura a principios de la década de 1980 detectada por la estadística de prueba $QLR$), se usan datos de 1981:Q1 a 2012:Q4. Las dos ecuaciones del modelo son

\begin{align*}
 GDPGR_t =& \, \beta_{10} + \beta_{11} GDPGR_{t-1} + \beta_{12} GDPGR_{t-2} + \gamma_{11} TSpread_{t-1} + \gamma_{12} TSpread_{t-2} + u_{1t}, \\
 TSpread_t =& \, \beta_{20} + \beta_{21} GDPGR_{t-1} + \beta_{22} GDPGR_{t-2} + \gamma_{21} TSpread_{t-1} + \gamma_{22} TSpread_{t-2} + u_{2t}.
\end{align*}

El conjunto de datos **us_macro_quarterly.xlsx** se proporciona en el sitio web complementario a @stock2015 y se puede descargar [aquí](http://wps.aw.com/aw_stock_ie_3/178/45691/11696965.cw/index.html). Contiene datos trimestrales sobre el PIB real de EE. UU. (es decir, ajustado a la inflación) de 1947 a 2004. Se comienza importando el conjunto de datos y aplicando un formato (ya se trabajó con este conjunto de datos en el Capítulo \@ref(IRSTP), por lo que se pueden omitir estos pasos si ya se han cargado los datos en el entorno de trabajo).

```{r, 887, eval=FALSE}
# cargar el conjunto de datos macroeconómicos de EE. UU.
USMacroSWQ <- read_xlsx("data/us_macro_quarterly.xlsx",
                         sheet = 1,
                         col_types = c("text", rep("numeric", 9)))

# establecer los nombres de las columnas
colnames(USMacroSWQ) <- c("Date", "GDPC96", "JAPAN_IP", "PCECTPI", "GS10", 
                          "GS1", "TB3MS", "UNRATE", "EXUSUK", "CPIAUCSL")

# formatear la columna de fecha
USMacroSWQ$Date <- as.yearqtr(USMacroSWQ$Date, format = "%Y:0%q")

# definir el PIB como objeto ts
GDP <- ts(USMacroSWQ$GDPC96,
          start = c(1957, 1), 
          end = c(2013, 4), 
          frequency = 4)

# definir el crecimiento del PIB como un objeto ts
GDPGrowth <- ts(400*log(GDP[-1]/GDP[-length(GDP)]),
                start = c(1957, 2), 
                end = c(2013, 4), 
                frequency = 4)

# tasa de interés de las letras del Tesoro a 3 meses como objeto 'ts'
TB3MS <- ts(USMacroSWQ$TB3MS,
            start = c(1957, 1), 
            end = c(2013, 4), 
            frequency = 4)

# tasa de interés de los bonos del Tesoro a 10 años como objeto 'ts'
TB10YS <- ts(USMacroSWQ$GS10, 
              start = c(1957, 1), 
              end = c(2013, 4), 
              frequency = 4)

# generar la serie diferencial por plazo
TSpread <- TB10YS - TB3MS
```

```{r, 888, echo=F, purl=F}
library(xts)
# cargar datos macroeconómicos de EE. UU.
USMacroSWQ <- read_xlsx("data/us_macro_quarterly.xlsx",
                         sheet = 1,
                         col_types = c("text", rep("numeric", 9)))

# establecer nombres de columna
colnames(USMacroSWQ) <- c("Date", "GDPC96", "JAPAN_IP", "PCECTPI", "GS10", "GS1", "TB3MS", "UNRATE", "EXUSUK", "CPIAUCSL")

# formatear columna de fecha
USMacroSWQ$Date <- as.yearqtr(USMacroSWQ$Date, format = "%Y:0%q")

# eludir error
GDP <- xts(USMacroSWQ$GDPC96, USMacroSWQ$Date)["1960::2013"]
GDPGrowth <- xts(400 * log(GDP/lag(GDP)))
GDP <- ts(GDP,
          start = c(1960, 1), 
          end = c(2013, 4), 
          frequency = 4)

GDPGrowth <- ts(GDPGrowth,
                start = c(1960, 1), 
                end = c(2013, 4), 
                frequency = 4)

# tasa de interés de las letras del Tesoro a 3 meses como objeto 'ts'
TB3MS <- ts(USMacroSWQ$TB3MS,
            start = c(1957, 1), 
            end = c(2013, 4), 
            frequency = 4)

# tasa de interés de las letras del Tesoro a 10 meses como objeto 'ts'
TB10YS <- ts(USMacroSWQ$GS10, 
              start = c(1957, 1), 
              end = c(2013, 4), 
              frequency = 4)

# generar la serie de diferencial por plazo
TSpread <- TB10YS - TB3MS
```

Se estiman ambas ecuaciones por separado por MCO y se usa **coeftest()** para obtener errores estándar robustos.

```{r, 889}
# estimar ambas ecuaciones usando 'dynlm()'
VAR_EQ1 <- dynlm(GDPGrowth ~ L(GDPGrowth, 1:2) + L(TSpread, 1:2), 
                 start = c(1981, 1), 
                 end = c(2012, 4))

VAR_EQ2 <- dynlm(TSpread ~ L(GDPGrowth, 1:2) + L(TSpread, 1:2),
                 start = c(1981, 1),
                 end = c(2012, 4))

# cambiar el nombre de los regresores para una mejor legibilidad
names(VAR_EQ1$coefficients) <- c("Intercept", 
                                 "Growth_t-1", 
                                 "Growth_t-2", 
                                 "TSpread_t-1", 
                                 "TSpread_t-2")
names(VAR_EQ2$coefficients) <- names(VAR_EQ1$coefficients)

# resúmenes robustos de coeficientes
coeftest(VAR_EQ1, vcov. = sandwich)
coeftest(VAR_EQ2, vcov. = sandwich)
```

Se termina con los siguientes resultados:

\begin{align*}
 GDPGR_t =& \, \underset{(0.46)}{0.52} + \underset{(0.11)}{0.29} GDPGR_{t-1} + \underset{(0.09)}{0.22} GDPGR_{t-2} -\underset{(0.36)}{0.90} TSpread_{t-1} + \underset{(0.39)}{1.33} TSpread_{t-2} \\
 TSpread_t =& \, \underset{(0.12)}{0.46} + \underset{(0.02)}{0.01} GDPGR_{t-1} -\underset{(0.03)}{0.06} GDPGR_{t-2} + \underset{(0.10)}{1.06} TSpread_{t-1} -\underset{(0.11)}{0.22} TSpread_{t-2} 
\end{align*}

La función **VAR()** se puede utilizar para obtener las mismas estimaciones de coeficientes que se presentaron anteriormente, ya que también se aplica MCO por ecuación.

```{r, 890}
# configurar datos para la estimación usando `VAR()`
VAR_data <- window(ts.union(GDPGrowth, TSpread), start = c(1980, 3), end = c(2012, 4))

# estimar los coeficientes del modelo usando `VAR()`
VAR_est <- VAR(y = VAR_data, p = 2)
VAR_est
```

**VAR()** devuelve una **lista** de objetos **lm** que se pueden pasar a las funciones habituales, por ejemplo **summary()** y, por lo tanto, es sencillo obtener estadísticas del modelo para el ecuaciones individuales.

```{r, 891}
# obtener el R^2 ajustado de la salida de 'VAR()'
summary(VAR_est$varresult$GDPGrowth)$adj.r.squared
summary(VAR_est$varresult$TSpread)$adj.r.squared
```

Se pueden utilizar los objetos del modelo individual para realizar pruebas de causalidad de Granger.

```{r, 892}
# Pruebas de causalidad de Granger:

# probar si el diferencial de plazo no tiene poder para explicar el crecimiento del PIB
linearHypothesis(VAR_EQ1, 
                 hypothesis.matrix = c("TSpread_t-1", "TSpread_t-2"),
                 vcov. = sandwich)

# probar si el crecimiento del PIB no tiene poder para explicar el diferencial de plazo
linearHypothesis(VAR_EQ2, 
                 hypothesis.matrix = c("Growth_t-1", "Growth_t-2"),
                 vcov. = sandwich)
```

Ambas pruebas de causalidad de Granger se rechazan al nivel de $5\%$. Esto es evidencia a favor de la conjetura de que el diferencial por plazo tiene poder para explicar el crecimiento del PIB y viceversa.

#### Pronósticos multivariados iterados usando un VAR iterado {-}

La idea de un pronóstico iterado para el período $T + 2$ basado en observaciones hasta el período $T$ es utilizar el pronóstico de un período adelantado como paso intermedio; es decir, el pronóstico para el período $T + 1$ se usa como una observación al predecir el nivel de una serie para el período $T + 2$. Esto se puede generalizar a un pronóstico de $h$ para el período futuro en el que todos los períodos intermedios entre $T$ y $T + h$ deben pronosticarse, ya que se utilizan como observaciones en el proceso. Los pronósticos iterados de múltiples períodos se resumen en el Concepto clave 16.2.

```{r, 893, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC16.2">
<h3 class = "right"> Concepto clave <br> 16.2 </h3>          
<h3 class = "left"> Pronósticos de varios períodos iterados </h3>
<p>

Los pasos para un *pronóstico de RA iterado de varios períodos* son:

1. Estimar el modelo AR ($p$) usando MCO y calcular el pronóstico de un período adelantado.

2. Utilizar el pronóstico de un período adelantado para obtener el pronóstico de dos períodos adelante.

3. Continuar iterando para obtener pronósticos más lejanos en el futuro.

Un *pronóstico VAR iterado de múltiples períodos* se realiza de la siguiente manera:

1. Estimar el modelo VAR ($p$) usando MCO por ecuación y calcular el pronóstico de un período adelantado para *todas* las variables en el VAR.

2. Utilizar los pronósticos de un período por delante para obtener los pronósticos de dos períodos adelante.

3. Continuar iterando para obtener pronósticos de todas las variables en el VAR en el futuro.

</p>
</div>
')
```

```{r, 894, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Pronósticos de varios períodos iterados]{16.2}

Los pasos para un \\textit{pronóstico de RA iterado de varios períodos} son:\\newline

\\begin{enumerate}
\\item Estimar el modelo AR ($p$) usando MCO y calcular el pronóstico de un período adelantado.
\\item Utilizar el pronóstico de un período adelantado para obtener el pronóstico de dos períodos adelante.
\\item Continuar iterando para obtener pronósticos más lejanos en el futuro.
\\end{enumerate}\\vspace{0.5cm}

Un \\textit{pronóstico VAR iterado de múltiples períodos} se realiza de la siguiente manera:\\newline

\\begin{enumerate}
\\item Estimar el modelo VAR ($p$) usando MCO por ecuación y calcular el pronóstico de un período adelantado para \\textit{todas} las variables en el VAR.
\\item Utilizar los pronósticos de un período por delante para obtener los pronósticos de dos períodos adelante.
\\item Continuar iterando para obtener pronósticos de todas las variables en el VAR en el futuro.
\\end{enumerate}
\\end{keyconcepts}
')
```

Dado que un VAR modela todas las variables usando rezagos de las otras variables respectivas, se necesitan calcular pronósticos para *todas* las variables. Puede resultar engorroso hacerlo cuando el VAR es grande, pero afortunadamente existen funciones **R** que facilitan esto. Por ejemplo, la función **predict()** se puede utilizar para obtener pronósticos multivariados iterados para modelos VAR estimados por la función **VAR()**.

El siguiente fragmento de código muestra cómo calcular pronósticos iterados para el crecimiento del PIB y el diferencial por plazo hasta el período 2015:Q1, que es $h = 10$, utilizando el objeto modelo **VAR_est**.

```{r, 895}
# calcular pronósticos iterados para el crecimiento del PIB y el diferencial por plazo para los próximos 10 trimestres
forecasts <- predict(VAR_est)
forecasts
```

Esto revela que el pronóstico de dos trimestres de crecimiento del PIB en 2013:Q2 utilizando datos hasta 2012:Q4 es $1.69$. Para el mismo período, el pronóstico de VAR iterado para el diferencial de plazo es $1.88$.

Las matrices devueltas por `predict(VAR_est)` también incluyen intervalos de predicción de $95\%$ (sin embargo, la función no se ajusta para la autocorrelación o heterocedasticidad de los errores).

También se pueden trazar los pronósticos iterados para ambas variables llamando a **plot()** en la salida de `predict(VAR_est)`.

```{r, 896, fig.align='center', fig.height=7}
# visualizar los pronósticos iterados
plot(forecasts)
```

#### Pronósticos directos para múltiples períodos {-}

Un pronóstico directo de múltiples períodos utiliza un modelo en el que los predictores se retrasan de manera adecuada, de modo que las observaciones disponibles se puedan usar *directamente* para realizar el pronóstico. La idea de la predicción directa de varios períodos se resume en el Concepto clave 16.3.

```{r, 897, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC16.3">
<h3 class = "right"> Concepto clave <br> 16.3 </h3>          
<h3 class = "left"> Pronósticos directos para múltiples períodos </h3>
<p>

Un *pronóstico directo para períodos múltiples* que pronostica períodos de $h$ en el futuro utilizando un modelo de $Y_t$ y un predictor adicional $X_t$ con retrasos de $p$ se realiza estimando primero:

\\begin{align*}
  Y_t =& \\, \\delta_0 + \\delta_1 Y_{t-h} + \\dots + \\delta_{p} Y_{t-p-h+1} + \\delta_{p+1} X_{t-h} \\\\
  +& \\dots + \\delta_{2p} Y_{t-p-h+1} + u_t,  
\\end{align*}

que luego se usa para calcular el pronóstico de $Y_{T + h}$ basado en observaciones durante el período $T$.

</p>
</div>
')
```

```{r, 898, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Direct Multiperiod Forecasts]{16.3}

Un \\textit{pronóstico directo para períodos múltiples} que pronostica períodos de $h$ en el futuro utilizando un modelo de $Y_t$ y un predictor adicional $X_t$ con retrasos de $p$ se realiza estimando primero:

\\begin{align*}
  Y_t =& \\, \\delta_0 + \\delta_1 Y_{t-h} + \\dots + \\delta_{p} Y_{t-p-h+1} + \\delta_{p+1} X_{t-h} \\\\
  +& \\dots + \\delta_{2p} Y_{t-p-h+1} + u_t,  
\\end{align*}

que luego se usa para calcular el pronóstico de $Y_{T + h}$ basado en observaciones durante el período $T$.

\\end{keyconcepts}
')
```

Por ejemplo, para obtener previsiones de dos trimestres de crecimiento del PIB y el margen de plazo, primero se estiman las ecuaciones:

\begin{align*}
 GDPGR_t =& \, \beta_{10} + \beta_{11} GDPGR_{t-2} + \beta_{12} GDPGR_{t-3} + \gamma_{11} TSpread_{t-2} + \gamma_{12} TSpread_{t-3} + u_{1t}, \\
 TSpread_t =& \, \beta_{20} + \beta_{21} GDPGR_{t-2} + \beta_{22} GDPGR_{t-3} + \gamma_{21} TSpread_{t-2} + \gamma_{22} TSpread_{t-3} + u_{2t}
\end{align*}

y luego se sustituyen los valores de $GDPGR_{2012:Q4}$, $GDPGR_{2012:Q3}$, $TSpread_{2012:Q4}$ y $TSpread_{2012:Q3}$ en ambas ecuaciones. Esto se hace fácilmente de forma manual.

```{r, 899}
# modelos de estimación para pronósticos directos a dos trimestres
VAR_EQ1_direct <- dynlm(GDPGrowth ~ L(GDPGrowth, 2:3) + L(TSpread, 2:3), 
                        start = c(1981, 1), end = c(2012, 4))

VAR_EQ2_direct <- dynlm(TSpread ~ L(GDPGrowth, 2:3) + L(TSpread, 2:3), 
                        start = c(1981, 1), end = c(2012, 4))

# calcular pronósticos directos a dos trimestres
coef(VAR_EQ1_direct) %*% c(1, # intercepto
                           window(GDPGrowth, start = c(2012, 3), end = c(2012, 4)), 
                           window(TSpread, start = c(2012, 3), end = c(2012, 4)))

coef(VAR_EQ2_direct) %*% c(1, # intercepto
                           window(GDPGrowth, start = c(2012, 3), end = c(2012, 4)), 
                           window(TSpread, start = c(2012, 3), end = c(2012, 4)))
```

Los economistas aplicados a menudo usan el método iterado, ya que estos pronósticos son más confiables en términos de $MSFE$, siempre que el modelo de un período adelantado se especifique correctamente. Si este no es el caso, por ejemplo, porque se cree que una ecuación en un VAR está mal especificada, puede ser beneficioso usar pronósticos directos, ya que el método iterado estará sesgado y, por lo tanto, tendrá un $MSFE$ más alto que el método directo. Consulte el Capítulo \@ref(OIPRUDFGLS) para obtener una discusión más detallada sobre las ventajas y desventajas de ambos métodos.

## Órdenes de integración y prueba de raíz unitaria DF-GLS {#OIPRUDFGLS}

Algunas series de tiempo económicas tienen tendencias más suaves que las variables que pueden describirse mediante modelos de recorridos aleatorios. Una forma de modelar estas series de tiempo es $$\Delta Y_t = \beta_0 + \Delta Y_{t-1} + u_t,$$ donde $u_t$ es un término de error no correlacionado en serie. Este modelo establece que la primera diferencia de una serie es un paseo aleatorio. En consecuencia, la serie de segundas diferencias de $Y_t$ es estacionaria. El Concepto clave 16.4 resume la notación.

```{r, 900, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC16.4">
<h3 class = "right"> Concepto clave <br> 16.4 </h3>          
<h3 class = "left"> Órdenes de integración, diferenciación y estacionariedad </h3>
<p>

+ Cuando una serie de tiempo $Y_t$ tiene una raíz autorregresiva unitaria, $Y_t$ se integra de orden uno. Esto a menudo se denota por $Y_t \\sim I(1)$. Simplemente se dice que $Y_t$ es $I(1)$. Si $Y_t$ es $I(1)$, su primera diferencia $\\Delta Y_t$ es estacionaria.

+ $Y_t$ es $I(2)$ cuando $Y_t$ necesita diferenciarse dos veces para obtener una serie estacionaria. Usando la notación presentada aquí, si $Y_t$ es $I(2)$, su primera diferencia $\\Delta Y_t$ es $I(1)$ y su segunda diferencia $\\Delta^2 Y_t$ es estacionaria. $Y_t$ es $I(d)$ cuando $Y_t$ debe diferenciarse $d$ veces para obtener una serie estacionaria.

+ Cuando $Y_t$ es estacionario, se integra del orden $0$ por lo que $Y_t$ es $I(0)$.

Es bastante fácil obtener diferencias de series de tiempo en <tt>R</tt>. Por ejemplo, la función <tt>diff()</tt> devuelve diferencias adecuadamente rezagadas e iteradas de vectores numéricos, matrices y objetos de series de tiempo de la clase <tt>ts</tt>.

</p>
</div>
')
```

```{r, 901, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Órdenes de integración, diferenciación y estacionariedad]{16.4}
\\begin{itemize}
\\item Cuando una serie de tiempo $Y_t$ tiene una raíz autorregresiva unitaria, $Y_t$ se integra de orden uno. Esto a menudo se denota por $Y_t \\sim I(1)$. Simplemente se dice que $Y_t$ es $I(1)$. Si $Y_t$ es $I(1)$, su primera diferencia $\\Delta Y_t$ es estacionaria.
\\item $Y_t$ es $I(2)$ cuando $Y_t$ necesita diferenciarse dos veces para obtener una serie estacionaria. Usando la notación presentada aquí, si $Y_t$ es $I(2)$, su primera diferencia $\\Delta Y_t$ es $I(1)$ y su segunda diferencia $\\Delta^2 Y_t$ es estacionaria. $Y_t$ es $I(d)$ cuando $Y_t$ debe diferenciarse $d$ veces para obtener una serie estacionaria.
\\item Cuando $Y_t$ es estacionario, se integra del orden $0$ por lo que $Y_t$ es $I(0)$.
\\end{itemize}\\vspace{0.5cm}

Es bastante fácil obtener diferencias de series de tiempo en \\texttt{R}. Por ejemplo, la función \\texttt{diff()} devuelve diferencias adecuadamente rezagadas e iteradas de vectores numéricos, matrices y objetos de series de tiempo de la clase \\texttt{ts}.

\\end{keyconcepts}
')
```

Se toma el nivel de precios de los EE. UU. Medido por el *Índice de precios de gastos de consumo personal* como ejemplo.

```{r, 902, fig.align='center'}
# definir el objeto ts del índice de precios PCE de EE. UU.
PCECTPI <- ts(log(USMacroSWQ$PCECTPI), 
              start = c(1957, 1), 
              end = c(2012, 4), 
              freq = 4)

# graficar logaritmo del índice de precios de PCE
plot(log(PCECTPI),
     main = "Logaritmo del índice de precios PCE de Estados Unidos",
     ylab = "Logaritmo",
     col = "steelblue", 
     lwd = 2)
```

El logaritmo del nivel de precios tiene una tendencia que varía suavemente. Esto es típico de una serie $I(2)$. Si el nivel de precios es realmente $I(2)$, las primeras diferencias de esta serie deberían ser $I(1)$. Dado que se está considerando el logaritmo del nivel de precios, se obtienen tasas de crecimiento tomando las primeras diferencias. Por lo tanto, la serie de niveles de precios diferenciados es la serie de tasas de inflación trimestrales. Esto se hace rápidamente en **R** usando la función **Delt()** del paquete **quantmod**. Como se explica en el Capítulo \@ref(DSTCS), al multiplicar las tasas de inflación trimestrales por $400$ se obtiene la tasa de inflación trimestral, medida en puntos porcentuales a una tasa anual.

```{r, 903, fig.align='center'}
# graficar la inflación de precios del PCE de EE. UU.
plot(400 * Delt(PCECTPI),
     main = "United States PCE Price Index",
     ylab = "Porcentaje por año",
     col = "steelblue", 
     lwd = 2)

# agregar una línea discontinua en y = 0
abline(0, 0, lty = 2)
```

La tasa de inflación se comporta de manera mucho más errática que el gráfico uniforme del logaritmo del índice de precios del PCE.

#### La prueba DF-GLS para una raíz unitaria {-}

La prueba DF-GLS para una raíz unitaria ha sido desarrollada por @elliott1996 y tiene mayor potencia que la prueba ADF cuando la raíz autorregresiva es grande, pero menor que uno; es decir, el DF-GLS tiene una mayor probabilidad de rechazar el falso nulo de una tendencia estocástica cuando los datos de la muestra provienen de una serie de tiempo que está cerca de integrarse.

La idea de la prueba DF-GLS es probar una raíz unitaria autorregresiva en la serie sin tendencia, mediante la cual las estimaciones de GLS de los componentes deterministas se utilizan para obtener la versión sin tendencia de la serie original.

Una función que realiza la prueba DF-GLS se implementa en el paquete **urca** (este paquete es una dependencia del paquete **vars**, por lo que ya debería estar cargado si se adjunta **vars**). La función que calcula la estadística de prueba es **ur.ers**.

```{r, 904}
# prueba DF-GLS para raíz unitaria en el PIB
summary(ur.ers(log(window(GDP, start = c(1962, 1), end = c(2012, 4))),
        model = "trend", 
        lag.max = 2))
```

El resumen de la prueba muestra que el estadístico de la prueba es de aproximadamente $-1.2$. El valor crítico de $10\% $ para la prueba DF-GLS es $-2.57$. Sin embargo, este **no** es el valor crítico apropiado para la prueba ADF cuando se incluyen una intersección y una tendencia temporal en la regresión de Dickey-Fuller: ¡Las distribuciones asintóticas de ambas estadísticas de prueba difieren y también sus valores críticos!

La prueba es del lado izquierdo, por lo que no se puede rechazar la hipótesis nula de que la inflación de EE. UU. no es estacionaria, utilizando la prueba DF-GLS.

## Cointegración

```{r, 905, eval = my_output == "html", results='asis', echo=F, purl=F}
cat('
<div class = "keyconcept" id="KC16.5">
<h3 class = "right"> Concepto clave <br> 16.5 </h3>          
<h3 class = "left"> Cointegration </h3>
<p>

Cuando $X_t$ y $Y_t$ son $I(1)$ y si existe un $\\theta$ tal que $Y_t - \\theta X_t$ es $I(0)$, $X_t$ y $Y_t$ están cointegrados. Dicho de otra manera, la cointegración de $X_t$ y $Y_t$ significa que $X_t$ y $Y_t$ tienen la misma tendencia estocástica o una común y que esta tendencia puede eliminarse tomando una diferencia específica de la serie de modo que la serie resultante sea estacionario.

Las funciones <tt>R</tt> para el análisis de cointegración se implementan en el paquete <tt>urca</tt>.

</p>
</div>
')
```

```{r, 906, eval = my_output == "latex", results='asis', echo=F, purl=F}
cat('\\begin{keyconcepts}[Cointegration]{16.5}

Cuando $X_t$ y $Y_t$ son $I(1)$ y si existe un $\\theta$ tal que $Y_t - \\theta X_t$ es $I(0)$, $X_t$ y $Y_t$ están cointegrados. Dicho de otra manera, la cointegración de $X_t$ y $Y_t$ significa que $X_t$ y $Y_t$ tienen la misma tendencia estocástica o una común y que esta tendencia puede eliminarse tomando una diferencia específica de la serie de modo que la serie resultante sea estacionario.\\newline

Las funciones \\texttt{R}  para el análisis de cointegración se implementan en el paquete \\texttt{urca}.

\\end{keyconcepts}
')
```

Como ejemplo, se reconsidera la relación entre las tasas de interés a corto y largo plazo con el ejemplo de las letras del Tesoro de EE. UU. a 3 meses, los bonos del Tesoro de EE. UU. a 10 años y el diferencial en sus tasas de interés que se han introducido en el Capítulo \@ref(PGMPI). El siguiente fragmento de código muestra cómo crear el gráfico.

```{r, 907, fig.align='center'}
# graficar ambas series de interés
plot(merge(as.zoo(TB3MS), as.zoo(TB10YS)), 
     plot.type = "single", 
     lty = c(2, 1),
     lwd = 2,
     xlab = "Fecha",
     ylab = "Porcentaje por año",
     ylim = c(-5, 17),
     main = "Tasas de interés")

# agregar la serie de diferencial por plazo
lines(as.zoo(TSpread), 
     col = "steelblue",
     lwd = 2,
     xlab = "Fecha",
     ylab = "Porcentaje por año",
     main = "Diferencial por plazo")

# sombrear el diferencial por plazo
polygon(c(time(TB3MS), rev(time(TB3MS))), 
        c(TB10YS, rev(TB3MS)),
        col = alpha("steelblue", alpha = 0.3),
        border = NA)

# agregar línea horizontal agregar 0
abline(0, 0)

# agregar una leyenda
legend("topright", 
       legend = c("TB3MS", "TB10YS", "Diferencial por plazo"),
       col = c("black", "black", "steelblue"),
       lwd = c(2, 2, 2),
       lty = c(2, 1, 1))
```

El gráfico sugiere que las tasas de interés a largo y a corto plazo están cointegradas: Ambas series de intereses parecen tener el mismo comportamiento a largo plazo. Comparten una tendencia estocástica común. El diferencial de plazo, que se obtiene tomando la diferencia entre las tasas de interés de largo y corto plazo, parece estacionario. De hecho, la teoría de expectativas de la estructura de términos sugiere que el coeficiente de cointegración $\theta$ es 1. Esto es consistente con el resultado visual.

#### Pruebas de cointegración {-}

Siguiendo el Concepto clave 16.5, parece natural construir una prueba para la cointegración de dos series de la siguiente manera: Si dos series $X_t$ y $Y_t$ están cointegradas, la serie obtenida tomando la diferencia $Y_t - \theta X_t$ debe ser estacionario. Si las series no están cointegradas, $Y_t - \theta X_t$ no es estacionaria. Esta es una suposición que se puede probar mediante una prueba de raíz unitaria. Se tiene que distinguir entre dos casos:

1. **$\theta$ es conocido.**

    El conocimiento de $\theta$ permite calcular las diferencias $z_t = Y_t - \theta X_t$ para que las pruebas de raíz unitaria de Dickey-Fuller y DF-GLS se puedan aplicar a $z_t$. Para estas pruebas, los valores críticos son los valores críticos de la prueba ADF o DF-GLS.

2. **$\theta$ es desconocido.**

    Si se desconoce $\theta$, debe estimarse antes de que se pueda aplicar la prueba de raíz unitaria. Esto se hace estimando la regresión $$Y_t = \alpha + \theta X_t + z_t$$ usando MCO (esto se conoce como la regresión de la primera etapa). Luego, se usa una prueba de Dickey-Fuller para probar la hipótesis de que $z_t$ es una serie no estacionaria. Esto se conoce como prueba Engle-Granger Augmented Dickey-Fuller para cointegración (o *prueba EG-ADF*) después de @engle1987. Los valores críticos para esta prueba son especiales, ya que la distribución nula asociada no es normal y depende del número de variables $I(1)$ utilizadas como regresores en la regresión de la primera etapa. Cuando solo existen dos variables presuntamente cointegradas (y, por lo tanto, se usa una sola variable $I(1)$ en la regresión de MCO de la primera etapa), los valores críticos para los niveles $10\% $, $5\%$ y $1\%$ son $-3.12$, $-3.41$ y $-3.96$.
    
#### Aplicación a las tasas de interés {-}

Como se mencionó anteriormente, la teoría de la estructura temporal sugiere que las tasas de interés a largo y corto plazo están cointegradas con un coeficiente de cointegración de $\theta = 1$. En la sección anterior se ha visto que existe evidencia visual de esta conjetura, ya que el diferencial de las tasas de interés a 10 años y a 3 meses parece estacionario.

Se continua usando pruebas formales (la ADF y la prueba DF-GLS) para ver si las series de tasas de interés individuales están integradas y si su diferencia es estacionaria (por ahora, se asume que se conoce $\theta = 1$). Ambos se hacen convenientemente usando las funciones **ur.df()** para el cálculo de la prueba ADF y **ur.ers()** para realizar la prueba DF-GLS. Se usan datos desde 1962:Q1 hasta 2012:Q4 y se emplean modelos que incluyen un término de deriva. Se establece el orden de retraso máximo en $6$ y se usa $AIC$ para seleccionar la longitud de retraso óptima.

```{r, 908}
# prueba de no estacionariedad de letras del tesoro a 3 meses usando la prueba ADF
ur.df(window(TB3MS, c(1962, 1), c(2012, 4)), 
      lags = 6, 
      selectlags = "AIC", 
      type = "drift")

# prueba de no estacionariedad de bonos del tesoro a 10 años usando la prueba ADF
ur.df(window(TB10YS, c(1962, 1), c(2012, 4)), 
      lags = 6, 
      selectlags = "AIC", 
      type = "drift")

# prueba de no estacionariedad de letras del tesoro a 3 meses usando la prueba DF-GLS
ur.ers(window(TB3MS, c(1962, 1), c(2012, 4)),
       model = "constant", 
       lag.max = 6)

# prueba de no estacionariedad de bonos del tesoro a 10 años utilizando la prueba DF-GLS
ur.ers(window(TB10YS, c(1962, 1), c(2012, 4)),
       model = "constant", 
       lag.max = 6)
```

El valor crítico correspondiente de $10\%$ para ambas pruebas es $-2.57$, por lo que no se puede rechazar las hipótesis nulas de no estacionario para ninguna de las series, incluso en el nivel de significancia de $10\%$.^[Nota: **ur.df()** informa dos estadísticas de prueba cuando existe una desviación en la regresión ADF. El primero de los cuales (el que interesa aquí) es el estadístico $t$ para la prueba de que el coeficiente en el primer rezago de la serie es 0. El segundo es el estadístico $t$ para la prueba de hipótesis que el término de deriva es igual a $0$.] Se concluye que es plausible modelar ambas series de tasas de interés como $I(1)$.

A continuación, se aplica la prueba ADF y DF-GLS para probar la no estacionariedad de la serie de márgenes de plazo, lo que significa que se prueba la no cointegración de las tasas de interés a corto y largo plazo.

```{r, 909}
# probar si el diferencial de plazo es estacionario (cointegración de tasas de interés) usando ADF
ur.df(window(TB10YS, c(1962, 1), c(2012, 4)) - window(TB3MS, c(1962, 1), c(2012 ,4)), 
      lags = 6, 
      selectlags = "AIC", 
      type = "drift")

# probar si el diferencial de plazo es estacionario (cointegración de las tasas de interés) utilizando la prueba DF-GLS
ur.ers(window(TB10YS, c(1962, 1), c(2012, 4)) - window(TB3MS, c(1962, 1), c(2012, 4)),
       model = "constant", 
       lag.max = 6)
```

La Tabla \@ref(tab:spreadcoint) resume los resultados.

| Series         | Estadístic0 de prueba ADF | Estadístico de prueba DF-GLS |
|----------------|:-------------------------:|:----------------------------:|
| TB3MS          |          $-2.10$          |            $-1.80$           |
| TB10YS         |          $-1.01$          |            $-0.94$           |
| TB10YS - TB3MS |          $-3.93$          |            $-3.86$           |

Table: (\#tab:spreadcoint) Estadísticos de prueba ADF y DF-GLS para series de tipos de interés

Ambas pruebas rechazan la hipótesis de no estacionariedad de la serie de diferenciales por plazo en el nivel de significancia de $1\%$, lo que constituye una fuerte evidencia a favor de la hipótesis de que el diferencial de plazos es estacionario, lo que implica la cointegración de las tasas de interés de largo y corto plazo.

Dado que la teoría sugiere que $\theta=1$, no existe necesidad de estimar $\theta$, por lo que no es necesario utilizar la prueba EG-ADF que permite que $\theta$ sea desconocido. Sin embargo, dado que es instructivo hacerlo, se calcula este estadístico de prueba. La regresión MCO de la primera etapa es $$TB10YS_t = \beta_0 + \beta_1 TB3MS_t + z_t.$$

```{r, 910}
# estimar la regresión de la primera etapa de la prueba EG-ADF
FS_EGADF <- dynlm(window(TB10YS, c(1962, 1), c(2012, 4)) ~ window(TB3MS, c(1962, 1), c(2012, 4)))
FS_EGADF
```

Así se tiene que: 

\begin{align*}
  \widehat{TB10YS}_t = 2.46 + 0.81 \cdot TB3MS_t,
\end{align*}

donde $\widehat{\theta} = 0.81$. A continuación, se toma la serie residual $\{\widehat{z_t}\}$ y se calcula el estadístico de prueba ADF.

```{r, 911}
# calcular los residuos
z_hat <- resid(FS_EGADF)

# calcular el estadístico de prueba ADF
ur.df(z_hat, lags = 6, type = "none", selectlags = "AIC")
```

El estadístico de prueba es $-3.19$, que es menor que el valor crítico de $10\%$, pero mayor que el valor crítico de $5\%$. Por lo tanto, la hipótesis nula de no cointegración puede rechazarse en el nivel de $10\%$, pero no en el nivel de $5\%$. Esto indica una potencia más baja de la prueba EG-ADF debido a la estimación de $\theta$: Cuando $\theta = 1$ es el valor correcto, se espera la potencia de la prueba ADF para una raíz unitaria en la serie de residuos \widehat{z} = TB10YS - TB3MS$ mayor que cuando se usa una estimación de $\widehat{\theta}$.

##### Un modelo de corrección de errores vectoriales para $TB10YS_t$ y $TB3MS$ {-}

Si dos $I(1)$ series de tiempo $X_t$ y $Y_t$ están cointegradas, sus diferencias son estacionarias y se pueden modelar en un VAR que se aumenta con el regresor $Y_{t-1} - \theta X_{t-1}$. Esto se denomina *modelo de corrección de errores vectoriales* (VECM) y $Y_{t} - \theta X_{t}$ se denomina *término de corrección de errores*. Los valores rezagados del término de corrección de errores son útiles para predecir $\Delta X_t$ y/o $\Delta Y_t$.

Se puede utilizar un VECM para modelar las dos tasas de interés consideradas en las secciones anteriores. Se especifica el VECM para incluir dos rezagos de ambas series como regresores y elegir $\theta = 1$, como sugiere la teoría (ver arriba).

```{r, 912}
TB10YS <- window(TB10YS, c(1962, 1), c(2012 ,4))
TB3MS <- window(TB3MS, c(1962, 1), c(2012, 4))

# configurar el término de corrección de errores
VECM_ECT <- TB10YS - TB3MS

# estimar ambas ecuaciones del VECM usando 'dynlm()'
VECM_EQ1 <- dynlm(d(TB10YS) ~ L(d(TB3MS), 1:2) + L(d(TB10YS), 1:2) + L(VECM_ECT))
VECM_EQ2 <- dynlm(d(TB3MS) ~ L(d(TB3MS), 1:2) + L(d(TB10YS), 1:2) + L(VECM_ECT))

# cambiar el nombre de los regresores para una mejor legibilidad
names(VECM_EQ1$coefficients) <- c("Intercept", "D_TB3MS_l1", "D_TB3MS_l2",
                                  "D_TB10YS_l1", "D_TB10YS_l2", "ect_l1")
names(VECM_EQ2$coefficients) <- names(VECM_EQ1$coefficients)

# resúmenes de coeficientes utilizando errores estándar de HAC
coeftest(VECM_EQ1, vcov. = NeweyWest(VECM_EQ1, prewhite = F, adjust = T))
coeftest(VECM_EQ2, vcov. = NeweyWest(VECM_EQ2, prewhite = F, adjust = T))
```

Por tanto, las dos ecuaciones estimadas del VECM son

\begin{align*}
 \widehat{\Delta TB3MS}_t =& \, -\underset{(0.11)}{0.06} + \underset{(0.11)}{0.24} \Delta TB3MS_{t-1} -\underset{(0.15)}{0.16} \Delta TB3MS_{t-2} \\ &+ \underset{(0.13)}{0.11} \Delta TB10YS_{t-1} -\underset{(0.11)}{0.15} \Delta TB10YS_{t-2} + \underset{(0.05)}{0.03} ECT_{t-1} \\
 \widehat{\Delta TB10YS}_t =& \, \underset{(0.06)}{0.12} -\underset{(0.07)}{0.00} \Delta TB3MS_{t-1} -\underset{(0.04)}{0.07} \Delta TB3MS_{t-2} \\ &+ \underset{(0.10)}{0.23} \Delta TB10YS_{t-1} -\underset{(0.07)}{0.07} \Delta TB10YS_{t-2} -\underset{(0.03)}{0.09} ECT_{t-1}.
\end{align*}

El resultado producido por **coeftest()** muestra que existe poca evidencia de que los valores rezagados de la serie de intereses diferenciados sean útiles para la predicción. Este hallazgo es más pronunciado para la ecuación de la serie diferenciada de la tasa de la letra del tesoro a 3 meses, donde el término de corrección de errores (el diferencial de plazo rezagado) no es significativamente diferente de cero en ningún nivel común de significancia. Sin embargo, para la tasa diferenciada de los bonos del tesoro a 10 años, el término de corrección de errores es estadísticamente significativo a $1\%$ con una estimación de $-0.09$. Esto se puede interpretar de la siguiente manera: Aunque ambas tasas de interés no son estacionarias, su relación de conintegración permite predecir el *cambio* en la tasa de los bonos del tesoro a 10 años utilizando el VECM. En particular, la estimación negativa del coeficiente del término de corrección de errores indica que habrá un cambio negativo en la tasa de los bonos del tesoro a 10 años del próximo período cuando la tasa de los bonos del tesoro a 10 años sea inusualmente alta en relación con la tasa del tesoro a 3 meses en el período actual.

## Agrupación de volatilidad y heterocedasticidad condicional autorregresiva

Las series de tiempo financieras suelen presentar un comportamiento que se conoce como *agrupamiento de volatilidad*: La volatilidad cambia con el tiempo y su grado muestra una tendencia a persistir; es decir, existen periodos de baja volatilidad y periodos donde la volatilidad es alta. Los econometristas llaman a esto *heterocedasticidad condicional autorregresiva*. La heterocedasticidad condicional es una propiedad interesante porque puede explotarse para pronosticar la varianza de períodos futuros.

Como ejemplo, se consideran los cambios diarios en el índice bursátil Whilshire 5000. Los datos están disponibles para su descarga en la [Base de datos económicos de la Reserva Federal](https://fred.stlouisfed.org/series/WILL5000INDFC). Para mantener la coherencia con el libro, se descargaron los datos desde el 29 de diciembre de 1989 hasta el 31 de diciembre de 2013 (es necesario elegir este intervalo de tiempo algo mayor, ya que más adelante se trabajará con los cambios diarios de la serie).

El siguiente fragmento de código muestra cómo formatear los datos y crear un gráfico a partir de ellos.

```{r, 913}
# importar datos del índice Wilshire 5000
W5000 <- read.csv2("data/Wilshire5000.csv", 
                   stringsAsFactors = F, 
                   header = T, 
                   sep = ",", 
                   na.strings = ".")

# transformar las columnas
W5000$DATE <- as.Date(W5000$DATE)
W5000$WILL5000INDFC <- as.numeric(W5000$WILL5000INDFC)

# eliminar NAs
W5000 <- na.omit(W5000)

# calcular cambios porcentuales diarios
W5000_PC <- data.frame("Date" = W5000$DATE, 
                       "Value" = as.numeric(Delt(W5000$WILL5000INDFC) * 100))
W5000_PC <- na.omit(W5000_PC)
```

```{r pcw5000, 914, fig.align='center', fig.cap="Rendimientos porcentuales diarios en el índice Wilshire 5000"}
# graficar cambios porcentuales
plot(W5000_PC, 
     ylab = "Porcentaje", 
     main = "Cambios porcentuales diarios",
     type="l", 
     col = "steelblue", 
     lwd = 0.5)

# agregar una línea horizontal en y = 0
abline(0, 0)
```

La serie de cambios porcentuales diarios en el índice de Wilshire parece fluctuar aleatoriamente alrededor de cero, lo que significa que existe poca autocorrelación. Esto se confirma mediante un gráfico de la función de autocorrelación de la muestra.

```{r acfw5000, 915, fig.align='center', fig.cap="Autocorrelación en los cambios diarios de precios del índice W5000"}
# graficar la autocorrelación de la muestra de los cambios porcentuales diarios
acf(W5000_PC$Value, main = "Serie 5000 de Wilshire")
```

En la Figura \@ref(fig:acfw5000) se ve que las autocorrelaciones son bastante débiles, por lo que es difícil predecir los resultados futuros utilizando, por ejemplo, un modelo AR.

Sin embargo, existe evidencia visual en \@ref(fig:pcw5000) de que la serie de retornos exhibe heterocedasticidad condicional, ya que se observan agrupaciones de volatilidad. Para algunas aplicaciones, es útil medir y pronosticar estos patrones. Esto se puede hacer utilizando modelos que asumen que la volatilidad puede describirse mediante un proceso autorregresivo.

### Modelos ARCH y GARCH {-}

Considere $$Y_t = \beta_0 + \beta_1 Y_{t-1} + \gamma_1 X_{t-1} + u_t,$$ un modelo de regresión ADL($1$,$1$). El econométrico Robert @engle1982 propuso modelar $\sigma^2_t = Var(u_t | u_{t-1},u_{t-2},\ldots)$, la varianza condicional del error $u_t$ dado su pasado, por un modelo de retraso distribuido de orden $p$,

\begin{align}
 \sigma^2_t = \alpha_0 + \alpha_1 u_{t-1}^2 + \alpha_2 u_{t-2}^2 + \dots + \alpha_p u_{t-p}^2, (\#eq:archmodel)
\end{align}

llamado modelo *de heterocedasticidad condicional autorregresiva* (ARCH) de orden $p$, o ARCH corto ($p$).^[Aunque se presenta el modelo ARCH como un componente en un modelo ADL($1$,$1$), se puede usar para modelar el término de error condicional de media cero de cualquier modelo de serie temporal.] Se supone $\alpha_0>0$ y $\alpha_1,\ldots,\alpha_p\geq0$ para asegurar una varianza positiva $\sigma_t^2>0$. La idea general se desprende de la estructura del modelo: Los coeficientes positivos $\alpha_0,\alpha_1,\dots,\alpha_p$ implican que los grandes errores cuadrados recientes conducen a una gran varianza y, por lo tanto, a grandes errores cuadrados en el período actual.

El modelo ARCH generalizado (GARCH), desarrollado por Tim @bollerslev1986, es una extensión del modelo ARCH, donde se permite que $\sigma^2_t$ dependa de sus propios retrasos y de los retrasos del término de error al cuadrado. El modelo GARCH ($p$,$q$) viene dado por

\begin{align}
 \sigma^2_t = \alpha_0 + \alpha_1 u_{t-1}^2 + \alpha_2 u_{t-2}^2 + \dots + \alpha_p u_{t-p}^2 + \phi_1 \sigma^2_{t-1} + \dots + \phi_p \sigma^2_{t-q}. (\#eq:garchmodel)
\end{align}

El modelo GARCH es un modelo ADL ($p$,$q$) y, por lo tanto, puede proporcionar parametrizaciones más parsimoniosas que el modelo ARCH.

### Aplicación a la volatilidad del precio de las acciones {-}

Las estimaciones de máxima verosimilitud de los modelos ARCH y GARCH son eficientes y tienen distribuciones normales en muestras grandes, de modo que se pueden aplicar los métodos habituales para realizar inferencias sobre los parámetros desconocidos. El paquete **fGarch** en **R** consiste en una colección de funciones para analizar y modelar el comportamiento heterocedástico en modelos de series de tiempo. La función **garchFit()** es algo sofisticada, dado que permite diferentes especificaciones del procedimiento de optimización, diferentes distribuciones de errores y mucho más (use `?GarchFit` para una descripción detallada de los argumentos). En particular, los errores estándar reportados por **garchFit()** son robustos.

El modelo GARCH ($1$,$1$) de cambios diarios en el índice Wilshire 5000 que se estimó viene dado por:

\begin{align}
  R_t =& \, \beta_0 + u_t \ , \ u_t \sim \mathcal{N}(0,\sigma^2_t), \\
  \sigma^2_t =& \, \alpha_0 + \alpha_1 u_{t-1}^2 + \phi_1 \sigma_{t-1}^2 (\#eq:w5000g11)
\end{align}

donde $R_t$ es el cambio porcentual en el período $t$. $\beta_0$, $\alpha_0$, $\alpha_1$ y $\phi_1$ son coeficientes desconocidos y $u_t$ es un término de error con media condicional cero. No se incluyen predictores rezagados en la ecuación de $R_t$ porque los cambios diarios en el índice Wilshire 5000 reflejan rendimientos bursátiles diarios que son esencialmente impredecibles. Se debe tener en cuenta que se supone que $u_t$ tiene una distribución normal y la varianza $\sigma^2_t$ depende de $t$, ya que sigue la recursividad GARCH ($1$,$1$) \@ref(eq:w5000g11).

Es sencillo estimar este modelo usando **garchFit()**.

```{r, 916}
# estimar el modelo GARCH (1,1) de cambios porcentuales diarios
GARCH_Wilshire <- garchFit(data = W5000_PC$Value, trace = F)
```

Se obtiene:

\begin{align}
  \widehat{R}_t =& \, \underset{(0.010)}{0.068}, (\#eq:w5000g11est1) \\
  \widehat{\sigma}^2_t =& \, \underset{(0.002)}{0.011} + \underset{(0.007)}{0.081} u_{t-1}^2 + \underset{(0.008)}{0.909} \sigma_{t-1}^2, (\#eq:w5000g11est2)
\end{align}

por lo que los coeficientes en $u_{t-1}^2$ y $\sigma^2_{t-1}$ son estadísticamente significativos en cualquier nivel común de significancia. Se puede demostrar que la persistencia de movimientos en $\sigma^2_t$ está determinada por la suma de ambos coeficientes, que aquí es $0.99$. Esto indica que los movimientos en la varianza condicional son muy persistentes, lo que implica períodos prolongados de alta volatilidad, lo que es consistente con la evidencia visual de la agrupación de volatilidad presentada anteriormente.

La varianza condicional estimada $\widehat{\sigma}^2_t$ se puede calcular conectando los residuos de \@ref(eq:w5000g11est1) en la ecuación \@ref(eq:w5000g11est2). Esto lo realiza automáticamente **garchFit()**, por lo que para obtener las desviaciones estándar condicionales estimadas $\widehat{\sigma}_t$ solo se tienen que leer los valores de **GARCH_Wilshire** agregando $\textit{@sigma.t}$.

Usando $\widehat{\sigma}_t$ se grafican bandas de $\pm$ una desviación estándar condicional junto con las desviaciones de la serie de cambios porcentuales diarios en el índice Wilshire 5000 de su media. El siguiente fragmento de código genera el gráfico.

```{r, 917, fig.align='center'}
# calcular las desviaciones de los cambios porcentuales de su media
dev_mean_W5000_PC <- W5000_PC$Value - GARCH_Wilshire@fit$coef[1]

# graficar la desviación de los cambios porcentuales de la media
plot(W5000_PC$Date, dev_mean_W5000_PC, 
     type = "l", 
     col = "steelblue",
     ylab = "Porcentaje", 
     xlab = "Fecha",
     main = "Bandas estimadas de + - una desviación estándar condicional",
     lwd = 0.2)

# agregar una línea horizontal en y = 0
abline(0, 0)

# agregar bandas de confianza GARCH(1,1) (una desviación estándar) al gráfico
lines(W5000_PC$Date, 
      GARCH_Wilshire@fit$coef[1] + GARCH_Wilshire@sigma.t, 
      col = "darkred", 
      lwd = 0.5)

lines(W5000_PC$Date, 
      GARCH_Wilshire@fit$coef[1] - GARCH_Wilshire@sigma.t, 
      col = "darkred", 
      lwd = 0.5)
```

Las bandas de las desviaciones estándar condicionales estimadas siguen bastante bien la heterocedasticidad observada en la serie de cambios diarios del índice Wilshire 5000. Esto es útil para cuantificar la volatilidad variable en el tiempo y el riesgo resultante para los inversores que poseen acciones resumidas por el índice. Además, este modelo GARCH también se puede utilizar para producir intervalos de pronóstico cuyos anchos dependen de la volatilidad de los períodos más recientes.

### Resumen {-}

+ Se ha discutido cómo las autorregresiones vectoriales se estiman convenientemente y se usan para pronosticar en **R** mediante funciones del paquete **vars**.

+ El paquete **urca** proporciona métodos avanzados para el análisis de la raíz unitaria y cointegración como las pruebas DF-GLS y EG-ADF. En una aplicación, se ha encontrado evidencia de que las tasas de interés a 3 meses y 10 años tienen una tendencia estocástica común (es decir, están cointegradas) y, por lo tanto, pueden modelarse utilizando un modelo de corrección de errores vectoriales.

+ Además, se ha introducido el concepto de agrupamiento de volatilidad y demostrado cómo se puede emplear la función **garchFit()** del paquete **fGarch** para estimar un modelo GARCH($1$,$1$) de la heterocedasticidad condicional inherente a los rendimientos del índice bursátil Wilshire 5000.