GARCH 波动率模型

波动率会聚集

日收益率几乎无法预测，但它们的幅度可以：今天的大幅波动使明天再出现大幅波动的可能性上升，平静的日子也成串出现。这种波动率聚集是金融中最确凿的事实之一，它推翻了简单 VaR 背后波动率恒定的假设。GARCH（广义自回归条件异方差）直接对条件方差建模，使风险和期权价格能跟随市场的情绪。

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

GARCH(1,1) 递推

模型有三个参数：var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]。omega 是下限，alpha 是对昨日冲击的反应，beta 是昨日方差的持续性。由于每一天都依赖前一天，所以要从初始方差出发用循环构建这个序列。日度股票的典型值是 alpha 约 0.05 到 0.1，beta 约 0.85 到 0.95：冲击衰减得很慢。

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

稳定水平与预测

若 alpha + beta < 1，方差均值回归，稳定在长期水平 omega / (1 - alpha - beta)；若等于或大于 1 则永远不会稳定。预测时，明天的方差由今天的收益和方差决定；再往后，预期平方收益等于方差，所以每一步是 omega + (alpha + beta) * 前一步：预测从今天的水平向长期水平滑行，速度由持续性 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])])

选择参数：似然函数

参数是未知的，要用数据通过最大似然来估计。若收益率服从方差为 var[t] 的正态分布，则一天的对数似然是 -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t])。把它们加总，得到给 (omega, alpha, beta) 的一组取值打分的一个数：最佳参数使观测到的收益率最可能出现。优化器（例如 scipy.optimize.minimize）会尝试很多组并保留最高的。

# 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

• 波动率会聚集：GARCH(1,1) 用循环构建方差，var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]。
• 当 alpha + beta < 1 时回归到 omega / (1 - alpha - beta)；预测使用 omega + (alpha + beta) * 前一步。
• 参数通过最大似然 sum(-0.5 * (log(2 * pi * var) + r**2 / var)) 选取；对尾部使用厚尾冲击。

下一节：布朗运动，期权定价模型背后的连续时间随机游走。

实验：评分内容

• 用 garch_variance 运行递推25 分
• 用 long_run_variance 找到它的稳定水平25 分
• 用 garch_forecast 预测未来几天25 分
• 用 garch_loglik 给拟合打分25 分