7. GARCH(p,q) y la estrategia ARIMA+GARCH sobre el S&P 500

Nivel 9: Series Temporales - ARIMA y GARCH

Pregunta natural: si aplicamos un proceso AR(p) a la varianza, ¿por qué no un MA(q)? ¿O un ARMA(p,q)? Esa es la motivación del ARCH generalizado, o GARCH, que se usa masivamente en la industria financiera porque muchos precios de activos son condicionalmente heterocedásticos.

Un GARCH(p,q) es un ARMA aplicado a la varianza de una serie: tiene un término autorregresivo y uno de media móvil. La parte AR(p) modela la varianza de los residuos al cuadrado; la parte MA(q) modela la varianza del proceso.

# GARCH(p,q):
# e[t] = sigma[t] * w[t],  con w[t] ruido blanco de media 0 y varianza 1
#
# sigma2[t] = alpha_0
#           + suma_i alpha_i * e[t-i]**2      (parte ARCH: choques pasados)
#           + suma_j beta_j  * sigma2[t-j]    (parte GARCH: varianza pasada)
#
# GARCH(1,1):
# sigma2[t] = alpha_0 + alpha_1*e[t-1]**2 + beta_1*sigma2[t-1]

Condición de estabilidad: alpha_1 + beta_1 debe ser menor que 1, o el modelo es inestable (la varianza incondicional no existe y las previsiones divergen). Con retornos diarios de índices se estiman valores como alpha_1 cerca de 0.08 y beta_1 cerca de 0.90: la suma se acerca mucho a 1, lo que significa volatilidad muy persistente. Un beta_1 alto dice que la volatilidad de hoy se parece a la de ayer; un alpha_1 alto, que reacciona con fuerza a los choques recientes.

Simulación y ajuste con el paquete arch:

np.random.seed(2)
a0, a1, b1 = 0.2, 0.5, 0.3

n = 10000
w = np.random.normal(size=n)
eps = np.zeros_like(w)
sigsq = np.zeros_like(w)

for i in range(1, n):
    sigsq[i] = a0 + a1 * (eps[i-1]**2) + b1 * sigsq[i-1]
    eps[i] = w[i] * np.sqrt(sigsq[i])

tsplot(eps, lags=30)      # ruido blanco aparente
tsplot(eps**2, lags=30)   # estructura clara en el cuadrado

from arch import arch_model     # pip install arch
print(arch_model(eps).fit(update_freq=5).summary())

Aplicación a una serie financiera. El procedimiento: recorrer combinaciones ARIMA(p,d,q), elegir los órdenes del GARCH según el ARIMA de menor AIC, ajustar el GARCH y examinar residuos y residuos al cuadrado.

Sobre retornos del S&P 500 el mejor modelo resulta ser ARIMA(4,0,3), con d cero. Los residuos parecen ruido blanco, pero el cuadrado de los residuos muestra autocorrelación clara. Ese es el disparador para el GARCH:

# la distribucion t de Student suele ajustar mejor las colas de los retornos
am = arch_model(model.resid, p=order[0], o=order[1], q=order[2],
                dist='StudentsT')
res = am.fit(update_freq=5, disp='off')
tsplot(res.resid, lags=30)
tsplot(res.resid**2, lags=30)

Si ahora tanto los residuos como sus cuadrados parecen ruido blanco discreto, el ajuste es bueno.

Estrategia ARIMA+GARCH sobre el S&P 500. Cada día n se ajustan el ARIMA optimo y después el GARCH sobre los retornos logarítmicos de los k días anteriores; se predice el retorno del día siguiente; si la predicción es positiva se va largo y si es negativa, corto, al cierre del día n; si coincide con la del día anterior no se hace nada.

windowLength = 252
foreLength = len(lrets) - windowLength

signal = 0 * lrets[-foreLength:]
for d in range(foreLength):
    TS = lrets[(1+d):(windowLength+d)]     # ventana movil
    res_tup = _get_best_model(TS)          # rejilla ARIMA por AIC, d fijado a 0
    order, model = res_tup[1], res_tup[2]
    am = arch_model(model.resid, p=order[0], o=order[1], q=order[2],
                    dist='StudentsT')
    res = am.fit(update_freq=5, disp='off')
    out = res.forecast(horizon=1, start=None, align='origin')
    signal.iloc[d] = np.sign(model.forecast(steps=1).iloc[0])   # la direccion la da la media condicional del ARIMA; el GARCH solo modela la varianza

# comparacion con comprar y mantener
rets = pd.DataFrame(index=signal.index)
rets['Buy and Hold'] = lrets[-foreLength:]
rets['Strategy'] = signal * rets['Buy and Hold']
eqCurves = rets.cumsum() + 1

El aviso del profesor es literal: el backtest no tiene en cuenta comisiones ni deslizamiento, así que el rendimiento real sería inferior al que ves aquí. Alternar largo y corto a diario paga muchísimas comisiones, y ajustar decenas de ARIMA más un GARCH por cada día hace que el indicador tarde muchísimo en generarse.

Un último detalle de ojo crítico: en el código original aparece q_ = order[1], que copia la d del ARIMA en el orden q del GARCH. Es una errata: debería ser order[2]. Además el parámetro o de arch_model es el orden de asimetría (efecto apalancamiento), no la d del ARIMA. Leer el código ajeno con esta desconfianza es parte del oficio.