El modelo GARCH: cómo predecir la volatilidad (con Python)

Si miras la gráfica de cualquier índice bursátil, verás que los días tranquilos vienen en rachas y los días agitados, también. Tras un desplome, el mercado no vuelve a la calma de un día para otro: sigue nervioso durante semanas. A esto se le llama agrupamiento de la volatilidad, y es uno de los hechos más sólidos de las finanzas.

El problema es que casi todo lo que hemos visto hasta ahora (el Monte Carlo, Black-Scholes o el VaR paramétrico) supone una volatilidad constante. El modelo GARCH es la forma estándar de dejar que la volatilidad cambie con el tiempo, y de preverla.

La idea: la volatilidad de mañana depende de la de hoy

Robert Engle propuso en 1982 el modelo ARCH, que le valió el Nobel de 2003, y Tim Bollerslev lo generalizó en 1986 al GARCH. La versión más usada, el GARCH(1,1), dice que la varianza de mañana es una mezcla de tres cosas:

\[ \sigma_t^2 = \omega + \alpha\, r_{t-1}^2 + \beta\, \sigma_{t-1}^2 \]

• \( \omega \): un nivel base, que fija la varianza de largo plazo.
• \( \alpha\, r_{t-1}^2 \): la reacción al movimiento de ayer. Un día de −5 % dispara la volatilidad esperada.
• \( \beta\, \sigma_{t-1}^2 \): la memoria. La volatilidad de ayer se arrastra hoy.

Dos números resumen el modelo. La persistencia \( \alpha + \beta \), que en acciones suele estar entre 0,95 y 0,99: los shocks de volatilidad tardan mucho en disiparse. Y la varianza de largo plazo, hacia la que la volatilidad siempre acaba volviendo:

\[ \sigma_L^2 = \frac{\omega}{1 - \alpha - \beta} \]

Simular un GARCH para entenderlo

La mejor forma de entender un modelo es simularlo. Usaremos parámetros típicos de un índice de acciones:

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm

def simular_garch(n, omega, alpha, beta, seed=0):
    rng = np.random.default_rng(seed)
    r = np.zeros(n)
    var = np.zeros(n)
    var[0] = omega / (1 - alpha - beta)              # arranca en la varianza de largo plazo
    for t in range(n):
        if t > 0:
            var[t] = omega + alpha * r[t - 1] ** 2 + beta * var[t - 1]
        r[t] = np.sqrt(var[t]) * rng.standard_normal()
    return r, np.sqrt(var)

omega, alpha, beta = 2e-6, 0.09, 0.89
r, vol = simular_garch(4000, omega, alpha, beta, seed=1)
print(f"Volatilidad anual de largo plazo: {np.sqrt(omega / (1 - alpha - beta) * 252):.1%}")
print(f"Volatilidad condicional: mínimo {vol.min() * np.sqrt(252):.1%}, máximo {vol.max() * np.sqrt(252):.1%}")

Volatilidad anual de largo plazo: 15.9% Volatilidad condicional: mínimo 8.6%, máximo 41.9%

Con una volatilidad media del 16 %, el modelo genera por sí solo periodos de calma con un 9 % y tormentas por encima del 40 %, sin que haya que programar ninguna crisis:

La huella del agrupamiento

Hay una forma sencilla de detectarlo en cualquier serie: las rentabilidades casi no se autocorrelacionan (no puedes predecir si mañana sube o baja), pero sus cuadrados sí (puedes predecir si mañana se moverá mucho o poco).

def autocorr(x, k=1):
    x = x - x.mean()
    return (x[:-k] * x[k:]).sum() / (x * x).sum()

print(f"Autocorrelación de r:  {autocorr(r):.3f}")
print(f"Autocorrelación de r²: {autocorr(r ** 2):.3f}")

Autocorrelación de r: 0.015 Autocorrelación de r²: 0.208

Haz la prueba con los precios reales de cualquier acción o índice: el resultado se parece mucho.

Estimar el GARCH por máxima verosimilitud

En la práctica tenemos las rentabilidades y queremos los parámetros. Si las rentabilidades son normales con varianza \( \sigma_t^2 \), la log-verosimilitud es:

\[ \ln \mathcal{L} = -\frac{1}{2} \sum_t \left[ \ln\!\left(2\pi\sigma_t^2\right) + \frac{r_t^2}{\sigma_t^2} \right] \]

Basta con calcular la varianza de cada día con la fórmula del GARCH y buscar los \( \omega, \alpha, \beta \) que maximizan esa expresión, respetando que sean positivos y que \( \alpha + \beta < 1 \):

def varianzas(params, r):
    omega, alpha, beta = params
    var = np.empty_like(r)
    var[0] = r.var()
    for t in range(1, len(r)):
        var[t] = omega + alpha * r[t - 1] ** 2 + beta * var[t - 1]
    return var

def menos_log_verosimilitud(params, r):
    omega, alpha, beta = params
    if omega <= 0 or alpha < 0 or beta < 0 or alpha + beta >= 1:
        return 1e10                                  # fuera de la región válida
    var = varianzas(params, r)
    return 0.5 * np.sum(np.log(2 * np.pi * var) + r ** 2 / var)

inicio = [r.var() * 0.05, 0.05, 0.90]
res = minimize(menos_log_verosimilitud, inicio, args=(r,), method="Nelder-Mead",
               options={"maxiter": 5000, "xatol": 1e-10, "fatol": 1e-10})
w, a, b = res.x
print(f"omega = {w:.2e}   alpha = {a:.3f}   beta = {b:.3f}   (reales: 2.00e-06, 0.090, 0.890)")
print(f"Persistencia alpha + beta = {a + b:.3f}")

omega = 2.19e-06 alpha = 0.097 beta = 0.882 (reales: 2.00e-06, 0.090, 0.890) Persistencia alpha + beta = 0.979

El estimador recupera los parámetros reales con bastante precisión. Con datos reales se usa la librería arch de Python, que hace esto mismo con errores estándar, otras distribuciones (como la t de Student) y variantes del modelo, pero es importante entender qué hay debajo.

Prever la volatilidad

La gran utilidad del GARCH es que da una previsión de la volatilidad a cualquier plazo. La varianza esperada dentro de \( h \) días vuelve hacia la de largo plazo a un ritmo marcado por la persistencia:

\[ \mathbb{E}\left[\sigma_{t+h}^2\right] = \sigma_L^2 + (\alpha + \beta)^{h-1}\left(\sigma_{t+1}^2 - \sigma_L^2\right) \]

var_hoy = varianzas(res.x, r)[-1]
var_manana = w + a * r[-1] ** 2 + b * var_hoy
var_largo = w / (1 - a - b)
for h in [1, 5, 21, 63, 252]:
    var_h = var_largo + (a + b) ** (h - 1) * (var_manana - var_largo)
    print(f"Previsión a {h:>3} días: volatilidad anual {np.sqrt(var_h * 252):.1%}")

print(f"Vida media de un shock de volatilidad: {np.log(0.5) / np.log(a + b):.0f} días")

Previsión a 1 días: volatilidad anual 26.8% Previsión a 5 días: volatilidad anual 26.1% Previsión a 21 días: volatilidad anual 23.7% Previsión a 63 días: volatilidad anual 19.7% Previsión a 252 días: volatilidad anual 16.4% Vida media de un shock de volatilidad: 33 días

Al final de la serie simulada el mercado está agitado (27 % anual), y el modelo prevé que la calma vuelva poco a poco: la mitad del exceso de volatilidad tarda unos 33 días de mercado, casi mes y medio, en disiparse.

Para qué se usa el GARCH

• VaR que reacciona. Un VaR con volatilidad histórica de varios años ignora que hoy el mercado está nervioso. Con la volatilidad GARCH de mañana, el VaR sube justo cuando debe:

sigma_manana = np.sqrt(var_manana)
print(f"VaR 99 % de mañana con GARCH:                {-norm.ppf(0.01) * sigma_manana:.2%}")
print(f"VaR 99 % con la volatilidad de toda la serie: {-norm.ppf(0.01) * r.std():.2%}")

VaR 99 % de mañana con GARCH: 3.92% VaR 99 % con la volatilidad de toda la serie: 2.33%

• Tamaño de las posiciones. Muchos fondos ajustan el tamaño de cada posición para que su riesgo sea constante: menos exposición cuando el GARCH prevé tormenta (volatility targeting).
• Valoración de opciones. La volatilidad que entra en la fórmula de una opción es la esperada durante su vida; el GARCH da una estimación de esa estructura temporal.
• Backtesting del riesgo. Las excepciones del VaR dejan de llegar en racha, que era el síntoma que vimos en el artículo del VaR.

Límites y extensiones

• Asimetría. En acciones, las caídas suben la volatilidad más que las subidas del mismo tamaño (el efecto apalancamiento). Los modelos GJR-GARCH y EGARCH lo recogen.
• Colas gruesas. Incluso condicionando a la volatilidad, las rentabilidades tienen colas más gruesas que la normal. Estimar con una t de Student suele mejorar el ajuste.
• Varios activos. Para carteras hacen falta correlaciones que también cambian con el tiempo (modelos DCC), y ahí la complejidad crece rápido.

Resumen

• La volatilidad se agrupa: las rentabilidades no se predicen, pero su tamaño sí.
• El GARCH(1,1) combina un nivel base, la reacción al último shock y la memoria de la volatilidad anterior.
• Se estima por máxima verosimilitud y da previsiones que vuelven a la media de largo plazo.
• Sirve para un VaR que reacciona, para dimensionar posiciones y para valorar opciones.

Ya sabes modelar el riesgo de un activo en el tiempo. El siguiente paso de la ruta es una estrategia que explota la relación entre dos activos: pairs trading y cointegración.

Fuentes

• R. F. Engle (1982), «Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation», Econometrica.
• T. Bollerslev (1986), «Generalized Autoregressive Conditional Heteroskedasticity», Journal of Econometrics.