量化交易中文教材

第 23 章 波动率与相关系数的估计

对应原书:Hull《Options, Futures, and Other Derivatives》第 9 版(Global Edition)第 23 章 Estimating Volatilities and Correlations。

与第 06 册的分工:第 06 册第 03a 章(Tsay)从时间序列建模的角度讲 ARCH/GARCH:ARCH 效应检验、ARMA 表示、峰度、IGARCH、GARCH-M 和模型定阶。本章从风险管理和衍生品定价的角度讲同一组模型:怎样每天更新一个波动率数字、怎样用它预测期权存续期内的平均波动率、怎样一致地更新整个协方差矩阵来算 VaR。两章的公式可以互相对照:Hull 的 \(u_{n-1}\) 就是 Tsay 的 \(a_{t-1}\),Hull 的 \(\sigma_n^2\) 就是 Tsay 的 \(\sigma_t^2\)。第 05 册第 17c 章则从计量经济学角度(持续性、多步预测与推断)讲同一组模型。

学习目标

读完本章,你应当能够:

  1. 说清楚日波动率监控中的三处简化(简单收益、零均值、除以 \(m\)),以及为什么要给近期观测更大的权重。
  2. 用 EWMA 和 GARCH(1,1) 做单步更新,会把 GARCH 参数换算成长期方差 \(V_L\)、衰减率和均值回归速度,并知道 EWMA 是 GARCH 的特例。
  3. 写出 GARCH(1,1) 的对数似然目标函数,用数值优化估计参数,知道方差锁定(variance targeting)和多初值的作用,并用 \(u_i^2/\sigma_i^2\) 的 Ljung–Box 统计量检验模型。
  4. 用 \(E[\sigma_{n+t}^2]=V_L+(\alpha+\beta)^t(\sigma_n^2-V_L)\) 预测未来方差,算出给 \(T\) 天期权定价的波动率期限结构,以及瞬时波动率变化对不同期限的传导。
  5. 用 EWMA/GARCH 更新协方差与相关系数,理解协方差矩阵必须半正定,并知道"方差和协方差必须用同一方法计算"这条实务纪律。

读前导读

这一章在解决什么问题。 第 22 章的 VaR 公式里有两个输入被当作已知:每个资产的日波动率和资产间的相关系数。本章讲这两个数每天怎么算。你在 CFA 里算过历史波动率:取一段收益,求样本标准差。它的问题是对每一天一视同仁,而市场明显有"波动率聚集":大波动后面跟着大波动,平静后面跟着平静。2008 年 9 月雷曼倒闭后的那几天,用过去两年等权算出的波动率显然低估了当下的风险。

本章的三种模型给出了递进的答案。EWMA 让越近的数据权重越大,权重按固定比例逐日衰减;GARCH(1,1) 在 EWMA 的基础上再加一个"长期平均波动率"的锚,让波动率在偏离后有回归的趋势。你在 CFA/FRM 里见过 EWMA 和 GARCH 的更新公式,本章补上三块:参数怎样用极大似然法从数据中估出来;GARCH 怎样预测未来几个月的平均波动率,从而给期权定价;协方差矩阵怎样一致地更新,才能保证算出的组合方差不为负。

需要先想起来的数学。

  • 几何级数。 \(1+\lambda+\lambda^2+\cdots=\frac{1}{1-\lambda}\)(\(|\lambda|<1\))。这是 EWMA 权重之和为 1 的原因,也和永续年金现值公式是同一个式子。例:\(\lambda=0.94\) 时,\((1-\lambda)\sum\lambda^{i-1}=0.06\times\frac{1}{0.06}=1\)。见 第 00 册第 04 章 级数与收敛。
  • 对数与求导求极值。 极大似然要对似然函数取对数(把连乘变成连加),再对参数求导令其为零。用到 \(\frac{d}{dv}\ln v=\frac1v\)、\(\frac{d}{dv}\frac1v=-\frac1{v^2}\)。见 第 00 册第 02 章 导数与泰勒展开。
  • 条件期望。 \(E[u_n^2\mid\text{第 }n-1\text{ 天的信息}]=\sigma_n^2\):在已知昨天信息的前提下,今天收益平方的期望就是模型给出的方差。"条件"二字的意思是"站在某一天、只用当时已知的信息"。见 第 00 册第 07 章 概率中的分析工具。
  • 定积分求平均值。 一个函数在 \([0,T]\) 上的平均值是 \(\frac1T\int_0^Tf(t)\,dt\)。指数函数的积分 \(\int_0^Te^{-at}dt=\frac{1-e^{-aT}}{a}\)。这用在 23.6.2 节求期权存续期内的平均方差。见 第 00 册第 03 章 积分。
  • 半正定矩阵。 对任意权重向量 \(w\) 都有 \(w^{\mathsf T}\Omega w\ge0\),即"任何组合的方差都不为负"。判别方法之一是所有特征值都 \(\ge0\)。见 第 00 册第 06 章 线性代数速成。

怎么读这一章。 23.1–23.4 节是核心,要能手算 EWMA 和 GARCH 的单步更新(例 23.1、23.2),并理解"GARCH = EWMA + 一份长期方差权重"。23.3.3 节的连续时间方程第一次可以只看结论"GARCH 有均值回归"。23.5 节极大似然:23.5.2 的推导要跟着走一遍,它解释了为什么除以 \(m\);23.5.3–23.5.5 的数值技巧和 Ljung–Box 检验可以先浏览。23.6 节预测公式 (23.13) 是期权定价最需要的结果,推导值得读;23.6.3 的传导公式可只看表 23.4 的结论。23.7 节的协方差更新和半正定条件、23.8 节的四指数例子直接连回第 22 章,建议必读。


23.0 为什么要估计"当前"波动率

BSM 模型(本册第 15 章)把波动率当常数,但真实市场的波动率时高时低,而且不能直接观测。估计波动率有两个用途,它们关心的时间尺度不同:

  1. 风险度量。用模型构建法算 VaR(本册第 22 章)时,持有期通常只有 1 天或 10 天,我们关心的是当前的波动率和相关系数。
  2. 衍生品定价。给一个 3 个月期权定价,需要的是期权整个存续期内的平均波动率,这就要预测波动率未来怎样演变。

本章的模型(EWMA、ARCH、GARCH)都在做同一件事:用过去的收益平方的加权平均跟踪时变的方差,近期的观测权重大。区别在于权重怎么分配,以及要不要给一个长期平均水平留一份权重。


23.1 估计波动率

23.1.1 记号

本章统一使用下列记号:

  • \(\sigma_n\):在第 \(n-1\) 天末估计的第 \(n\) 天的波动率;\(\sigma_n^2\) 称为方差率(variance rate)。
  • \(S_i\):第 \(i\) 天末市场变量的值。
  • \(u_i\):第 \(i\) 天的收益。

"在第 \(n-1\) 天末估计第 \(n\) 天"这个约定非常重要:它保证了估计只用到当时已知的信息。做回测时,第 \(n\) 天的仓位只能用 \(\sigma_n\)(而不能用包含第 \(n\) 天收益的估计),否则就有前视偏差。

23.1.2 从标准估计到监控公式

第 15.4 节给出过波动率的标准无偏估计。用最近 \(m\) 天的连续复利收益 \(u_i=\ln(S_i/S_{i-1})\):

\[\sigma_n^2=\frac{1}{m-1}\sum_{i=1}^{m}(u_{n-i}-\bar u)^2,\qquad \bar u=\frac1m\sum_{i=1}^m u_{n-i}\tag{23.1}\]

为了监控日波动率,风险管理中作三处简化:

  1. 把 \(u_i\) 改成百分比变化(简单收益):
    \[u_i=\frac{S_i-S_{i-1}}{S_{i-1}}\tag{23.2}\]
    这与第 22 章 VaR 中波动率的定义一致。
  2. 令 \(\bar u=0\)。一天的期望收益相对于日标准差非常小,忽略它几乎没有影响。
  3. 用 \(m\) 代替 \(m-1\)。这把无偏估计换成了极大似然估计(23.5 节会证明)。

三处修改几乎不改变数值,得到简洁的形式:

\[\sigma_n^2=\frac1m\sum_{i=1}^m u_{n-i}^2\tag{23.3}\]

注意下标含义:本章的下标 \(i\) 指同一变量的不同日期;而第 22 章 \(x_i\) 的下标指同一天的不同变量。

23.1.3 加权方案

式 (23.3) 对过去 \(m\) 个 \(u^2\) 一视同仁。如果目标是估计当前的波动率,越近的观测应当越重要:

\[\sigma_n^2=\sum_{i=1}^m\alpha_i u_{n-i}^2,\qquad \alpha_i>0,\quad \alpha_i<\alpha_j\ (i>j),\quad \sum_{i=1}^m\alpha_i=1\tag{23.4}\]

ARCH(\(m\)) 模型(autoregressive conditional heteroscedasticity,自回归条件异方差;Engle 1982)在此基础上再给一个长期平均方差率 \(V_L\) 一份权重 \(\gamma\):

\[\sigma_n^2=\gamma V_L+\sum_{i=1}^m\alpha_i u_{n-i}^2,\qquad \gamma+\sum_{i=1}^m\alpha_i=1\tag{23.5}\]

令 \(\omega=\gamma V_L\),写成估计时用的形式:

\[\sigma_n^2=\omega+\sum_{i=1}^m\alpha_i u_{n-i}^2\tag{23.6}\]

白话解释:\(\gamma V_L\) 这一项可以看作"长期平均水平的锚"。纯粹的加权平均 (23.4) 只看最近 \(m\) 天,近期平静时估计就很低,近期动荡时估计就很高。ARCH 规定每天的估计都必须保留一份 \(\gamma\) 的权重给长期平均方差 \(V_L\),这样估计不会离长期水平太远。这和信用评级中的"穿越周期"(through-the-cycle)思路相近:既看当前,也要锚定长期均值。 写成 \(\omega=\gamma V_L\) 只是为了估计方便:回归时 \(\omega\) 是一个常数项,估出 \(\omega\) 和各 \(\alpha_i\) 后,再由 \(\gamma=1-\sum\alpha_i\) 反推 \(V_L=\omega/\gamma\)。

"条件异方差"这个名字的含义是:给定过去的信息,第 \(n\) 天收益的方差 \(\sigma_n^2\) 是变化的。ARCH 检验与 ARCH 模型的统计性质见第 06 册第 03a 章 3.3–3.4 节。


23.2 指数加权移动平均模型(EWMA)

23.2.1 递推公式

EWMA(exponentially weighted moving average)是式 (23.4) 的特例:权重按指数递减,\(\alpha_{i+1}=\lambda\alpha_i\),\(0<\lambda<1\)。它可以化成一个极简的更新公式:

\[\sigma_n^2=\lambda\sigma_{n-1}^2+(1-\lambda)u_{n-1}^2\tag{23.7}\]

读法是:今天对明天方差的估计 = \(\lambda\) × 昨天的估计 + \((1-\lambda)\) × 今天的收益平方。

23.2.2 为什么是指数权重

把式 (23.7) 中的 \(\sigma_{n-1}^2\) 再用同一公式展开:

\[\sigma_n^2=\lambda\big[\lambda\sigma_{n-2}^2+(1-\lambda)u_{n-2}^2\big]+(1-\lambda)u_{n-1}^2=(1-\lambda)\big(u_{n-1}^2+\lambda u_{n-2}^2\big)+\lambda^2\sigma_{n-2}^2\]

继续代入 \(m\) 次:

\[\sigma_n^2=(1-\lambda)\sum_{i=1}^m\lambda^{i-1}u_{n-i}^2+\lambda^m\sigma_{n-m}^2\]

\(m\) 很大时末项可以忽略。所以 EWMA 正是式 (23.4) 中取 \(\alpha_i=(1-\lambda)\lambda^{i-1}\) 的情形:每个权重是前一个的 \(\lambda\) 倍,权重之和 \((1-\lambda)\sum\lambda^{i-1}\to1\)。

两个常用的换算:

  • 半衰期(half-life)\(h\):权重衰减到一半所需的天数,\(\lambda^h=0.5\),即 \(h=\ln0.5/\ln\lambda\)。\(\lambda=0.94\) 对应约 11.2 天,\(\lambda=0.97\) 对应约 22.8 天。商业风险模型常用"半衰期 \(h\) 天"来描述衰减,它与 \(\lambda=0.5^{1/h}\) 是同一回事。
  • 有效窗口:权重的"平均年龄"约为 \(1/(1-\lambda)\) 天,\(\lambda=0.94\) 时约 17 天。

推导拆解:两个换算各一步。 权重之和:\((1-\lambda)\sum_{i=1}^{\infty}\lambda^{i-1}=(1-\lambda)\cdot\frac{1}{1-\lambda}=1\),用的是几何级数求和,与永续年金 \(PV=\frac{C}{r}\) 的推导相同。 半衰期:要找 \(h\) 使第 \(h+1\) 天前的权重是最新权重的一半,\(\lambda^h=0.5\),两边取对数得 \(h\ln\lambda=\ln0.5\)。\(\lambda=0.94\):\(\ln0.94=-0.0619\),\(h=0.693/0.0619=11.2\)。 平均年龄:\(\sum_{i\ge1}i\cdot(1-\lambda)\lambda^{i-1}=\frac{1}{1-\lambda}\),这是几何分布的期望。它和债券的麦考利久期是同一个概念:以权重(现值占比)加权的平均时间。所以"EWMA 的有效窗口"可以理解为"权重的久期"。

23.2.3 例题

例 23.1 \(\lambda=0.90\),第 \(n-1\) 天的波动率估计为每天 1%,当天市场变量上涨 2%。于是 \(\sigma_{n-1}^2=0.0001\),\(u_{n-1}^2=0.0004\):

\[\sigma_n^2=0.9\times0.0001+0.1\times0.0004=0.00013,\qquad \sigma_n=\sqrt{0.00013}=1.14\%\text{/天}\]

直觉:昨天的估计说 \(u_{n-1}^2\) 的期望是 \(0.0001\),实际却是 \(0.0004\),于是估计上调;如果实际收益平方小于预期,估计就下调。

23.2.4 EWMA 的特点

  • 存储极小:只需保存当前的方差估计和最近一次价格。新价格到来,算出收益,用式 (23.7) 更新,旧数据即可丢弃。这让 EWMA 特别适合流式计算和成千上万个风险因子的日常更新。
  • 跟踪变化:前一天波动大,\(u_{n-1}^2\) 就大,估计随之上升。
  • \(\lambda\) 控制反应速度:\(\lambda\) 小,最新观测权重大,估计本身跳得厉害;\(\lambda\) 接近 1,估计平滑但对新信息反应慢。

RiskMetrics(JPMorgan 1994 年公开)用 \(\lambda=0.94\) 更新日波动率。选择依据是:在多种市场变量上,这个 \(\lambda\) 给出的方差预测最接近已实现方差率(realized variance,这里定义为其后 25 天 \(u_i^2\) 的等权平均;原书习题 23.19 让你复现这一做法)。


23.3 GARCH(1,1) 模型

23.3.1 定义

GARCH(1,1) 由 Bollerslev(1986)提出。它与 EWMA 的关系,就像式 (23.5) 与式 (23.4) 的关系——多给长期平均方差率 \(V_L\) 一份权重:

\[\sigma_n^2=\gamma V_L+\alpha u_{n-1}^2+\beta\sigma_{n-1}^2,\qquad \gamma+\alpha+\beta=1\tag{23.8}\]
  • EWMA 是它的特例:\(\gamma=0\),\(\alpha=1-\lambda\),\(\beta=\lambda\)。
  • "(1,1)"表示 \(\sigma_n^2\) 依赖最近 1 个 \(u^2\) 和最近 1 个方差估计。一般的 GARCH(\(p,q\)) 用最近 \(p\) 个 \(u^2\) 和最近 \(q\) 个方差估计(下标约定在不同书中可能对调,见第 06 册第 03a 章)。实务中 GARCH(1,1) 最常用。
  • 对股票,波动率在下跌时倾向于上升(杠杆效应,见第 20 章),这时可以用让 \(\sigma_n\) 依赖 \(u_{n-1}\) 符号的非对称变体,如 EGARCH(Nelson 1990)及 Engle & Ng(1993)讨论的模型。

估计时令 \(\omega=\gamma V_L\):

\[\sigma_n^2=\omega+\alpha u_{n-1}^2+\beta\sigma_{n-1}^2\tag{23.9}\]

估出 \(\omega,\alpha,\beta\) 之后,\(\gamma=1-\alpha-\beta\),\(V_L=\omega/\gamma\)。平稳性条件是 \(\alpha+\beta<1\);否则长期方差的权重 \(\gamma\) 为负,模型没有稳定的长期水平。

例 23.2 某日数据估出

\[\sigma_n^2=0.000002+0.13u_{n-1}^2+0.86\sigma_{n-1}^2\]

则 \(\gamma=1-0.13-0.86=0.01\),\(V_L=0.000002/0.01=0.0002\),长期日波动率 \(\sqrt{0.0002}=1.4\%\)。若第 \(n-1\) 天的估计为 \(\sigma_{n-1}=1.6\%\)(\(\sigma_{n-1}^2=0.000256\)),当天下跌 1%(\(u_{n-1}^2=0.0001\)):

\[\sigma_n^2=0.000002+0.13\times0.0001+0.86\times0.000256=0.00023516,\qquad \sigma_n=1.53\%\text{/天}\]

23.3.2 权重结构

把式 (23.9) 中的 \(\sigma_{n-1}^2,\sigma_{n-2}^2,\dots\) 反复代入,会发现 \(u_{n-i}^2\) 的权重是 \(\alpha\beta^{i-1}\),以速率 \(\beta\) 指数衰减。所以 \(\beta\) 可以理解为衰减率(decay rate),作用和 EWMA 的 \(\lambda\) 相同。例如 \(\beta=0.9\) 时,\(u_{n-2}^2\) 的重要性是 \(u_{n-1}^2\) 的 90%,\(u_{n-3}^2\) 是 81%。

一句话概括:GARCH(1,1) = 指数衰减的权重 + 一份给长期平均方差的权重。

推导拆解:权重 \(\alpha\beta^{i-1}\) 是怎么来的。 把 \(\sigma_{n-1}^2=\omega+\alpha u_{n-2}^2+\beta\sigma_{n-2}^2\) 代入 (23.9):\(\sigma_n^2=\omega+\alpha u_{n-1}^2+\beta\omega+\alpha\beta u_{n-2}^2+\beta^2\sigma_{n-2}^2\)。 再代入 \(\sigma_{n-2}^2\):多出 \(\beta^2\omega+\alpha\beta^2u_{n-3}^2+\beta^3\sigma_{n-3}^2\)。 规律是:\(u_{n-i}^2\) 的系数为 \(\alpha\beta^{i-1}\);常数项累积为 \(\omega(1+\beta+\beta^2+\cdots)\to\frac{\omega}{1-\beta}\);末项 \(\beta^k\sigma_{n-k}^2\) 随 \(k\) 增大趋于 0。 所以 GARCH 的结构和 EWMA 展开后一模一样,只是多了一个常数项。这个常数项就是"长期方差那一份权重"在展开后的样子。

23.3.3 均值回归

GARCH(1,1) 的一个关键性质是方差均值回归(mean reversion)。可以证明,在以天为时间单位时,它对应于方差 \(V\) 的连续时间过程

\[dV=a(V_L-V)\,dt+\xi V\,dz,\qquad a=1-\alpha-\beta,\quad \xi=\alpha\sqrt2\]

推导思路(原书习题 23.14):由式 (23.8),

\[\sigma_n^2-\sigma_{n-1}^2=(1-\alpha-\beta)(V_L-\sigma_{n-1}^2)+\alpha\,(u_{n-1}^2-\sigma_{n-1}^2)\]

第一项是确定性的漂移:方差高于 \(V_L\) 时被拉下来,低于时被推上去,速度为 \(a=1-\alpha-\beta\)。第二项是随机冲击:给定 \(\sigma_{n-1}\),\(u_{n-1}\) 条件正态,\(u_{n-1}^2\) 的均值是 \(\sigma_{n-1}^2\),方差是 \(E[u^4]-\sigma^4=3\sigma^4-\sigma^4=2\sigma_{n-1}^4\),所以 \(\alpha(u_{n-1}^2-\sigma_{n-1}^2)\) 的标准差是 \(\alpha\sqrt2\,\sigma_{n-1}^2\),这正是 \(\xi V\,dz\) 项。

这是一个随机波动率模型,第 27a 章会进一步讨论这类模型在期权定价中的作用。

金融直觉:\(dV=a(V_L-V)dt+\xi V\,dz\) 的漂移项和利率模型中的 Vasicek 模型 \(dr=a(b-r)dt+\sigma dz\) 结构相同:变量偏离长期水平时,被一股与偏离量成正比的力拉回。\(a\) 是回归速度,S&P 500 的估计中 \(a=1-0.9935=0.0065\)/天,对应的"偏离半衰期"约为 \(\ln2/0.0065\approx107\) 个交易日,约 5 个月。也就是说,一次波动率冲击大约 5 个月消退一半,这与经验上"危机后的高波动通常持续数月"吻合。 符号说明:\(dz\) 是维纳过程的增量,表示一小段时间内的标准化随机冲击,在 \(\Delta t\) 内约为 \(\epsilon\sqrt{\Delta t}\)(\(\epsilon\) 为标准正态)。


23.4 模型选择

实践中方差率确实有均值回归倾向:高波动不会永远持续,低波动迟早会被打破。GARCH(1,1) 有均值回归,EWMA 没有,所以 GARCH(1,1) 在理论上更合理。

两者的关系是:\(\omega=0\) 时 GARCH(1,1) 退化为 EWMA。如果数据拟合出的最优 \(\omega\) 为负,GARCH(1,1) 就不稳定,这时应改用 EWMA。实务中两者并存:EWMA 简单稳健、参数少(只有一个 \(\lambda\),且可以对所有风险因子统一取值),常用于大规模协方差矩阵;GARCH 用于需要预测较长期限波动率的场合,例如给期权定价。


23.5 极大似然估计

23.5.1 基本思想

极大似然法(maximum likelihood method)选择参数,使已观测到的数据出现的概率(似然,likelihood)最大。一般理论见第 03 册(Wasserman)第 09 章 9.3–9.7 节,这里只讲在波动率模型中怎样用。

一个简单例子:随机抽 10 只股票,发现 1 只某天下跌、9 只没有下跌。设下跌概率为 \(p\),似然为 \(p(1-p)^9\)。对 \(p\) 求导令其为零得 \(\hat p=0.1\),与直觉一致。

23.5.2 估计常数方差

设 \(X\sim N(0,v)\),观测到 \(u_1,\dots,u_m\)。每个观测的似然是正态密度,联合似然为

\[\prod_{i=1}^m\frac{1}{\sqrt{2\pi v}}\exp\Big(-\frac{u_i^2}{2v}\Big)\tag{23.10}\]

取对数并去掉常数和常数乘子,等价于最大化

\[\sum_{i=1}^m\Big(-\ln v-\frac{u_i^2}{v}\Big)\tag{23.11}\]

即 \(-m\ln v-\sum u_i^2/v\)。对 \(v\) 求导:\(-m/v+\sum u_i^2/v^2=0\),得 \(\hat v=\frac1m\sum u_i^2\)。这就证明了 23.1 节第 3 处简化:除以 \(m\) 得到的是极大似然估计。

推导拆解:从 (23.10) 到 (23.11) 的每一步。 第一步,取对数,连乘变连加:\(\ln\prod_i(\cdots)=\sum_i\left[-\frac12\ln(2\pi)-\frac12\ln v-\frac{u_i^2}{2v}\right]\)。取对数不改变最大值的位置,因为对数是单调递增函数。 第二步,\(-\frac12\ln(2\pi)\) 与 \(v\) 无关,去掉;剩下每项都有因子 \(\frac12\),乘以 2 也不改变最大值位置。于是得到 \(\sum_i(-\ln v-u_i^2/v)\)。 第三步,对 \(v\) 求导:\(-\ln v\) 的导数是 \(-\frac1v\),共 \(m\) 项;\(-\frac{u_i^2}{v}\) 的导数是 \(+\frac{u_i^2}{v^2}\)。令和为零,两边乘 \(v^2\) 得 \(-mv+\sum u_i^2=0\)。 直觉:极大似然找的是"让观测到的这组收益看起来最不意外"的方差。\(v\) 太小,大收益的概率密度极低(\(-u^2/v\) 项惩罚);\(v\) 太大,所有收益的概率密度都被摊薄(\(-\ln v\) 项惩罚)。两者平衡点正好是 \(u^2\) 的平均值。

23.5.3 估计 EWMA 与 GARCH(1,1) 的参数

令 \(v_i=\sigma_i^2\) 为模型对第 \(i\) 天方差的估计,并假设给定方差时 \(u_i\) 条件正态。最优参数最大化

\[\sum_{i=1}^m\Big(-\ln v_i-\frac{u_i^2}{v_i}\Big)\tag{23.12}\]

它与式 (23.11) 形式完全相同,只是常数 \(v\) 换成了随时间变化的 \(v_i\)。由于 \(v_i\) 通过递推依赖参数,没有解析解,需要数值搜索。数值优化方法见第 04 册。

原书表 23.1 的实例(S&P 500,2005-07-18 至 2010-08-13,共 1,279 天)。表的各列是日期、天数 \(i\)、\(S_i\)、\(u_i=(S_i-S_{i-1})/S_{i-1}\)、\(v_i\)、\(-\ln v_i-u_i^2/v_i\)。第 3 天用 \(v_3=u_2^2\) 初始化,之后用式 (23.9) 递推,目标是让最后一列之和最大。示例行:第 3 天 \(u=0.004759\),\(v=0.00004531\),似然项 9.5022。结果:

模型 参数估计 目标函数最大值
GARCH(1,1) \(\omega=0.0000013465\),\(\alpha=0.083394\),\(\beta=0.910116\) 10,228.2349
GARCH(1,1) + 方差锁定 \(\alpha=0.08445\),\(\beta=0.9101\)(\(V_L\) 取样本方差 0.0002412) 10,228.1941
EWMA \(\lambda=0.9374\) 10,192.5104

由 GARCH 估计,\(V_L=\omega/(1-\alpha-\beta)=0.0000013465/0.006490=0.0002075\),长期日波动率 \(\sqrt{0.0002075}=1.4404\%\)。原书图 23.1、23.2 画出了 S&P 500 指数及其 GARCH(1,1) 日波动率:多数时候日波动率低于 2%,信用危机期间高达每天 5%(同期 VIX 也极高,见第 15.11 节)。

方差锁定(variance targeting;Engle & Mezrich 1996):令 \(V_L\) 等于样本方差(或其他合理值),则 \(\omega=V_L(1-\alpha-\beta)\),只需估计 \(\alpha,\beta\) 两个参数。本例中目标函数只比完全估计略低,但估计通常更稳健——长期方差这个最难估的量直接取自数据。

EWMA 只有一个参数,目标函数明显低于 GARCH:差了约 36 个对数似然单位,说明"长期方差的那一份权重"在这组数据上很有价值。

白话解释:"差 36 个单位"有多大,可以用似然比检验来衡量。注意 (23.12) 的目标函数在推导时乘了 2(见上一个讲解框),它等于"2 × 对数似然"加一个常数,所以两模型目标函数之差就直接是似然比统计量 \(2\Delta\ln L\approx36\)。EWMA 是 GARCH 加了两个约束(\(\omega=0\)、\(\alpha+\beta=1\))的特例,远大于自由度 2 的 \(\chi^2\) 分布 5% 临界值 5.99。所以统计上可以明确拒绝"EWMA 足够好"。相比之下,方差锁定只比完整 GARCH 低 0.04,说明把 \(V_L\) 固定为样本方差几乎没有损失。

23.5.4 数值实现的技巧

  • 参数缩放。优化器在各参数量级相近时表现最好。原书用 Excel Solver 时,让单元格分别存 \(\omega\times10^5\)、\(10\alpha\) 和 \(\beta\),再换算回来算似然。用 Python 时同样适用;也可以像 arch 包那样把收益乘以 100 再估计。
  • 多个初值。似然面可能有局部极大值,应从几组不同初值出发,取最好的结果。
  • 初始化。原书表 23.1 用第一个收益平方初始化方差,23.8 节则用总体方差初始化;两者对最终估计影响很小,但对前几十天的方差估计有影响。
  • 约束。\(\omega>0\)、\(\alpha,\beta\ge0\)、\(\alpha+\beta<1\)。

23.5.5 模型好不好

GARCH 的前提是波动率聚集:\(u_i^2\) 大时,其后的 \(u_{i+1}^2,u_{i+2}^2,\dots\) 也倾向于大,即 \(u_i^2\) 有正自相关。好的 GARCH 模型应当把这种自相关"解释掉"——标准化后的序列 \(u_i^2/\sigma_i^2\) 应该几乎没有自相关。

原书表 23.2(S&P 500,滞后 1–15):\(u_i^2\) 的自相关全为正,在 0.121–0.431 之间(如滞后 1 为 0.183,滞后 2 为 0.385,滞后 11 为 0.431);\(u_i^2/\sigma_i^2\) 的自相关有正有负,绝对值都很小(滞后 1 为 0.063,滞后 10 为 0.083,其余多在 0.04 以下)。

用 Ljung–Box 统计量(Ljung & Box 1978)做正式检验。序列有 \(m\) 个观测,\(\eta_k\) 是滞后 \(k\) 的自相关,\(K\) 是考虑的滞后数:

\[Q=m\sum_{k=1}^K w_k\eta_k^2,\qquad w_k=\frac{m+2}{m-k}\]

在"零自相关"的原假设下 \(Q\) 近似服从自由度为 \(K\) 的 \(\chi^2\) 分布。\(K=15\) 时 95% 临界值约为 25。本例 \(u_i^2\) 的 \(Q\approx1{,}566\)(强自相关),\(u_i^2/\sigma_i^2\) 的 \(Q=21.7\),低于 25,说明 GARCH 基本消除了自相关。


23.6 用 GARCH(1,1) 预测未来波动率

23.6.1 期望方差的路径

由 \(\sigma_n^2=(1-\alpha-\beta)V_L+\alpha u_{n-1}^2+\beta\sigma_{n-1}^2\),两边减去 \(V_L\):

\[\sigma_n^2-V_L=\alpha(u_{n-1}^2-V_L)+\beta(\sigma_{n-1}^2-V_L)\]

同样的关系在未来第 \(n+t\) 天成立:\(\sigma_{n+t}^2-V_L=\alpha(u_{n+t-1}^2-V_L)+\beta(\sigma_{n+t-1}^2-V_L)\)。取期望并利用 \(E[u_{n+t-1}^2]=E[\sigma_{n+t-1}^2]\),得到 \(E[\sigma_{n+t}^2-V_L]=(\alpha+\beta)E[\sigma_{n+t-1}^2-V_L]\)。反复使用:

\[E[\sigma_{n+t}^2]=V_L+(\alpha+\beta)^t(\sigma_n^2-V_L)\tag{23.13}\]

推导拆解:关键一步"\(E[u_{n+t-1}^2]=E[\sigma_{n+t-1}^2]\)"为什么成立。 模型假设 \(u_k=\sigma_k\epsilon_k\),\(\epsilon_k\) 是与过去无关的标准正态。站在第 \(k-1\) 天,\(\sigma_k\) 已知,所以条件期望 \(E[u_k^2\mid\text{第 }k-1\text{ 天信息}]=\sigma_k^2E[\epsilon_k^2]=\sigma_k^2\)。 站在今天看更远的未来,\(\sigma_k\) 本身还不知道,再对它取期望,就得到 \(E[u_k^2]=E[\sigma_k^2]\)(这一步叫"重期望法则":先在较近的时点取条件期望,再在今天取期望)。 代入后,\(\alpha(\cdot)\) 和 \(\beta(\cdot)\) 两项合并成 \((\alpha+\beta)E[\sigma_{n+t-1}^2-V_L]\)。偏离量每过一天乘以 \((\alpha+\beta)\),这就是一个公比为 \(\alpha+\beta\) 的几何数列,\(t\) 天后为 \((\alpha+\beta)^t\) 倍。 金融类比:和 CFA 里"均值回归的 AR(1) 过程的多步预测 \(x_{t+h}=\mu+\phi^h(x_t-\mu)\)"完全相同,\(\alpha+\beta\) 扮演 \(\phi\) 的角色。

这是用第 \(n-1\) 天末的信息预测第 \(n+t\) 天的方差。三种情况:

  • EWMA(\(\alpha+\beta=1\)):期望未来方差等于当前方差,没有均值回归。这也是第 06 册第 03a 章 3.6 节把 RiskMetrics 称为 IGARCH 的原因——它是"单位根"的 GARCH。
  • \(\alpha+\beta<1\):末项随 \(t\) 衰减,预测趋向 \(V_L\),回归速度由 \(1-\alpha-\beta\) 决定(原书图 23.3:当前方差高于或低于 \(V_L\) 时,期望路径分别从上方或下方趋近 \(V_L\))。
  • \(\alpha+\beta>1\):长期方差权重为负,过程"均值逃离"(mean fleeing)。

数值例 S&P 500 估计结果 \(\alpha+\beta=0.9935\),\(V_L=0.0002075\)。设当前方差 0.0003(日波动 1.732%):

  • 10 天后期望方差 \(0.0002075+0.9935^{10}(0.0003-0.0002075)=0.0002942\),日波动 1.72%,仍明显高于长期的 1.44%;
  • 500 天后期望方差 \(0.0002110\),日波动 1.45%,已接近长期水平。

23.6.2 波动率期限结构

在第 \(n\) 天,记 \(V(t)=E[\sigma_{n+t}^2]\),令 \(a=\ln\dfrac{1}{\alpha+\beta}\),式 (23.13) 可写成连续形式:

\[V(t)=V_L+e^{-at}\,[V(0)-V_L]\]

\(V(t)\) 是 \(t\) 天后的瞬时方差率。期权定价需要的是存续期内的平均方差率(第 27a 章会说明,在波动率随机但与股价无关时,用平均方差定价是合理的近似):

\[\frac1T\int_0^T V(t)\,dt=V_L+\frac{1-e^{-aT}}{aT}\,[V(0)-V_L]\]

\(T\) 越大,结果越接近 \(V_L\)。定义 \(\sigma(T)\) 为给 \(T\) 天期权定价所用的年化波动率,每年按 252 个交易日:

\[\sigma(T)^2=252\Big(V_L+\frac{1-e^{-aT}}{aT}\,[V(0)-V_L]\Big)\tag{23.14}\]

推导拆解:从离散到连续,再到平均值。 第一步,\((\alpha+\beta)^t=e^{t\ln(\alpha+\beta)}=e^{-at}\),其中 \(a=\ln\frac{1}{\alpha+\beta}>0\)。这只是把幂写成指数,与"离散复利换成连续复利"的做法相同。 第二步,求平均:\(\frac1T\int_0^T\left[V_L+e^{-at}(V(0)-V_L)\right]dt\)。常数 \(V_L\) 的平均还是 \(V_L\);\(\int_0^Te^{-at}dt=\frac{1-e^{-aT}}{a}\),除以 \(T\) 得 \(\frac{1-e^{-aT}}{aT}\)。 第三步,乘 252 把日方差年化,就是 (23.14)。 系数 \(\frac{1-e^{-aT}}{aT}\) 在 \(T\to0\) 时趋于 1(短期期权几乎完全用当前方差),在 \(T\) 很大时趋于 0(长期期权几乎完全用长期方差)。你可能在利率期限结构的 Nelson–Siegel 模型里见过同一个函数,它在那里也描述"短端到长端的过渡"。

当前波动率高于长期水平时,期限结构向下倾斜;低于长期水平时向上倾斜。GARCH 算出的期限结构一般不等于市场隐含波动率的期限结构,但常被用来预测隐含波动率期限结构怎样响应波动率的变化。

原书表 23.3 S&P 500 参数,\(a=\ln(1/0.99351)=0.006511\),\(V(0)=0.0003\),\(T\) 以天计:

期权期限(天) 10 30 50 100 500
期权波动率(年化 %) 27.36 27.10 26.87 26.35 24.32

23.6.3 波动率变化如何传导到各期限

把式 (23.14) 改写为 \(\sigma(T)^2=252\Big(V_L+\dfrac{1-e^{-aT}}{aT}\Big(\dfrac{\sigma(0)^2}{252}-V_L\Big)\Big)\),对 \(\sigma(0)\) 求导:\(2\sigma(T)\,d\sigma(T)=\dfrac{1-e^{-aT}}{aT}\cdot2\sigma(0)\,d\sigma(0)\),所以

\[\Delta\sigma(T)\approx\frac{1-e^{-aT}}{aT}\cdot\frac{\sigma(0)}{\sigma(T)}\,\Delta\sigma(0)\tag{23.15}\]

原书表 23.4 \(V(0)=0.0003\),\(\sigma(0)=\sqrt{252}\times\sqrt{0.0003}=27.50\%\),瞬时波动率上升 100 个基点到 28.50%:

期权期限(天) 10 30 50 100 500
波动率上升(%) 0.97 0.92 0.87 0.77 0.33

实务意义:许多金融机构据此计算账簿对波动率的暴露。算 vega 时,不假设所有隐含波动率一律平移 1%,而是让冲击随期限递减:10 天期权用 0.97%,30 天用 0.92%,以此类推。这比平行移动更接近市场实际——短期隐含波动率的变动通常比长期大得多。


23.7 相关系数

23.7.1 协方差是基本变量

两个变量 \(X,Y\) 的相关系数为

\[\rho=\frac{\text{cov}(X,Y)}{\sigma_X\sigma_Y},\qquad \text{cov}(X,Y)=E[(X-\mu_X)(Y-\mu_Y)]\]

相关系数更直观,但在分析中真正被更新的基本变量是协方差——正如前面 EWMA/GARCH 更新的是方差而不是波动率。

设 \(x_i,y_i\) 为 \(X,Y\) 第 \(i\) 天的百分比变化,\(\sigma_{x,n},\sigma_{y,n}\) 为第 \(n\) 天的日波动率估计,\(\text{cov}_n\) 为协方差估计,则相关系数估计为 \(\text{cov}_n/(\sigma_{x,n}\sigma_{y,n})\)。

等权、零均值时:\(\sigma_{x,n}^2=\frac1m\sum x_{n-i}^2\),\(\sigma_{y,n}^2=\frac1m\sum y_{n-i}^2\),以及

\[\text{cov}_n=\frac1m\sum_{i=1}^m x_{n-i}y_{n-i}\tag{23.16}\]

EWMA 版本:

\[\text{cov}_n=\lambda\,\text{cov}_{n-1}+(1-\lambda)x_{n-1}y_{n-1}\]

23.7.2 例题

例 23.3 \(\lambda=0.95\),第 \(n-1\) 天的相关系数估计为 0.6,\(X,Y\) 的日波动率估计为 1% 和 2%。于是协方差为 \(0.6\times0.01\times0.02=0.00012\)。当天 \(X,Y\) 分别变动 0.5% 和 2.5%:

\[\sigma_{x,n}^2=0.95\times0.01^2+0.05\times0.005^2=0.00009625\]
\[\sigma_{y,n}^2=0.95\times0.02^2+0.05\times0.025^2=0.00041125\]
\[\text{cov}_n=0.95\times0.00012+0.05\times0.005\times0.025=0.00012025\]

新的波动率为 0.981% 和 2.028%,新相关系数为 \(0.00012025/(0.00981\times0.02028)=0.6044\)。

23.7.3 GARCH 协方差与多元 GARCH

GARCH(1,1) 的协方差版本为

\[\text{cov}_n=\omega+\alpha x_{n-1}y_{n-1}+\beta\,\text{cov}_{n-1}\]

长期平均协方差为 \(\omega/(1-\alpha-\beta)\)。与式 (23.13)、(23.14) 类比,可以预测未来协方差和期权存续期内的平均协方差。更一般地,多元 GARCH(multivariate GARCH)以一致的方式更新整个协方差矩阵;其模型族(VEC、BEKK、DCC 等)属于 Tsay 原书第 10 章多元波动率模型的内容,见第 06 册第 10 章。

23.7.4 协方差矩阵的一致性条件

\(N\) 个变量的方差–协方差矩阵 \(\Omega\) 必须满足

\[w^{\mathsf T}\Omega w\ge0\quad\text{对所有 }N\times1\text{ 向量 }w\tag{23.17}\]

即 \(\Omega\) 半正定(positive semidefinite)。理由很简单:\(w^{\mathsf T}\Omega w\) 是组合 \(w_1x_1+\dots+w_Nx_N\) 的方差,方差不可能为负。矩阵半正定的判别和性质见第 01 册第 07a 章。

反例 相关矩阵

\[\begin{pmatrix}1&0&0.9\\0&1&0.9\\0.9&0.9&1\end{pmatrix}\]

变量 1、2 都与变量 3 高度相关,却彼此不相关,这本身就很可疑。取 \(w=(1,1,-1)^{\mathsf T}\):\(w^{\mathsf T}\Omega w=3+2(0)-2(0.9)-2(0.9)=3-3.6=-0.6<0\),所以它不是合法的相关矩阵(它的最小特征值为 \(-0.27\))。对 \(3\times3\) 相关矩阵,内部一致的条件为

\[\rho_{12}^2+\rho_{13}^2+\rho_{23}^2-2\rho_{12}\rho_{13}\rho_{23}\le1\]

实务纪律:要保证半正定,方差和协方差必须用一致的方法计算。方差用最近 \(m\) 天等权,协方差也必须用同样的 \(m\) 天等权;方差用 \(\lambda=0.94\) 的 EWMA,协方差也必须用 \(\lambda=0.94\)。原因是:如果每个元素都是同一组"外积" \(x_ix_i^{\mathsf T}\) 的非负加权平均,那么整个矩阵就是半正定矩阵的非负组合,自然半正定;一旦方差和协方差用不同的权重,这个保证就没有了。

推导拆解:为什么"同一组外积的非负加权平均"一定半正定。 设第 \(i\) 天的收益向量为 \(x_i\)(各资产当天收益组成的列向量),外积 \(x_ix_i^{\mathsf T}\) 是一个矩阵,其 \((j,k)\) 元素是资产 \(j\) 与 \(k\) 当天收益的乘积。EWMA 协方差矩阵就是 \(\Omega=\sum_i c_ix_ix_i^{\mathsf T}\),\(c_i\ge0\) 是权重。 对任意权重向量 \(w\):\(w^{\mathsf T}\Omega w=\sum_ic_i(w^{\mathsf T}x_i)(x_i^{\mathsf T}w)=\sum_ic_i(w^{\mathsf T}x_i)^2\ge0\)。 中间的 \(w^{\mathsf T}x_i\) 就是组合 \(w\) 在第 \(i\) 天的收益,平方后非负,再乘非负权重相加,结果一定非负。 金融含义:这样算出的"组合方差"本质上就是"组合历史收益平方的加权平均",当然不会为负。如果方差用 \(\lambda=0.94\)、协方差用 \(\lambda=0.97\),矩阵就不再能写成这种形式,上面的论证失效。


23.8 例:四指数组合的 EWMA VaR

延续第 22.2 节的例子:2008 年 9 月 25 日的组合包含 DJIA 400 万美元、FTSE 100 300 万美元、CAC 40 100 万美元、Nikkei 225 200 万美元,用截至该日的 500 天日收益。

等权估计。相关矩阵(原书表 23.5):DJIA–FTSE 0.489,DJIA–CAC 0.496,DJIA–Nikkei 0.062,FTSE–CAC 0.918,FTSE–Nikkei 0.201,CAC–Nikkei 0.211。协方差矩阵(表 23.6)对角线为 0.0001227、0.0002010、0.0001950、0.0001909。由式 (22.3),组合损失(单位:千美元)的方差为 8,761.833,标准差 93.60,一日 99% VaR \(=2.33\times93.60=217.757\),即 217,757 美元(对比历史模拟法的 253,385 美元)。

EWMA(\(\lambda=0.94\))。协方差矩阵(表 23.7)对角线为 0.0004801、0.0010314、0.0009535、0.0002541。组合损失方差为 40,995.765,标准差 202.474,一日 99% VaR \(=2.33\times202.474=471.025\),即 471,025 美元,是等权结果的两倍多。

原因:多头组合的标准差随各资产的波动率和相关系数上升而上升,而 2008 年 9 月前夕两者都在升高。

日波动率(%) DJIA FTSE 100 CAC 40 Nikkei 225
等权 1.11 1.42 1.40 1.38
EWMA 2.19 3.21 3.09 1.59

EWMA 相关系数(表 23.9):DJIA–FTSE 0.611,DJIA–CAC 0.629,DJIA–Nikkei 0.113,FTSE–CAC 0.971,FTSE–Nikkei 0.409,CAC–Nikkei 0.342,普遍高于等权估计。相关性在不利市况下倾向上升,这是风险管理中最重要的经验事实之一:分散化恰恰在最需要的时候失效。

金融直觉:把这个例子与第 22 章连起来看。同一个组合、同一天,历史模拟给出 25.3 万,等权方差-协方差给出 21.8 万,EWMA 方差-协方差给出 47.1 万。差异主要来自"用多旧的数据代表明天":等权把两年前的平静日子和上周的动荡日子同等看待;EWMA 的有效窗口只有约 17 天,几乎只看危机中的最近几周。哪个"对"取决于用途:监管资本希望稳定、不随市场剧烈摆动,日内风控则希望及时反应。巴塞尔 2.5 同时要求计算普通 VaR 和压力 VaR,正是为了兼顾两方面。

这里 EWMA 的初始方差设为总体方差(而不是像表 23.1 那样用第一个收益平方);两种初始化对最终结果影响很小。


量化实战

本章内容在量化中的用途

  1. 风险模型。EWMA 协方差是日常风控、组合 VaR 和压力监控的标准工具。商业风险模型常用"半衰期 \(h\) 天"描述衰减,等价于 \(\lambda=0.5^{1/h}\);它们也常对波动率和相关使用不同的半衰期(波动率用较短的半衰期以快速反应,相关用较长的半衰期以保持稳定)。从 23.7.4 节可知这样做会失去"自动半正定"的保证,需要事后修正(如对相关矩阵做特征值截断或求最近相关矩阵),或者先分别估计波动率与相关矩阵,再用 \(D R D\) 的形式组合——只要相关矩阵 \(R\) 本身是由同一组标准化收益的外积一致估计的,\(DRD\) 就仍然半正定。
  2. 波动率目标与风险平价。按 \(\sigma_n\) 缩放仓位(目标波动率 / 预测波动率)是 CTA、风险平价和许多因子组合的常规做法。\(\lambda\) 的选择是"反应速度 vs 换手"的权衡:\(\lambda\) 小,风险控制更及时,但仓位换手更大。
  3. 信号标准化。动量、反转、价差等信号常除以 EWMA 波动率再做横截面比较或设阈值;\(u_i/\sigma_i\) 也是检验"信号残差里是否还有可利用结构"的起点(Ljung–Box)。
  4. 期权定价与波动率交易。式 (23.14) 给出模型的波动率期限结构,可以与隐含波动率期限结构比较,寻找相对价值;式 (23.15) 是按期限缩放 vega 冲击的依据,用于把一个期权簿的 vega 风险合并成一个数字。
  5. 回测纪律。严格使用"第 \(n-1\) 天末估计第 \(n\) 天"的约定;GARCH 参数若用全样本估计再回放,就引入了前视偏差,应该用滚动或扩展窗口重新估计。

Python 示例:从参数估计到 VaR

下面的代码分五步:(1) 用原书表 23.1 的 S&P 500 参数模拟 2,500 天 GARCH(1,1) 收益;(2) 自己写式 (23.12) 的对数似然,分别估计完整 GARCH、方差锁定 GARCH 和 EWMA,并与 arch 包对照;(3) 用 Ljung–Box 检验 \(u^2\) 与 \(u^2/\sigma^2\);(4) 用式 (23.14)、(23.15) 算波动率期限结构和冲击传导;(5) 对一个四资产组合比较等权与 EWMA 协方差下的 VaR,并演示方差、协方差混用不同 \(\lambda\) 会破坏半正定。

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
from statsmodels.stats.diagnostic import acorr_ljungbox
from arch import arch_model

rng = np.random.default_rng(23)

# ---------- 1) 用 Hull 表 23.1 的参数模拟 GARCH(1,1) 日收益 ----------
w0, a0, b0 = 1.3465e-6, 0.083394, 0.910116
n = 2500
u = np.empty(n); v = np.empty(n)
v[0] = w0 / (1 - a0 - b0)
for t in range(n):
    if t > 0:
        v[t] = w0 + a0 * u[t-1]**2 + b0 * v[t-1]
    u[t] = np.sqrt(v[t]) * rng.standard_normal()

# ---------- 2) 极大似然:最大化 sum(-ln v_i - u_i^2 / v_i),式 (23.12) ----------
def garch_var(params, u, v0):
    w, a, b = params
    v = np.empty_like(u); v[0] = v0
    for t in range(1, len(u)):
        v[t] = w + a * u[t-1]**2 + b * v[t-1]
    return v

def loglik(params, u, v0):
    v = garch_var(params, u, v0)
    return np.sum(-np.log(v[1:]) - u[1:]**2 / v[1:])   # 第 1 天只用于初始化

v0 = u[:20].var()          # 用前 20 天样本方差初始化(原书表 23.1 用第一个 u^2,影响很小)
# 参数缩放:优化 (w*1e5, 10a, b),与原书 Solver 技巧相同
def neg_garch(x):
    w, a, b = x[0] / 1e5, x[1] / 10, x[2]
    if w <= 0 or a < 0 or b < 0 or a + b >= 1: return 1e10
    return -loglik((w, a, b), u, v0)
best = min((minimize(neg_garch, x0, method="Nelder-Mead", options={"xatol":1e-8,"fatol":1e-8,"maxiter":5000})
            for x0 in ([0.1, 0.5, 0.9], [0.5, 1.0, 0.85], [0.05, 0.3, 0.95])), key=lambda r: r.fun)
w, a, b = best.x[0] / 1e5, best.x[1] / 10, best.x[2]
VL = w / (1 - a - b)
print(f"GARCH MLE: omega={w:.3e} alpha={a:.4f} beta={b:.4f} a+b={a+b:.4f} "
      f"长期日波动率={np.sqrt(VL):.3%} 目标函数={-best.fun:.2f}")

# 方差锁定:VL=样本方差,只估 alpha、beta
VL_s = u.var()
def neg_vt(x):
    a_, b_ = x
    if a_ < 0 or b_ < 0 or a_ + b_ >= 1: return 1e10
    return -loglik((VL_s * (1 - a_ - b_), a_, b_), u, v0)
vt = minimize(neg_vt, [0.08, 0.9], method="Nelder-Mead")
print(f"方差锁定:  alpha={vt.x[0]:.4f} beta={vt.x[1]:.4f} 目标函数={-vt.fun:.2f}")

# EWMA:omega=0, alpha=1-lambda, beta=lambda,只估 lambda
def neg_ewma(lam):
    lam = float(lam[0])
    if not 0 < lam < 1: return 1e10
    return -loglik((0.0, 1 - lam, lam), u, v0)
ew = minimize(neg_ewma, [0.94], method="Nelder-Mead")
print(f"EWMA:      lambda={ew.x[0]:.4f} 目标函数={-ew.fun:.2f}")

# 与 arch 包对照(收益放大 100 倍以改善数值条件)
res = arch_model(100 * u, mean="Zero", vol="GARCH", p=1, q=1).fit(disp="off")
print("arch 包:  ", {k: round(float(x), 4) for k, x in res.params.items()},
      f"-> omega 折回日收益尺度 {res.params['omega']/1e4:.3e}")

# ---------- 3) 模型检验:u^2 与 u^2/sigma^2 的 Ljung-Box(15) ----------
vhat = garch_var((w, a, b), u, v0)
q_raw = acorr_ljungbox(u[1:]**2, lags=[15])["lb_stat"].iloc[0]
q_std = acorr_ljungbox(u[1:]**2 / vhat[1:], lags=[15])["lb_stat"].iloc[0]
print(f"Ljung-Box Q(15): u^2 = {q_raw:.1f}, u^2/sigma^2 = {q_std:.1f} (95% 临界值 25.0)")

# ---------- 4) 波动率期限结构 (23.14) 与按期限缩放的 vega 冲击 (23.15) ----------
V0 = vhat[-1]
aa = np.log(1 / (a + b))
print(f"当前日波动率 {np.sqrt(V0):.3%}, 长期 {np.sqrt(VL):.3%}")
for T in [10, 30, 50, 100, 250, 500]:
    k = (1 - np.exp(-aa * T)) / (aa * T)
    sT = np.sqrt(252 * (VL + k * (V0 - VL)))
    s0 = np.sqrt(252 * V0)
    print(f"  T={T:>3}天  sigma(T)={sT:.2%}  sigma(0)上升1个百分点时 sigma(T)上升 {k*s0/sT:.2f} 个百分点")

# ---------- 5) 多资产 EWMA 协方差、VaR 与一致性(半正定) ----------
N, m = 4, 500
C_calm = 0.0001 * np.array([[1.2,.6,.6,.1],[.6,2.,1.8,.3],[.6,1.8,2.,.3],[.1,.3,.3,1.9]])
X = rng.multivariate_normal(np.zeros(N), C_calm, size=m)
X[-40:] *= 2.0                          # 最后 40 天进入高波动期
pos = np.array([4, 3, 1, 2]) * 1e3      # 头寸,单位:千美元

def ewma_cov(X, lam):
    S = np.cov(X.T, bias=True)          # 用全样本(零均值近似)协方差初始化
    for x in X:
        S = lam * S + (1 - lam) * np.outer(x, x)
    return S

for name, S in [("等权", X.T @ X / m), ("EWMA λ=0.94", ewma_cov(X, 0.94))]:
    sd = np.sqrt(pos @ S @ pos)
    print(f"{name:12s} 组合日标准差 {sd:7.2f} 千美元, 1日99% VaR {norm.ppf(0.99)*sd:7.2f} 千美元,"
          f" 最小特征值 {np.linalg.eigvalsh(S).min():.2e}")

# 方差和协方差用不同 λ 会破坏一致性:方差用 0.99、协方差用 0.80
S_v, S_c = ewma_cov(X, 0.99), ewma_cov(X, 0.80)
S_mix = S_c.copy(); np.fill_diagonal(S_mix, np.diag(S_v))
d = np.sqrt(np.diag(S_mix)); print("混用 λ 的相关矩阵:\n", np.round(S_mix / np.outer(d, d), 3))
print("混用 λ 的最小特征值:", f"{np.linalg.eigvalsh(S_mix).min():.2e}")
O = np.array([[1, 0, .9], [0, 1, .9], [.9, .9, 1]]); wv = np.array([1, 1, -1])
print("原书反例 w'Ωw =", wv @ O @ wv, " 特征值:", np.round(np.linalg.eigvalsh(O), 4))

关键输出:

GARCH MLE: omega=1.819e-06 alpha=0.0827 beta=0.9051 a+b=0.9878 长期日波动率=1.219% 目标函数=20143.51
方差锁定:  alpha=0.0808 beta=0.9053 目标函数=20143.35
EWMA:      lambda=0.9270 目标函数=20117.46
arch 包:   {'omega': 0.0182, 'alpha[1]': 0.0827, 'beta[1]': 0.905} -> omega 折回日收益尺度 1.820e-06
Ljung-Box Q(15): u^2 = 766.8, u^2/sigma^2 = 10.2 (95% 临界值 25.0)
当前日波动率 1.284%, 长期 1.219%
  T= 10天  sigma(T)=20.33%  sigma(0)上升1个百分点时 sigma(T)上升 0.94 个百分点
  T= 30天  sigma(T)=20.22%  sigma(0)上升1个百分点时 sigma(T)上升 0.84 个百分点
  T= 50天  sigma(T)=20.13%  sigma(0)上升1个百分点时 sigma(T)上升 0.76 个百分点
  T=100天  sigma(T)=19.96%  sigma(0)上升1个百分点时 sigma(T)上升 0.59 个百分点
  T=250天  sigma(T)=19.68%  sigma(0)上升1个百分点时 sigma(T)上升 0.32 个百分点
  T=500天  sigma(T)=19.53%  sigma(0)上升1个百分点时 sigma(T)上升 0.17 个百分点
等权           组合日标准差  105.88 千美元, 1日99% VaR  246.32 千美元, 最小特征值 2.44e-05
EWMA λ=0.94  组合日标准差  218.40 千美元, 1日99% VaR  508.08 千美元, 最小特征值 7.33e-05
混用 λ 的相关矩阵:
 [[1.    1.781 2.074 0.197]
 [1.781 1.    2.139 0.147]
 [2.074 2.139 1.    0.08 ]
 [0.197 0.147 0.08  1.   ]]
混用 λ 的最小特征值: -6.15e-04
原书反例 w'Ωw = -0.6000000000000001  特征值: [-0.2728  1.      2.2728]

读输出时注意几点:

  • 自写 MLE 与 arch 包的结果一致。真实参数 \(\alpha+\beta=0.9935\),2,500 天样本估出 0.9878:持续性参数 \(\alpha+\beta\) 和长期方差 \(V_L\) 是最难估准的量,这正是方差锁定有吸引力的原因。
  • 方差锁定的目标函数只比完整 GARCH 低 0.16,EWMA 低约 26,与原书 S&P 500 的结论模式相同。
  • \(u^2\) 的 Ljung–Box \(Q\) 远超临界值,标准化后降到 10.2,模型把波动率聚集解释掉了。
  • 期限结构中,当前波动率略高于长期水平,所以 \(\sigma(T)\) 随期限缓慢下降;瞬时波动率的冲击对 10 天期权几乎完全传导(0.94),对 500 天期权只传导 0.17。由于本例 \(\alpha+\beta\) 比原书小,衰减比表 23.4 快。
  • 最后 40 天波动翻倍后,EWMA VaR 是等权的两倍多,与 23.8 节的四指数例子同理。
  • 方差用慢的 \(\lambda\)、协方差用快的 \(\lambda\),在波动突然上升时协方差涨得比方差快,"相关系数"超过 1,矩阵出现负特征值。这样的矩阵送进均值–方差优化器,会被利用出"负方差"的虚假对冲组合。

本章小结

BSM 等主流模型假设波动率恒定,但实际上波动率随机变化且不能直接观测。日波动率监控把方差率写成过去收益平方的加权平均,近期观测权重更大。EWMA 的权重按 \(\lambda\) 指数衰减,只需保存一个状态,RiskMetrics 取 \(\lambda=0.94\);GARCH(1,1) 在指数衰减之外给长期平均方差 \(V_L\) 一份权重,因而具有均值回归,可以预测未来方差路径和期权存续期的平均波动率,EWMA 是它 \(\omega=0\) 的特例。参数用极大似然估计,目标是最大化 \(\sum(-\ln v_i-u_i^2/v_i)\),方差锁定可以减少一个参数;模型好坏看它能否消除 \(u_i^2\) 的自相关(Ljung–Box)。由 GARCH 得到的波动率期限结构和"冲击随期限递减"的规律被用于按期限缩放 vega。每个方差模型都有对应的协方差模型,可用来更新 VaR 所需的整个协方差矩阵,但方差与协方差必须用一致的方法计算,才能保证矩阵半正定。危机前夕波动率和相关性同时上升,使 EWMA VaR 远高于等权估计。

概念 / 公式 内容
监控用方差估计 \(\sigma_n^2=\frac1m\sum_{i=1}^m u_{n-i}^2\),\(u_i\) 为简单收益,零均值,除以 \(m\)(即 MLE)
ARCH(\(m\)) \(\sigma_n^2=\omega+\sum_{i=1}^m\alpha_iu_{n-i}^2\),\(\omega=\gamma V_L\)
EWMA \(\sigma_n^2=\lambda\sigma_{n-1}^2+(1-\lambda)u_{n-1}^2\),权重 \((1-\lambda)\lambda^{i-1}\),半衰期 \(\ln0.5/\ln\lambda\)
GARCH(1,1) \(\sigma_n^2=\omega+\alpha u_{n-1}^2+\beta\sigma_{n-1}^2\),\(V_L=\omega/(1-\alpha-\beta)\),需 \(\alpha+\beta<1\)
连续时间对应 \(dV=(1-\alpha-\beta)(V_L-V)dt+\alpha\sqrt2\,V\,dz\)(时间以天计)
对数似然 最大化 \(\sum_i(-\ln v_i-u_i^2/v_i)\);方差锁定 \(\omega=V_L(1-\alpha-\beta)\)
模型检验 \(u_i^2/\sigma_i^2\) 的 Ljung–Box:\(Q=m\sum_k\frac{m+2}{m-k}\eta_k^2\),\(K=15\) 时临界值约 25
方差预测 \(E[\sigma_{n+t}^2]=V_L+(\alpha+\beta)^t(\sigma_n^2-V_L)\)
期权波动率期限结构 \(\sigma(T)^2=252\big(V_L+\frac{1-e^{-aT}}{aT}[V(0)-V_L]\big)\),\(a=\ln\frac1{\alpha+\beta}\)
冲击传导 \(\Delta\sigma(T)\approx\frac{1-e^{-aT}}{aT}\frac{\sigma(0)}{\sigma(T)}\Delta\sigma(0)\)
EWMA 协方差 \(\text{cov}_n=\lambda\,\text{cov}_{n-1}+(1-\lambda)x_{n-1}y_{n-1}\)
一致性条件 \(w^{\mathsf T}\Omega w\ge0\)(半正定);方差、协方差用同一方法;\(3\times3\):\(\rho_{12}^2+\rho_{13}^2+\rho_{23}^2-2\rho_{12}\rho_{13}\rho_{23}\le1\)

练习

基础题

  1. \(\lambda=0.94\),昨天的日波动率估计为 1.5%,今天收益为 −3%。用 EWMA 求明天的日波动率估计(参考原书习题 23.3、23.7)。 提示:\(\sigma^2=0.94\times0.015^2+0.06\times0.03^2=0.0002655\),\(\sigma=1.629\%\)。注意收益符号不影响 EWMA。

  2. GARCH(1,1) 为 \(\sigma_n^2=0.000003+0.04u_{n-1}^2+0.94\sigma_{n-1}^2\)。求长期平均方差、长期日波动率和年化波动率(按 252 天),并写出它对应的均值回归方程(参考原书习题 23.10)。 提示:\(\gamma=0.02\),\(V_L=0.00015\),日波动率 1.225%,年化 19.44%;\(dV=0.02(0.00015-V)dt+0.04\sqrt2\,V\,dz\)。

  3. 求 \(\lambda=0.94\) 与 \(\lambda=0.97\) 的 EWMA 半衰期。如果某风险模型说"波动率半衰期 42 天",对应的 \(\lambda\) 是多少? 提示:11.2 天、22.8 天;\(\lambda=0.5^{1/42}=0.9836\)。

  4. \(\lambda=0.94\),\(X,Y\) 的日波动率估计为 1% 和 1.2%,相关系数 0.5。今天 \(X\) 涨 1.5%,\(Y\) 跌 0.5%。更新两者的波动率和相关系数(参考原书习题 23.9、23.11)。 提示:\(\sigma_x^2=0.0001075\),\(\sigma_y^2=0.00013686\),\(\text{cov}=0.94\times0.00006+0.06\times(-0.000075)=0.0000519\);新波动率 1.037%、1.170%,新相关系数 0.428。

  5. 某分析师给出三个资产的相关系数 \(\rho_{12}=0.8\),\(\rho_{13}=0.8\),\(\rho_{23}=0.2\)。这个矩阵内部一致吗? 提示:\(0.64+0.64+0.04-2\times0.8\times0.8\times0.2=1.064>1\),不一致(最小特征值约 −0.036)。直觉:1 与 2、1 与 3 都高度相关,2 与 3 就不可能只有 0.2。

  6. 解释 EWMA 与 GARCH(1,1) 的区别,并说明为什么 GARCH(1,1) 的最优 \(\omega\) 为负时应改用 EWMA(原书习题 23.1、23.2)。 提示:GARCH 多一份给 \(V_L\) 的权重,因而均值回归;\(\omega<0\) 意味着 \(\gamma V_L<0\),模型不稳定。

进阶题

  1. 用第 2 题的 GARCH 模型,设当前日波动率为 1.8%。(a) 求 20 天后期望日波动率;(b) 求 30 天和 252 天期权定价所用的年化波动率;(c) 若瞬时年化波动率上升 1 个百分点,这两个期权的波动率各上升多少(参考原书习题 23.10、23.20)? 提示:(a) \(0.00015+0.98^{20}\times0.000174=0.0002662\),日波动率 1.631%;(b) \(a=\ln(1/0.98)=0.0202\),\(\sigma(30)=26.59\%\),\(\sigma(252)=21.53\%\);(c) \(\sigma(0)=28.57\%\),分别上升约 0.81 和 0.26 个百分点。

  2. 证明 GARCH(1,1) 与连续时间模型 \(dV=a(V_L-V)dt+\xi V\,dz\)(\(a=1-\alpha-\beta\),\(\xi=\alpha\sqrt2\),时间以天计)等价,并把它换算成时间以年计、每年 252 个交易日的形式(原书习题 23.14)。 提示:见 23.3.3 节推导。时间以年计时 \(dt_{\text{年}}=dt_{\text{天}}/252\):漂移速度变为 \(252a\);\(dz\) 的方差按时间缩放,\(\xi\) 变为 \(\xi\sqrt{252}\);若把 \(V\) 也换成年方差率,\(V_L\) 乘以 252,方程形式不变。

  3. 说明为什么"方差用等权 500 天、协方差用 EWMA"会可能得到非半正定矩阵,而"全部用同一个 \(\lambda\) 的 EWMA"一定半正定。 提示:同一个 \(\lambda\) 时 \(\Omega=\sum_i c_i\,x_ix_i^{\mathsf T}\)(\(c_i\ge0\)),每个 \(x_ix_i^{\mathsf T}\) 半正定,非负组合仍半正定;混用权重后没有这种结构。参考代码第 (5) 步。

  4. 编程:在代码第 (1) 步生成的序列上,复现 RiskMetrics 选 \(\lambda\) 的方法——对一组 \(\lambda\) 计算 \(\sum_i(v_i-\beta_i)^2\),其中 \(\beta_i\) 为第 \(i\) 天起之后 25 天 \(u^2\) 的等权平均,找出使之最小的 \(\lambda\),与 MLE 得到的 \(\lambda\) 比较(原书习题 23.19)。 提示:两种准则衡量的东西不同——MLE 看一步预测的似然,RiskMetrics 准则看对未来 25 天平均方差的预测,后者通常偏好更大的 \(\lambda\)(更平滑)。

原书推荐习题:23.3、23.8(EWMA 与 GARCH 单步更新);23.9、23.11(协方差与相关系数递推更新);23.10、23.20(长期波动率、期限结构、冲击传导,覆盖 23.6 节全部公式);23.14(GARCH 与连续时间随机波动率模型等价,衔接第 27a 章);23.19、23.22(用真实数据做已实现方差拟合与 MLE,最接近量化实务)。其他:23.4、23.6(参数变化的影响);23.12、23.13(汇率换算后指数的波动率与相关系数:\(Z=XY\) 时百分比变化近似相加);23.15、23.16、23.21(用作者网站数据重算四指数 VaR)。


原书对照

本章小节 原书小节 PDF 页码 原书页码
23.0 为什么要估计当前波动率 Chapter 23 引言 544 543
23.1 估计波动率 23.1 Estimating Volatility 544–546 543–545
23.2 EWMA 23.2 The Exponentially Weighted Moving Average Model 546–547 545–546
23.3 GARCH(1,1) 23.3 The GARCH(1,1) Model 548–549 547–548
23.4 模型选择 23.4 Choosing between the Models 549–550 548–549
23.5 极大似然估计 23.5 Maximum Likelihood Methods 550–555 549–554
23.6 用 GARCH 预测未来波动率 23.6 Using GARCH(1,1) to Forecast Future Volatility 555–558 554–557
23.7 相关系数 23.7 Correlations 558–560 557–559
23.8 四指数组合的 EWMA VaR 23.8 Application of EWMA to Four-Index Example 561–562 560–561
小结与习题 Summary, Further Reading, Practice Questions 563–566 562–565

延伸阅读:Bollerslev (1986) GARCH 原始论文;Engle (1982) ARCH 原始论文;Engle & Mezrich (1995, 1996) 方差锁定与 GARCH 预测;Engle & Ng (1993) 非对称消息冲击;Cumby, Figlewski & Hasbrook (1993) EGARCH 预测;Noh, Engle & Kane (1994) 用 GARCH 预测波动率为期权定价。时间序列视角的完整处理见第 06 册第 03a、03b 章。