Modelos de volatilidade GARCH

A volatilidade se agrupa

Os retornos diários são quase imprevisíveis, mas seu tamanho não: um grande movimento hoje torna outro grande amanhã mais provável, e os dias calmos vêm em sequência. Esse agrupamento da volatilidade é um dos fatos mais bem estabelecidos em finanças e quebra a suposição de volatilidade constante por trás do VaR simples. O GARCH (heterocedasticidade condicional autorregressiva generalizada) modela diretamente a variância condicional, de modo que o risco e os preços de opções acompanham o humor do 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

A recursão GARCH(1,1)

O modelo tem três parâmetros: var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]. omega é um piso, alpha a reação ao choque de ontem e beta a persistência da variância de ontem. Como cada dia depende do anterior, a série é construída com um laço, a partir de uma variância inicial. Em ações diárias o típico é alpha em torno de 0,05 a 0,1 e beta perto de 0,85 a 0,95: um choque se dissipa devagar.

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

Onde se estabiliza, e a previsão

Se alpha + beta < 1, a variância reverte à média e se estabiliza no nível de longo prazo omega / (1 - alpha - beta); se for 1 ou mais, nunca o faz. Nas previsões, a variância de amanhã é conhecida a partir do retorno e da variância de hoje. Além disso, o retorno ao quadrado esperado iguala a variância, então cada passo é omega + (alpha + beta) * anterior: a previsão desliza do nível de hoje para o de longo prazo, com velocidade dada pela persistência 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])])

Escolher os parâmetros: a verossimilhança

Os parâmetros não são conhecidos; são estimados nos dados por máxima verossimilhança. Se os retornos são normais com variância var[t], a log-verossimilhança de um dia é -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t]). Somá-la dá um número que avalia uma escolha de (omega, alpha, beta): os melhores parâmetros tornam os retornos observados mais prováveis. Um otimizador (por exemplo scipy.optimize.minimize) tenta muitas trincas e fica com a melhor.

# 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

• A volatilidade se agrupa: o GARCH(1,1) constrói a variância com um laço, var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1].
• Com alpha + beta < 1 reverte a omega / (1 - alpha - beta); as previsões usam omega + (alpha + beta) * anterior.
• Os parâmetros são escolhidos por máxima verossimilhança, sum(-0.5 * (log(2 * pi * var) + r**2 / var)); use choques de caudas pesadas para as caudas.

A seguir: o movimento browniano, o passeio aleatório em tempo contínuo que sustenta os modelos de precificação de opções.

Lab: o que é avaliado

• Rodar a recursão com garch_variance25 pts
• Achar onde ela se estabiliza com long_run_variance25 pts
• Prever os dias seguintes com garch_forecast25 pts
• Pontuar o ajuste com garch_loglik25 pts