Modelos de volatilidad GARCH

La volatilidad se agrupa

Los rendimientos diarios son casi impredecibles, pero su tamaño no: un gran movimiento hoy hace más probable otro grande mañana, y los días tranquilos vienen en rachas. Este agrupamiento de la volatilidad es uno de los hechos mejor establecidos en finanzas y rompe el supuesto de volatilidad constante tras el VaR simple. GARCH (heterocedasticidad condicional autorregresiva generalizada) modela directamente la varianza condicional, de modo que el riesgo y los precios de opciones pueden seguir el humor del mercado.

import numpy as np

rng = np.random.default_rng(4)
n = 3000
omega, alpha, beta = 2e-6, 0.10, 0.85
var = np.empty(n)
returns = np.empty(n)
var[0] = omega / (1 - alpha - beta)
returns[0] = np.sqrt(var[0]) * rng.normal()
for t in range(1, n):
    var[t] = omega + alpha * returns[t - 1] ** 2 + beta * var[t - 1]
    returns[t] = np.sqrt(var[t]) * rng.normal()

# Clustering: the size of a return today predicts the size tomorrow, though the sign is a coin flip
print(round(float(np.corrcoef(returns[:-1], returns[1:])[0, 1]), 3))                 # about 0
print(round(float(np.corrcoef(np.abs(returns[:-1]), np.abs(returns[1:]))[0, 1]), 3))  # clearly positive

La recursión GARCH(1,1)

El modelo tiene tres parámetros: var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]. omega es un suelo, alpha la reacción al shock de ayer y beta la persistencia de la varianza de ayer. Como cada día depende del anterior, la serie se construye con un bucle, partiendo de una varianza inicial. En acciones diarias lo típico es alpha alrededor de 0,05 a 0,1 y beta cerca de 0,85 a 0,95: un shock se desvanece despacio.

import numpy as np

def garch_variance(returns, omega, alpha, beta, initial_var):
    var = np.empty(len(returns))
    var[0] = initial_var
    for t in range(1, len(returns)):
        var[t] = omega + alpha * returns[t - 1] ** 2 + beta * var[t - 1]
    return var

returns = np.array([0.01, -0.02, 0.015, 0.0, -0.03, 0.01, 0.0, 0.0])
var = garch_variance(returns, 1e-6, 0.08, 0.9, 1e-4)
print(np.round(np.sqrt(var) * 100, 3))             # daily volatility in %: it jumps after the -3% day, then fades

Dónde se estabiliza, y el pronóstico

Si alpha + beta < 1, la varianza revierte a la media y se estabiliza en el nivel de largo plazo omega / (1 - alpha - beta); si es 1 o más, nunca lo hace. Para los pronósticos, la varianza de mañana se conoce a partir del rendimiento y la varianza de hoy. Más allá, el rendimiento al cuadrado esperado iguala la varianza, así que cada paso es omega + (alpha + beta) * anterior: el pronóstico se desliza del nivel de hoy hacia el de largo plazo, a una velocidad fijada por la persistencia alpha + beta.

def long_run_variance(omega, alpha, beta):
    if alpha + beta >= 1:
        raise ValueError("alpha + beta must be below 1")
    return omega / (1 - alpha - beta)

def garch_forecast(omega, alpha, beta, last_var, last_return, steps):
    out = [omega + alpha * last_return ** 2 + beta * last_var]
    for _ in range(steps - 1):
        out.append(omega + (alpha + beta) * out[-1])
    return out

omega, alpha, beta = 1e-6, 0.08, 0.9
print(long_run_variance(omega, alpha, beta) ** 0.5)                        # about 0.00707, 0.71% a day
forecast = garch_forecast(omega, alpha, beta, 4e-4, -0.03, 60)             # after a -3% day in a volatile spell
print([round(v ** 0.5 * 100, 3) for v in (forecast[0], forecast[9], forecast[29], forecast[59])])

Elegir los parámetros: la verosimilitud

Los parámetros no se conocen; se estiman con los datos por máxima verosimilitud. Si los rendimientos son normales con varianza var[t], la log-verosimilitud de un día es -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t]). Sumarla da un número que puntúa una elección de (omega, alpha, beta): los mejores parámetros hacen más probables los rendimientos observados. Un optimizador (por ejemplo scipy.optimize.minimize) prueba muchas ternas y se queda con la mejor.

# Score three parameter choices on data simulated with alpha = 0.10, beta = 0.85
import numpy as np

rng = np.random.default_rng(4)
n = 3000
true = (2e-6, 0.10, 0.85)
var = np.empty(n); returns = np.empty(n)
var[0] = true[0] / (1 - true[1] - true[2]); returns[0] = np.sqrt(var[0]) * rng.normal()
for t in range(1, n):
    var[t] = true[0] + true[1] * returns[t - 1] ** 2 + true[2] * var[t - 1]
    returns[t] = np.sqrt(var[t]) * rng.normal()

def loglik(params):
    omega, alpha, beta = params
    v = np.empty(n)
    v[0] = returns.var()
    for t in range(1, n):
        v[t] = omega + alpha * returns[t - 1] ** 2 + beta * v[t - 1]
    return float(np.sum(-0.5 * (np.log(2 * np.pi * v) + returns ** 2 / v)))

for params in ((2e-6, 0.10, 0.85), (2e-6, 0.30, 0.50), (1e-5, 0.0, 0.0)):
    print(params, round(loglik(params), 1))                       # the true parameters score highest

• La volatilidad se agrupa: GARCH(1,1) construye la varianza con un bucle, var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1].
• Con alpha + beta < 1 revierte a omega / (1 - alpha - beta); los pronósticos usan omega + (alpha + beta) * anterior.
• Los parámetros se eligen por máxima verosimilitud, sum(-0.5 * (log(2 * pi * var) + r**2 / var)); use shocks de colas gruesas para las colas.

A continuación: el movimiento browniano, el paseo aleatorio en tiempo continuo que sustenta los modelos de valoración de opciones.

Lab: qué se califica

• Ejecutar la recursión con garch_variance25 pts
• Hallar dónde se estabiliza con long_run_variance25 pts
• Pronosticar los días siguientes con garch_forecast25 pts
• Puntuar el ajuste con garch_loglik25 pts