第 04 章 风险模型

风险模型给出未来一段时间资产收益的协方差矩阵 \(\boldsymbol\Sigma\)。它有三个用途：组合优化（要用 \(\boldsymbol\Sigma^{-1}\)）、风险预测（组合波动、VaR）、风险归因（收益与风险来自哪些因子）。难点在于股票数 \(N\) 往往与样本长度 \(T\) 同一量级，样本协方差的误差在求逆时被放大。本章的结论是：样本协方差在 \(q=N/T\) 不小时会让最小方差组合的真实风险高出 \(1/(1-q)\) 倍，而模型自己报告的风险却只有真实的 \((1-q)\) 倍左右；收缩、因子结构和随机矩阵去噪都是在给估计加结构，其中与真实结构相符的因子模型效果最好；波动和相关随时间变化，要用 EWMA 或 DCC 跟踪。

学习目标

1. 理解样本协方差的维度灾难：特征值散开、条件数爆炸、最小方差组合的风险被低估。
2. 会实现 Ledoit–Wolf 收缩，知道收缩强度的含义和它的局限。
3. 会构建统计因子模型（PCA）与基本面因子模型（截面回归），用 Woodbury 公式高效求逆，用偏差统计量检验风险预测。
4. 会用 EWMA 和 DCC 跟踪时变的波动与相关，并用偏差统计量和 QLIKE 评价预测。
5. 会用 Marchenko–Pastur 边界区分信号与噪声特征值，实现特征值截断去噪，知道它在真实数据中的失效情形。

读前导读

这一章在解决什么问题。 风险模型对应的是风控部门的日常：给每个组合算出"预期波动是多少、风险来自哪里、离风险预算还有多远"。你在 CFA 里学过组合方差 \(\sigma_p^2=\boldsymbol w'\boldsymbol\Sigma\boldsymbol w\)，也学过 VaR；那里协方差矩阵 \(\boldsymbol\Sigma\) 是"给定的"。本章要回答的是：\(\boldsymbol\Sigma\) 本身从哪里来、估得准不准、估不准会怎样。

对应到实际工作：

• 4.2 维度灾难：风控部门如果直接用过去一年 250 天的数据算 500 只股票的协方差矩阵，会得到一个"看起来很精确、其实大量是噪声"的矩阵。用它做优化，组合会集中押在"样本里碰巧风险很低"的方向上，真实风险远超报告值。这是"模型风险"的典型案例。
• 4.3 收缩：相当于给估计加一个"审慎调整"：不完全相信样本，往一个保守的基准拉一点。和信用评级里把个别违约率估计向行业平均"收缩"是同一种思路。
• 4.4 因子风险模型：就是 Barra、Axioma 这类商业风险模型的原理。风控报告里的"行业暴露、风格暴露、特质风险"拆分，就来自这个结构。风险预算按因子分配，也是基于它。
• 4.5 EWMA 与 DCC：RiskMetrics 用的就是 EWMA（\(\lambda=0.94\)）。它回答"危机来了，模型多快能反应过来"。
• 4.6 随机矩阵去噪：区分"真实的共同风险来源"和"样本噪声"的一种统计工具。

需要先想起来的数学。 本章是全册线性代数最密集的一章，下面几项是关键。

1. 矩阵求逆与线性方程组：\(\boldsymbol\Sigma^{-1}\boldsymbol 1\) 就是解方程组 \(\boldsymbol\Sigma\boldsymbol x=\boldsymbol 1\)。代码里用 np.linalg.solve 而不是先求逆再相乘，因为前者更快更稳。见 第 00 册第 06 章 线性代数速成。
2. 特征值与特征向量：对协方差矩阵，特征向量是一组互相正交的"组合方向"，对应的特征值是这个方向组合的方差。最大特征值通常对应"市场组合"，最小特征值对应"风险最低的对冲组合"。例：两只股票方差都是 1、相关 0.8，则特征值为 \(1.8\)（等权多头方向）和 \(0.2\)（一多一空方向）。
3. 条件数：最大特征值除以最小特征值。条件数越大，矩阵越接近"不可逆"，求逆时误差被放大得越厉害。上面的例子条件数为 9；相关 0.99 时变为 199。
4. 迹 \(\mathrm{tr}\) 与 Frobenius 范数 \(\|\cdot\|_F\)：迹是对角元之和，等于所有特征值之和（即总方差）。Frobenius 范数是"所有元素平方和再开方"，衡量两个矩阵整体差多少。
5. 正定矩阵：对任何非零 \(\boldsymbol w\) 都有 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w>0\)，即"任何组合的方差都为正"。合格的协方差矩阵至少要半正定。
6. 拉格朗日乘子：GMV 公式来自"在 \(\boldsymbol 1'\boldsymbol w=1\) 约束下最小化 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w\)"。见 第 00 册第 05 章 多元微积分与优化。

怎么读这一章。 必读 4.1、4.2（结论和两个近似式）、4.4（因子模型结构与偏差统计量）。4.3 的 LW 公式第一次只需理解"\(\delta\) 由'噪声有多大'和'离目标有多远'决定"。4.5 读 EWMA 和结论即可，DCC 的似然函数代码可以跳过。4.6 第一次可以只看去噪的三步和"共同波动会破坏 MP 边界"这条结论。

4.1 问题与动机

均值–方差组合的解含有 \(\boldsymbol\Sigma^{-1}\)（第 04 册第 12 章 12.10 节），全局最小方差组合（GMV）为

推导拆解：GMV 公式三步得到。目标是 \(\min_{\boldsymbol w}\ \boldsymbol w'\boldsymbol\Sigma\boldsymbol w\)，约束 \(\boldsymbol 1'\boldsymbol w=1\)（权重和为 1）。 第一步，写拉格朗日函数 \(\mathcal L=\boldsymbol w'\boldsymbol\Sigma\boldsymbol w-\lambda(\boldsymbol 1'\boldsymbol w-1)\)。 第二步，对 \(\boldsymbol w\) 求偏导并令其为零：\(2\boldsymbol\Sigma\boldsymbol w-\lambda\boldsymbol 1=\boldsymbol 0\)，所以 \(\boldsymbol w=\tfrac{\lambda}{2}\boldsymbol\Sigma^{-1}\boldsymbol 1\)。这里用了"二次型 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w\) 的梯度是 \(2\boldsymbol\Sigma\boldsymbol w\)"，相当于标量 \(\sigma^2w^2\) 求导得 \(2\sigma^2w\)。 第三步，用约束定出常数：\(\boldsymbol 1'\boldsymbol w=1\) 推出 \(\tfrac\lambda2=1/(\boldsymbol 1'\boldsymbol\Sigma^{-1}\boldsymbol 1)\)，代回即得公式。 关键观察：权重完全由 \(\boldsymbol\Sigma^{-1}\) 决定，协方差估计中的任何误差都会经过"求逆"传到权重上。这就是为什么 GMV 是检验协方差估计质量的"试金石"。

它不需要预期收益，是检验协方差估计质量最干净的试金石：给定估计 \(\hat{\boldsymbol\Sigma}\)，算出 \(\hat{\boldsymbol w}\)，再用真实的 \(\boldsymbol\Sigma\) 计算它的真实方差 \(\hat{\boldsymbol w}'\boldsymbol\Sigma\hat{\boldsymbol w}\)，与理论最优 \(\boldsymbol w^{*\prime}\boldsymbol\Sigma\boldsymbol w^*\) 比较。合成数据的好处是真协方差已知，可以直接做这种比较；在真实数据中只能用样本外实现的波动近似。

本章的所有合成数据都有因子结构和厚尾（多元 \(t\) 或独立 \(t\) 新息），第 4.5 节还加入了 GARCH 波动和危机时的相关跃升。它们没有模拟的是：资产的进入和退出、因子结构本身的漂移、流动性危机中协方差的突变。

本章代码按顺序在同一会话中运行：4.2 节的代码定义了 make_true_cov、gmv 等函数和随机数生成器 rng，后面各段都会用到；4.6 节的代码还要用到 4.4 节定义的 pca_model。

4.2 样本协方差的维度灾难

4.2.1 特征值偏差

\(\boldsymbol\Sigma\) 有 \(N(N+1)/2\) 个参数，\(T\) 个观测只提供 \(NT\) 个数。\(q=N/T\) 不小时，即使真协方差是单位阵，样本特征值也会散开到 Marchenko–Pastur（MP）区间

内。直观上，最大的样本特征值是"在所有方向中挑出方差最大的那个"，必然偏高；最小的必然偏低。第 01 册第 04a 章 4a.5.3 节的 Ky Fan 迹极值定理给出严格版本：前 \(k\) 个特征值之和是矩阵的凸函数，由 Jensen 不等式，样本前 \(k\) 个特征值之和的期望不小于真值；最小特征值之和反之。\(q\ge1\) 时样本协方差奇异，有 \(N-T+1\) 个零特征值。

4.2.2 最小方差组合：风险被系统性低估

GMV 的权重由 \(\hat{\boldsymbol\Sigma}^{-1}\) 决定，而 \(\hat{\boldsymbol\Sigma}^{-1}\) 的最大特征值来自 \(\hat{\boldsymbol\Sigma}\) 偏低的最小特征值。优化器会把权重押在"样本看起来风险很小、实际上并不小"的方向上。正态假设下有两个有用的近似：

两者之比约为 \((1-q)^{-2}\)：\(q=0.5\) 时，模型以为风险只有最优的一半，实际上是最优的两倍。

白话解释：为什么样本特征值会"散开"？用一个比喻：你从 200 只股票里挑一个"历史波动最大的组合"，你总能挑出一个比真实情况更极端的，因为你在大量方向中挑选，挑中的往往是"恰好样本里波动偏大"的那个方向。最小特征值同理，是"恰好样本里看起来最稳"的方向。\(q=N/T\) 越大（股票越多、样本越短），可挑的方向越多、每个方向的估计越粗，散开越严重。 这对 GMV 是致命的：优化器的任务就是找风险最小的方向，于是它会精确地押在那些"样本里看起来最稳、实际并不稳"的方向上。这就像审计中"管理层倾向于选择最有利的估计"：最优化本身就是一种系统性的偏向选择，它会主动利用估计误差。

金融直觉：表中 \(T=300\)、\(q=0.67\) 那一行：模型报告的风险只有最优的 0.32 倍，真实风险却是最优的 3 倍，两者相差近 10 倍。如果风控部门用这个模型的输出设定风险限额，限额会被持续突破，而模型自己毫无察觉。这就是"模型风险"在协方差估计上的具体形态。

import numpy as np
import pandas as pd

rng = np.random.default_rng(404)

# ---------- (a) 真协方差为单位阵：样本特征值也会散开 ----------
N, T = 200, 250
q = N / T
X = rng.standard_normal((T, N))
S = np.cov(X, rowvar=False)
ev = np.linalg.eigvalsh(S)
lo, hi = (1 - np.sqrt(q))**2, (1 + np.sqrt(q))**2
print(f'N={N}, T={T}, q=N/T={q:.2f}；真特征值全为 1')
print(f'样本特征值范围 [{ev.min():.3f}, {ev.max():.3f}]，Marchenko–Pastur 边界 [{lo:.3f}, {hi:.3f}]')
print(f'样本协方差条件数 {ev.max()/ev.min():.0f}（真值为 1）')

# ---------- (b) 真协方差为因子结构：最小方差组合的"预测风险"与"真实风险" ----------
def make_true_cov(N, k=4, seed=1):
    g = np.random.default_rng(seed)
    B = g.normal(0, 1, (N, k)) * np.array([1.0, 0.6, 0.5, 0.4]) ; B[:, 0] = g.normal(1, 0.3, N)
    F = np.diag([0.012, 0.006, 0.005, 0.004])**2                    # 日频因子波动
    D = np.diag(g.uniform(0.01, 0.03, N)**2)                         # 特质日波动 1%–3%
    return B @ F @ B.T + D

def gmv(S):
    w = np.linalg.solve(S, np.ones(len(S)))
    return w / w.sum()

Sig = make_true_cov(N)
w_opt = gmv(Sig); v_opt = w_opt @ Sig @ w_opt
L = np.linalg.cholesky(Sig)
rows = []
for T in [300, 400, 600, 1000, 2500]:
    pred, real = [], []
    for rep in range(20):
        Xr = rng.standard_normal((T, N)) @ L.T
        Sh = np.cov(Xr, rowvar=False)
        w = gmv(Sh)
        pred.append(w @ Sh @ w / v_opt); real.append(w @ Sig @ w / v_opt)
    q = N / T
    rows.append([T, q, np.mean(pred), 1 - q, np.mean(real), 1 / (1 - q)])
print(pd.DataFrame(rows, columns=['T', 'q=N/T', '预测/最优', '≈1-q', '真实/最优', '≈1/(1-q)']).round(3).to_string(index=False))

输出：

N=200, T=250, q=N/T=0.80；真特征值全为 1
样本特征值范围 [0.014, 3.443]，Marchenko–Pastur 边界 [0.011, 3.589]
样本协方差条件数 249（真值为 1）
   T  q=N/T  预测/最优  ≈1-q  真实/最优  ≈1/(1-q)
 300  0.667  0.318 0.333  3.027     3.000
 400  0.500  0.506 0.500  2.028     2.000
 600  0.333  0.660 0.667  1.498     1.500
1000  0.200  0.788 0.800  1.238     1.250
2500  0.080  0.909 0.920  1.090     1.087

\(q=0.8\) 时，真特征值全为 1，样本特征值却散布在 0.014 到 3.44，条件数 249，与 MP 边界 \([0.011,3.589]\) 吻合。GMV 的模拟结果与两个近似式几乎完全一致。实务含义：用 250 个交易日估计 200 只股票的样本协方差，优化出来的组合真实风险约为最优的 5 倍，而模型会告诉你它很安全。

4.3 收缩估计

4.3.1 Ledoit–Wolf

收缩估计把样本协方差 \(\boldsymbol S\) 向一个结构化目标 \(\boldsymbol F\) 拉近：\(\hat{\boldsymbol\Sigma}=\delta\boldsymbol F+(1-\delta)\boldsymbol S\)。这是第 03 册第 12 章 12.7 节 Stein 悖论在矩阵上的版本：牺牲一点偏差，换取方差的大幅下降。Ledoit 与 Wolf（2004）以 \(\boldsymbol F=m\boldsymbol I\)（\(m=\mathrm{tr}\boldsymbol S/N\)）为目标，在 Frobenius 损失下求最优 \(\delta\)，并给出一个不需要知道真值的相合估计：

\(d^2\) 衡量样本协方差离目标多远，\(\bar b^2\) 衡量样本协方差自身的估计误差。误差越大、离目标越近，就收缩得越多。收缩到 \(m\boldsymbol I\) 不改变特征向量，只把所有特征值向均值 \(m\) 拉近，因此条件数一定下降。

推导拆解：逐项看公式在算什么。 \(m=\mathrm{tr}\boldsymbol S/N\) 是所有股票方差的平均值，目标 \(m\boldsymbol I\) 的意思是"假设每只股票方差都等于平均水平、彼此不相关"。这是一个极端简化、但非常稳定的估计。 \(d^2=\tfrac1N\|\boldsymbol S-m\boldsymbol I\|_F^2\)：把 \(\boldsymbol S\) 和目标逐元素相减、平方、求和，再除以 \(N\)，衡量"样本矩阵离目标有多远"。 \(\bar b^2\)：每一天的观测 \(\boldsymbol x_t\) 自己就能形成一个"单日协方差" \(\boldsymbol x_t\boldsymbol x_t'\)。\(\boldsymbol S\) 是这 \(T\) 个单日矩阵的平均。单日矩阵围绕平均值的离散程度，除以 \(T\) 后（再除以 \(T\) 是因为平均值的方差 = 单个的方差 / \(T\)），就是 \(\boldsymbol S\) 本身的估计误差。这和"均值的标准误 \(=s/\sqrt T\)"是同一个思路，只是对象从数字换成了矩阵。 \(\hat\delta=\bar b^2/d^2\)（截在 1 以内）：噪声占"离目标距离"的比例。噪声越大，越不该信样本，越往目标靠。 数值直觉：代码输出 \(\delta=0.097\)，即最终估计 = 9.7% 目标 + 90.3% 样本。

from sklearn.covariance import LedoitWolf, OAS

def lw_identity(X):
    """Ledoit–Wolf (2004)：向 m·I 收缩，m = tr(S)/N。返回 (估计, 收缩强度 δ)。"""
    T, N = X.shape
    Xc = X - X.mean(0)
    S = Xc.T @ Xc / T
    m = np.trace(S) / N
    d2 = np.sum((S - m * np.eye(N))**2) / N                       # 样本协方差离目标多远
    b2_bar = sum(np.sum((np.outer(x, x) - S)**2) for x in Xc) / T**2 / N   # S 自身的估计误差
    b2 = min(b2_bar, d2)
    delta = b2 / d2
    return delta * m * np.eye(N) + (1 - delta) * S, delta

def multivariate_t(T, L, nu=5):
    """厚尾收益：多元 t，协方差为 L L'。"""
    Z = rng.standard_normal((T, L.shape[0])) @ L.T
    W = rng.chisquare(nu, (T, 1)) / nu
    return Z / np.sqrt(W) * np.sqrt((nu - 2) / nu)

X = multivariate_t(250, L)
S_lw, delta = lw_identity(X)
sk = LedoitWolf().fit(X)
print(f'手写 LW 收缩强度 δ={delta:.4f}，sklearn δ={sk.shrinkage_:.4f}，两者矩阵最大差 {np.abs(S_lw - sk.covariance_).max():.2e}')

def evaluate(est_fn, T, reps=20):
    fro, ratio, cond = [], [], []
    for _ in range(reps):
        X = multivariate_t(T, L)
        Sh = est_fn(X)
        fro.append(np.linalg.norm(Sh - Sig) / np.linalg.norm(Sig))
        w = gmv(Sh); ratio.append(w @ Sig @ w / v_opt)
        e = np.linalg.eigvalsh(Sh); cond.append(e.max() / e.min())
    return np.mean(fro), np.mean(ratio), np.median(cond)

ests = {'样本协方差': lambda X: np.cov(X, rowvar=False),
        'LW→单位阵': lambda X: LedoitWolf().fit(X).covariance_,
        'OAS': lambda X: OAS().fit(X).covariance_}
rows = []
for T in [250, 500, 1000]:
    for name, f in ests.items():
        rows.append([T, name, *evaluate(f, T)])
print(pd.DataFrame(rows, columns=['T', '估计量', '相对 Frobenius 误差', 'GMV 真实/最优', '条件数(中位数)'])
      .round(3).to_string(index=False))
print(f'真协方差条件数 {np.linalg.cond(Sig):.0f}')

输出：

手写 LW 收缩强度 δ=0.0970，sklearn δ=0.0970，两者矩阵最大差 2.17e-19
   T    估计量  相对 Frobenius 误差  GMV 真实/最优  条件数(中位数)
 250  样本协方差            0.393      5.334 13345.043
 250 LW→单位阵            0.347      1.802   368.498
 250    OAS            0.366      2.052   776.356
 500  样本协方差            0.255      1.836  1240.124
 500 LW→单位阵            0.261      1.520   439.512
 500    OAS            0.252      1.611   681.787
1000  样本协方差            0.229      1.364   590.111
1000 LW→单位阵            0.197      1.287   416.596
1000    OAS            0.217      1.344   532.949
真协方差条件数 299

手写实现与 sklearn 完全一致。\(T=250\) 时，LW 把 GMV 的真实风险从最优的 5.3 倍降到 1.8 倍，条件数从一万多降到三四百。但也要看到它的局限：

• Frobenius 误差并不总是改善（\(T=500\) 时 LW 略差于样本协方差）。Frobenius 损失由大特征值主导，而 GMV 关心的是小特征值；两个目标不一致。
• 收缩到单位阵会压低市场方向的方差。真实股票协方差的第一特征值（市场）远大于均值 \(m\)，收缩会低估所有多头组合的风险（下一节的随机多头组合偏差统计量 1.12 就是这个效应）。
• 因此实务中更常用结构化的目标：单指数模型、常相关矩阵，或者直接用因子模型（Ledoit 与 Wolf 在其他论文中给出了这些目标的收缩强度公式）。

4.4 因子风险模型

4.4.1 结构

因子模型（第 06 册第 09 章 9.1 节）把协方差写成

\(\boldsymbol X\) 为 \(N\times k\) 暴露矩阵，\(\boldsymbol F\) 为 \(k\times k\) 因子协方差，\(\boldsymbol D\) 为对角特质方差。参数从 \(N(N+1)/2\) 降到 \(Nk+k(k+1)/2+N\)。两类做法：

• 统计因子模型：\(\boldsymbol X\) 与 \(\boldsymbol F\) 都由 PCA 从收益中估计（第 06 册第 09 章 9.4 节）。不需要外部数据，但因子没有经济含义，个数需要选择，因子的符号和顺序可能随估计窗口跳动。
• 基本面因子模型（BARRA 类）：\(\boldsymbol X\) 是已知的行业哑变量和风格暴露，每天做截面回归得到因子收益 \(\hat{\boldsymbol f}_t\)，\(\hat{\boldsymbol F}\) 是 \(\hat{\boldsymbol f}_t\) 的协方差，\(\hat{\boldsymbol D}\) 来自回归残差（第 06 册第 09 章 9.3.1 节）。因子可解释，可以直接做风险归因；暴露当天就能更新，对新股、风格突变反应快。

金融直觉：这个公式就是 CFA 里单指数模型 \(\sigma_i^2=\beta_i^2\sigma_m^2+\sigma_{\varepsilon,i}^2\)、\(\sigma_{ij}=\beta_i\beta_j\sigma_m^2\) 的多因子版本。\(\boldsymbol X\) 的每一行是一只股票对各因子的"β"（属于哪个行业、市值风格暴露多少），\(\boldsymbol F\) 是因子之间的协方差，\(\boldsymbol D\) 是每只股票自己的特质风险。两只股票之间的协方差只能通过共同因子产生，特质部分互不相关。 参数数量的对比很直观：500 只股票，样本协方差要估 \(500\times501/2=125{,}250\) 个参数；11 个因子的模型只要 \(500\times11+66+500\approx6{,}066\) 个，而且暴露 \(\boldsymbol X\) 在基本面模型里是已知的，真正要从收益估计的只有 \(\boldsymbol F\) 的 66 个和 \(\boldsymbol D\) 的 500 个。这就是"给估计加结构"换来稳定性的本质。 风控报告里的"风险归因"也来自这里：组合方差 \(\boldsymbol w'\boldsymbol\Sigma\boldsymbol w=(\boldsymbol X'\boldsymbol w)'\boldsymbol F(\boldsymbol X'\boldsymbol w)+\boldsymbol w'\boldsymbol D\boldsymbol w\)，前一项是"因子风险"（\(\boldsymbol X'\boldsymbol w\) 就是组合对各因子的暴露），后一项是"特质风险"。

4.4.2 Woodbury：不形成 \(N\times N\) 矩阵

优化器要的是 \(\boldsymbol\Sigma^{-1}\boldsymbol b\) 而不是 \(\boldsymbol\Sigma^{-1}\) 本身。由 Sherman–Morrison–Woodbury 公式（第 01 册第 00 章 0.7 节）

只需要解一个 \(k\times k\) 的系统，计算量从 \(O(N^3)\) 降到 \(O(Nk^2)\)，内存从 \(O(N^2)\) 降到 \(O(Nk)\)。

白话解释：Woodbury 公式的作用是"把大矩阵求逆换成小矩阵求逆"。右边所有需要求逆的东西只有两个：\(\boldsymbol D^{-1}\)（对角矩阵，求逆就是每个对角元取倒数，几乎不花时间）和中间那个 \(k\times k\) 矩阵（\(k=11\)，极小）。\(O(N^3)\) 是"计算量大约与 \(N^3\) 成正比"的记号：\(N=3000\) 时 \(N^3=2.7\times10^{10}\)，而 \(Nk^2=3000\times121\approx3.6\times10^5\)，差五个数量级。 对照代码：XtDi = Xb.T * Dinv1 就是 \(\boldsymbol X'\boldsymbol D^{-1}\)（每一列除以对应的特质方差）；inner 是 \(\boldsymbol F^{-1}+\boldsymbol X'\boldsymbol D^{-1}\boldsymbol X\)；最后一行把公式右乘向量 \(\boldsymbol 1\)，全程没有出现 \(3000\times3000\) 的矩阵。 验证方法：把公式右边乘上 \(\boldsymbol\Sigma=\boldsymbol X\boldsymbol F\boldsymbol X'+\boldsymbol D\)，展开后各项相消得到单位阵（练习 3）。

4.4.3 偏差统计量

风险模型最常用的检验是偏差统计量：对一个组合，计算样本外每期收益除以模型预测波动 \(z_t=r_{p,t}/\hat\sigma_{p,t}\)，再求 \(z_t\) 的标准差 \(b\)。\(b\approx1\) 表示预测无偏，\(b>1\) 表示低估风险。正态下 \(b\) 的抽样误差约为 \(\sqrt{1/(2T)}\)，\(T=250\) 时 95% 区间约为 \([0.91,1.09]\)；厚尾会让区间更宽。检验要对两类组合做：随机组合（检验模型的整体校准），以及用该模型优化出的组合（检验模型在优化器最"喜欢"的方向上是否可靠）。后者才是真正的压力测试。

金融直觉：偏差统计量就是风险模型的"回溯测试"，与 VaR 回溯测试同一思路。如果模型每天说"今天组合波动 1%"，那么实际收益除以 1% 应当像一个标准差为 1 的序列。\(b=1.5\) 表示实际波动是预测的 1.5 倍，风险被低估了三分之一。 抽样误差 \(\sqrt{1/(2T)}\) 的来源：正态样本的标准差估计量，其标准误约为 \(\sigma/\sqrt{2T}\)。\(T=250\) 时为 \(0.045\)，\(\pm1.96\) 倍即 \([0.91,1.09]\)。 为什么要特别检验"优化出的组合"？随机组合是"平均方向"，模型在平均方向上通常不太差；优化器却专门找模型认为风险最低的方向，那恰恰是估计误差最集中的地方。类比：审计抽样如果只随机抽，可能查不出问题；针对管理层"最有动机操纵"的科目做专项测试，才更能暴露缺陷。

import time

# ---------- 真实世界：8 个行业 + 3 个风格因子，暴露已知（基本面因子模型的设定） ----------
N, K, P = 500, 8, 3
g = np.random.default_rng(42)
Xexp = np.column_stack([np.eye(K)[g.integers(0, K, N)], g.standard_normal((N, P))])   # N×11 暴露
C = np.full((K + P, K + P), 0.1); C[:K, :K] = 0.6; np.fill_diagonal(C, 1.0)  # 行业间相关 0.6（共同的市场成分）
sd_f = np.r_[np.full(K, 0.011), np.full(P, 0.004)]                          # 行业日波动 1.1%，风格 0.4%
F_true = np.outer(sd_f, sd_f) * C
sd_eps = g.uniform(0.01, 0.03, N)
Sig2 = Xexp @ F_true @ Xexp.T + np.diag(sd_eps**2)
w_opt2 = gmv(Sig2); v_opt2 = w_opt2 @ Sig2 @ w_opt2
Lf = np.linalg.cholesky(F_true)

def simulate(T):
    nu = 5; sc = np.sqrt((nu - 2) / nu)
    f = (rng.standard_normal((T, K + P)) @ Lf.T) / np.sqrt(rng.chisquare(nu, (T, 1)) / nu) * sc
    e = rng.standard_t(nu, (T, N)) * sc * sd_eps
    return f @ Xexp.T + e

def fundamental(Rw):
    """逐日截面 OLS 求因子收益（第 06 册 9.3.1 节），再拼出 Σ = X F X' + D。"""
    fhat = np.linalg.lstsq(Xexp, Rw.T, rcond=None)[0].T                  # T×11
    resid = Rw - fhat @ Xexp.T
    return Xexp @ np.cov(fhat, rowvar=False) @ Xexp.T + np.diag(resid.var(0))

def pca_model(Rw, k):
    """统计因子模型：样本协方差前 k 个主成分 + 对角特质方差。"""
    S = np.cov(Rw, rowvar=False)
    val, vec = np.linalg.eigh(S)
    B = vec[:, -k:] * np.sqrt(val[-k:])
    return B @ B.T + np.diag(np.maximum(np.diag(S - B @ B.T), 1e-8))

def bias_stat(w, Sh, Rtest):
    """偏差统计量：样本外 r_p / σ_pred 的标准差。≈1 表示风险预测无偏，>1 表示低估风险。"""
    return np.std(Rtest @ w / np.sqrt(w @ Sh @ w))

Ttrain = 250
ests = {'LW→单位阵': lambda R: LedoitWolf().fit(R).covariance_,
        'PCA k=11': lambda R: pca_model(R, 11),
        '基本面因子': fundamental}
res = {k: [] for k in ests}
for rep in range(10):
    Rtr, Rte = simulate(Ttrain), simulate(250)
    wr = rng.dirichlet(np.ones(N), 50)                                    # 50 个随机多头组合
    for k, f in ests.items():
        Sh = f(Rtr); w = gmv(Sh)
        res[k].append([w @ Sig2 @ w / v_opt2, bias_stat(w, Sh, Rte),
                       np.mean([bias_stat(x, Sh, Rte) for x in wr])])
print(f'N={N}，训练窗口 T={Ttrain}（样本协方差秩 ≤ {Ttrain-1}，奇异，无法直接求 GMV）')
print(pd.DataFrame({k: np.mean(v, 0) for k, v in res.items()},
                   index=['GMV 真实/最优', 'GMV 偏差统计量', '随机组合偏差统计量']).T.round(3).to_string())

# ---------- Woodbury：不形成 N×N 矩阵也能求 Σ^{-1}1 ----------
Nbig = 3000
Xb = np.column_stack([np.eye(K)[g.integers(0, K, Nbig)], g.standard_normal((Nbig, P))])
Db = g.uniform(0.01, 0.03, Nbig)**2
t0 = time.perf_counter()
Sb = Xb @ F_true @ Xb.T + np.diag(Db)
x_direct = np.linalg.solve(Sb, np.ones(Nbig)); t_direct = time.perf_counter() - t0
t0 = time.perf_counter()
Dinv1 = 1 / Db; XtDi = Xb.T * Dinv1                                       # X' D^{-1}
inner = np.linalg.inv(F_true) + XtDi @ Xb                                  # k×k
x_wood = Dinv1 - XtDi.T @ np.linalg.solve(inner, XtDi @ np.ones(Nbig))     # Σ^{-1}1
t_wood = time.perf_counter() - t0
print(f'\nN={Nbig}：直接求解 {t_direct*1000:.0f} ms，Woodbury {t_wood*1000:.1f} ms，'
      f'相对差 {np.abs(x_wood - x_direct).max() / np.abs(x_direct).max():.1e}')

输出：

N=500，训练窗口 T=250（样本协方差秩 ≤ 249，奇异，无法直接求 GMV）
          GMV 真实/最优  GMV 偏差统计量  随机组合偏差统计量
LW→单位阵        1.516      2.900      1.123
PCA k=11      1.308      2.110      1.000
基本面因子         1.070      1.071      1.000

N=3000：直接求解 87 ms，Woodbury 0.1 ms，相对差 6.5e-13

\(N=500\)、\(T=250\) 时样本协方差奇异，无法直接使用。三种可用的估计中：

• LW 优化出的 GMV 真实风险是最优的 1.52 倍，而模型预测的波动只有实际的三分之一左右（偏差统计量 2.90）。随机多头组合的偏差统计量 1.12，说明市场方向也被低估了。
• PCA 统计模型好一些，但优化组合的偏差统计量仍为 2.11：11 个主成分中混入了噪声方向，而特质方差 \(\boldsymbol D\) 在那些被优化器利用的方向上偏低。
• 基本面模型的结构与真实世界一致，GMV 只比最优差 7%，偏差统计量 1.07，在抽样误差范围内。

白话解释：这段代码的设计是"三方比武、同一考卷"。先造一个真实世界：500 只股票、8 个行业 + 3 个风格因子，暴露已知；因子收益和特质收益都是厚尾的。然后每次模拟 250 天训练数据和 250 天测试数据，让三种方法各自估一个协方差矩阵，用它求 GMV，再用三个尺子打分：(1) 用真协方差算这个 GMV 的真实风险，与理论最优比；(2) 在测试数据上算 GMV 的偏差统计量；(3) 在测试数据上算 50 个随机多头组合的平均偏差统计量。 fundamental 函数是 Barra 方法的最简版：每天用已知暴露 \(\boldsymbol X\) 对 500 只股票的收益做一次截面回归，得到当天 11 个因子收益；250 天的因子收益算协方差得 \(\hat{\boldsymbol F}\)；回归残差的方差得 \(\hat{\boldsymbol D}\)。pca_model 不知道暴露，从样本协方差的前 \(k\) 个特征向量"猜"出暴露。两者的差别正是"有没有外部信息（行业分类、风格暴露）"。

这个对比的前提是基本面模型的结构是对的。真实世界里暴露矩阵也是近似的：行业分类不完美、风格因子遗漏、暴露有测量误差。商业风险模型会在此基础上做多项修正，例如对因子协方差做特征值调整、对特质方差做贝叶斯收缩、把日频协方差换算到月频时做 Newey–West 式的序列相关调整。Woodbury 一节显示，\(N=3000\) 时直接构造并求解需要近 90 毫秒（耗时因机器而异），Woodbury 只需零点几毫秒，结果相同到 \(10^{-12}\)；在每日重新优化几千只股票、或在回测中重复上千次时，这个差别是决定性的。

4.5 时变波动与相关：EWMA 与 DCC

前面假设 \(\boldsymbol\Sigma\) 不随时间变化。真实收益有波动聚集，相关性在危机中上升，用过去一年等权估计的协方差会在危机开始时严重低估风险。两种常用的跟踪方法：

EWMA（第 06 册第 10 章 10.2 节）：\(\hat{\boldsymbol\Sigma}_t=\lambda\hat{\boldsymbol\Sigma}_{t-1}+(1-\lambda)\boldsymbol r_{t-1}\boldsymbol r_{t-1}'\)。半衰期 \(\ln0.5/\ln\lambda\)，\(\lambda=0.94\) 约 11 天，\(\lambda=0.97\) 约 23 天。一个参数，计算快，自动半正定；但所有元素共用一个衰减速度，且没有均值回复。

DCC（第 06 册第 10 章 10.6 节）：先对每个资产拟合单变量 GARCH，得到标准化残差 \(\boldsymbol\epsilon_t\)，再让相关矩阵按

演化。波动和相关分开建模，波动有均值回复；两步估计时第二步只有 \((a,b)\) 两个参数。

评价预测用两个指标：偏差统计量，以及 QLIKE 损失 \(\frac1T\sum_t\big(\ln\hat\sigma_t^2+r_t^2/\hat\sigma_t^2\big)\)。QLIKE 是正态似然的负值（差一个常数），对低估风险的惩罚比高估更重，而且在 \(r_t^2\) 只是真实方差的噪声代理时仍能正确排序预测模型。

白话解释：EWMA 的递推式读作"今天的风险估计 = 94% 昨天的估计 + 6% 昨天实际发生的波动"。\(\boldsymbol r\boldsymbol r'\) 是单日收益的外积，对角线是各资产收益的平方（单日方差的粗略估计），非对角线是两两乘积（单日协方差的粗略估计）。展开递推，相当于对过去每天的 \(\boldsymbol r\boldsymbol r'\) 做加权平均，权重按 \(\lambda\) 的幂次几何递减：昨天 \((1-\lambda)\)，前天 \((1-\lambda)\lambda\)，依此类推（几何级数，见 第 00 册第 04 章）。这正是 RiskMetrics 的做法。 DCC 多了两点：一是先把每个资产的波动用 GARCH 单独建模，再研究"标准化之后"的相关，二是相关有长期均值 \(\bar{\boldsymbol Q}\)，危机过后会回归。公式里的 \(\mathrm{diag}(\boldsymbol Q_t)^{-1/2}\boldsymbol Q_t\,\mathrm{diag}(\boldsymbol Q_t)^{-1/2}\) 只是把 \(\boldsymbol Q_t\) 的对角线归一为 1，即"协方差除以两个标准差得到相关系数"的矩阵写法。

推导拆解：QLIKE 为什么对低估惩罚更重？固定真实方差 \(\sigma^2\)，看单期期望损失 \(\ln\hat\sigma^2+\sigma^2/\hat\sigma^2\) 作为 \(\hat\sigma^2\) 的函数。对 \(\hat\sigma^2\) 求导：\(1/\hat\sigma^2-\sigma^2/\hat\sigma^4=0\)，得 \(\hat\sigma^2=\sigma^2\)，即预测正确时损失最小。数值例：设 \(\sigma^2=1\)，预测成 0.5 时损失 \(\ln0.5+2=1.31\)；预测成 2 时损失 \(\ln2+0.5=1.19\)；最优值为 1。同样差一倍，低估的损失更大。这和风控的偏好一致：低估风险比高估风险更危险。

from arch import arch_model
from scipy.optimize import minimize

# ---------- 3 个资产（可理解为 3 个因子）的日收益：GARCH 波动 + 两状态相关（平时 0.2，危机 0.7） ----------
T, k = 3000, 3
state = np.zeros(T, int)
for t in range(1, T):                                   # 马尔可夫状态：危机平均持续 50 天
    p_stay = 0.98 if state[t-1] else 0.995
    state[t] = state[t-1] if rng.random() < p_stay else 1 - state[t-1]
C0 = np.full((k, k), 0.2); np.fill_diagonal(C0, 1); C1 = np.full((k, k), 0.7); np.fill_diagonal(C1, 1)
Lc = [np.linalg.cholesky(C0), np.linalg.cholesky(C1)]
om, al, be = np.array([0.02, 0.03, 0.015]), np.array([0.08, 0.10, 0.06]), np.array([0.90, 0.88, 0.92])
h = om / (1 - al - be); r = np.zeros((T, k)); u = np.zeros(k)
for t in range(T):
    if t: h = om + al * u**2 + be * h                    # GARCH 部分由"平时口径"的冲击 u 驱动
    u = np.sqrt(h) * (Lc[state[t]] @ (rng.standard_t(6, k) / np.sqrt(1.5)))
    r[t] = u * np.sqrt(1 + 1.0 * state[t])               # 危机时方差再翻倍；收益单位：%
split = 1500                                             # 前 1500 天估计参数，后 1500 天样本外评估

# ---------- 预测器：都只用 t-1 及之前的信息预测第 t 天的协方差 ----------
def ewma_cov(r, lam):
    H = np.empty((len(r), k, k)); H[0] = np.cov(r[:250], rowvar=False)
    for t in range(1, len(r)):
        H[t] = lam * H[t-1] + (1 - lam) * np.outer(r[t-1], r[t-1])
    return H

def rolling_cov(r, win=250):
    H = np.full((len(r), k, k), np.nan)
    for t in range(win, len(r)):
        H[t] = np.cov(r[t-win:t], rowvar=False)
    return H

def dcc_cov(r, split):
    sig, eps = np.empty_like(r), np.empty_like(r)
    for j in range(k):                                   # 第一步：逐个 GARCH(1,1)，参数只用训练段估计
        fit = arch_model(r[:split, j], vol='GARCH', p=1, q=1, dist='t').fit(disp='off')
        full = arch_model(r[:, j], vol='GARCH', p=1, q=1, dist='t').fix(fit.params)
        sig[:, j] = full.conditional_volatility           # σ_t 只依赖 t-1 及之前的数据
        eps[:, j] = (r[:, j] - fit.params['mu']) / sig[:, j]
    Qbar = np.cov(eps[:split], rowvar=False)
    def corr_path(a, b, e):
        Q = Qbar.copy(); Rs = np.empty((len(e), k, k))
        for t in range(len(e)):
            if t: Q = (1 - a - b) * Qbar + a * np.outer(e[t-1], e[t-1]) + b * Q
            d = 1 / np.sqrt(np.diag(Q)); Rs[t] = Q * np.outer(d, d)
        return Rs
    def nll(p):                                          # 第二步：相关部分的负对数似然（第 06 册 10.6 节）
        a, b = p
        if a < 0 or b < 0 or a + b >= 0.999:
            return 1e10
        Rs = corr_path(a, b, eps[:split])
        _, logdet = np.linalg.slogdet(Rs)
        quad = np.einsum('ti,tij,tj->t', eps[:split], np.linalg.inv(Rs), eps[:split])
        return 0.5 * np.sum(logdet + quad)
    a, b = minimize(nll, [0.03, 0.95], method='Nelder-Mead').x
    Rs = corr_path(a, b, eps)
    return sig[:, :, None] * Rs * sig[:, None, :], (a, b)

H_dcc, (a_hat, b_hat) = dcc_cov(r, split)
preds = {'滚动 250 日': rolling_cov(r), 'EWMA λ=0.97': ewma_cov(r, 0.97),
         'EWMA λ=0.94': ewma_cov(r, 0.94), 'DCC-GARCH': H_dcc}
print(f'DCC 估计：a={a_hat:.3f}, b={b_hat:.3f}；样本外危机天数占比 {state[split:].mean():.1%}')

ports = {'等权多头': np.ones(k) / k, '多空 (1,-1,0)': np.array([1., -1, 0])}
rows = []
for pn, w in ports.items():
    rp = r[split:] @ w
    for name, H in preds.items():
        v = np.einsum('i,tij,j->t', w, H[split:], w)
        z = rp / np.sqrt(v)
        rows.append([pn, name, z.std(), np.mean(np.log(v) + rp**2 / v), np.mean(np.abs(z) > 2.576)])
print(pd.DataFrame(rows, columns=['组合', '模型', '偏差统计量', 'QLIKE(越小越好)', '|z|>2.58 频率(应≈1%)'])
      .round(3).to_string(index=False))

输出：

DCC 估计：a=0.030, b=0.941；样本外危机天数占比 18.5%
         组合          模型  偏差统计量  QLIKE(越小越好)  |z|>2.58 频率(应≈1%)
       等权多头    滚动 250 日  1.095        0.570              0.034
       等权多头 EWMA λ=0.97  1.037        0.409              0.022
       等权多头 EWMA λ=0.94  1.049        0.403              0.022
       等权多头   DCC-GARCH  0.952        0.406              0.014
多空 (1,-1,0)    滚动 250 日  1.073        1.751              0.031
多空 (1,-1,0) EWMA λ=0.97  1.030        1.511              0.023
多空 (1,-1,0) EWMA λ=0.94  1.050        1.494              0.022
多空 (1,-1,0)   DCC-GARCH  1.012        1.472              0.019

结论：

• 滚动 250 日的 QLIKE 明显最差，偏差统计量 1.07–1.10，\(|z|>2.58\) 的频率是名义 1% 的 3 倍以上。它对危机的反应要等几个月，而危机平均只持续 50 天。
• EWMA 大幅改善，两种 \(\lambda\) 差别不大。
• DCC 的 QLIKE 与 EWMA 相近或略好，极端偏离频率最接近名义水平（1.4%–1.9%）；等权多头的偏差统计量 0.95，略为高估风险，因为 GARCH 的均值回复在危机结束后更快地把波动拉回长期水平。

在实际的风险模型中，EWMA/DCC 一般用在因子协方差 \(\boldsymbol F\)（\(k\) 只有几十）上，而不是 \(N\times N\) 的个股协方差上；特质方差单独用 EWMA 或 GARCH。波动的半衰期通常比相关的半衰期短，例如波动用 1–3 个月、相关用 6–12 个月，因为相关的估计噪声更大。

4.6 随机矩阵去噪

如果收益只是独立噪声，样本相关矩阵的特征值应落在 MP 区间内。Laloux 等（1999）发现，美国股票相关矩阵的绝大多数特征值都在 MP 区间之内，只有少数几个（市场和若干行业）明显超出。这启发了一种去噪方法：

1. 计算样本相关矩阵的特征分解，噪声方差取 \(\sigma^2=1-\lambda_{\max}/N\)（扣除市场模式解释的部分），\(\lambda_+=\sigma^2(1+\sqrt q)^2\)；
2. 把不超过 \(\lambda_+\) 的特征值全部替换为它们的均值（保持迹不变），特征向量不动；
3. 重建相关矩阵，把对角线恢复为 1，再乘回样本标准差得到协方差。

这本质上是"保留 \(k\) 个信号方向、其余方向等方差"的统计因子模型，\(k\) 由 MP 边界自动决定。

白话解释：去噪的思路是"先画一条噪声线，线以上的当真，线以下的一视同仁"。MP 边界 \(\lambda_+\) 是"如果完全没有共同因子，样本特征值最多能大到哪里"。超出它的特征值（例如市场、几个大行业）被认为是真实的风险来源，原样保留；低于它的特征值大小参差不齐，但这些差异主要是样本噪声，于是全部换成它们的平均值。这样做不改变总方差（迹不变），但消除了"有些方向看起来风险特别低"的假象，GMV 也就不会再押在这些方向上。 为什么 \(\sigma^2=1-\lambda_{\max}/N\)？相关矩阵的迹是 \(N\)（对角线全是 1），即总方差为 \(N\)。最大特征值代表的市场模式已经解释了 \(\lambda_{\max}\)，剩下平均到每个方向的"噪声方差"约为 \((N-\lambda_{\max})/N\)。 数值例：本节的模拟中 \(q=0.6\)，\(\lambda_{\max}\approx85\)，\(N=300\)，\(\sigma^2\approx0.72\)，\(\lambda_+\approx0.72\times(1+\sqrt{0.6})^2\approx2.26\)，与输出一致。

# ---------- 随机矩阵去噪：把落在 Marchenko–Pastur 噪声带内的特征值"压平" ----------
def rmt_clip(R):
    """输入收益矩阵 T×N。对样本相关矩阵做特征值截断：
    λ+ 以下的特征值替换为它们的均值（保持迹不变），再恢复对角线为 1，最后乘回样本标准差。"""
    T, N = R.shape
    sd = R.std(0, ddof=1)
    C = np.corrcoef(R, rowvar=False)
    val, vec = np.linalg.eigh(C)
    q = N / T
    sigma2 = 1 - val[-1] / N                              # Laloux 等：扣除市场模式解释的方差
    lam_plus = sigma2 * (1 + np.sqrt(q))**2
    noise = val <= lam_plus
    val_c = val.copy(); val_c[noise] = val[noise].mean()
    Cc = vec @ np.diag(val_c) @ vec.T
    d = 1 / np.sqrt(np.diag(Cc)); Cc = Cc * np.outer(d, d)
    return Cc * np.outer(sd, sd), (~noise).sum(), lam_plus

N, T = 300, 500
Sig3 = make_true_cov(N, seed=7)
L3 = np.linalg.cholesky(Sig3); w3 = gmv(Sig3); v3 = w3 @ Sig3 @ w3
sd3 = np.sqrt(np.diag(Sig3)); Ctrue = Sig3 / np.outer(sd3, sd3)
ev_true = np.linalg.eigvalsh(Ctrue)
print(f'真相关矩阵：最大 6 个特征值 {np.round(ev_true[-6:][::-1], 2).tolist()}')

ests = {'样本协方差': lambda R: np.cov(R, rowvar=False),
        'LW→单位阵': lambda R: LedoitWolf().fit(R).covariance_,
        'RMT 截断': lambda R: rmt_clip(R)[0],
        'PCA k=4': lambda R: pca_model(R, 4)}
def indep_t(T, L, nu=5):
    """厚尾但没有"共同波动"：各分量独立的 t 新息再做线性组合。"""
    return rng.standard_t(nu, (T, L.shape[0])) * np.sqrt((nu - 2) / nu) @ L.T

rows = {k: [] for k in ests}; n_sig = []
for rep in range(10):
    R = indep_t(T, L3)
    n_sig.append(int(rmt_clip(R)[1]))
    for k, f in ests.items():
        Sh = f(R); sdh = np.sqrt(np.diag(Sh))
        w = gmv(Sh)
        rows[k].append([np.linalg.norm(Sh / np.outer(sdh, sdh) - Ctrue) / np.linalg.norm(Ctrue),
                        w @ Sig3 @ w / v3, np.linalg.cond(Sh)])
R = indep_t(T, L3)
_, ns, lp = rmt_clip(R)
ev = np.linalg.eigvalsh(np.corrcoef(R, rowvar=False))
print(f'q=N/T={N/T:.2f}，λ+={lp:.2f}；10 次模拟中超出 λ+ 的特征值个数：{n_sig}（真因子数 4）')
print(f'样本最大 6 个特征值 {np.round(ev[-6:][::-1], 2).tolist()}；最小特征值 {ev[0]:.3f}（真值 {ev_true[0]:.3f}）')
Rc = multivariate_t(T, L3)                                # 有共同波动（所有股票同一天一起放大）
evc = np.linalg.eigvalsh(np.corrcoef(Rc, rowvar=False))
print(f'若收益有共同波动：超出 λ+ 的个数 {int(rmt_clip(Rc)[1])}，第 5、6 大特征值 {np.round(evc[-6:-4][::-1], 2).tolist()}')
print(pd.DataFrame({k: np.mean(v, 0) for k, v in rows.items()},
                   index=['相关矩阵相对误差', 'GMV 真实/最优', '条件数']).T.round(3).to_string())

输出：

真相关矩阵：最大 6 个特征值 [87.19, 7.6, 3.86, 2.17, 0.97, 0.95]
q=N/T=0.60，λ+=2.26；10 次模拟中超出 λ+ 的特征值个数：[4, 5, 4, 4, 4, 4, 4, 4, 4, 4]（真因子数 4）
样本最大 6 个特征值 [84.8, 8.92, 4.45, 2.66, 2.16, 2.11]；最小特征值 0.033（真值 0.255）
若收益有共同波动：超出 λ+ 的个数 16，第 5、6 大特征值 [3.48, 3.11]
         相关矩阵相对误差  GMV 真实/最优       条件数
样本协方差       0.144      2.453  3582.222
LW→单位阵      0.142      1.894  1541.461
RMT 截断      0.143      1.149   481.757
PCA k=4     0.110      1.077   546.782

在噪声相互独立的设定下，MP 边界几乎总是准确识别出 4 个真因子。去噪后 GMV 的真实风险从最优的 2.45 倍降到 1.15 倍，远好于 LW（1.89 倍）。用正确个数的 PCA 因子模型（\(k=4\)）表现最好，因为它还用对角矩阵描述了异质的特质方差。三种估计的相关矩阵 Frobenius 误差几乎相同，再次说明 Frobenius 误差与组合优化的目标不一致。

一个重要的反例在倒数第二行：当收益存在共同波动（所有股票同一天一起放大，这正是真实市场的常态），噪声特征值会被推出 MP 区间，超出 \(\lambda_+\) 的特征值从 4 个变成 16 个。MP 定律假设各观测独立同分布，共同波动破坏了这个假设。实务中的对策是先用 EWMA 或 GARCH 波动对收益做标准化，再做特征值分析；或者不依赖 MP 边界，而用样本外的 GMV 风险或 QLIKE 选择 \(k\)。

金融直觉：为什么"所有股票同一天一起放大"会骗过 MP 边界？代码里 multivariate_t 让每一天所有股票共用一个随机缩放因子（危机日全体放大），而 indep_t 让每只股票各自抽厚尾冲击。前者更像真实市场：恐慌日所有股票都剧烈波动。这种"共同的波动状态"本身就在制造股票之间的相关，于是样本里冒出许多额外的"大特征值"，看起来像真实因子，其实只是波动状态的变化。这提醒我们：统计工具的假设（这里是"各天独立同分布"）在金融数据里经常不成立，套用前要先检查。

4.7 常见陷阱与检查清单

• [ ] 报告 \(q=N/T\)；\(q>0.1\) 时不要直接对样本协方差求逆。
• [ ] 用优化器输出的组合做偏差统计量检验，而不只是随机组合。
• [ ] 不用 Frobenius 误差作为选择风险模型的唯一标准；以样本外 GMV 风险或 QLIKE 为主。
• [ ] 收缩到单位阵会低估市场方向风险；多头组合优先考虑结构化目标或因子模型。
• [ ] 因子模型用 Woodbury 求解，不形成 \(N\times N\) 矩阵；检查 \(\boldsymbol D\) 的对角元为正。
• [ ] PCA 因子数和 RMT 边界在有共同波动时会失真，先做波动标准化。
• [ ] 波动与相关使用不同的半衰期；危机期单独检查偏差统计量与极端偏离频率。
• [ ] 风险预测与回测使用同一时点信息，协方差估计窗口不包含预测期。

本章小结

样本协方差在 \(q=N/T\) 不可忽略时，特征值向两端散开，最小方差组合的真实风险约为最优的 \(1/(1-q)\) 倍，模型预测的风险却只有最优的 \((1-q)\) 倍，本章的模拟与这两个近似式几乎完全吻合。Ledoit–Wolf 收缩到单位阵能显著降低条件数和组合风险，但会压低市场方向的方差；与真实结构一致的基本面因子模型在本章的模拟中把 GMV 风险降到最优的 1.07 倍、偏差统计量接近 1，并可用 Woodbury 公式以 \(O(Nk^2)\) 的代价求解。时变方面，滚动等权协方差对危机反应过慢，EWMA 和 DCC 的 QLIKE 明显更低，DCC 的极端偏离频率最接近名义水平。随机矩阵去噪用 MP 边界自动选因子数，在噪声独立时效果接近正确设定的 PCA 模型，但共同波动会让噪声特征值越过边界，使用前应先做波动标准化。

练习

1. 在 4.2 节中把收益改为多元 \(t\)（自由度 4），重做 GMV 实验。两个近似式还成立吗？偏离的方向是什么？
2. 实现 Ledoit–Wolf 向"常相关矩阵"收缩的版本（目标的对角元为样本方差，非对角元为平均样本相关乘以两个标准差），在 4.4 节的设定下与收缩到单位阵比较随机多头组合的偏差统计量。 提示：只需把 \(d^2\) 中的目标换掉，\(\bar b^2\) 的近似可以沿用。
3. 证明 4.4.2 节的 Woodbury 公式，并写出 GMV 权重 \(\boldsymbol\Sigma^{-1}\boldsymbol 1\) 只用 \(O(Nk^2)\) 运算的完整算法。
4. 在 4.4 节的基本面模型中故意遗漏一个风格因子，观察随机组合与 GMV 组合的偏差统计量如何变化。哪一种组合更早暴露模型缺陷？
5. 在 4.5 节中把 EWMA 改为"波动用 \(\lambda=0.94\)、相关用 \(\lambda=0.99\)"的双半衰期版本，与单一 \(\lambda\) 比较 QLIKE。
6. 在 4.6 节有共同波动的数据上，先用每日横截面收益的平方和的 EWMA 估计"共同波动"，把收益除以它，再做 RMT 截断，检查超出 \(\lambda_+\) 的特征值个数是否回到 4 附近。

延伸阅读

• 第 01 册第 00 章 0.7 节（Sherman–Morrison–Woodbury 公式）与第 04a 章（Hermitian 矩阵的特征值不等式）。
• 第 03 册第 12 章 12.7 节（Stein 悖论）与第 14 章（多元模型）。
• 第 04 册附录 A.3 节（因子模型协方差的快速求解与病态）。
• 第 06 册第 09 章（主成分分析与因子模型）、第 10 章（多元波动率模型：EWMA、DCC）。
• 第 08 册第 23 章（波动率与相关系数的估计）。
• Ledoit, O., & Wolf, M. (2004). A Well-Conditioned Estimator for Large-Dimensional Covariance Matrices. Journal of Multivariate Analysis, 88(2), 365–411.
• Laloux, L., Cizeau, P., Bouchaud, J.-P., & Potters, M. (1999). Noise Dressing of Financial Correlation Matrices. Physical Review Letters, 83(7), 1467–1470.
• Engle, R. (2002). Dynamic Conditional Correlation: A Simple Class of Multivariate GARCH Models. Journal of Business & Economic Statistics, 20(3), 339–350.
• Grinold, R. C., & Kahn, R. N. (2000). Active Portfolio Management (2nd ed.). McGraw-Hill. 关于多因子风险模型的章节。