第 10 章 多元波动率模型

本章对应 Tsay 原书第 10 章。第 03a、03b 章的 GARCH 类模型描述了单个资产波动率的时变性；本章把它推广到多个资产，研究条件协方差矩阵 \(\boldsymbol\Sigma_t\) 如何随时间演化。组合风险、动态对冲比、多资产 VaR、风险平价与最小方差组合都需要 \(\boldsymbol\Sigma_t\) 的预测；而"危机时相关性上升、分散化失效"这一现象，也只有在时变相关的框架里才能被正确刻画。本章的重点是 EWMA、DCC 和 Tsay 的 Cholesky 参数化。第 08 册第 23 章用 EWMA 与 GARCH(1,1) 逐日更新协方差和相关系数，是本章 10.2 节的实务版本，可对照阅读。

学习目标

1. 理解多元波动率建模的两大难点——维数灾难与正定约束，能计算各类模型的参数个数。
2. 掌握 EWMA 协方差估计（RiskMetrics），会用极大似然估计衰减因子 \(\lambda\)。
3. 了解 DVEC 与 BEKK 模型的形式、优缺点。
4. 掌握两种重参数化：\(\boldsymbol\Sigma_t=\boldsymbol D_t\boldsymbol\rho_t\boldsymbol D_t\) 与 Cholesky 分解 \(\boldsymbol\Sigma_t=\boldsymbol L_t\boldsymbol G_t\boldsymbol L_t'\)，理解 Cholesky 参数的回归含义（\(q_{21,t}\) 就是时变 β）。
5. 掌握常相关（CCC）模型、时变相关模型和动态条件相关（DCC）模型，会用两步法估计 DCC。
6. 会用多元波动率模型计算组合 VaR、时变对冲比，了解高维时的序贯 Cholesky 建模和因子–波动率模型。

读前导读

这一章在解决什么问题。 单变量 GARCH 告诉你"明天这只股票的波动率是多少"；做组合风险、VaR 或对冲，你还需要"明天这几只资产之间的协方差和相关系数是多少"。本章就是把 GARCH 从一个数（方差）推广到一个矩阵（协方差矩阵 \(\boldsymbol\Sigma_t\)）。你在 FRM/CFA 里见过 RiskMetrics 的 EWMA 波动率和组合 VaR 公式 \(\sqrt{\boldsymbol w'\boldsymbol\Sigma\boldsymbol w}\)；本章讲的是 \(\boldsymbol\Sigma\) 应该怎样每天更新，以及为什么"危机时相关性上升"必须用时变相关模型才能抓住。

推广的难点只有两个：矩阵元素太多（维数灾难），以及每天的矩阵都必须"合法"——任何组合的方差都要为正（正定约束）。本章所有模型都是在这两个难点之间找平衡：EWMA 最省参数；DVEC、BEKK 是"直接推广"，参数很多；CCC、DCC 把"波动率"和"相关系数"拆开分别建模；Cholesky 参数化把协方差矩阵变成一串回归——其中的回归系数就是你熟悉的时变 β 和最小方差对冲比 \(h^*=\mathrm{Cov}(\Delta S,\Delta F)/\mathrm{Var}(\Delta F)\)。

需要先想起来的数学。

• 正定矩阵。 对任意非零向量 \(\boldsymbol w\) 都有 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w>0\)，即"任何组合的方差都为正"。半正定允许等于 0。两个有用的事实：正定加半正定仍正定；\(\boldsymbol A\boldsymbol S\boldsymbol A'\) 在 \(\boldsymbol S\) 半正定时也半正定（因为 \(\boldsymbol w'\boldsymbol A\boldsymbol S\boldsymbol A'\boldsymbol w=(\boldsymbol A'\boldsymbol w)'\boldsymbol S(\boldsymbol A'\boldsymbol w)\ge0\)）。见 第 00 册第 06 章 线性代数速成。
• Cholesky（LDL）分解。 正定矩阵可写成"单位下三角 × 正对角 × 单位上三角"。2×2 时 \(\begin{bmatrix}4&3\\3&9\end{bmatrix}=\begin{bmatrix}1&0\\0.75&1\end{bmatrix}\begin{bmatrix}4&0\\0&6.75\end{bmatrix}\begin{bmatrix}1&0.75\\0&1\end{bmatrix}\)。本章会说明它其实就是一次回归。见第 00 册第 06 章。
• 多元正态的对数似然。 单变量是 \(-\frac12[\ln\sigma^2+a^2/\sigma^2]\)（省略常数）；多元版本把 \(\sigma^2\) 换成行列式 \(|\boldsymbol\Sigma|\)、把 \(a^2/\sigma^2\) 换成二次型 \(\boldsymbol a'\boldsymbol\Sigma^{-1}\boldsymbol a\)。估计所有多元 GARCH 都是最大化这个式子。见第 00 册第 06 章与 第 00 册第 07 章 概率中的分析工具。
• 条件期望与条件方差。 \(F_{t-1}\) 表示"截至 \(t-1\) 的全部信息"，\(\mathrm{Cov}(\boldsymbol a_t|F_{t-1})\) 就是"站在昨天收盘时看，今天收益的协方差矩阵"。见第 00 册第 07 章。
• 几个记号。 \(\odot\) 是逐元素相乘（Hadamard 积），如 \(\begin{bmatrix}1&2\\3&4\end{bmatrix}\odot\begin{bmatrix}5&6\\7&8\end{bmatrix}=\begin{bmatrix}5&12\\21&32\end{bmatrix}\)；\(\mathrm{diag}\{\cdot\}\) 是以括号内各数为对角元的对角阵；\(\boldsymbol\Sigma^{1/2}\) 是满足 \(\boldsymbol\Sigma^{1/2}(\boldsymbol\Sigma^{1/2})'=\boldsymbol\Sigma\) 的矩阵"平方根"（可以取 Cholesky 因子）。

怎么读这一章。 核心必读：10.1、10.2（EWMA）、10.4（两种重参数化，尤其 10.4.2 Cholesky 的回归解释）、10.6（DCC 与两步估计）、10.9（多资产 VaR）。10.3 的 DVEC 和 BEKK 只需记住形式和优缺点，看懂参数个数表即可。10.5 的 CCC 要读，时变相关的两种方法可以先看结论（"不同模型的相关预测差别很大"）。10.7、10.8 适合处理高维问题的读者，10.10 多元 t 第一次可以跳过公式。建议顺序：10.1 → 10.2 → 10.4 → 10.6 → 10.9 → 示例一、二 → 其余各节。

10.1 问题设定

记 \(\boldsymbol r_t=\boldsymbol\mu_t+\boldsymbol a_t\)，\(\boldsymbol\mu_t=E(\boldsymbol r_t|F_{t-1})\) 是条件均值，通常用第 08a 章的 VARMA（可加外生变量）刻画：

多元波动率指新息的条件协方差矩阵 \(\boldsymbol\Sigma_t=\mathrm{Cov}(\boldsymbol a_t|F_{t-1})\)。建模就是刻画 \(\boldsymbol\Sigma_t\) 的时间演化。两大难点：

• 维数灾难：\(\boldsymbol\Sigma_t\) 有 \(k(k+1)/2\) 个不同元素（\(k=5\) 时 15 个，\(k=100\) 时 5050 个），每个元素都要一条动态方程。
• 正定约束：每个时点的 \(\boldsymbol\Sigma_t\) 都必须正定，否则组合方差可能为负。

本章介绍几类相对简单、实用的模型，特别是允许时变相关系数的模型——它们可以用来估计时变的市场 β 和对冲比。

10.2 指数加权估计（EWMA）

最朴素的估计是等权样本协方差 \(\hat{\boldsymbol\Sigma}=\frac1{t-1}\sum_{j=1}^{t-1}\boldsymbol a_j\boldsymbol a_j'\)。为强调近期信息，用指数平滑：

权重和为 1。\(t\) 足够大时化为递推

这就是 EWMA 协方差。它只有一个参数，自动正定（只要初值正定），计算极快，是业界最常用的基准（RiskMetrics 日频取 \(\lambda=0.94\)）。若 \(\boldsymbol a_t|F_{t-1}\sim N(\boldsymbol 0,\boldsymbol\Sigma_t)\)，可以用极大似然联合估计 \(\lambda\) 和均值参数：

用递推得到的 \(\hat{\boldsymbol\Sigma}_t\) 代替 \(\boldsymbol\Sigma_t\)。

推导拆解：(10.2) 怎么变成递推式。 第一步，\(t\) 很大时 \(\lambda^{t-1}\approx0\)，归一化常数 \(\frac{1-\lambda}{1-\lambda^{t-1}}\approx1-\lambda\)，于是 \(\hat{\boldsymbol\Sigma}_t\approx(1-\lambda)\sum_{j\ge1}\lambda^{j-1}\boldsymbol a_{t-j}\boldsymbol a_{t-j}'\)。权重 \((1-\lambda)\lambda^{j-1}\) 是几何级数，加起来等于 1。 第二步，把 \(j=1\) 这一项单独拿出来：\((1-\lambda)\boldsymbol a_{t-1}\boldsymbol a_{t-1}'+\lambda\cdot\left[(1-\lambda)\sum_{j\ge2}\lambda^{j-2}\boldsymbol a_{t-j}\boldsymbol a_{t-j}'\right]\)。方括号里正是 \(\hat{\boldsymbol\Sigma}_{t-1}\)。 第三步，似然函数是多元正态对数密度逐日相加：\(\ln|\boldsymbol\Sigma_t|\) 惩罚"把协方差估得太大"，二次型 \(\boldsymbol a_t'\boldsymbol\Sigma_t^{-1}\boldsymbol a_t\) 惩罚"实际冲击相对估计太大"。哪个 \(\lambda\) 让两者平衡得最好，就选哪个。

金融直觉：\(\lambda=0.94\) 时，昨天的冲击权重 6%，一个月（约 21 个交易日）前的冲击权重只剩 \(0.06\times0.94^{20}\approx1.7\%\)；权重衰减到一半大约需要 \(\ln0.5/\ln0.94\approx11\) 天。所以 EWMA 本质上是"有效窗口约一两个月、越近越重要"的滚动协方差。它没有常数项把预测拉回长期均值，这就是正文说的"IGARCH 结构"。

例 10.1（恒生指数与日经 225）。 2006-01-04 至 2008-12-30 的日对数收益（%），只取两市同时开市的日子，713 个观测。单变量 GARCH(1,1)：

两者都接近 IGARCH（\(\alpha+\beta\approx1\)），这是次贷危机推高波动所致。二元 EWMA 的估计为 \(\hat\lambda=1-0.0695\approx0.9305\)，处于实务常见范围。EWMA 波动率比单变量 GARCH 更平滑，但形态相似。

EWMA 的局限也很明显：所有元素共用一个 \(\lambda\)；没有均值回复（\(\alpha+\beta=1\) 的 IGARCH 结构），长期预测停留在当前水平，不会回到无条件水平。

10.3 多元 GARCH 模型

综述见 Bauwens, Laurent & Rombouts (2004)。

10.3.1 对角 VEC（DVEC）模型

Bollerslev, Engle & Wooldridge (1988) 把 EWMA 推广为

\(\boldsymbol A_i,\boldsymbol B_j\) 对称，\(\odot\) 是 Hadamard（逐元素）积。二元 DVEC(1,1) 写开就是逐元素的 GARCH(1,1)：

EWMA 是其特例：\(\boldsymbol A_0=\boldsymbol 0\)，\(\boldsymbol A_1=(1-\lambda)\boldsymbol J\)，\(\boldsymbol B_1=\lambda\boldsymbol J\)（\(\boldsymbol J\) 为全 1 矩阵）。优点是简单；缺点是不保证 \(\boldsymbol\Sigma_t\) 正定，且不允许波动率之间的动态依赖（波动溢出）。

例 10.2（辉瑞与默克月收益，1965–2008，528 个观测）。

拟合的时变相关在 0.37 到 0.83 之间变动。（原书正文报告的 PFE 平方残差 \(Q(12)=12.35(0.42)\) 与软件输出不符，输出中 PFE 平方残差为 22.08(0.037)，12.35 是 MRK 标准化残差的统计量。）

10.3.2 BEKK 模型

Engle & Kroner (1995)：

\(\boldsymbol A\) 下三角。每一项都是"矩阵 × 半正定矩阵 × 转置"的形式，所以只要 \(\boldsymbol A\boldsymbol A'\) 正定，\(\boldsymbol\Sigma_t\) 就正定；\(\boldsymbol A_i,\boldsymbol B_j\) 的非对角元允许波动溢出。缺点：参数对滞后冲击和波动的影响没有直接解释（要做矩阵乘法才看得出）；参数个数 \(k^2(m+s)+k(k+1)/2\) 增长很快，经验上许多估计不显著。

推导拆解：为什么 DVEC 不保证正定、BEKK 保证？ BEKK：任取组合 \(\boldsymbol w\)，\(\boldsymbol w'\boldsymbol A_1(\boldsymbol a\boldsymbol a')\boldsymbol A_1'\boldsymbol w=(\boldsymbol w'\boldsymbol A_1\boldsymbol a)^2\ge0\)，是一个数的平方；\(\boldsymbol w'\boldsymbol B_1\boldsymbol\Sigma_{t-1}\boldsymbol B_1'\boldsymbol w=(\boldsymbol B_1'\boldsymbol w)'\boldsymbol\Sigma_{t-1}(\boldsymbol B_1'\boldsymbol w)\ge0\)，是另一个组合的方差。再加上严格为正的 \(\boldsymbol w'\boldsymbol A\boldsymbol A'\boldsymbol w\)，总和必为正。 DVEC：逐元素相乘没有这种"平方"结构。比如 \(\boldsymbol A_1\) 的非对角元比对角元大很多时，协方差的更新幅度会超过两个方差的更新幅度，隐含的相关系数可能超过 1，矩阵就不再正定。

金融直觉："波动溢出"指一个市场昨天的大波动会推高另一个市场今天的波动，比如美股暴跌后亚洲市场次日波动加大。DVEC 里 \(\sigma_{11,t}\) 只看 \(a_{1,t-1}^2\)，没有这个渠道；BEKK 通过矩阵乘法把各资产的冲击混在一起，能表达溢出，代价是参数多、难以直接解读。

例 10.3（同辉瑞/默克数据，BEKK(1,1)）。 ARCH 矩阵 \(\begin{bmatrix}0.213&0.063\\0.100&0.182\end{bmatrix}\)、GARCH 矩阵 \(\begin{bmatrix}0.909&-0.008\\-0.059&0.982\end{bmatrix}\)，非对角元均不显著（原书拟合方程中 GARCH 矩阵左上角印为 0.901，与参数表的 0.909 不一致）。与 DVEC 相比，BEKK 的时变相关波动更大。

参数个数比较（GARCH(1,1) 型）：

10.4 两种重参数化

利用 \(\boldsymbol\Sigma_t\) 的对称性，把它换成一组更易建模的量。

10.4.1 用相关系数

\(\boldsymbol\Sigma_t\) 的演化由 \(k\) 个条件方差和 \(k(k-1)/2\) 个条件相关系数决定，合起来记为 \(k(k+1)/2\) 维向量 \(\boldsymbol\Xi_t\)（\(k=2\) 时 \(\boldsymbol\Xi_t=(\sigma_{11,t},\sigma_{22,t},\rho_{21,t})'\)）。二元正态下对数似然为

优点是直接建模方差和相关；缺点是 \(k\ge3\) 时似然复杂，且要用约束保证 \(\boldsymbol\rho_t\) 正定。

10.4.2 Cholesky 分解

\(\boldsymbol\Sigma_t\) 正定，所以存在单位下三角阵 \(\boldsymbol L_t\) 和正对角阵 \(\boldsymbol G_t\) 使

二元情形。 \(\boldsymbol L_t=\begin{bmatrix}1&0\\q_{21,t}&1\end{bmatrix}\)，\(\boldsymbol G_t=\mathrm{diag}(g_{11,t},g_{22,t})\)，展开得

回归解释是这一参数化的精髓。考虑条件回归 \(a_{2t}=\beta a_{1t}+b_{2t}\)：最小二乘系数 \(\beta=\sigma_{21,t}/\sigma_{11,t}\)，残差方差 \(\sigma_{22,t}-\sigma_{21,t}^2/\sigma_{11,t}\)，且 \(b_{2t}\perp a_{1t}\)。所以

• \(g_{11,t}\) 是 \(a_{1t}\) 的条件方差；
• \(q_{21,t}\) 是 \(a_{2t}\) 对 \(a_{1t}\) 的条件回归系数——若 \(a_{1t}\) 是市场、\(a_{2t}\) 是个股，\(q_{21,t}\) 就是时变 β，也就是最小方差对冲比；
• \(g_{22,t}\) 是回归残差（特质）的条件方差。

Cholesky 分解等价于正交变换 \(b_{1t}=a_{1t}\)，\(b_{2t}=a_{2t}-q_{21,t}a_{1t}\)。

推导拆解：(10.13) 到 (10.14) 以及回归解释。 第一步，乘开 \(\boldsymbol L_t\boldsymbol G_t\boldsymbol L_t'\)：\(\begin{bmatrix}1&0\\q&1\end{bmatrix}\begin{bmatrix}g_{11}&0\\0&g_{22}\end{bmatrix}\begin{bmatrix}1&q\\0&1\end{bmatrix}=\begin{bmatrix}g_{11}&qg_{11}\\qg_{11}&q^2g_{11}+g_{22}\end{bmatrix}\)，逐元素对照 \(\boldsymbol\Sigma_t\) 就是 (10.13)。 第二步，从上往下解：第一个式子直接给 \(g_{11}\)；第二个除以 \(g_{11}\) 得 \(q\)；第三个移项得 \(g_{22}=\sigma_{22}-q^2\sigma_{11}=\sigma_{22}-\sigma_{21}^2/\sigma_{11}\)。 第三步，为什么这是回归：\(\mathrm{Cov}(b_{2t},a_{1t})=\mathrm{Cov}(a_{2t},a_{1t})-q\mathrm{Var}(a_{1t})=\sigma_{21}-q\sigma_{11}=0\)。所以 \(q=\sigma_{21}/\sigma_{11}\) 恰好是让残差与解释变量不相关的系数，即 OLS 的 β；\(g_{22}=\sigma_{22}(1-\rho^2)\) 就是"总方差 ×（1 − \(R^2\)）"，即残差方差。 数值例子（练习 4）：\(\sigma_{11}=4\)、\(\sigma_{22}=9\)、\(\sigma_{21}=3\)，得 \(q=0.75\)、\(g_{22}=9-9/4=6.75\)。读法：资产 1 涨 1%，资产 2 平均跟涨 0.75%；资产 2 的 9 个单位方差里，2.25 个来自资产 1，6.75 个是"特质"部分。

金融直觉：期货套保里的最小方差对冲比就是现货对期货回归的斜率。若 \(a_{1t}\) 是期货、\(a_{2t}\) 是现货，\(q_{21,t}\) 就是每天更新的最优对冲比，\(g_{22,t}\) 是对冲后剩余的基差风险。这一参数化让模型直接输出交易员关心的量。

一般情形：令 \(b_{1t}=a_{1t}\)，对 \(i=2,\dots,k\) 做逐次回归

即 \(\boldsymbol b_t=\boldsymbol L_t^{-1}\boldsymbol a_t\)，\(\mathrm{Cov}(\boldsymbol b_t)=\boldsymbol G_t\) 对角。由正交性，\(\sigma_{ii,t}=\sum_{v=1}^iq_{iv,t}^2g_{vv,t}\)，\(\sigma_{ij,t}=\sum_{v=1}^jq_{iv,t}q_{jv,t}g_{vv,t}\)（\(j<i\)，\(q_{vv,t}=1\)）。

似然大大简化。 \(|\boldsymbol L_t|=1\)，所以 \(|\boldsymbol\Sigma_t|=\prod_ig_{ii,t}\)；在正态假设下 \(\boldsymbol b_t\sim N(\boldsymbol 0,\boldsymbol G_t)\)，

不需要矩阵求逆，变成 \(k\) 个单变量似然之和。

推导拆解：(10.20) 从多元正态似然 \(-\frac12[\ln|\boldsymbol\Sigma_t|+\boldsymbol a_t'\boldsymbol\Sigma_t^{-1}\boldsymbol a_t]\) 来。 行列式：\(|\boldsymbol\Sigma_t|=|\boldsymbol L_t||\boldsymbol G_t||\boldsymbol L_t'|\)（乘积的行列式等于行列式的乘积）。三角阵的行列式是对角元之积，\(\boldsymbol L_t\) 对角全是 1，所以 \(|\boldsymbol L_t|=1\)，剩下 \(|\boldsymbol G_t|=\prod_ig_{ii,t}\)，取对数变成求和。 二次型：\(\boldsymbol\Sigma_t^{-1}=(\boldsymbol L_t')^{-1}\boldsymbol G_t^{-1}\boldsymbol L_t^{-1}\)，所以 \(\boldsymbol a_t'\boldsymbol\Sigma_t^{-1}\boldsymbol a_t=(\boldsymbol L_t^{-1}\boldsymbol a_t)'\boldsymbol G_t^{-1}(\boldsymbol L_t^{-1}\boldsymbol a_t)=\boldsymbol b_t'\boldsymbol G_t^{-1}\boldsymbol b_t=\sum_ib_{it}^2/g_{ii,t}\)，因为 \(\boldsymbol G_t\) 是对角阵。 两部分合起来，每个 \(i\) 贡献一项 \(\ln g_{ii,t}+b_{it}^2/g_{ii,t}\)，正好是单变量 GARCH 似然的形式。

三个优点：(1) 只要 \(g_{ii,t}>0\) 就正定，对 \(\ln g_{ii,t}\) 建模即可完全去掉约束；(2) 参数有清晰的回归解释；(3) 相关系数 \(\rho_{21,t}=q_{21,t}\sqrt{\sigma_{11,t}/\sigma_{22,t}}\)——即使 \(q_{21,t}\) 是常数，只要方差比时变，相关仍是时变的，这是它与相关系数参数化的主要区别。主要缺点：结果依赖分量排序（\(a_{1t}\) 不被变换），推断稍复杂。实践中通常把市场指数或最"外生"的资产排第一。

10.5 二元 GARCH 模型

多元 GARCH 用不含随机冲击的"确定性方程"描述 \(\boldsymbol\Xi_t\) 的演化。

10.5.1 常相关（CCC）模型

Bollerslev (1990) 假设 \(\rho_{21,t}\equiv\rho_{21}\)，只需对两个方差建模：

令 \(\boldsymbol\eta_t=\boldsymbol a_t^2-\boldsymbol\Xi_t^*\)，可写成 \(\boldsymbol a_t^2=\boldsymbol\alpha_0+(\boldsymbol\alpha_1+\boldsymbol\beta_1)\boldsymbol a_{t-1}^2+\boldsymbol\eta_t-\boldsymbol\beta_1\boldsymbol\eta_{t-1}\)，即平方新息的二元 ARMA(1,1)——单变量 GARCH(1,1) 的直接推广。由此：

推导拆解：这里 \(\boldsymbol\Xi_t^*=(\sigma_{11,t},\sigma_{22,t})'\) 是两个条件方差组成的向量，\(\boldsymbol a_t^2=(a_{1t}^2,a_{2t}^2)'\)。\(\boldsymbol\eta_t\) 是"实际平方冲击 − 条件方差"，即方差的预测误差，它的条件均值为 0。 第一步，在 (10.22) 中把 \(\boldsymbol\Xi_t^*\) 换成 \(\boldsymbol a_t^2-\boldsymbol\eta_t\)，\(\boldsymbol\Xi_{t-1}^*\) 换成 \(\boldsymbol a_{t-1}^2-\boldsymbol\eta_{t-1}\)： \(\boldsymbol a_t^2-\boldsymbol\eta_t=\boldsymbol\alpha_0+\boldsymbol\alpha_1\boldsymbol a_{t-1}^2+\boldsymbol\beta_1(\boldsymbol a_{t-1}^2-\boldsymbol\eta_{t-1})\)。 第二步，整理即得正文的式子。\(\boldsymbol a_{t-1}^2\) 的系数 \(\boldsymbol\alpha_1+\boldsymbol\beta_1\) 起 AR 作用，\(-\boldsymbol\beta_1\boldsymbol\eta_{t-1}\) 起 MA 作用。 第三步，无条件方差：两边取期望，平稳时 \(E(\boldsymbol a_t^2)=E(\boldsymbol a_{t-1}^2)\)，\(E(\boldsymbol\eta)=\boldsymbol 0\)，解得 \((\boldsymbol I-\boldsymbol\alpha_1-\boldsymbol\beta_1)^{-1}\boldsymbol\alpha_0\)，对应单变量的 \(\omega/(1-\alpha-\beta)\)。 验算例 10.4（对角 CCC）：恒生 \(0.079/(1-0.145-0.833)=0.079/0.022\approx3.6\)，即日波动率约 1.9%。

1. 若 \(\boldsymbol\alpha_1+\boldsymbol\beta_1\) 的特征值都在 (0,1) 内，\(\boldsymbol a_t^2\) 弱平稳，无条件方差为 \((\boldsymbol I-\boldsymbol\alpha_1-\boldsymbol\beta_1)^{-1}\boldsymbol\alpha_0\)，协方差 \(\rho_{21}\sigma_1\sigma_2\)。
2. \(\alpha_{12}=\beta_{12}=0\) 时 \(a_{1t}\) 的波动不依赖 \(a_{2t}\) 的过去波动，反之亦然；\(\boldsymbol\alpha_1,\boldsymbol\beta_1\) 都对角时退化为两个独立的单变量 GARCH。
3. 预测：\(\boldsymbol\Xi_h^*(1)=\boldsymbol\alpha_0+\boldsymbol\alpha_1\boldsymbol a_h^2+\boldsymbol\beta_1\boldsymbol\Xi_h^*\)，\(\boldsymbol\Xi_h^*(\ell)=\boldsymbol\alpha_0+(\boldsymbol\alpha_1+\boldsymbol\beta_1)\boldsymbol\Xi_h^*(\ell-1)\)；协方差预测为 \(\hat\rho_{21}\sqrt{\sigma_{11,h}(\ell)\sigma_{22,h}(\ell)}\)。

例 10.4（恒生/日经，对角 CCC）。 \(\sigma_{11,t}=0.079+0.145a_{1,t-1}^2+0.833\sigma_{11,t-1}\)，\(\sigma_{22,t}=0.054+0.105a_{2,t-1}^2+0.875\sigma_{22,t-1}\)，常相关 0.668。标准化残差的 \(Q_2(4)=17.29\)、\(Q_2(12)=48.21\)，拟合合理。

例 10.5（IBM 与 S&P 500 月收益，1926–1999）。

常相关 0.614（标准误 0.020）。非对角元显著，两收益的波动率之间存在反馈关系。

10.5.2 时变相关模型

常相关假设常被数据拒绝：IBM 与 S&P 500 的 120 个月滚动相关随时间明显变化，且近年下降。Tse (2000) 给出了检验常相关的 LM 统计量。

方法一：直接建模相关系数。 相关为正时用 logistic 变换保证 \(\rho\in(0,1)\)：

驱动项是上期标准化新息的乘积（"实现相关"），这是相关系数的 GARCH(1,1)。符号未知时用 Fisher 变换 \(q_t=\ln\frac{1+\rho_t}{1-\rho_t}\)。

白话解释：相关系数必须落在 \((-1,1)\) 内，直接对 \(\rho_t\) 写线性方程，某一天就可能算出 1.2。做法是让一个可以取任意实数的 \(q_t\) 去跟随线性方程，再用一个"压缩函数"把它映射到允许的区间。logistic 函数 \(e^q/(1+e^q)\) 把整条实数轴压到 \((0,1)\)：\(q=0\) 时 \(\rho=0.5\)，\(q=2\) 时 \(\rho\approx0.88\)，\(q\to-\infty\) 时 \(\rho\to0\)。Fisher 变换反解出来是 \(\rho=(e^q-1)/(e^q+1)\)，把实数轴压到 \((-1,1)\)。这和信用模型里用 logit 保证违约概率在 0 到 1 之间是同一个技巧。 驱动项 \(a_{1}a_{2}/\sqrt{\sigma_{11}\sigma_{22}}\) 是两个标准化冲击的乘积：两个市场昨天同向大幅波动时它为大的正数，把今天的相关往上推。

例 10.5 中估计为 \(q_t=-2.024+3.983\rho_{t-1}+0.088\frac{a_{1,t-1}a_{2,t-1}}{\sqrt{\sigma_{11,t-1}\sigma_{22,t-1}}}\)，全部高度显著。拟合相关平均 0.612（接近常相关 0.614），但随时间波动、近年变小；最大对数似然从 −3691.21 升至 −3679.64，显著改进。1 步预测：常相关模型 \(\hat{\boldsymbol\Sigma}(1)=\begin{bmatrix}71.09&21.83\\21.83&17.79\end{bmatrix}\)；时变相关模型 \(\begin{bmatrix}75.15&23.48\\23.48&24.70\end{bmatrix}\)，相关预测 0.545。

方法二：Cholesky 参数化。 二元时 \(\boldsymbol\Xi_t=(g_{11,t},g_{22,t},q_{21,t})'\)：

\(b_{1t}\) 是单变量 GARCH，\(b_{2t}\) 是二元 GARCH，\(q_{21,t}\) 自相关并以 \(a_{2,t-1}\) 为解释变量，似然用 (10.20)。例 10.5 的估计：

时变相关 \(\rho_t=q_{21,t}\sqrt{g_{11,t}}/\sqrt{g_{22,t}+q_{21,t}^2g_{11,t}}\)（10.30）。相关有很强的持续性（0.9915），路径比方法一更平滑，同样显示下降趋势；两种时变相关模型的最大似然都约为 −3672。但 1 步相关预测只有 0.203（\(\hat{\boldsymbol\Sigma}(1)=\begin{bmatrix}73.45&7.34\\7.34&17.87\end{bmatrix}\)），远低于前两个模型的 0.55 左右。不同多元波动率模型对相关的预测可以差很多，而相关预测直接决定组合风险和对冲比——模型选择本身就是一种风险。

10.6 动态条件相关（DCC）模型

基于 \(\boldsymbol\Sigma_t=\boldsymbol D_t\boldsymbol\rho_t\boldsymbol D_t\)，先用单变量 GARCH 刻画每个 \(\sigma_{ii,t}\)，再对相关矩阵 \(\boldsymbol\rho_t\) 建立极简的动态方程——所有相关对共用两个参数。

Tse & Tsui (2002)：

\(\bar{\boldsymbol\rho}\) 为单位对角正定阵（常取样本相关矩阵），\(\boldsymbol\psi_{t-1}\) 是用最近 \(m\) 期标准化新息算出的局部样本相关矩阵。\(\theta_i\ge0\)、\(\theta_1+\theta_2<1\) 时 \(\boldsymbol\rho_t\) 是正定矩阵的凸组合，自动正定。

Engle (2002)：

\(\epsilon_{it}=a_{it}/\sqrt{\sigma_{ii,t}}\) 是标准化新息，\(\bar{\boldsymbol Q}\) 是其无条件协方差，\(\boldsymbol J_t=\mathrm{diag}\{q_{11,t}^{-1/2},\dots,q_{kk,t}^{-1/2}\}\) 把 \(\boldsymbol Q_t\) 归一化为相关矩阵。\(\boldsymbol Q_t\) 本身就是标准化新息协方差的"GARCH(1,1)"。

两步估计（本教材补充，Engle 2002）。正态下 \(-2\ln L\) 可以拆成两部分：

于是第一步逐个拟合单变量 GARCH 得到 \(\boldsymbol\epsilon_t\)，第二步固定 \(\boldsymbol\epsilon_t\)、用样本协方差代替 \(\bar{\boldsymbol Q}\)，只对 \((\theta_1,\theta_2)\) 最大化相关部分。两步估计相合但非完全有效，好处是无论 \(k\) 多大，第二步都只有两个参数。

推导拆解：分解式怎么来的。正态下每期 \(-2\ell_t=k\ln2\pi+\ln|\boldsymbol\Sigma_t|+\boldsymbol a_t'\boldsymbol\Sigma_t^{-1}\boldsymbol a_t\)。 行列式：\(|\boldsymbol\Sigma_t|=|\boldsymbol D_t|\,|\boldsymbol\rho_t|\,|\boldsymbol D_t|\)，取对数得 \(2\ln|\boldsymbol D_t|+\ln|\boldsymbol\rho_t|\)。 二次型：\(\boldsymbol\Sigma_t^{-1}=\boldsymbol D_t^{-1}\boldsymbol\rho_t^{-1}\boldsymbol D_t^{-1}\)，而 \(\boldsymbol D_t^{-1}\boldsymbol a_t=\boldsymbol\epsilon_t\)，所以 \(\boldsymbol a_t'\boldsymbol\Sigma_t^{-1}\boldsymbol a_t=\boldsymbol\epsilon_t'\boldsymbol\rho_t^{-1}\boldsymbol\epsilon_t\)。 技巧在于"加一项再减一项"：\(\boldsymbol\epsilon_t'\boldsymbol\rho_t^{-1}\boldsymbol\epsilon_t=\boldsymbol\epsilon_t'\boldsymbol\epsilon_t+(\boldsymbol\epsilon_t'\boldsymbol\rho_t^{-1}\boldsymbol\epsilon_t-\boldsymbol\epsilon_t'\boldsymbol\epsilon_t)\)，其中 \(\boldsymbol\epsilon_t'\boldsymbol\epsilon_t=\boldsymbol a_t'\boldsymbol D_t^{-2}\boldsymbol a_t=\sum_ia_{it}^2/\sigma_{ii,t}\)。把只含方差的项归到第一组，就得到"\(k\) 个单变量 GARCH 似然之和"；剩下的只含 \(\boldsymbol\rho_t\)。 为什么"相合但非有效"：第二步把第一步的估计当成真值，忽略了第一步的估计误差，标准误需要修正；但两步都瞄准正确的参数，样本足够大时都会收敛到真值。

白话解释：\(\boldsymbol Q_t\) 的递推保证它正定，但对角元不一定等于 1，所以要用 \(\boldsymbol J_t\boldsymbol Q_t\boldsymbol J_t\) 归一化——相当于把协方差除以两边标准差得到相关系数，\(\rho_{ij,t}=q_{ij,t}/\sqrt{q_{ii,t}q_{jj,t}}\)。\(\theta_1\) 决定相关对新冲击的反应速度，\(\theta_2\) 决定持续性，\(\theta_1+\theta_2\) 越接近 1，相关偏离长期水平后回归得越慢。

共同缺点：\(\theta_1,\theta_2\) 是标量，所有相关对的动态完全相同，\(k\) 大时难以自圆其说。

Tsay (2006) 的推广：(1) 标准化新息服从多元 Student-t（10.8 节）；(2) 边际波动率含杠杆效应：

\(\boldsymbol\Lambda_i\) 对角，\(\boldsymbol A_{t-1}=\mathrm{diag}\{a_{1,t-1},\dots,a_{k,t-1}\}\)，\(\boldsymbol L_{t-1}\) 的对角元在 \(a_{i,t-1}<0\) 时等于 \(a_{i,t-1}\)、否则为 0；相关方程为 \(\boldsymbol\rho_t=(1-\theta_1-\theta_2)\hat{\boldsymbol\rho}+\theta_1\boldsymbol\psi_{t-1}+\theta_2\boldsymbol\rho_{t-1}\)（10.32）。

例 10.6（美元/欧元、日元/美元汇率与 IBM、Dell 日收益）。 1999–2004 年，1496 个观测，四个序列都有厚尾（超额峰度 2.0–6.2）。全模型对数似然 −9175.80；限制 \(\boldsymbol\Lambda_1=\lambda\boldsymbol I\) 后 −9176.62；最终限制模型（两只股票共用 ARCH 系数、两种汇率共用常数与 ARCH 系数）−9177.44，只用 9 个参数刻画四维收益的时变波动与相关：\(\lambda=0.9603\)，ARCH 系数 0.0248（汇率）与 0.0372（股票），多元 t 自由度 \(v=7.918\)，相关方程系数 0.9809 与 0.0137。

记号说明：原书表中把 0.98 记为 \(\theta_1\)，但按 (10.32) 的写法 \(\theta_1\) 乘的是局部相关 \(\boldsymbol\psi_{t-1}\)。从数值看，0.98 应是滞后相关矩阵 \(\boldsymbol\rho_{t-1}\) 的系数（高持续性），0.014 是 \(\boldsymbol\psi_{t-1}\) 的系数，这与 Engle 型 DCC 的典型估计（新息系数小、持续性系数接近 1）一致。

结论：LR 检验不能拒绝最终限制模型；两只股票 \(0.0372+0.9603=0.9975\approx1\)，呈 IGARCH；与 69 日滚动窗口相比，模型波动率对大冲击反应更快、相关更平滑。加入杠杆项后对数似然升至 −9169.04，杠杆效应只对股票显著（\(\boldsymbol\Lambda_3=\mathrm{diag}\{0,0,0.0159,0.0114\}\)），LR=15.16，p=0.0005。

10.7 高维波动率模型

10.7.1 序贯 Cholesky 建模

Cholesky 变换有嵌套结构：\(b_{it}\) 只依赖 \(b_{jt}\)（\(j<i\)），因此可以一个一个地加入资产：

1. 选最关心的市场指数（或股票），建单变量波动率模型；
2. 加入第二个序列，做正交变换，建二元模型（以第 1 步估计为初值）；
3. 加入第三个序列，同样处理（以二元估计为初值）；
4. 重复直到包含所有序列，每步都做模型检验。

经验表明这能大幅简化高维建模、显著减少计算时间。

例 10.7（S&P 500、Cisco、Intel 日收益，1991–1999，2275 个观测）。 排序为 (SP5, CSCO, INTC)。均值方程显示：指数不依赖两只个股的过去收益，Cisco 依赖指数滞后 2 期，Intel 依赖指数滞后 1 期——与第 08a 章 IBM 的结论一致。最终三维模型中，\(q_{ij,t}\) 都有很高的持续性（AR(1) 系数 0.82–0.98），例如

特征：指数的波动不依赖个股的过去波动，经 Cholesky 逆变换后个股的波动依赖市场的过去波动；收益波动大时相关系数上升，与"金融危机期间市场间相关上升"的国际实证一致。

10.7.2 通向多元随机波动率

令 \(\boldsymbol v_t=(\ln g_{11,t},\ln g_{22,t},\ln g_{33,t})'\)，\(\boldsymbol q_t=(q_{21,t},q_{31,t},q_{32,t})'\)，确定性模型 \(\boldsymbol v_t=\boldsymbol c_1+\boldsymbol\beta_1\boldsymbol v_{t-1}\)、\(\boldsymbol q_t=\boldsymbol c_2+\boldsymbol\beta_2\boldsymbol q_{t-1}\) 加上噪声项，就得到简单的多元随机波动率（SV）模型，需要用 MCMC 估计（第 12a、12b 章；Chib, Nardari & Shephard 1999）。

10.8 因子–波动率模型

另一条降维路线是用第 09 章的因子：

1. 对新息 \(\boldsymbol a_t\) 做 PCA，选出解释大部分变异的前几个主成分；
2. 对这些主成分建（少数几个）波动率模型；
3. 把每个 \(a_{it}\) 的波动率表示为主成分波动率的函数。

例 10.8（IBM 与 S&P 500 月收益）。 第一主成分解释约 82.5% 的方差，\(x_t=0.796r_{1t}+0.605r_{2t}\)。对它拟合 GARCH：\(\sigma_t^2=3.834+0.110a_{t-1}^2+0.825\sigma_{t-1}^2\)。再以 \(\sigma_t^2\) 为共同波动率因子：

相关方程与 10.5.2 节方法一几乎相同。均值方程无序列相关，但平方残差 \(Q_2^*(8)=61.95\)（p=0.0004）显示高阶条件异方差没处理好——单一因子只解释 82.5% 的方差。总体是参数更少的合理近似。原书用两步估计（先估因子波动再视为已知），简单但未必有效；若共同因子已知，可联合估计。（原书一处把系数 0.796 误印为 0.769。）

10.9 应用：多资产 VaR

投资者持有 Cisco 与 Intel 各 100 万美元多头，用日收益数据末端的 1 步预测和 5% 分位数计算日 VaR。组合 VaR（第 07a 章）为

金融直觉：这就是两资产组合标准差公式 \(\sigma_p^2=w_1^2\sigma_1^2+w_2^2\sigma_2^2+2\rho w_1w_2\sigma_1\sigma_2\) 两边乘以分位数 \(z_\alpha\) 后的形式：均值为零的正态下，\(\mathrm{VaR}_i=z_\alpha\times\)头寸\(_i\times\sigma_i\)，组合 VaR 同理，于是 \(z_\alpha\) 可以整体提出根号。\(\rho=1\) 时组合 VaR 等于简单加总（没有分散化），\(\rho<1\) 时小于加总。所以相关预测从 0.475 降到 0.345，抵消了 Cholesky 模型较高的单资产 VaR。

以第一行为例：\(\hat r_1=0.626\)，\(q_1=0.626-1.65\sqrt{4.152}=-2.736\%\)，\(\mathrm{VaR}_1=\$27{,}360\)；同理 \(\mathrm{VaR}_2=\$38{,}840\)，代入公式得 $57,117（本教材已核算三行组合 VaR）。三种方法结果相近，差距约 1100 美元。注意：单资产 VaR 中包含了非零均值，而上述合成公式严格成立需要均值为零（或对组合整体计算分位数），所以这里只是近似。Cholesky 模型的两个单资产 VaR 之和最大，但相关预测低得多（0.345），组合 VaR 反而居中——相关预测对组合风险的影响不亚于方差预测。

10.10 多元 Student-t 分布

多元高斯新息常无法刻画收益的厚尾。\(k\) 维、自由度 \(v\)、均值 0、尺度阵 \(\boldsymbol I\) 的多元 t 密度为

各分量方差为 \(v/(v-2)\)，标准化 \(\boldsymbol\epsilon_t=\sqrt{(v-2)/v}\,\boldsymbol x\) 使协方差为 \(\boldsymbol I\)。在波动率模型中 \(\boldsymbol a_t=\boldsymbol\Sigma_t^{1/2}\boldsymbol\epsilon_t\)，

用 Cholesky 分解时 \(|\boldsymbol\Sigma_t|^{1/2}=\prod_jg_{jj,t}^{1/2}\)，\(\boldsymbol a_t'\boldsymbol\Sigma_t^{-1}\boldsymbol a_t=\sum_jb_{jt}^2/g_{jj,t}\)，无需矩阵求逆。但要注意：多元 t 下 \(b_{jt}\) 不相关但不独立，似然不能像正态那样拆成单变量似然的乘积。

白话解释：多元 t 可以理解为"多元正态再乘上一个所有分量共用的随机放大倍数"——市场平静的日子倍数小，恐慌的日子倍数大。因为这个倍数是共用的，一个分量出现极端值时，其他分量也更可能出现极端值，所以即使它们的相关系数为 0，也不独立。\(\Gamma(\cdot)\) 是伽马函数（阶乘向实数的推广，\(\Gamma(n)=(n-1)!\)），在这里只起归一化常数的作用，让密度积分为 1。 对风险管理的含义：多元 t 下资产会"一起出极端值"，这正是危机时分散化失效的统计表现；自由度 \(v\) 越小，这种尾部共振越强。例 10.6 估计 \(v\approx7.9\)，尾部明显厚于正态。

估计说明（原书附录）。 原书用 RATS 实现例 10.5 的常相关、logistic 时变相关和 Cholesky 模型，以及例 10.7 的三维 Cholesky 模型：用 nonlin 声明参数、frml 写方差递推和对数似然、maximize(method=bhhh) 求解。本章的 Python 示例用 scipy.optimize 实现相同的思路。

量化实战

应用场景

风险管理。 组合 VaR/ES 的核心输入是 \(\boldsymbol\Sigma_t\) 的 1 步（或多步）预测。EWMA 是行业标准的基准；DCC 允许均值回复和时变相关，是更精细的选择。压力测试应使用危机情景下的条件相关，而不是长期平均相关——10.7 节"波动大时相关上升"意味着分散化恰在最需要的时候失效。

组合优化与资产配置。 最小方差、风险平价、均值–方差组合的权重都由 \(\boldsymbol\Sigma_t\) 决定，用 DCC 预测替代滚动样本协方差可以更快适应市场状态。高维时（数百只股票）一般先用第 09 章的因子模型把问题降到几十个因子，再对因子协方差用 EWMA 或 DCC，特质方差用单变量 GARCH——这就是商业风险模型的做法。

动态对冲与时变 β。 Cholesky 参数中的 \(q_{21,t}=\sigma_{21,t}/\sigma_{11,t}\) 就是时变回归系数：个股对市场的时变 β、期货最优对冲比（最小方差对冲比 \(h^*=\mathrm{Cov}(\Delta S,\Delta F)/\mathrm{Var}(\Delta F)\)，见第 08 册）、配对交易中动态的对冲比 \(\gamma_t\)。

工程注意。 排序依赖性：把市场指数放第一位；数值稳定性：Cholesky 似然无需求逆，对 \(\ln g\) 建模免去约束；估计初值：用序贯法逐个加入资产。

Python 示例一：DCC 两步估计、模型比较与组合 VaR

模拟 3 资产 DCC-GARCH(1,1) 日收益（真值 \(\theta_1=0.04\)、\(\theta_2=0.94\)），比较 EWMA、CCC、DCC 和 60 日滚动窗口对真实相关路径的估计精度，最后计算下一日的组合 VaR。

import numpy as np
from scipy.optimize import minimize, minimize_scalar
from arch import arch_model

rng = np.random.default_rng(10)

# ---------- 1. 模拟 3 资产 DCC-GARCH(1,1) 日收益（%），Engle (2002) 设定 ----------
T, k = 2500, 3
omega = np.array([0.02, 0.03, 0.05]); alp = np.array([0.08, 0.10, 0.06]); bet = np.array([0.90, 0.87, 0.92])
Qbar = np.array([[1.0, 0.5, 0.3], [0.5, 1.0, 0.6], [0.3, 0.6, 1.0]])
th1, th2 = 0.04, 0.94                                     # 真值：Q_t = (1-θ1-θ2)Q̄ + θ1 εε' + θ2 Q_{t-1}
h = omega / (1 - alp - bet); Q = Qbar.copy(); eps_prev = np.zeros(k); a_prev = np.zeros(k)
A = np.zeros((T, k)); Rtrue = np.zeros((T, k, k))
for t in range(T):
    h = omega + alp * a_prev**2 + bet * h
    Q = (1 - th1 - th2) * Qbar + th1 * np.outer(eps_prev, eps_prev) + th2 * Q
    J = np.diag(1 / np.sqrt(np.diag(Q))); Rt = J @ Q @ J
    eps = np.linalg.cholesky(Rt) @ rng.standard_normal(k)
    a = np.sqrt(h) * eps
    A[t], Rtrue[t] = a, Rt
    eps_prev, a_prev = eps, a

# ---------- 2. EWMA 协方差：RiskMetrics λ=0.94 与极大似然估计 λ ----------
def ewma_path(A, lam):
    S = np.cov(A[:50].T); out = np.zeros((len(A), k, k))
    for t in range(len(A)):
        out[t] = S
        S = (1 - lam) * np.outer(A[t], A[t]) + lam * S
    return out
def ewma_nll(lam):
    S = ewma_path(A, lam)
    sign, logdet = np.linalg.slogdet(S)
    quad = np.einsum("ti,tij,tj->t", A, np.linalg.inv(S), A)
    return 0.5 * (logdet + quad)[50:].sum()
lam_hat = minimize_scalar(ewma_nll, bounds=(0.85, 0.999), method="bounded").x
print(f"EWMA 极大似然 λ = {lam_hat:.4f}")

# ---------- 3. DCC 两步估计：先逐个 GARCH(1,1)，再估相关方程 ----------
H = np.zeros((T, k)); E = np.zeros((T, k))
for i in range(k):
    res = arch_model(A[:, i], mean="Zero", vol="GARCH", p=1, q=1).fit(disp="off")
    H[:, i] = res.conditional_volatility**2
    E[:, i] = A[:, i] / res.conditional_volatility
    w, a1, b1 = res.params[["omega", "alpha[1]", "beta[1]"]]
    print(f"资产{i+1} GARCH: ω={w:.3f} α={a1:.3f} β={b1:.3f}  (真值 {omega[i]}, {alp[i]}, {bet[i]})")

Qb = np.cov(E.T, bias=True)
def dcc_path(theta, E):
    t1, t2 = theta; Qt = Qb.copy(); R = np.zeros((len(E), k, k))
    for t in range(len(E)):
        d = 1 / np.sqrt(np.diag(Qt)); R[t] = Qt * np.outer(d, d)
        Qt = (1 - t1 - t2) * Qb + t1 * np.outer(E[t], E[t]) + t2 * Qt
    return R
def dcc_nll(theta):
    if theta[0] < 0 or theta[1] < 0 or theta.sum() >= 0.999:
        return 1e10
    R = dcc_path(theta, E)
    _, logdet = np.linalg.slogdet(R)
    quad = np.einsum("ti,tij,tj->t", E, np.linalg.inv(R), E)
    return 0.5 * (logdet + quad - (E**2).sum(1)).sum()   # 似然的相关部分
fit = minimize(dcc_nll, x0=[0.03, 0.90], method="Nelder-Mead")
print(f"DCC: θ1(εε')={fit.x[0]:.4f}, θ2(Q_t-1)={fit.x[1]:.4f}   (真值 {th1}, {th2})")

# 常相关（CCC）是 θ1=θ2=0 的特例：似然比检验
R_ccc = np.corrcoef(E.T)
_, ld = np.linalg.slogdet(R_ccc)
nll_ccc = 0.5 * (T * ld + np.einsum("ti,ij,tj->", E, np.linalg.inv(R_ccc), E) - (E**2).sum())
print(f"LR(CCC vs DCC) = {2 * (nll_ccc - fit.fun):.1f}  （χ²_2 的 1% 临界值 9.21）")

# ---------- 4. 相关路径的估计精度：DCC / CCC / EWMA / 60 日滚动 ----------
R_dcc = dcc_path(fit.x, E)
S_ew = ewma_path(A, 0.94)
R_ew = S_ew / np.sqrt(np.einsum("tii,tjj->tij", S_ew, S_ew))
R_roll = np.array([np.corrcoef(A[max(0, t - 60):t].T) if t >= 60 else R_ccc for t in range(T)])
iu = np.triu_indices(k, 1)
for name, Rh in (("DCC", R_dcc), ("CCC", np.broadcast_to(R_ccc, Rtrue.shape)),
                 ("EWMA 0.94", R_ew), ("60日滚动", R_roll)):
    err = (Rh[:, iu[0], iu[1]] - Rtrue[:, iu[0], iu[1]])[100:]
    print(f"{name:9s}: 相关系数 RMSE = {np.sqrt((err**2).mean()):.4f}")
rt = Rtrue[:, 1, 2]
print(f"真实 ρ23 的波动范围: {rt.min():.2f} ~ {rt.max():.2f}")

# ---------- 5. 下一日协方差预测与组合 VaR（每个资产多头 100 万） ----------
hT = np.array([arch_model(A[:, i], mean="Zero", vol="GARCH", p=1, q=1).fit(disp="off")
               .forecast(horizon=1).variance.values[-1, 0] for i in range(k)])
t1, t2 = fit.x
QT = Qb.copy()
for t in range(T):
    QT = (1 - t1 - t2) * Qb + t1 * np.outer(E[t], E[t]) + t2 * QT
d = 1 / np.sqrt(np.diag(QT)); RT = QT * np.outer(d, d)
Sig = np.sqrt(np.outer(hT, hT)) * RT
pos = np.array([1e6, 1e6, 1e6])
VaR_i = 1.645 * np.sqrt(hT) / 100 * pos
VaR_p = 1.645 * np.sqrt(pos @ Sig @ pos) / 100
print("单资产 VaR:", VaR_i.round(0), " 简单加总:", VaR_i.sum().round(0))
print(f"DCC 组合 VaR = {VaR_p:,.0f}；  用 CCC 相关 = "
      f"{1.645*np.sqrt(pos @ (np.sqrt(np.outer(hT,hT))*R_ccc) @ pos)/100:,.0f}")

输出：

EWMA 极大似然 λ = 0.9593
资产1 GARCH: ω=0.021 α=0.083 β=0.895  (真值 0.02, 0.08, 0.9)
资产2 GARCH: ω=0.034 α=0.106 β=0.861  (真值 0.03, 0.1, 0.87)
资产3 GARCH: ω=0.067 α=0.063 β=0.909  (真值 0.05, 0.06, 0.92)
DCC: θ1(εε')=0.0326, θ2(Q_t-1)=0.9545   (真值 0.04, 0.94)
LR(CCC vs DCC) = 313.5  （χ²_2 的 1% 临界值 9.21）
DCC      : 相关系数 RMSE = 0.0190
CCC      : 相关系数 RMSE = 0.1489
EWMA 0.94: 相关系数 RMSE = 0.0780
60日滚动    : 相关系数 RMSE = 0.1085
真实 ρ23 的波动范围: 0.12 ~ 0.81
单资产 VaR: [17042. 19202. 23228.]  简单加总: 59472.0
DCC 组合 VaR = 51,296；  用 CCC 相关 = 47,679

• EWMA 的 ML 估计 \(\hat\lambda=0.959\)，比 RiskMetrics 的 0.94 更"慢"——与例 10.1 一样，\(\lambda\) 应根据数据估计，0.94 只是经验值。
• 两步 DCC 恢复了真实参数（0.033/0.955 对 0.04/0.94），持续性 \(\theta_1+\theta_2=0.987\) 与真值 0.98 接近。
• CCC 被强烈拒绝（LR=313.5）。严格地说，\(\theta_1=0\) 时 \(\theta_2\) 不可识别，LR 统计量的分布非标准，这里的卡方临界值只作参考；但统计量如此之大，结论不受影响。
• 相关路径精度：DCC 的 RMSE（0.019）远低于 EWMA（0.078）、60 日滚动（0.109）和 CCC（0.149）。真实 \(\rho_{23}\) 在 0.12 到 0.81 之间波动，常相关假设会严重误估某些时段的组合风险。
• 组合 VaR：三个单资产 VaR 简单加总为 59,472，DCC 组合 VaR 为 51,296（分散化收益），而用常相关算出 47,679——当前条件相关高于长期平均，CCC 会低估约 7% 的风险。

Python 示例二：Cholesky 时变 β 与动态对冲

模拟市场与个股：个股的真实 β 在 0.6–1.7 之间缓慢变动，并在第 2000 天上跳 0.3（结构转变）。用四种方法估计时变 β 并用于对冲：静态 β、250 日滚动 OLS、EWMA 协方差比、Tsay 的 Cholesky GARCH。Cholesky 模型的 \(q_{21,t}\) 方程用两种驱动项：原书 (10.28) 的 \(a_{2,t-1}\)，以及本教材的一个变体——用"β 的意外" \(a_{1,t-1}b_{2,t-1}/g_{11,t-1}\) 驱动（它是对数似然对 \(q_{21}\) 的得分方向，作用类似递推最小二乘的修正量）。

import numpy as np
from scipy.optimize import minimize

rng = np.random.default_rng(4)

# ---------- 1. 模拟：市场 GARCH(1,1) + 个股，beta 缓慢时变（结构转变） ----------
T = 3000
beta_true = 1.0 + 0.4 * np.sin(2 * np.pi * np.arange(T) / 1500) + np.where(np.arange(T) > 2000, 0.3, 0.0)
m = np.zeros(T); s = np.zeros(T); h = 1.0; g = 1.5; m_prev = 0.0; u_prev = 0.0
for t in range(T):
    h = 0.02 + 0.08 * m_prev**2 + 0.90 * h            # 市场条件方差
    g = 0.05 + 0.05 * u_prev**2 + 0.92 * g            # 个股特质条件方差
    m[t] = np.sqrt(h) * rng.standard_normal()
    u = np.sqrt(g) * rng.standard_normal()
    s[t] = beta_true[t] * m[t] + u
    m_prev, u_prev = m[t], u

# ---------- 2. Tsay 的 Cholesky 时变相关 GARCH，式 (10.28) 的简化版 ----------
#  b1 = a1 = m；b2 = a2 - q21_t a1
#  g11_t = α10 + α11 b1^2 + β11 g11；  g22_t = α20 + α22 b2^2 + β22 g22
#  q21 方程两种驱动项：
#    "tsay"  : q21_t = γ0 + γ1 q21_{t-1} + γ2 a2_{t-1}                 （原书 10.28）
#    "score" : q21_t = γ0 + γ1 q21_{t-1} + γ2 a1_{t-1} b2_{t-1}/g11_{t-1}（本教材的变体：
#              驱动项是 beta 的"意外"，类似递推最小二乘的修正量）
def filt(p, a1, a2, driver):
    a10, a11, b11, c0, c1, c2, a20, a22, b22 = p
    n = len(a1); g11 = np.var(a1); g22 = np.var(a2); q = c0 / (1 - c1)
    G11, G22, Q, B2 = (np.zeros(n) for _ in range(4))
    for t in range(n):
        if t > 0:
            g11 = a10 + a11 * a1[t-1]**2 + b11 * g11
            x = a2[t-1] if driver == "tsay" else a1[t-1] * B2[t-1] / G11[t-1]
            q = c0 + c1 * q + c2 * x
            g22 = a20 + a22 * B2[t-1]**2 + b22 * g22
        G11[t], Q[t], G22[t] = g11, q, g22
        B2[t] = a2[t] - q * a1[t]
    return G11, G22, Q, B2
def nll(p, driver):
    if (np.array(p)[[0, 1, 2, 6, 7, 8]] < 0).any() or p[1] + p[2] >= 1 or p[7] + p[8] >= 1 or abs(p[4]) >= 1:
        return 1e10
    G11, G22, Q, B2 = filt(p, m, s, driver)
    return 0.5 * (np.log(G11) + m**2 / G11 + np.log(G22) + B2**2 / G22).sum()   # 式 (10.20)
p0 = [0.05, 0.05, 0.9, 0.05, 0.95, 0.0, 0.1, 0.05, 0.9]
names = ["α10", "α11", "β11", "γ0", "γ1", "γ2", "α20", "α22", "β22"]
q21 = {}
for drv in ("tsay", "score"):
    fit = minimize(nll, p0, args=(drv,), method="Nelder-Mead",
                   options={"maxiter": 20000, "maxfev": 20000, "xatol": 1e-6, "fatol": 1e-6})
    fit = minimize(nll, fit.x, args=(drv,), method="Nelder-Mead", options={"maxiter": 20000, "maxfev": 20000})
    print(f"[{drv}] -logL={fit.fun:.1f}", {n: round(float(v), 4) for n, v in zip(names, fit.x)})
    G11, G22, q21[drv], B2 = filt(fit.x, m, s, drv)

# ---------- 3. 其他时变 beta 估计：EWMA 协方差比、250 日滚动 OLS ----------
lam = 0.97; c11 = np.var(m[:50]); c21 = np.cov(m[:50], s[:50])[0, 1]
beta_ew = np.zeros(T)
for t in range(T):
    beta_ew[t] = c21 / c11                           # 用 t-1 及以前信息
    c11 = lam * c11 + (1 - lam) * m[t]**2
    c21 = lam * c21 + (1 - lam) * m[t] * s[t]
beta_roll = np.full(T, np.nan)
for t in range(250, T):
    x, y = m[t-250:t], s[t-250:t]
    beta_roll[t] = (x @ y) / (x @ x)
beta_static = (m[:250] @ s[:250]) / (m[:250] @ m[:250])

# ---------- 4. 比较：beta 跟踪误差与对冲效果（对冲后收益方差） ----------
idx = np.arange(250, T)
print(f"{'方法':16s}{'beta RMSE':>10s}{'对冲后方差':>12s}")
for name, b in (("静态(前250日)", np.full(T, beta_static)), ("250日滚动OLS", beta_roll),
                ("EWMA λ=0.97", beta_ew), ("Cholesky(10.28)", q21["tsay"]),
                ("Cholesky 变体", q21["score"])):
    rmse = np.sqrt(((b[idx] - beta_true[idx])**2).mean())
    hedged = s[idx] - b[idx] * m[idx]
    print(f"{name:16s}{rmse:10.3f}{hedged.var():12.3f}")
print(f"不对冲方差 {s[idx].var():.3f}；完美对冲(真 beta)方差 {(s[idx]-beta_true[idx]*m[idx]).var():.3f}")
q = q21["score"]
rho = q * np.sqrt(G11) / np.sqrt(G22 + q**2 * G11)   # 式 (10.30)
print(f"Cholesky 模型隐含相关系数范围 {rho[idx].min():.2f} ~ {rho[idx].max():.2f}")

输出：

[tsay] -logL=3348.3 {'α10': 0.0364, 'α11': 0.1051, 'β11': 0.8521, 'γ0': 0.0311, 'γ1': 0.9727, 'γ2': -0.0058, 'α20': 0.0355, 'α22': 0.0359, 'β22': 0.9433}
[score] -logL=3322.8 {'α10': 0.0356, 'α11': 0.1034, 'β11': 0.8546, 'γ0': 0.0035, 'γ1': 0.9969, 'γ2': 0.0104, 'α20': 0.0402, 'α22': 0.0386, 'β22': 0.9375}
方法               beta RMSE       对冲后方差
静态(前250日)            0.290       1.737
250日滚动OLS            0.165       1.709
EWMA λ=0.97          0.187       1.717
Cholesky(10.28)      0.291       1.735
Cholesky 变体          0.150       1.701
不对冲方差 2.752；完美对冲(真 beta)方差 1.691
Cholesky 模型隐含相关系数范围 0.26 ~ 0.89

• 原书设定 (10.28) 的 \(q_{21,t}\) 以 \(a_{2,t-1}\) 为驱动，而 \(a_{2,t-1}\) 本身几乎不含"β 变了"的信息，所以估计出 \(\gamma_2\approx0\)，\(q_{21,t}\) 退化为常数，跟踪效果与静态 β 相同。这并不说明 Cholesky 参数化不好，而是说明驱动项要选对：在 IBM/S&P 500 这类数据上它够用，在 β 发生系统性漂移时就不够。
• 变体以 \(a_{1,t-1}b_{2,t-1}/g_{11,t-1}\) 驱动，参数个数相同，但 \(-\ln L\) 降低 25.5，β 的跟踪误差（0.150）在四种方法中最小，对冲后方差（1.701）最接近用真实 β 的理论下限（1.691）。
• 不对冲时方差 2.75，任何合理的对冲都去掉约 38% 的方差；方法之间的差别相对较小，但在大资金、高杠杆的对冲账户中，这一差别会直接变成损益。
• Cholesky 模型隐含的相关系数在 0.26–0.89 之间变化——β 和方差比同时时变，相关随之时变（10.4.2 节的第三个优点）。

本章小结

多元波动率模型刻画条件协方差矩阵 \(\boldsymbol\Sigma_t\)，难点在于维数灾难和正定约束。EWMA 只有一个参数，是实务基准；DVEC 逐元素 GARCH，简单但不保证正定、没有溢出；BEKK 保证正定、允许溢出，但参数多且难解释。两种重参数化把问题变得可控：\(\boldsymbol D_t\boldsymbol\rho_t\boldsymbol D_t\) 把方差和相关分开，催生了 CCC、logistic 时变相关和 DCC；Cholesky 分解 \(\boldsymbol L_t\boldsymbol G_t\boldsymbol L_t'\) 等价于逐次回归正交化，参数是时变回归系数和残差方差，似然可拆分、无约束，但依赖排序，并可序贯推广到高维。DCC 用两步法估计，第二步无论维数多高都只有两个参数，是最常用的多元波动率模型。实证规律：波动大时相关上升；个股波动受市场过去波动影响，反之不然；股票日波动接近 IGARCH 且有杠杆效应；不同模型的相关预测可能差别很大，直接影响 VaR 和对冲比。

练习

基础

1. 证明 EWMA 递推式中，若初值 \(\hat{\boldsymbol\Sigma}_1\) 正定，则所有 \(\hat{\boldsymbol\Sigma}_t\) 正定。 提示：正定矩阵与半正定矩阵的正系数组合仍正定。
2. 写出 DVEC(1,1) 退化为 EWMA 的参数条件，并说明为什么 DVEC 一般不保证正定。
3. 由 (10.13) 推导 (10.14)，并说明 \(g_{22,t}\) 的回归含义。
4. 给定 \(\sigma_{11}=4\)、\(\sigma_{22}=9\)、\(\sigma_{21}=3\)，求 Cholesky 参数 \(g_{11},q_{21},g_{22}\) 与相关系数；若改变排序，\(q_{21}\) 变为多少？ 提示：\(g_{11}=4\)，\(q_{21}=0.75\)，\(g_{22}=6.75\)，\(\rho=0.5\)；换序后 \(q_{21}=3/9=1/3\)。
5. 用 10.9 节表中数据核算第三行的组合 VaR。 提示：\(\sqrt{30504^2+39512^2+2\times0.345\times30504\times39512}\approx57{,}648\)。

进阶

6. 证明 DCC 的 \(-2\ln L\) 可以分解为"波动部分 + 相关部分"，并说明为什么两步估计相合但非有效。
7. 推导常相关二元 GARCH(1,1) 的无条件方差 \((\boldsymbol I-\boldsymbol\alpha_1-\boldsymbol\beta_1)^{-1}\boldsymbol\alpha_0\)，并给出 \(\ell\) 步方差预测的闭式表达。
8. 修改示例一：把新息改为多元 t（\(v=6\)），比较高斯两步 DCC 与把相关部分似然换成多元 t 的估计结果；计算 1% 组合 VaR 并比较。
9. 修改示例一：用 50 个资产（\(\bar{\boldsymbol Q}\) 随机生成），比较 DCC 两步法与 BEKK 的可行性（参数个数、估计时间）。
10. 在示例二中把 EWMA 的 \(\lambda\) 改为用对冲后方差最小化来选择（在形成期上），再在后半段样本外评估；与 Cholesky 变体比较。

原书推荐习题：10.10（三种方法计算多元 VaR 并比较，直接对应风险管理实务）、10.6（三维 DCC，\(\boldsymbol\rho\) 取样本相关阵）、10.7 + 10.8（同一数据比较 logistic 时变相关与 Cholesky 时变相关）、10.1（EWMA 估计 \(\lambda\)）、10.2、10.3（DVEC 与 BEKK）。

原书对照

（原书印刷页码约等于 PDF 页码减 20。）