Capítulo 6 Modelos Univariados y Multivariados de Volatilidad

6.1 Motivación

La exposición de los modelos univariados de este capítulo –los hechos estilizados que los motivan, la especificación ARCH y GARCH y sus extensiones asimétricas– sigue a Brooks (2019, cap. 9); la parte multivariada sigue a Lütkepohl (2005, cap. 16).

En este capítulo discutiremos los modelos de Heterocedasticidad Condicional Autorregresiva (ARCH, por sus siglas en inglés) y los modelos de Heterocedasticidad Condicional Autorregresiva Generalizados (GARCH, por sus siglas en inglés), los cuales tienen la característica de modelar situaciones como las que ilustra la Figura 6.2. Es decir:

Existen zonas donde la variación de los datos es mayor y zonas donde la variación es más estable–a estas situaciones se les conoce como de variabilidad por clúster–, y
los datos corresponden a información de alta frecuencia.

Iniciemos por las bibliotecas necesarias:

library(expm)
library(Matrix)
library(ggplot2)
library(quantmod)
library(moments)
library(dynlm)
library(broom)
library(FinTS)
library(lubridate)
library(forecast)
library(readxl)
library(MASS)
library(rugarch)
library(tsbox)
library(MTS)
library(rmgarch)
library(Rcpp)

Para el análisis de temas financieros existe un paquete de mucha utilidad llamado quantmod. En primer lugar, este paquete permite acceder a datos financieros de un modo muy simple: es posible descargar series financieras desde Yahoo Finance, la FRED (Federal Reserve Economic Data), entre otras fuentes. Por otro lado, también permite realizar gráficos de uso común en finanzas con unas cuantas líneas de código. Usaremos datos de la cotización del bitcoin (ticker: “BTC-USD”) descargados de Yahoo Finance, con corte al 7 de julio de 2026; la serie está almacenada en la carpeta BD/ del repositorio para garantizar la reproducibilidad de los resultados que se discuten en el texto. Con el mismo procedimiento y la misma fecha de corte se congeló la serie del Ethereum (BD/ETH_USD.RData), que utilizaremos en el ejemplo multivariado del final del capítulo.

## Datos congelados (corte: 2026-07-07) para reproducibilidad.
## Quien desee replicar el ejercicio con datos vivos puede descomentar
## las siguientes líneas para descargar la serie actualizada desde
## Yahoo Finance (los resultados numéricos cambiarán con cada descarga):
#
# options("getSymbols.warning4.0" = FALSE)
# BTC <- getSymbols("BTC-USD", src = "yahoo", auto.assign = FALSE)
# BTC <- na.omit(BTC)
# save(BTC, file = "BD/BTC_USD.RData")

load("BD/BTC_USD.RData")

chartSeries(BTC,TA='addBBands();
                    addBBands(draw="p");
                    addVo();
                    addMACD()',# subset='2021',
                theme="white")

head(BTC)

##            BTC-USD.Open BTC-USD.High BTC-USD.Low BTC-USD.Close
## 2014-09-17      465.864      468.174     452.422       457.334
## 2014-09-18      456.860      456.860     413.104       424.440
## 2014-09-19      424.103      427.835     384.532       394.796
## 2014-09-20      394.673      423.296     389.883       408.904
## 2014-09-21      408.085      412.426     393.181       398.821
## 2014-09-22      399.100      406.916     397.130       402.152
##            BTC-USD.Volume BTC-USD.Adjusted
## 2014-09-17       21056800          457.334
## 2014-09-18       34483200          424.440
## 2014-09-19       37919700          394.796
## 2014-09-20       36863600          408.904
## 2014-09-21       26580100          398.821
## 2014-09-22       24127600          402.152

tail(BTC)

##            BTC-USD.Open BTC-USD.High BTC-USD.Low BTC-USD.Close
## 2026-07-02     60004.77     62117.88    59532.32      61485.30
## 2026-07-03     61492.07     62879.01    61176.72      62544.20
## 2026-07-04     62545.13     63398.41    62287.59      63088.30
## 2026-07-05     63089.06     63935.85    62413.99      63547.88
## 2026-07-06     63551.02     64597.57    61275.83      63995.02
## 2026-07-07     63994.60     64257.63    62623.93      63297.39
##            BTC-USD.Volume BTC-USD.Adjusted
## 2026-07-02    40109297349         61485.30
## 2026-07-03    26131813598         62544.20
## 2026-07-04    18608397613         63088.30
## 2026-07-05    18267466154         63547.88
## 2026-07-06    36552222303         63995.02
## 2026-07-07    31026493556         63297.39

Para fines del ejercicio de este capítulo, usaremos el precio ajustado del activo. Esto nos servirá para calcular el rendimiento diario, o puesto en lenguaje de series temporales podemos decir que usaremos la serie en diferencias logarítmicas.

plot(BTC$`BTC-USD.Adjusted`)

Figura 6.1: Evolución del precio del Bitcoin

Una de las preguntas relevantes al observar la serie en diferencias, es si podríamos afirmar que esta serie cumple con el supuesto de homocedasticidad. Para ello, la Figura 6.2 muestra que las variaciones en el precio del Bitcoin muestran un escenario en el que no se cumple dicho supuesto.

logret <- ts(diff(log(BTC$`BTC-USD.Adjusted`))[-1])

plot(logret)

Figura 6.2: Evolución del rendimiento (diferencias logarítmicas) del Bitcoin

6.1.1 Valor en Riesgo (VaR)

Utilicemos como ejemplo el Valor en Riesgo (VaR, Value at Risk). Se trata del cuantil $\alpha$ de la distribución de rendimientos: la pérdida que no se excede con probabilidad $1 - \alpha$. Supongamos un $\alpha = 0.05$; de esta forma, la Figura 6.3 ilustra la región de la distribución que consideraríamos como el VaR.

alpha <- 0.05

# VaR histórico (no condicional): cuantil alpha de los rendimientos
VaR <- quantile(logret, alpha)

# Déficit esperado (Expected Shortfall): pérdida promedio condicional
# a que se rebase el VaR
ES <- mean(logret[logret < VaR])

knitr::kable(
  data.frame(Medida = c("VaR (5\\%)", "Déficit esperado (5\\%)"),
             Valor  = round(100*c(as.numeric(VaR), ES), 4)),
  col.names = c("Medida", "Rendimiento diario (\\%)"),
  caption   = "Valor en Riesgo y déficit esperado históricos al 5\\% para los rendimientos diarios del Bitcoin.",
  align     = c("l", "r"),
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.1: Valor en Riesgo y déficit esperado históricos al 5% para los rendimientos diarios del Bitcoin.
Medida	Rendimiento diario (%)
VaR (5%)	-5.4010
Déficit esperado (5%)	-8.5707

Es decir, en el 5% de los días con peor desempeño el Bitcoin perdió al menos 5.4% de su valor, y la pérdida promedio en esos días –el déficit esperado o expected shortfall– fue de 8.57%. La distinción es relevante en la administración de riesgos: el VaR indica el umbral de pérdida, mientras que el déficit esperado informa sobre la magnitud de las pérdidas cuando ese umbral se rebasa.

df_ret <- data.frame(logret = as.numeric(logret))

ggplot(df_ret, aes(x = logret)) +
  geom_histogram(fill = 'lightblue', color = 'white', bins = 30) +
  geom_histogram(data = subset(df_ret, logret < as.numeric(VaR)),
                 fill = 'red', color = 'white', bins = 30) +
  geom_vline(xintercept = as.numeric(VaR), linetype = 'dashed') +
  labs(x = 'Rendimiento diario (diferencias logarítmicas)',
       y = 'Frecuencia')

$Histograma de los rendimientos diarios del Bitcoin. En rojo, el 5\% de los días con las mayores pérdidas: la frontera de esa región es el Valor en Riesgo$

Figura 6.3: Histograma de los rendimientos diarios del Bitcoin. En rojo, el 5% de los días con las mayores pérdidas: la frontera de esa región es el Valor en Riesgo

Ahora bien, una de las preguntas que nos podemos hacer es si los rendimientos del Bitcoin se aproximan a una distribución normal. Para ello, la Figura 6.4 ilustra esta comparación, de la cual podemos observar que una prueba de normalidad rechaza esa hipótesis–el estadístico Jarque-Bera indica que la serie de rendimientos no tiene una distribución normal–.

set.seed(1234)

# Distribución normal de referencia: misma media y desviación estándar
normal_dist <- rnorm(100000, mean(logret), sd(logret))

# VaR al 5% que implicaría el supuesto de normalidad
VaR_n <- quantile(normal_dist, 0.05)

ggplot() +
  geom_density(aes(x = as.numeric(logret), colour = 'Rendimientos observados'),
               linewidth = 0.8) +
  geom_density(aes(x = normal_dist, colour = 'Normal de referencia'),
               linewidth = 0.8) +
  scale_colour_manual(values = c('Rendimientos observados' = 'darkred',
                                 'Normal de referencia'    = 'darkblue')) +
  labs(x = 'Rendimiento diario', y = 'Densidad') +
  theme(legend.position = 'bottom', legend.title = element_blank())

Figura 6.4: Densidad estimada de los rendimientos del Bitcoin frente a una distribución normal con la misma media y varianza

La comparación es reveladora: la densidad de los rendimientos observados es más apuntada en el centro y tiene colas más gruesas que la normal. La consecuencia práctica se aprecia al comparar cuantiles. Al 5%, el VaR que implica el supuesto de normalidad (5.64%) es incluso algo más severo que el histórico (5.4%), porque en la zona central la normal tiene más masa que la distribución empírica. Pero al adentrarnos en la cola la situación se invierte: al 1% el cuantil histórico es 10.54% frente a 7.96% bajo normalidad. Es decir, el supuesto de normalidad subestima justamente el riesgo de los eventos extremos, que es el que más importa en la administración de riesgos.

vector_ret <- as.vector(logret)

##Kurtosis
round(kurtosis(vector_ret),2)

## [1] 14.86

##Sesgo
round(skewness(vector_ret),2)

## [1] -0.71

6.1.2 Prueba de normalidad

Formalicemos la observación anterior con la prueba de Jarque-Bera, cuya hipótesis nula es que la serie se distribuye como una normal, es decir, $H_0: K = S = 0$, donde $K$ es el exceso de curtosis y $S$ el sesgo:

jarque.test(vector_ret)

## 
##  Jarque-Bera Normality Test
## 
## data:  vector_ret
## JB = 25617, p-value < 0.00000000000000022
## alternative hypothesis: greater

6.2 Modelos ARCH y GARCH Univariados

Para plantear el modelo, supongamos–por simplicidad–que hemos construido y estimado un modelo AR(1). Es decir, asumamos que el proceso subyacente para la media condicional está dada por: \[\begin{equation} X_t = a_0 + a_1 X_{t-1} + U_t \end{equation}\]

Donde $| a_1 |< 1$ para garantizar la convergencia del proceso en el largo plazo, en el cual: \[\begin{eqnarray*} \mathbb{E}[X_t] & = & \frac{a_0 }{1 - a_1} = \mu \\ Var[X_t] & = & \frac{\sigma^2}{1 - a_1^2} \end{eqnarray*}\]

Ahora, supongamos que este tipo de modelos pueden ser extendidos y generalizados a un modelo ARMA(p, q), que incluya otras variables exógenas. Denotemos a $\mathbf{Z}_t$ como el conjunto que incluye los componentes AR, MA y variables exógenas que pueden explicar a $X_t$ de forma que el proceso estará dado por: \[\begin{equation} X_t = \mathbf{Z}_t \boldsymbol{\beta} + U_t \end{equation}\]

Donde $U_t$ es un proceso estacionario que representa el error asociado a un proceso ARMA(p, q) y donde siguen siendo válidos los supuestos: \[\begin{eqnarray*} \mathbb{E}[U_t] & = & 0 \\ Var[U_t] & = & \sigma^2 \end{eqnarray*}\]

No obstante, en este caso podemos suponer que existe autocorrelación en el término de error al cuadrado que puede ser capturada por un proceso similar a uno de medias móviles (MA) dado por: \[\begin{equation} U_t^2 = \omega + \alpha_1 U_{t-1}^2 + \alpha_2 U_{t-2}^2 + \ldots + \alpha_q U_{t-q}^2 + \nu_t \end{equation}\]

Donde $\nu_t$ es un ruido blanco y $U_{t-i} = X_{t-i} - \mathbf{Z}_{t-i} \boldsymbol{\beta}$, $i = 1, 2 ,$. Si bien los procesos son estacionarios por los supuestos antes enunciados, la varianza condicional estará dada por: \[\begin{eqnarray*} \sigma^2_{t | t-1} & = & Var[ U_t | \Omega_{t-1} ] \\ & = & \mathbb{E}[ U^2_t | \Omega_{t-1} ] \end{eqnarray*}\]

Donde $\Omega_{t-1} = \{U_{t-1}, U_{t-2}, \ldots \}$ es el conjunto de toda la información pasada de $U_t$ y observada hasta el momento $t-1$, por lo que: \[\begin{equation*} U_t | \Omega_{t-1} \sim \mathbb{D}(0, \sigma^2_{t | t-1}) \end{equation*}\]

Así, de forma similar a un proceso MA(q) podemos decir que la varianza condicional tendrá efectos ARCH de orden $q$ (ARCH(q)) cuando: \[\begin{equation} \sigma^2_{t | t-1} = \omega + \alpha_1 U_{t-1}^2 + \alpha_2 U_{t-2}^2 + \ldots + \alpha_q U_{t-q}^2 \tag{6.1} \end{equation}\]

Donde $\mathbb{E}[\nu_t] = 0$, $\omega > 0$ y $\alpha_i \geq 0$, para $i = 1, 2, \ldots, q$. Estas condiciones son necesarias para garantizar que la varianza sea positiva. En general, la varianza condicional se expresa de la forma $\sigma^2_{t | t-1}$, no obstante, para facilitar la notación, nos referiremos en cada caso a esta simplemente como $\sigma^2_{t}$.

Podemos generalizar esta situación si asumimos a la varianza condicional como dependiente de los valores de la varianza rezagados, es decir, como si fuera un proceso AR de orden $p$ para la varianza y juntándolo con la ecuación (6.1). Bollerslev (1986) y Taylor (1986) generalizaron el problema de heterocedasticidad condicional. El modelo se conoce como GARCH(p, q), el cual se especifica como: \[\begin{eqnarray} \sigma^2_t & = & \omega + \alpha_1 U_{t-1}^2 + \alpha_2 U_{t-2}^2 + \ldots + \alpha_q U_{t-q}^2 \\ \nonumber & & + \beta_1 \sigma^2_{t-1} + \beta_2 \sigma^2_{t-2} + \ldots + \beta_p \sigma^2_{t-p} \tag{6.2} \end{eqnarray}\]

Donde las condiciones de no negatividad son que $\omega > 0$, $\alpha_i \geq 0$ para $i = 1, 2, \ldots, q$, y $\beta_j \geq 0$ para $j = 1, 2, \ldots, p$. Además, otra condición de convergencia es que: \[\begin{equation*} \alpha_1 + \ldots + \alpha_q + \beta_1 + \ldots + \beta_p < 1 \end{equation*}\]

Usando el operador rezago $L$ en la ecuación (6.2) podemos obtener: \[\begin{equation} \sigma^2_t = \omega + \alpha(L) U_t^2 + \beta(L) \sigma^2_t \tag{6.3} \end{equation}\]

De donde podemos establecer: \[\begin{equation} \sigma^2_t = \frac{\omega}{1 - \beta(L)} + \frac{\alpha(L)}{1 - \beta(L)} U_t^2 \end{equation}\]

Por lo que la ecuación (6.2) del GARCH(p, q) representa un ARCH($\infty$): \[\begin{equation} \sigma^2_t = \frac{\omega}{1 - \beta_1 - \beta_2 - \ldots - \beta_p} + \sum_{i = 1}^\infty c_i U_{t-i}^2 \end{equation}\]

donde $c_i$ son los coeficientes de la expansión en serie de potencias $\frac{\alpha(L)}{1 - \beta(L)} = \sum_{i=1}^\infty c_i L^i$.

6.2.1 Propiedades: varianza no condicional y colas pesadas

Vale la pena detenerse en las propiedades del caso más simple, el ARCH(1), porque ilustran por qué esta familia de modelos reproduce los hechos estilizados de las series financieras que documentamos en la sección de motivación. Para el ARCH(1), $\sigma^2_t = \omega + \alpha_1 U^2_{t-1}$, y aplicando esperanzas a ambos lados, la varianza no condicional del proceso es constante:

\[\begin{equation} \bar{\sigma}^2 = \mathbb{E}[U_t^2] = \omega + \alpha_1 \mathbb{E}[U_{t-1}^2] \quad \Longrightarrow \quad \bar{\sigma}^2 = \frac{\omega}{1 - \alpha_1}, \tag{6.4} \end{equation}\]

siempre que $\alpha_1 < 1$. Es decir, el proceso $U_t$ es estacionario y homocedástico sin condicionar, aunque su varianza condicional cambie periodo a periodo: la heterocedasticidad es un fenómeno condicional a la información pasada. Esto explica el agrupamiento de volatilidad: un choque grande en $_{t-1}$ eleva $\sigma^2_t$ y hace más probables choques grandes en $_t$, pero en el largo plazo la varianza revierte a $\bar{\sigma}^2$.

Más aún, si el término estandarizado $U_t / \sigma_t$ es normal, la curtosis no condicional del ARCH(1) es:

\[\begin{equation} \kappa = \frac{\mathbb{E}[U_t^4]}{(\mathbb{E}[U_t^2])^2} = \frac{3 (1 - \alpha_1^2)}{1 - 3 \alpha_1^2} > 3, \quad \text{si} \quad 3 \alpha_1^2 < 1. \tag{6.5} \end{equation}\]

Es decir, aun cuando las innovaciones condicionales sean normales, la distribución no condicional de $U_t$ tiene colas más pesadas que la normal. Esta es exactamente la característica que encontramos en los rendimientos del Bitcoin al inicio del capítulo (curtosis muestral muy superior a 3), y es una de las razones del éxito empírico de los modelos ARCH/GARCH en finanzas.

6.2.2 Estimación por Máxima Verosimilitud

Los modelos GARCH no pueden estimarse por mínimos cuadrados ordinarios, pues la varianza condicional $\sigma^2_t$ es una función no lineal de los parámetros y no es observable. La estrategia estándar es la estimación por máxima verosimilitud. Si suponemos que $U_t | \Omega_{t-1} \sim N(0, \sigma_t^2)$, la función log-verosimilitud condicional de la muestra es:

\[\begin{equation} ln L(\theta) = -\frac{1}{2} \sum_{t=1}^{T} \left[ ln(2\pi) + ln \, \sigma^2_t(\theta) + \frac{U_t^2(\theta)}{\sigma^2_t(\theta)} \right], \tag{6.6} \end{equation}\]

donde $\theta$ agrupa los parámetros de la media condicional ($\boldsymbol{\beta}$) y de la varianza condicional ($\omega, \alpha_i, \beta_j$). Para cada valor candidato de $\theta$, la recursión de la ecuación (6.2) genera la trayectoria completa de $\sigma^2_t$, con lo que se evalúa (6.6); un optimizador numérico busca el máximo. Éste es el procedimiento que ejecuta internamente la función ugarchfit() del paquete rugarch que utilizaremos en los ejemplos.

Dos observaciones prácticas. Primera: si la distribución condicional verdadera no es normal, el estimador anterior se interpreta como de cuasi-máxima verosimilitud (QML); sigue siendo consistente, aunque los errores estándar deben corregirse. Segunda: dado que las series financieras suelen exhibir colas pesadas incluso condicionalmente, es común especificar una distribución $t$ de Student para las innovaciones (en rugarch, con el argumento distribution.model = "std").

6.2.3 Ejemplo ARCH(1)

Hasta ahora, las distribuciones utilizadas para medir el Valor en Riesgo de un activo (Bitcoin, en este caso) asumen que no existe correlación serial en los retornos diarios del activo. Observemos un par de gráficas de la función de autocorrelación para corroborar este hecho–ver la Figura 6.5–.

acf(logret)

Figura 6.5: Función de autocorrelación de los rendimientos del Bitcoin

La idea de clusterización de volatilidad asume que períodos de alta volatilidad serán seguidos por alta volatilidad y viceversa. Por esta razón, la función de autocorrelación útil para saber si existen clústers de volatilidad es utilizando el valor absoluto, ya que lo que importa es saber si la serie está autocorrelacionada en la magnitud de los movimientos. La Figura 6.6 muestra el comportamiento de los rendimientos en valor absoluto y la Figura 6.7 su función de autocorrelación.

plot(abs(logret))

Figura 6.6: Evolución de los rendimientos del Bitcoin en valor absoluto

acf(abs(logret))

Figura 6.7: Función de autocorrelación de los rendimientos del Bitcoin en valor absoluto

Otra manera de corroborar esta idea es volviendo IID nuestra serie de datos y observar que de este modo se pierde la autocorrelación serial, lo que refuerza la idea de que en esta serie existen clusters de volatilidad.

logret_random <- sample(as.vector(logret), size =  length(logret), replace = FALSE)

acf(abs(logret_random))

Figura 6.8: Función de autocorrelación del valor absoluto de los rendimientos del Bitcoin reordenados aleatoriamente

par(mfrow = c(1,2))
plot(logret)
plot(logret_random, type = 'l')

Figura 6.9: Rendimientos del Bitcoin: serie original vs. serie reordenada aleatoriamente

Ahora busquemos una prueba formal. La prueba de efectos ARCH de Engle (1982) es una prueba de multiplicadores de Lagrange: se estima la media condicional de la serie, se recuperan los residuales $\hat{U}_t$ y se regresan sus cuadrados sobre sus propios rezagos, \[\hat{U}^2_t = \gamma_0 + \gamma_1 \hat{U}^2_{t-1} + \ldots + \gamma_q \hat{U}^2_{t-q} + \nu_t,\] para contrastar la hipótesis nula $H_0 : \gamma_1 = \ldots = \gamma_q = 0$ (ausencia de efectos ARCH) mediante la estadística $T \cdot R^2$, que se distribuye $\chi^2_{(q)}$. El Cuadro 6.2 reporta el resultado de esa regresión auxiliar y de la prueba para los rendimientos del Bitcoin:

# Media condicional (una constante, por simplicidad) y residuales al cuadrado
logret_mean <- dynlm(logret ~ 1)

ehatsq <- ts(resid(logret_mean)^2)

# Regresión auxiliar de la prueba de efectos ARCH
ARCH_m <- dynlm(ehatsq ~ L(ehatsq))

# Prueba formal de efectos ARCH (Engle, 1982)
arch_lm <- ArchTest(logret, lags = 1, demean = TRUE)

knitr::kable(
  data.frame(
    Concepto = c("Coeficiente de $\\hat{U}^2_{t-1}$",
                 "Estadístico $t$ del coeficiente",
                 "$R^2$ de la regresión auxiliar",
                 "Estadístico ARCH-LM ($T \\cdot R^2$)",
                 "Grados de libertad",
                 "Valor $p$"),
    Valor = c(round(coef(summary(ARCH_m))[2, 1], 4),
              round(coef(summary(ARCH_m))[2, 3], 4),
              round(summary(ARCH_m)$r.squared, 4),
              round(as.numeric(arch_lm$statistic), 4),
              as.numeric(arch_lm$parameter),
              signif(arch_lm$p.value, 4))),
  col.names = c("Concepto", "Valor"),
  caption   = "Prueba de efectos ARCH sobre los rendimientos del Bitcoin: regresión auxiliar y estadística ARCH-LM con un rezago.",
  align     = c("l", "r"),
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.2: Prueba de efectos ARCH sobre los rendimientos del Bitcoin: regresión auxiliar y estadística ARCH-LM con un rezago.
Concepto	Valor
Coeficiente de $\hat{U}^2_{t-1}$	0.1352
Estadístico $t$ del coeficiente	8.9595
$R^2$ de la regresión auxiliar	0.0183
Estadístico ARCH-LM ($T \cdot R^2$)	78.8401
Grados de libertad	1.0000
Valor $p$	0.0000

El coeficiente del residual al cuadrado rezagado es positivo y significativo, y la estadística ARCH-LM rechaza contundentemente la hipótesis nula de ausencia de efectos ARCH. Se justifica, por lo tanto, modelar explícitamente la varianza condicional.

La prueba puede repetirse para diferentes especificaciones de ARCH, sin embargo, para efectos ilustrativos usaremos un ARCH(1). Estimemos un ARCH(1), considerando la siguiente especificación:

\[\begin{eqnarray*} Y_t & = & \mu + U_t, \quad U_t = \sqrt{h_t}\,\varepsilon_t \\ h_t & = & \omega + \alpha_1 U_{t-1}^2 \\ \varepsilon_t & \sim & iid(0,1) \end{eqnarray*}\]

Donde, para no multiplicar la notación, escribimos $h_t \equiv \sigma^2_{t|t-1}$ y $U_t = \sqrt{h_t}\,\varepsilon_t$ es la innovación de la ecuación de la media: nótese que lo que entra en la ecuación de la varianza es el cuadrado de esa innovación, $U^2_{t-1}$, y no el del choque estandarizado $\varepsilon_{t-1}$, cuya varianza es uno por construcción.

Considerando que la media es una AR(2), los resultados de esta estimación se muestran en el Cuadro 6.3. En la estimación empleamos una distribución $t$ de Student estandarizada para las innovaciones (argumento distribution.model = "std" de rugarch), en línea con las colas pesadas de los rendimientos que documentamos al inicio del capítulo.

library(rugarch)

model.spec = ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1, 0)),
                        mean.model     = list(armaOrder = c(2, 0)),
                        distribution.model = "std")

arch.fit = ugarchfit(spec = model.spec, data = logret, solver = 'solnp')

df_arch <- as.data.frame(arch.fit@fit$matcoef)
df_arch$Signif <- ifelse(df_arch[, 4] < 0.001, "***",
                  ifelse(df_arch[, 4] < 0.01,  "**",
                  ifelse(df_arch[, 4] < 0.05,  "*",
                  ifelse(df_arch[, 4] < 0.1,   ".", ""))))
knitr::kable(
  df_arch,
  col.names = c("Estimado", "Error Est.", "Estad. $t$", "Prob. $(>|t|)$", "Signif."),
  caption   = "Estimación del modelo ARCH(1) con media AR(2).",
  align     = c("r", "r", "r", "r", "c"),
  digits    = 4,
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.3: Estimación del modelo ARCH(1) con media AR(2).
	Estimado	Error Est.	Estad. $t$	Prob. $(>\|t\|)$	Signif.
mu	0.0014	0.0003	4.2651	0.0000	***
ar1	-0.0465	0.0152	-3.0598	0.0022	**
ar2	-0.0042	0.0125	-0.3386	0.7349
omega	0.0016	0.0003	4.8146	0.0000	***
alpha1	0.6852	0.1717	3.9900	0.0001	***
shape	2.3936	0.1110	21.5583	0.0000	***

6.2.4 Ejemplo GARCH(1,0)

Estimemos ahora el caso complementario: un modelo en el que la varianza condicional depende únicamente de su propio rezago y no de los errores al cuadrado. Con la convención $GARCH(p, q)$ de la ecuación (6.2) –donde $p$ es el número de rezagos de $\sigma^2_t$ y $q$ el de rezagos de $U^2_t$– se trata de un $GARCH(1, 0)$: \[\begin{eqnarray*} Y_t & = & \mu + U_t, \quad U_t = \sqrt{h_t}\,\varepsilon_t \\ h_t & = & \omega + \beta_1 h_{t-1} \\ \varepsilon_t & \sim & iid(0,1) \end{eqnarray*}\]

Conviene advertir aquí una fuente frecuente de confusión: el argumento garchOrder de ugarchspec() recibe los órdenes en el orden inverso al de nuestra notación, es decir, garchOrder = c(q, p). Por eso el $GARCH(1,0)$ se especifica con garchOrder = c(0, 1).

Considerando que la media es una AR(2), los resultados de esta estimación se muestran en el Cuadro 6.4.

model.spec = ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(0, 1)),
                        mean.model     = list(armaOrder = c(2, 0)),
                        distribution.model = "std")

fit.garch.n = ugarchfit(spec = model.spec, data = logret, solver = "solnp")

df_garch01 <- as.data.frame(fit.garch.n@fit$matcoef)
df_garch01$Signif <- ifelse(df_garch01[, 4] < 0.001, "***",
                     ifelse(df_garch01[, 4] < 0.01,  "**",
                     ifelse(df_garch01[, 4] < 0.05,  "*",
                     ifelse(df_garch01[, 4] < 0.1,   ".", ""))))
knitr::kable(
  df_garch01,
  col.names = c("Estimado", "Error Est.", "Estad. $t$", "Prob. $(>|t|)$", "Signif."),
  caption   = "Estimación del modelo GARCH(1,0) con media AR(2).",
  align     = c("r", "r", "r", "r", "c"),
  digits    = 4,
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.4: Estimación del modelo GARCH(1,0) con media AR(2).
	Estimado	Error Est.	Estad. $t$	Prob. $(>\|t\|)$	Signif.
mu	0.0014	0.0003	4.1675	0.0000	***
ar1	-0.0572	0.0129	-4.4362	0.0000	***
ar2	-0.0040	0.0124	-0.3197	0.7492
omega	0.0000	0.0000	59.1393	0.0000	***
beta1	0.9976	0.0000	21619.8572	0.0000	***
shape	2.5751	0.0283	91.0808	0.0000	***

6.2.5 Selección GARCH(p,q) óptimo

¿Cómo seleccionamos el orden adecuado para un GARCH(p,q)? Una respuesta la ofrecen los criterios de información, partiendo de la especificación general:

$Y_t = \mu + U_t, \quad U_t = \sqrt{h_t}\,\varepsilon_t$

$h_t = \omega + \sum_{i=1}^p \beta_i h_{t-i} + \sum_{j=1}^q \alpha_j U^2_{t-j}$

$\varepsilon_t \sim iid(0,1)$

6.2.5.1 Criterios de información

ic <- infocriteria(fit.garch.n)
knitr::kable(
  data.frame(Criterio = rownames(ic), Valor = round(as.numeric(ic), 6)),
  caption  = "Criterios de información del modelo GARCH(1,0).",
  align    = c("l", "r"),
  booktabs = TRUE
)

Cuadro 6.5: Criterios de información del modelo GARCH(1,0).
Criterio	Valor
Akaike	-4.164250
Bayes	-4.155385
Shibata	-4.164253
Hannan-Quinn	-4.161119

6.2.5.2 Selección del modelo óptimo

Podemos hacer una búsqueda del mejor modelo de entre varios que probemos en un espectro de hasta un GARCH(4,4), manteniendo para la media condicional la misma especificación AR(2) empleada en las estimaciones anteriores. Los resultados son los reportados en el siguiente Cuadro.

source("lag_opt_GARCH.R")
crit_garch <- as.data.frame(Lag_Opt_GARCH(logret, 4, 4, arma = c(2, 0)))
knitr::kable(
  crit_garch,
  col.names = c("$q$", "$p$", "AIC", "Óptimo"),
  caption   = "Criterios de información AIC en función de $p$ y $q$ para GARCH($p$,$q$).",
  align     = c("c", "c", "r", "c"),
  digits    = 5,
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.6: Criterios de información AIC en función de $p$ y $q$ para GARCH($p$,$q$).
$q$	$p$	AIC	Óptimo
1	1	-4.29787	0
1	2	-4.29896	0
1	3	-4.30017	1
1	4	-4.29971	0
2	1	-4.29741	0
2	2	-4.29849	0
2	3	-4.29989	0
2	4	-4.29937	0
3	1	-4.29672	0
3	2	-4.29780	0
3	3	-4.29943	0
3	4	-4.29890	0
4	1	-4.29633	0
4	2	-4.29762	0
4	3	-4.29882	0
4	4	-4.29844	0

6.2.5.3 Estimación de modelo óptimo

De esta forma, con los datos al corte señalado, el modelo óptimo es un GARCH(3,1). Los resultados de la estimación se muestran en el Cuadro 6.7.

model.spec = ugarchspec(variance.model = list(model = 'sGARCH', garchOrder = c(1, 3)),
                        mean.model     = list(armaOrder = c(2, 0)),
                        distribution.model = "std")

model.fit = ugarchfit(spec = model.spec, data = logret, solver = 'solnp')

df_opt <- as.data.frame(model.fit@fit$matcoef)
df_opt$Signif <- ifelse(df_opt[, 4] < 0.001, "***",
                 ifelse(df_opt[, 4] < 0.01,  "**",
                 ifelse(df_opt[, 4] < 0.05,  "*",
                 ifelse(df_opt[, 4] < 0.1,   ".", ""))))
knitr::kable(
  df_opt,
  col.names = c("Estimado", "Error Est.", "Estad. $t$", "Prob. $(>|t|)$", "Signif."),
  caption   = "Estimación del modelo GARCH(3,1) óptimo con media AR(2).",
  align     = c("r", "r", "r", "r", "c"),
  digits    = 4,
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.7: Estimación del modelo GARCH(3,1) óptimo con media AR(2).
	Estimado	Error Est.	Estad. $t$	Prob. $(>\|t\|)$	Signif.
mu	0.0011	0.0003	3.8379	0.0001	***
ar1	-0.0451	0.0145	-3.1025	0.0019	**
ar2	-0.0031	0.0131	-0.2366	0.8129
omega	0.0000	0.0000	3.2907	0.0010	***
alpha1	0.1829	0.0228	8.0070	0.0000	***
beta1	0.3148	0.1230	2.5594	0.0105
beta2	0.1457	0.1159	1.2572	0.2087
beta3	0.3556	0.0861	4.1300	0.0000	***
shape	3.2339	0.1473	21.9486	0.0000	***

6.2.5.4 Pronósticos con el modelo GARCH óptimo

Antes de pasar al código, conviene precisar qué pronostica un GARCH. Para un GARCH(1,1), la varianza condicional esperada a $h$ periodos satisface la recursión $\mathbb{E}_t[\sigma^2_{t+h}] = \omega + (\alpha_1 + \beta_1) \, \mathbb{E}_t[\sigma^2_{t+h-1}]$, cuya solución es:

\[\begin{equation} \mathbb{E}_t[\sigma^2_{t+h}] = \bar{\sigma}^2 + (\alpha_1 + \beta_1)^{h-1} \left( \sigma^2_{t+1} - \bar{\sigma}^2 \right), \quad \text{con} \quad \bar{\sigma}^2 = \frac{\omega}{1 - \alpha_1 - \beta_1}. \tag{6.7} \end{equation}\]

Es decir, el pronóstico de la varianza revierte geométricamente hacia la varianza no condicional $\bar{\sigma}^2$ a una velocidad gobernada por la persistencia $\alpha_1 + \beta_1$: cuanto más cercana a 1 sea la suma, más lentamente se disipa un episodio de alta volatilidad. La misma lógica se extiende al GARCH(p, q) con persistencia $\sum_i \alpha_i + \sum_j \beta_j$.

Para realizar pronósticos con la estimación de un GARCH, utilizando el paquete rugarch, es necesario utilizar la función ugarchforecast(). Emplearemos el modelo GARCH(3,1) con media AR(2) que resultó óptimo en la sección anterior (objeto model.fit):

spec = getspec(model.fit)

setfixed(spec) <- as.list(coef(model.fit))

Esta función precisa como argumentos nuestra estimación del modelo GARCH, con una modificación en la manera en que se presentan los coeficientes, realizada en la última línea del código anterior y que llamamos spec. n.ahead es el número de periodos que vamos a pronosticar, n.roll señala el número de pronósticos móviles que utilizaremos, en caso de que haya más información para realizar el pronóstico. Finalmente damos como input nuestro set de datos y como producto obtendremos el pronostico de Sigma tanto como de la serie.

forecast = ugarchforecast(spec, n.ahead = 12, n.roll = 0, logret)

sigma(forecast)

##      4311-01-01
## T+1  0.02073343
## T+2  0.02168883
## T+3  0.02195810
## T+4  0.02215006
## T+5  0.02260453
## T+6  0.02294560
## T+7  0.02324236
## T+8  0.02359079
## T+9  0.02392004
## T+10 0.02423233
## T+11 0.02455141
## T+12 0.02486508

fitted(forecast)

##       4311-01-01
## T+1  0.001665495
## T+2  0.001152099
## T+3  0.001136026
## T+4  0.001138346
## T+5  0.001138291
## T+6  0.001138286
## T+7  0.001138287
## T+8  0.001138287
## T+9  0.001138287
## T+10 0.001138287
## T+11 0.001138287
## T+12 0.001138287

forecast

## 
## *------------------------------------*
## *       GARCH Model Forecast         *
## *------------------------------------*
## Model: sGARCH
## Horizon: 12
## Roll Steps: 0
## Out of Sample: 0
## 
## 0-roll forecast [T0=]:
##        Series   Sigma
## T+1  0.001665 0.02073
## T+2  0.001152 0.02169
## T+3  0.001136 0.02196
## T+4  0.001138 0.02215
## T+5  0.001138 0.02260
## T+6  0.001138 0.02295
## T+7  0.001138 0.02324
## T+8  0.001138 0.02359
## T+9  0.001138 0.02392
## T+10 0.001138 0.02423
## T+11 0.001138 0.02455
## T+12 0.001138 0.02487

6.3 Extensiones: asimetría y GARCH en media

El GARCH(p, q) trata de forma simétrica a los choques: en la ecuación (6.2) solo importa la magnitud $U^2_{t-1}$, no su signo. Sin embargo, un hecho estilizado adicional de los mercados financieros es el efecto apalancamiento (leverage effect): las caídas de precios suelen incrementar la volatilidad más que las alzas de la misma magnitud. Existen dos extensiones clásicas para capturar esta asimetría.

El modelo TGARCH o GJR-GARCH de Glosten, Jagannathan y Runkle (1993) añade un término que solo se activa cuando el choque es negativo:

\[\begin{equation} \sigma^2_t = \omega + \alpha_1 U^2_{t-1} + \gamma \, U^2_{t-1} \, \mathbb{I}(U_{t-1} < 0) + \beta_1 \sigma^2_{t-1}, \tag{6.8} \end{equation}\]

donde $\mathbb{I}(\cdot)$ es la función indicadora. Si $\gamma > 0$, un choque negativo incrementa la varianza condicional en $\alpha_1 + \gamma$, mientras que uno positivo solo en $\alpha_1$.

El modelo EGARCH de Nelson (1991) especifica el logaritmo de la varianza condicional en función de los choques estandarizados $z_t = U_t / \sigma_t$:

\[\begin{equation} ln \, \sigma^2_t = \omega + \beta_1 \, ln \, \sigma^2_{t-1} + \alpha_1 \left( |z_{t-1}| - \mathbb{E}|z_{t-1}| \right) + \gamma \, z_{t-1}, \tag{6.9} \end{equation}\]

con dos ventajas: al modelar $ln \, \sigma^2_t$, no se requieren restricciones de no negatividad sobre los parámetros, y el término $\gamma \, z_{t-1}$ captura la asimetría de forma directa (un $\gamma < 0$ implica efecto apalancamiento).

Por último, el modelo GARCH en media (GARCH-M) de Engle, Lilien y Robins (1987) permite que la varianza condicional afecte a la media condicional:

\[\begin{equation} X_t = \mathbf{Z}_t \boldsymbol{\beta} + \lambda \sigma^2_t + U_t, \tag{6.10} \end{equation}\]

donde $\lambda$ se interpreta como el premio al riesgo: si $\lambda > 0$, los periodos de mayor volatilidad esperada se compensan con mayores rendimientos esperados, en línea con la teoría de portafolios. Todas estas variantes están disponibles en rugarch: basta cambiar el argumento variance.model = list(model = "gjrGARCH") o "eGARCH" en ugarchspec(), o añadir archm = TRUE en la especificación de la media para el GARCH-M.

6.4 Modelos ARCH y GARCH Multivariados

De forma similar a los modelos univariados, los modelos multivariados de heterocedasticidad condicional asumen una estructura de la media condicional. En este caso, descrita por un VAR(p) cuyo proceso estocástico $\mathbf{X}$ es estacionario de dimensión $k$. De esta forma, la expresión reducida del modelo o el proceso VAR(p) estará dado por: \[\begin{equation} \mathbf{X}_t = \boldsymbol{\delta} + \mathbf{A_1} \mathbf{X}_{t-1} + \mathbf{A_2} \mathbf{X}_{t-2} + \ldots + \mathbf{A_p} \mathbf{X}_{t-p} + \mathbf{U}_{t} \end{equation}\]

Donde cada uno de las $\mathbf{A_i}$, $i = 1, 2, \ldots, p$, son matrices cuadradas de dimensión $k$ y $\mathbf{U}_t$ representa un vector de dimensión $k \times 1$ con los residuales en el momento del tiempo $t$ que son un proceso puramente aleatorio. También se incorpora un vector de términos constantes denominado como $\boldsymbol{\delta}$, el cual es de dimensión $k \times 1$ –en este caso también es posible incorporar procesos determinísticos adicionales–.

Así, suponemos que el término de error tendrá estructura de vector: \[\begin{equation*} \mathbf{U}_t = \begin{bmatrix} U_{1t} \\ U_{2t} \\ \vdots \\ U_{Kt} \end{bmatrix} \end{equation*}\]

De forma que diremos que: \[\begin{equation*} \mathbf{U}_t | \Omega_{t-1} \sim (0, \Sigma_{t | t-1}) \end{equation*}\]

Dicho lo anterior, entonces, el modelo ARCH(q) multivariado será descrito por: \[\begin{equation} Vech(\Sigma_{t | t-1}) = \boldsymbol{\gamma}_0 + \Gamma_1 Vech(\mathbf{U}_{t-1} \mathbf{U}_{t-1}') + \ldots + \Gamma_q Vech(\mathbf{U}_{t-q} \mathbf{U}_{t-q}') \tag{6.11} \end{equation}\]

Donde $Vech$ es un operador que apila en un vector la parte triangular inferior (incluyendo la diagonal) de la matriz a la cual se le aplique, $\boldsymbol{\gamma}_0$ es un vector de constantes, $\Gamma_i$, $i = 1, 2, \ldots$ son matrices de coeficientes asociados a la estimación.

Para ilustrar la ecuación (6.11), tomemos un ejemplo de $K = 2$, de esta forma tenemos que un M-ARCH(1) será: \[\begin{equation*} \Sigma_{t | t-1} = \begin{bmatrix} \sigma^2_{1, t | t-1} & \sigma_{12, t | t-1} \\ \sigma_{21, t | t-1} & \sigma^2_{2, t | t-1} \end{bmatrix} = \begin{bmatrix} \sigma_{11, t} & \sigma_{12, t} \\ \sigma_{21, t} & \sigma_{22, t} \end{bmatrix} = \Sigma_{t} \end{equation*}\]

Donde hemos simplificado la notación de las varianzas y la condición de que están en función de $t-1$. Así, \[\begin{equation*} Vech(\Sigma_{t}) = Vech \begin{bmatrix} \sigma_{11, t} & \sigma_{12, t} \\ \sigma_{21, t} & \sigma_{22, t} \end{bmatrix} = \begin{bmatrix} \sigma_{11, t} \\ \sigma_{12, t} \\ \sigma_{22, t} \end{bmatrix} \end{equation*}\]

De esta forma, podemos establecer el modelo M-ARCH(1) con $K = 2$ será de la forma: \[\begin{equation*} \begin{bmatrix} \sigma_{11, t} \\ \sigma_{12, t} \\ \sigma_{22, t} \end{bmatrix} = \begin{bmatrix} \gamma_{10} \\ \gamma_{20} \\ \gamma_{30} \end{bmatrix} + \begin{bmatrix} \gamma_{11} & \gamma_{12} & \gamma_{13} \\ \gamma_{21} & \gamma_{22} & \gamma_{23} \\ \gamma_{31} & \gamma_{32} & \gamma_{33} \end{bmatrix} \begin{bmatrix} U^2_{1, t-1} \\ U_{1, t-1} U_{2, t-1} \\ U^2_{2, t-1} \end{bmatrix} \end{equation*}\]

Como notarán, este tipo de procedimientos implica la estimación de muchos parámetros. En esta circunstancia, se suelen estimar modelos restringidos para reducir el número de coeficientes estimados. Por ejemplo, podríamos querer estimar un caso como: \[\begin{equation*} \begin{bmatrix} \sigma_{11, t} \\ \sigma_{12, t} \\ \sigma_{22, t} \end{bmatrix} = \begin{bmatrix} \gamma_{10} \\ \gamma_{20} \\ \gamma_{30} \end{bmatrix} + \begin{bmatrix} \gamma_{11} & 0 & 0 \\ 0 & \gamma_{22} & 0 \\ 0 & 0 & \gamma_{33} \end{bmatrix} \begin{bmatrix} U^2_{1, t-1} \\ U_{1, t-1} U_{2, t-1} \\ U^2_{2, t-1} \end{bmatrix} \end{equation*}\]

Finalmente y de forma análoga al caso univariado, podemos plantear un modelo M-GARCH(p, q) como: \[\begin{eqnarray} Vech(\Sigma_{t | t-1}) & = & \boldsymbol{\gamma}_0 + \sum_{j = 1}^q \Gamma_j Vech(\mathbf{U}_{t-j} \mathbf{U}_{t-j}') \nonumber \\ & & + \sum_{m = 1}^p \mathbf{G}_m Vech(\Sigma_{t-m | t-m-1}) \tag{6.12} \end{eqnarray}\]

Donde cada una de las $\mathbf{G}_m$ es una matriz de coeficientes. Para ilustrar este caso, retomemos el ejemplo anterior, pero ahora para un modelo M-GARCH(1, 1) con $K = 2$ de forma que tendríamos: \[\begin{eqnarray*} \begin{bmatrix} \sigma_{11, t} \\ \sigma_{12, t} \\ \sigma_{22, t} \end{bmatrix} & = & \begin{bmatrix} \gamma_{10} \\ \gamma_{20} \\ \gamma_{30} \end{bmatrix} + \begin{bmatrix} \gamma_{11} & \gamma_{12} & \gamma_{13} \\ \gamma_{21} & \gamma_{22} & \gamma_{23} \\ \gamma_{31} & \gamma_{32} & \gamma_{33} \end{bmatrix} \begin{bmatrix} U^2_{1, t-1} \\ U_{1, t-1} U_{2, t-1} \\ U^2_{2, t-1} \end{bmatrix} \\ & & + \begin{bmatrix} g_{11} & g_{12} & g_{13} \\ g_{21} & g_{22} & g_{23} \\ g_{31} & g_{32} & g_{33} \end{bmatrix} \begin{bmatrix} \sigma_{11, t-1} \\ \sigma_{12, t-1} \\ \sigma_{22, t-1} \end{bmatrix} \end{eqnarray*}\]

6.4.1 Especificaciones restringidas de uso frecuente

La formulación $Vech$ anterior es la más general, pero también la menos práctica: con $K$ series y un modelo de orden $(1,1)$ se deben estimar del orden de $K^4/2$ parámetros –para $K = 5$, más de trescientos–, además de que hay que garantizar que la matriz $\Sigma_{t|t-1}$ resulte definida positiva en cada periodo. Por ello, en el trabajo aplicado se utilizan versiones restringidas. Las tres familias más comunes son:

VECH diagonal y BEKK. La versión diagonal impone que las matrices $\Gamma_j$ y $\mathbf{G}_m$ sean diagonales, de modo que cada varianza y cada covarianza dependen únicamente de su propio pasado –el caso ilustrado arriba–. La parametrización BEKK de Engle y Kroner (1995) escribe en cambio $\Sigma_t = \mathbf{C}'\mathbf{C} + \mathbf{A}' \mathbf{U}_{t-1} \mathbf{U}_{t-1}' \mathbf{A} + \mathbf{B}' \Sigma_{t-1} \mathbf{B}$, forma cuadrática que garantiza por construcción que $\Sigma_t$ sea definida positiva.
Correlación condicional constante (CCC). Bollerslev (1990) propone estimar un GARCH univariado para cada serie y suponer que la matriz de correlaciones $\mathbf{R}$ es constante, de modo que $\Sigma_t = \mathbf{D}_t \mathbf{R} \mathbf{D}_t$, donde $\mathbf{D}_t$ es la matriz diagonal de desviaciones estándar condicionales. El número de parámetros se reduce drásticamente, a costa de un supuesto fuerte: las correlaciones no cambian en el tiempo.
Correlación condicional dinámica (DCC). Engle (2002) relaja ese supuesto permitiendo que la matriz de correlaciones evolucione, $\Sigma_t = \mathbf{D}_t \mathbf{R}_t \mathbf{D}_t$, con $\mathbf{R}_t$ gobernada por una recursión tipo GARCH sobre los residuales estandarizados. Es hoy la especificación más utilizada en finanzas, porque combina parsimonia con la posibilidad de estudiar el contagio entre mercados. En R se estima con las funciones dccspec() y dccfit() del paquete rmgarch, que se apoyan en las estimaciones univariadas de rugarch. Es la especificación que aplicaremos en el ejemplo de la siguiente sección.

6.5 Pruebas para detectar efectos ARCH

La prueba que mostraremos es conocida como una ARCH-LM, la cual está basada en una regresión de los residuales estimados de un modelo VAR(p) o cualquier otra estimación que deseemos probar, con el objeto de determinar si existen efectos ARCH –esta prueba se puede simplificar para el caso univariado–.

Partamos de plantear: \[\begin{eqnarray} Vech(\hat{\mathbf{U}}_t \hat{\mathbf{U}}_t') & = & \mathbf{B}_0 + \mathbf{B}_1 Vech(\hat{\mathbf{U}}_{t-1} \hat{\mathbf{U}}_{t-1}') + \ldots \\ \nonumber & & + \mathbf{B}_q Vech(\hat{\mathbf{U}}_{t-q} \hat{\mathbf{U}}_{t-q}') + \varepsilon_t \tag{6.13} \end{eqnarray}\]

Dada la estimación en la ecuación (6.13), plantearemos la estructura de hipótesis dada por: \[\begin{eqnarray*} H_0 & : & \mathbf{B}_1 = \mathbf{B}_2 = \ldots = \mathbf{B}_q = \mathbf{0} \\ H_a & : & \text{al menos una } \mathbf{B}_i \neq \mathbf{0} \end{eqnarray*}\]

La estadística de prueba será determinada por: \[\begin{equation} LM_{M-ARCH} = \frac{1}{2} T K (K + 1) - T \cdot Traza \left( \hat{\Sigma}_{ARCH} \hat{\Sigma}^{-1}_{0} \right) \sim \chi^2_{[q K^2 (K + 1)^2 / 4]} \end{equation}\]

Donde la matriz $\hat{\Sigma}_{ARCH}$ se calcula de acuerdo con la ecuación (6.13) y la matriz $\hat{\Sigma}_{0}$ sin considerar una estructura dada para los errores.

6.6 Ejemplo: correlación dinámica entre Bitcoin y Ethereum

Apliquemos el modelo DCC a un sistema bivariado. A la serie del Bitcoin que hemos usado en todo el capítulo le agregamos la del Ethereum (ticker ETH-USD), congelada con la misma fecha de corte en BD/ETH_USD.RData. La pregunta de interés es de manejo de riesgo: ¿qué tan estable es la correlación entre ambos criptoactivos? Si fuera constante, una cartera que los combine tendría un beneficio de diversificación predecible; si se dispara justo en los episodios de turbulencia, la diversificación falla precisamente cuando más se necesita.

load("BD/BTC_USD.RData")
load("BD/ETH_USD.RData")

# Rendimientos diarios de ambos activos, sobre las fechas en común
ret2 <- na.omit(merge(diff(log(BTC$`BTC-USD.Adjusted`)),
                      diff(log(ETH$`ETH-USD.Adjusted`)),
                      join = "inner"))
colnames(ret2) <- c("BTC", "ETH")

La muestra conjunta tiene 3162 observaciones, de 10/11/2017 a 07/07/2026 –el Ethereum empieza a cotizar después que el Bitcoin, de modo que el periodo común es más corto–, y su correlación no condicional es de 0.796.

Antes de estimar conviene verificar que el sistema en efecto presenta efectos ARCH multivariados, con la prueba que planteamos en la sección anterior. La función MarchTest() del paquete MTS reporta varias versiones del contraste:

prueba_march <- MTS::MarchTest(as.matrix(ret2), lag = 5)

## Q(m) of squared series(LM test):  
## Test statistic:  530.8596  p-value:  0 
## Rank-based Test:  
## Test statistic:  892.6097  p-value:  0 
## Q_k(m) of squared series:  
## Test statistic:  158.2136  p-value:  0 
## Robust Test(5%) :  219.5226  p-value:  0

Las cuatro versiones rechazan la hipótesis nula de ausencia de efectos ARCH con valores $p$ nulos, de modo que modelar la matriz de varianzas y covarianzas condicional está plenamente justificado. Especificamos entonces un DCC(1,1) con márgenes $GARCH(1,1)$ y distribución $t$ de Student –coherentes con las colas pesadas que documentamos al inicio del capítulo– y una $t$ multivariada para el sistema:

esp_uni <- ugarchspec(variance.model = list(model = "sGARCH",
                                            garchOrder = c(1, 1)),
                      mean.model     = list(armaOrder = c(0, 0)),
                      distribution.model = "std")

esp_dcc <- dccspec(uspec = multispec(replicate(2, esp_uni)),
                   dccOrder = c(1, 1), distribution = "mvt")

dcc_fit <- dccfit(esp_dcc, data = ret2)

par_dcc <- dcc_fit@mfit$matcoef[c("[Joint]dcca1", "[Joint]dccb1",
                                  "[Joint]mshape"), ]
rownames(par_dcc) <- c("$\\alpha$ (DCC)", "$\\beta$ (DCC)",
                       "Grados de libertad")

knitr::kable(
  par_dcc,
  col.names = c("Estimado", "Error Est.", "Estad. $t$", "Prob. $(>|t|)$"),
  caption   = "Parámetros de la ecuación de correlación dinámica del modelo DCC-GARCH(1,1) para los rendimientos de Bitcoin y Ethereum.",
  align     = c("r", "r", "r", "r"),
  digits    = 4,
  booktabs  = TRUE,
  escape    = FALSE
)

Cuadro 6.8: Parámetros de la ecuación de correlación dinámica del modelo DCC-GARCH(1,1) para los rendimientos de Bitcoin y Ethereum.
	Estimado	Error Est.	Estad. $t$
$\alpha$ (DCC)	0.0491	0.0081	6.0412
$\beta$ (DCC)	0.9453	0.0103	91.9684
Grados de libertad	4.0000	0.7160	5.5870

Los dos parámetros de la ecuación de correlación del Cuadro 6.8 son significativos, lo que descarta de entrada la hipótesis de correlación constante: si $\alpha = \beta = 0$, el DCC colapsa al modelo CCC de Bollerslev (1990). Su suma, 0.9944, se interpreta igual que la persistencia de un GARCH univariado: la correlación se mueve, pero lo hace con mucha inercia, de modo que los desplazamientos son duraderos. Los grados de libertad estimados de la $t$ multivariada (4) confirman además que la distribución conjunta tiene colas mucho más pesadas que la normal.

rho_t <- xts(rcor(dcc_fit)[1, 2, ], order.by = index(ret2))

plot(zoo(rho_t), col = "darkblue", lwd = 1,
     xlab = "Tiempo", ylab = "Correlación condicional",
     main = "Correlación dinámica entre Bitcoin y Ethereum")
abline(h = cor(ret2)[1, 2], lty = 2, col = "darkred")

Correlación condicional entre los rendimientos diarios de Bitcoin y Ethereum estimada con un modelo DCC-GARCH(1,1). La línea punteada es la correlación no condicional de toda la muestra

Figura 6.10: Correlación condicional entre los rendimientos diarios de Bitcoin y Ethereum estimada con un modelo DCC-GARCH(1,1). La línea punteada es la correlación no condicional de toda la muestra

La Figura 6.10 muestra por qué el supuesto de correlación constante habría sido inadecuado. En los primeros meses de la muestra la correlación es baja e inestable –llega a caer a 0.014 en 11/12/2017, cuando el Ethereum era un activo aún incipiente–, y después se estabiliza en niveles altos. El punto relevante para la administración de riesgos es lo que ocurre en los episodios de estrés: durante el desplome de marzo de 2020 la correlación promedió 0.914, es decir, los dos activos cayeron prácticamente al unísono. Es el fenómeno que en la literatura financiera se conoce como contagio, y la razón por la que una cartera diversificada entre criptoactivos ofrece mucha menos protección de lo que sugeriría su correlación promedio.

6.7 Resumen del capítulo

En este capítulo introdujimos los modelos de heterocedasticidad condicional, diseñados para capturar la variabilidad agrupada (volatility clustering) que caracteriza a las series financieras de alta frecuencia.

Motivación. Las series de rendimientos financieros presentan dos hechos estilizados que los modelos ARMA estándar no capturan: (i) períodos de alta volatilidad seguidos de alta volatilidad y períodos de baja volatilidad seguidos de baja volatilidad, y (ii) colas pesadas en la distribución de los rendimientos. La prueba ARCH-LM de Engle (1982) contrasta formalmente si los residuales al cuadrado presentan autocorrelación, lo que indica efectos ARCH.

Modelos ARCH y GARCH univariados. El modelo ARCH($q$) de Engle (1982) especifica la varianza condicional como función de los cuadrados de los errores pasados: \[h_t = \omega + \alpha_1 U_{t-1}^2 + \ldots + \alpha_q U_{t-q}^2\] El modelo GARCH($p,q$) de Bollerslev (1986) añade términos de varianza condicional rezagada, logrando parsimonia frente al ARCH($\infty$) equivalente: \[h_t = \omega + \sum_{i=1}^q \alpha_i U_{t-i}^2 + \sum_{j=1}^p \beta_j h_{t-j}\] La condición de estacionariedad débil requiere $\sum_{i} \alpha_i + \sum_{j} \beta_j < 1$; bajo ella el proceso tiene varianza no condicional constante $\bar{\sigma}^2 = \omega / (1 - \sum_i \alpha_i - \sum_j \beta_j)$ y colas más pesadas que la normal, y los pronósticos de varianza revierten geométricamente a $\bar{\sigma}^2$ a la velocidad de la persistencia $\sum_i \alpha_i + \sum_j \beta_j$. La estimación se realiza por máxima verosimilitud (o cuasi-máxima verosimilitud) evaluando recursivamente la varianza condicional. La selección del orden $(p,q)$ se realiza con los criterios AIC, BIC, Shibata y Hannan-Quinn disponibles en la función infocriteria() del paquete rugarch. En la práctica, el GARCH(1,1) es suficiente para la mayoría de las series financieras.

Extensiones. El TGARCH/GJR-GARCH y el EGARCH capturan el efecto apalancamiento (los choques negativos elevan la volatilidad más que los positivos), mientras que el GARCH-M incorpora la varianza condicional en la ecuación de la media como premio al riesgo.

Valor en Riesgo (VaR). El VaR al nivel $\alpha$ se obtiene como el cuantil $\alpha$ de la distribución de rendimientos. Los modelos GARCH mejoran la estimación dinámica del VaR al actualizar la varianza condicional periodo a periodo.

Modelos ARCH y GARCH multivariados (M-GARCH). En sistemas de $K$ variables, el M-GARCH extiende el modelo univariado apilando las varianzas y covarianzas condicionales mediante el operador $Vech(\cdot)$. La proliferación de parámetros obliga a trabajar con versiones restringidas: VECH diagonal, BEKK, correlación condicional constante (CCC) y correlación condicional dinámica (DCC). La prueba M-ARCH-LM generaliza la prueba univariada y sigue una distribución $\chi^2$ bajo la hipótesis nula de no efectos ARCH. El ejemplo del par Bitcoin-Ethereum ilustra el resultado que motiva todo el enfoque: la correlación entre dos activos no es constante y tiende a elevarse en los episodios de estrés, justo cuando la diversificación debería proteger.

El capítulo siguiente extiende el análisis a datos en panel, permitiendo explotar tanto la dimensión temporal como la transversal de los datos.

6.8 Ejercicios

Los siguientes ejercicios combinan preguntas teóricas con aplicaciones en R. Los marcados como (Teórico) se resuelven analíticamente; los marcados como (Computacional) utilizan los paquetes quantmod, FinTS y rugarch empleados en este capítulo. Tenga presente que los ejercicios con datos descargados en vivo producirán resultados distintos según la fecha de descarga.

(Teórico) Considere un proceso ARCH(1): $\sigma^2_t = \omega + \alpha_1 U^2_{t-1}$, con $\omega > 0$ y $0 \leq \alpha_1 < 1$.
1. Defina $\nu_t = U^2_t - \sigma^2_t$ y demuestre que el cuadrado del error sigue un proceso AR(1): $U^2_t = \omega + \alpha_1 U^2_{t-1} + \nu_t$. Verifique que $\mathbb{E}[\nu_t | \Omega_{t-1}] = 0$.
2. A partir del inciso anterior, obtenga la varianza no condicional $\bar{\sigma}^2 = \omega / (1 - \alpha_1)$ y explique por qué se requiere $\alpha_1 < 1$.
3. Explique por qué las condiciones $\omega > 0$ y $\alpha_1 \geq 0$ son necesarias para que la varianza condicional esté bien definida.
4. Utilice la ecuación del inciso a. para explicar el fenómeno de agrupamiento de volatilidad (volatility clustering): ¿por qué un choque grande en $t-1$ hace más probable observar choques grandes en $t$?
(Teórico) Considere un proceso GARCH(1,1): $\sigma^2_t = \omega + \alpha_1 U^2_{t-1} + \beta_1 \sigma^2_{t-1}$, con $\alpha_1 + \beta_1 < 1$.
1. Demuestre que el pronóstico de la varianza condicional a $h$ pasos satisface $\mathbb{E}[\sigma^2_{t+h} | \Omega_t] - \bar{\sigma}^2 = (\alpha_1 + \beta_1)^{h-1} \left( \sigma^2_{t+1} - \bar{\sigma}^2 \right)$, donde $\bar{\sigma}^2 = \omega / (1 - \alpha_1 - \beta_1)$.
2. Concluya que los pronósticos de varianza revierten geométricamente a $\bar{\sigma}^2$ y que la velocidad de reversión depende de la persistencia $\alpha_1 + \beta_1$.
3. Defina la vida media de un choque de volatilidad como el horizonte $h^*$ tal que $(\alpha_1 + \beta_1)^{h^*} = 0.5$, es decir, $h^* = ln(0.5) / ln(\alpha_1 + \beta_1)$. Calcúlela para $\alpha_1 + \beta_1 = 0.95$ y para $\alpha_1 + \beta_1 = 0.99$, y comente la diferencia.
(Teórico) Sobre las extensiones asimétricas estudiadas en este capítulo:
1. Escriba la especificación de la varianza condicional del modelo GJR-GARCH(1,1), incluyendo la variable indicadora de choques negativos.
2. Explique qué signo se espera para el coeficiente de asimetría cuando existe efecto apalancamiento (leverage effect) en una serie de rendimientos accionarios y cuál es la intuición económica de dicho efecto.
3. Señale dos diferencias entre el GJR-GARCH y el EGARCH, incluyendo la razón por la cual el EGARCH no requiere restricciones de no negatividad sobre sus parámetros.
4. Explique en qué situación es preferible estimar el modelo con una distribución $t$ de Student en lugar de una normal, y cómo se relaciona esta elección con la curtosis de los rendimientos financieros.
(Computacional) Replique el análisis exploratorio de este capítulo con un activo distinto al Bitcoin, por ejemplo el Ethereum: descargue la serie ETH-USD con getSymbols() del paquete quantmod (o utilice cualquier otro activo financiero líquido de su interés).
1. Calcule los rendimientos diarios como diferencias logarítmicas del precio ajustado y grafique el precio y los rendimientos. ¿Observa episodios de agrupamiento de volatilidad?
2. Grafique el histograma de los rendimientos junto con la densidad normal de igual media y varianza, calcule la curtosis y aplique la prueba de Jarque-Bera. ¿Presenta la serie colas pesadas?
3. Estime un modelo AR para la media de los rendimientos y aplique la prueba de efectos ARCH con ArchTest() del paquete FinTS. ¿Se justifica modelar la varianza condicional?
(Computacional) Continúe con la serie de rendimientos del ejercicio anterior.
1. Estime un GARCH(1,1) con media AR(1) mediante ugarchspec() y ugarchfit(), primero con distribución normal (distribution.model = "norm") y después con $t$ de Student (distribution.model = "std"). Compare ambos ajustes con infocriteria().
2. Reporte la persistencia estimada $\hat{\alpha}_1 + \hat{\beta}_1$ y calcule la vida media de los choques de volatilidad definida en el Ejercicio 2.
3. Busque el orden GARCH($p$, $q$) óptimo con $p, q \leq 3$ comparando criterios de información (puede adaptar la función Lag_Opt_GARCH() utilizada en este capítulo). ¿Confirma la afirmación de que el GARCH(1,1) suele ser suficiente en la práctica?
(Computacional) Con el modelo GARCH óptimo del ejercicio anterior:
1. Genere un pronóstico de la varianza condicional a 30 días con ugarchforecast() y grafique la volatilidad pronosticada ($\hat{\sigma}_{t+h}$). ¿Se observa la reversión a la media derivada en el Ejercicio 2?
2. Calcule el VaR al 5% para cada día del horizonte de pronóstico utilizando la varianza condicional pronosticada y la distribución estimada del modelo.
3. Compare el VaR dinámico del inciso anterior con el VaR incondicional calculado como el cuantil histórico del 5% de los rendimientos. ¿Qué ventaja ofrece el enfoque GARCH para la administración de riesgos?
(Computacional) Extienda el ejemplo multivariado del capítulo, que estimó un DCC-GARCH(1,1) para el par Bitcoin-Ethereum.
1. Estime, sobre los mismos datos, un modelo de correlación condicional constante (CCC): basta con fijar dccOrder = c(0, 0) o, equivalentemente, comparar contra la correlación no condicional. Contraste ambos modelos con los criterios de información de infocriteria() y verifique que la evidencia favorece al DCC.
2. Grafique, junto con la correlación condicional, las dos volatilidades condicionales marginales (sigma(dcc_fit)). ¿Los episodios de correlación alta coinciden con los de volatilidad alta?
3. Calcule la varianza condicional de una cartera con pesos iguales, $\sigma^2_{p,t} = 0.25(\sigma^2_{1,t} + \sigma^2_{2,t} + 2\sigma_{12,t})$, y compárela con la que se obtendría suponiendo correlación constante. ¿En qué periodos el supuesto de correlación constante subestima el riesgo de la cartera?
4. Repita el ejercicio con un par de activos de naturaleza distinta –por ejemplo, un criptoactivo y un índice accionario como el S&P 500 (^GSPC)– y compare el nivel y la estabilidad de la correlación con los del par Bitcoin-Ethereum.

Concepto	Valor
Coeficiente de \(\hat{U}^2_{t-1}\)	0.1352
Estadístico \(t\) del coeficiente	8.9595
\(R^2\) de la regresión auxiliar	0.0183
Estadístico ARCH-LM (\(T \cdot R^2\))	78.8401
Grados de libertad	1.0000
Valor \(p\)	0.0000

	Estimado	Error Est.	Estad. \(t\)
\(\alpha\) (DCC)	0.0491	0.0081	6.0412
\(\beta\) (DCC)	0.9453	0.0103	91.9684
Grados de libertad	4.0000	0.7160	5.5870