稳健协方差 MCD 估计：用最小协方差行列式抗住极端值污染

先说结论：样本协方差矩阵是量化里用得最多、也最脆弱的估计量——它的崩溃点是 0，意味着单单几个极端值就能把整个相关结构拖到面目全非。 合成实验实测：两个资产真实相关 0.7，只往里塞 10% 的反向污染样本，经典皮尔逊相关直接崩到 −0.068（连符号都翻了），而 MCD（最小协方差行列式）稳健估计还原到 0.611。这个差距不是学术洁癖：把两套协方差分别拿去算最小方差组合，真实权重是 (1.04, −0.04)，经典协方差给出 (0.62, 0.38)——它把本该轻微做空的资产配成了 38% 的多头。MCD 给出 (0.95, 0.05)，基本还原。本文用纯 numpy 从零实现 Fast-MCD 的 C-step 迭代，并在最后拆穿『稳健估计=无脑抗噪』这个最常见的误解。

上图是全文的缩影：蓝点是正常样本（真实相关 0.7 的椭圆云），红叉是 10% 的反向污染。黄色虚线椭圆（经典协方差 2σ）被红叉硬生生扭成了近乎圆形甚至反向倾斜——相关结构被抹平了；绿色实线椭圆（MCD）几乎贴合黑色点线（真实结构），像没看见那些红叉一样。

一、样本协方差的原罪：崩溃点为零#

先厘清一个统计概念——崩溃点（breakdown point）：一个估计量在被多大比例的任意坏数据污染时，仍能给出有界的、有意义的结果。样本均值的崩溃点是 0：只要一个观测跑到无穷远，均值就跟着跑到无穷远。样本协方差同理。

这在量化里是实打实的伤害。协方差矩阵是三大件的共同地基：

• 组合优化：马科维茨最小方差、风险平价，全靠 Σ 求逆算权重
• 风险模型：VaR、成分风险贡献、压力测试
• 统计套利：配对/篮子的价差方差、协整残差的标准化

而金融数据天生带极端值：财报跳空、闪崩、流动性枯竭日的价格、数据源的错误 tick。这些点未必是”错误数据”，很多是真实但非典型的市场状态。问题在于：样本协方差对每个点的权重完全相同，一个 −8σ 的日子和一个 +0.1σ 的日子，在 np.cov 里贡献的地位一样重。极端值的平方项在协方差公式里被放大，于是几个点就足以主导整个估计。

看合成实验的数字：真实协方差对角线是 (1.0, 2.25)，经典估计给出 (2.22, 3.67)——方差被污染样本抬高了一倍多。相关系数从 0.7 崩到 −0.068，符号都反了。这不是”精度下降”，是”结论错误”。

二、MCD 的核心思想：找最紧致的那半个数据#

最小协方差行列式（Minimum Covariance Determinant, Rousseeuw 1984）的想法极其干脆：

与其让所有点平等地污染估计，不如只用一个”最干净”的子集来估协方差。

怎么定义”最干净”？用协方差矩阵的行列式。行列式 det(Σ) 是数据云”体积”的度量——数据越紧致，椭球体积越小，行列式越小。而极端值会把椭球撑大。所以：

在所有大小为 h 的样本子集中，找出使 det(Σ) 最小的那个子集，用它来估计均值和协方差。

极端值天然会撑大行列式，因此它们不会出现在”行列式最小”的那个子集里——它们被算法自动排除，而不是靠人工设阈值去剔。这就是 MCD 抗污染的本质。

h 的选择决定了崩溃点：h ≈ N/2 时崩溃点最高（可抗 50% 污染），但估计效率低；实务中常取 h ≈ 0.75N，在抗污染和统计效率间平衡。本文实验用 h = 0.75 × 400 = 300。

三、Fast-MCD 与 C-step：从组合爆炸到可计算#

暴力枚举所有 C(N, h) 个子集在计算上不可能（400 选 300 是天文数字）。Rousseeuw & Van Driessen (1999) 的 Fast-MCD 靠一个关键定理让它可算——C-step（Concentration step）：

给定当前子集算出的 (μ, Σ)，对全部 N 个点算马氏距离，取距离最小的 h 个点组成新子集。定理保证：新子集的 det(Σ) 一定不大于旧子集。于是从任意初始子集出发，反复 C-step，行列式单调下降直到收敛到局部最优。

import numpy as np

def mahalanobis_sq(X, mu, cov):
    d = X - mu
    inv = np.linalg.inv(cov)
    return np.einsum("ij,jk,ik->i", d, inv, d)

def cstep(X, subset_idx, h):
    """一次 concentration step：用当前子集估计 -> 取最近 h 个点"""
    mu = X[subset_idx].mean(0)
    cov = np.cov(X[subset_idx].T, bias=False)
    if np.linalg.det(cov) <= 1e-12:
        cov += 1e-6 * np.eye(X.shape[1])
    d2 = mahalanobis_sq(X, mu, cov)
    new_idx = np.argsort(d2)[:h]
    det = np.linalg.det(np.cov(X[new_idx].T, bias=False))
    return new_idx, det

单次 C-step 只能到局部最优，所以要多起点：随机抽 p+1 个点做初始估计，跑到收敛，记录行列式，重复几百次取全局最小。

def fast_mcd(X, h, n_trials=500, max_cstep=50):
    N, p = X.shape
    best_det, best_idx = np.inf, None
    for _ in range(n_trials):
        init = np.random.choice(N, p + 1, replace=False)
        mu0 = X[init].mean(0)
        cov0 = np.cov(X[init].T, bias=False)
        if np.linalg.det(cov0) <= 1e-12:
            cov0 += 1e-3 * np.eye(p)
        d2 = mahalanobis_sq(X, mu0, cov0)
        idx_c = np.argsort(d2)[:h]
        prev = np.inf
        for _ in range(max_cstep):
            idx_c, det = cstep(X, idx_c, h)
            if abs(prev - det) < 1e-10:
                break
            prev = det
        if det < best_det:
            best_det, best_idx = det, idx_c
    return best_idx

C-step 的单调性是 Fast-MCD 全部效率的来源：它把”在天文数字个子集里搜索”变成”从几百个起点各自梯度下降”，实测 500 起点在几秒内就能稳定命中全局最优子集。

四、一致性校正：别忘了这一步#

用 h 个”最紧致”的点估出的协方差，会系统性地偏小——因为你故意只留了内圈的点，把外圈正常的尾部也切掉了。直接拿这个 raw 估计去算马氏距离，阈值会算错。

修正办法是乘一个一致性因子 c，让 MCD 在纯高斯数据（无污染）下渐近无偏。标准做法基于卡方分布的中位数：

def chi2_median(p):
    # p 自由度卡方分布中位数的 Wilson-Hilferty 近似
    return p * (1 - 2 / (9 * p))**3

med_d2 = np.median(mahalanobis_sq(X, mu_mcd, cov_mcd_raw))
c = med_d2 / chi2_median(p)
cov_mcd = c * cov_mcd_raw   # 一致性校正后的 MCD 协方差

直觉：如果数据真是高斯的，理论上马氏距离平方的中位数应该等于 chi2_median(p)；实测中位数除以理论值就是缩放偏差，乘回去即校正。跳过这一步，你的稳健马氏距离阈值会系统性偏松，异常检测的准确率大打折扣。

五、实测：相关系数、异常检测、组合权重#

三个维度看 MCD vs 经典协方差在 10% 反向污染下的表现。

相关系数还原：

真实 0.700，经典 −0.068（符号翻转，彻底错误），MCD 0.611（有轻微低估但方向和量级都对）。0.611 vs 真值 0.7 的差距主要来自 h=0.75N 切掉了部分正常尾部——这是抗污染换来的必然代价，用更大的 h 会更接近真值但抗污染能力下降。

异常检测（距离-距离图）：

横轴是经典马氏距离，纵轴是 MCD 稳健马氏距离。经典距离（用被污染的协方差算）把污染点自己”拉进”了正常范围——它们参与了估计，自然显得不异常，只抓到 11/40。MCD 距离把 40 个污染点里的 27 个清晰甩到右上方超出阈值，误报仅 1 个。这就是稳健估计做异常检测的正确姿势：先用不受污染影响的估计量定义”正常”，再拿它去衡量每个点有多异常——而不是让嫌疑犯参与定义什么叫正常。

组合权重（最扎心的一环）：

把三套协方差分别代入最小方差组合公式 w = Σ⁻¹𝟙 / (𝟙ᵀΣ⁻¹𝟙)：

真实最优组合几乎满仓 A、轻微做空 B。经典协方差因为把相关结构算反了，给出 (0.62, 0.38)——把该做空的资产配成了 38% 的多头，这是一个方向性错误，不是量级偏差。真金白银的组合会因为几十个污染样本承受完全错误的暴露。MCD 给出 (0.95, 0.05)，虽不完美但没有方向性错误。

六、为什么最大化非高斯性能解出独立源——不，这是 MCD 不是 ICA#

（此处保持警醒：MCD 靠的是”最小体积子集”，与非高斯性无关。）MCD 的引擎是几何的——找体积最小的椭球；它不假设分布形态，只假设”大多数点来自同一个椭球分布，少数是污染”。这也划定了它的适用边界：如果污染比例超过 1−h/N，或者污染点本身聚成一个紧致的小簇（比自身正常数据还紧致），MCD 会把这个污染簇误当成主体。 这不是 bug，是任何基于”多数即正常”的稳健方法的共同软肋。

七、反面教材：稳健 ≠ 无脑抗噪#

MCD 很强，但有三个必须知道的陷阱，否则会误用。

陷阱一：h 选太大，照样被污染拖偏。 h 是抗污染能力的旋钮，不是越大越好也不是越小越好。如果你的数据污染比例是 15%，却设 h = 0.9N（只排除 10%），那剩下的 5% 污染仍会进入”干净子集”，估计照样被拉偏。h 的选择必须匹配你对污染比例的先验：宁可高估污染比例（h 小一点），也别让污染混进核心子集。 代价是 h 越小、统计效率越低、正常数据也被切掉越多（本文 0.611 vs 0.7 的低估就是这么来的）。

陷阱二：多起点不够，收敛到坏的局部最优。 C-step 只保证单调下降到局部最优，不保证全局。起点太少（比如只跑 10 次），很可能所有起点都落进被污染主导的盆地，输出一个”看起来收敛了但其实错的”结果。必须几百个随机起点，且初始子集要用 p+1 个点（最小的能定义协方差的子集）来保证多样性。

陷阱三：把 MCD 当成”检测异常就删掉”的预处理。 常见误用：跑 MCD 标出异常点 → 删掉 → 对剩下的数据用经典协方差。这在两种情况下危险：(1) 被标记的”异常”可能是真实的、信息量最大的市场状态（危机日的相关性飙升是真实的 regime，不是噪声），删掉等于把风险模型里最该关注的部分抹掉；(2) 硬删除引入选择偏差，剩余数据的尾部被人为截断，方差系统性偏小。MCD 的正确用法是产出稳健的 (μ, Σ) 本身，以及一个连续的异常度打分，而不是一刀切的删除决策。 是否剔除、如何降权，应该结合业务判断，而不是让算法替你决定哪天的市场”不算数”。

结语#

样本协方差的崩溃点为 0，这是量化风险模型里一个被长期低估的脆弱性。MCD 用一个几何上极干净的想法——找体积最小的核心子集——把崩溃点提到接近 50%，代价是牺牲一部分统计效率和一点计算量。本文的合成实验把这个价值量化到了组合权重层面：10% 的污染足以让经典协方差给出方向性错误的持仓，而 MCD 稳住了大方向。

但记住最后一节：稳健估计不是免死金牌。h 要匹配污染先验，起点要够多，更重要的是——不要把”统计上异常”简单等同于”应该删除”。金融数据里最异常的那些天，往往正是风险模型最该学会的那些天。

（本文所有数字来自纯 numpy 合成实验，随机种子固定，可复现。真实市场数据的污染结构比合成实验复杂，MCD 的表现会依赖污染是否满足”少数且分散”的假设。）