量化投资中的”假发现”陷阱
假设你测试了 1000 个量化因子，每次使用 \(\alpha = 0.05\) 的显著性水平：
- 即使所有因子都无效（\(H_0\) 全为真），期望的”显著”因子数 = \(1000 \times 0.05 = 50\) 个
- 至少出现 1 个假阳性的概率：
\[\text{FWER} = 1 - (1 - \alpha)^m = 1 - (1-0.05)^{1000} \approx 100\%\]
核心问题：当同时进行大量假设检验时，第一类错误（假阳性）发生的概率急剧膨胀。
为什么“显著”不等于“真的有用”
刚接触统计检验时，容易把“显著”直接理解成“这个结论已经很可靠”。实际上，至少要区分三层含义：
- 统计显著：结果不像纯噪声那样随机
- 经济显著：效应大小足以覆盖交易成本、冲击成本和风险暴露
- 样本外稳健：换一段数据、换一个市场后，结论仍大体成立
本章讨论的多重检验，主要是在修补第一层。当你一次测试很多候选因子时，原本看似显著的结果中会混入大量纯噪声。
为什么“单次 5%”不等于“整体 5%”
如果每次检验都使用 \(\alpha = 0.05\)，那么至少出现 1 次假阳性的概率会随着检验次数快速上升：
| 1 |
\(1 - (1-0.05)^1 = 5\%\) |
| 10 |
\(1 - (1-0.05)^{10} \approx 40\%\) |
| 100 |
\(1 - (1-0.05)^{100} \approx 99.4\%\) |
| 1000 |
\(1 - (1-0.05)^{1000} \approx 100\%\) |
所以，单次检验中的“5% 误报率”，绝不能直接理解为整轮筛选中的“5% 误报率”。
假设检验的基本框架回顾
| 不拒绝 \(H_0\) |
正确接受 (U) |
第二类错误 (T) |
| 拒绝 \(H_0\) |
第一类错误 (V) |
正确拒绝 (S) |
关键定义：
- 第一类错误（\(\alpha\)）：拒绝了一个真实的原假设（假阳性）
- 第二类错误（\(\beta\)）：未能拒绝一个虚假的原假设（假阴性）
- 检验效力 (Power)：\(1 - \beta\)，正确拒绝虚假原假设的概率
p 值不是“原假设为真”的概率
很多初学者最容易误解的是 p 值的含义。
p 值真正表示的是：
- 在原假设 \(H_0\) 为真的前提下
- 观察到当前样本，或更极端样本结果的概率
p 值并不表示：
- 原假设为真的概率
- 这个因子未来继续有效的概率
- 这个结果有多大的经济价值
因此，p 值小只能说明“这份样本与原假设不太相容”，并不自动推出“这个结论一定真实有效”。
两类错误的经济学含义
统计学的两类错误：
| 判定”有效” |
I类错误（假阳性） |
✅ 正确发现 |
| 判定”无效” |
✅ 正确排除 |
II类错误（假阴性） |
在量化投资中的代价：
- I类错误（假发现）：把无效因子当作Alpha → 亏钱
- II类错误（遗漏）：错过真实Alpha → 少赚钱
非对称损失：
- 在学术发表中，I类错误代价更大（错误结论难以撤回）
- 在量化投资中，I类错误直接导致资金损失
- 因此，多重检验校正在金融中格外重要
模拟实验：1000 个无效因子的检验
Code
import numpy as np # 导入numpy用于数值计算
import matplotlib.pyplot as plt # 导入绑图库
from scipy import stats # 导入统计模块
plt.rcParams['font.family'] = ['Source Han Serif SC', 'Noto Serif CJK SC', 'SimSun'] # 设置中文字体优先级
plt.rcParams['axes.unicode_minus'] = False # 解决负号显示问题
np.random.seed(42) # 设置随机种子保证可重现
num_tests = 1000 # 假设检验次数（1000个因子）
null_test_statistics = np.random.standard_normal(num_tests) # 从标准正态分布生成检验统计量（所有H0为真）
null_p_values = 2 * (1 - stats.norm.cdf(np.abs(null_test_statistics))) # 计算双侧p值
num_false_positives = np.sum(null_p_values < 0.05) # 统计p<0.05的假阳性个数
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5)) # 创建1行2列子图
# 左图：p值直方图
axes[0].hist(null_p_values, bins=50, color='#008080', alpha=0.7, edgecolor='white') # 绘制p值频率直方图
axes[0].axvline(x=0.05, color='#E3120B', linestyle='--', linewidth=2, label='α = 0.05') # 标注显著性阈值
axes[0].set_xlabel('p 值') # x轴标签
axes[0].set_ylabel('频次') # y轴标签
axes[0].set_title(f'p值分布（假阳性: {num_false_positives} 个）') # 标题（含假阳性数）
axes[0].legend() # 显示图例
# 右图：-log10(p)散点图
axes[1].scatter(range(num_tests), -np.log10(null_p_values), s=5, alpha=0.5, color='gray') # 绘制-log10(p)散点
axes[1].axhline(y=-np.log10(0.05), color='#E3120B', linestyle='--', linewidth=2, label='α = 0.05') # 显著性阈值线
axes[1].set_xlabel('因子编号') # x轴标签
axes[1].set_ylabel('$-\\log_{10}(p)$') # y轴标签
axes[1].set_title('所有因子都无效，但假阳性大量出现') # 标题
axes[1].legend() # 显示图例
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
一批检验为什么要被看成同一个 family
多重检验里的 family，可以理解为“应当一起结算错误预算的一组检验”。
常见例子：
- 对同一只股票测试 50 条均线策略
- 对 80 多个行业同时检验 Alpha
- 对 300 个候选因子做一轮统一筛选
只要这些检验共同服务于同一轮研究结论或投资决策，就应该被看成同一个 family，不能只盯着单个 p 值解释。
FWER 与 FDR：两种不同的错误控制策略
Bonferroni 的直觉是“把总犯错预算切成 m 份”
把总错误预算 \(\alpha\) 理解成一笔固定预算：
- 如果只做 1 次检验，这 0.05 可以全部给这 1 次
- 如果做 10 次检验，每次最多分到 0.005
- 如果做 100 次检验，每次最多分到 0.0005
| 5 |
0.0100 |
| 50 |
0.0010 |
| 500 |
0.0001 |
这能保护整体错误率，但代价是：检验越多，每次分到的“单次犯错额度”越小，真信号也会越来越难通过。
FWER控制的数学原理
FWER（族错误率）：
\[
\text{FWER} = P(\text{至少一个假阳性}) = 1 - (1-\alpha)^m
\]
当 \(m\) 很大时，FWER趋近于1
例子：\(\alpha = 0.05\)，测试 \(m = 100\) 个假设
\[
\text{FWER} = 1 - (1-0.05)^{100} \approx 0.994
\]
几乎必然出现假阳性！
FDR（假发现率）：
\[
\text{FDR} = E\left[\frac{\text{假阳性数}}{\text{所有拒绝数}}\right] = E\left[\frac{V}{R}\right]
\]
- FDR关注的是”在你宣称的发现中，有多少比例是假的”
- FDR \(\leq\) FWER，控制FDR更宽松但更实用
Bonferroni：最简单的 FWER 控制方法
核心思想：将显著性水平均匀分配给每一次检验。
\[\alpha_{\text{adjusted}} = \frac{\alpha}{m}\]
等价地，对 p 值进行调整：
\[p_{\text{adjusted}} = \min(m \times p_{\text{原始}}, \; 1)\]
优点：简单、直观、对任何依赖结构都有效
缺点：当 \(m\) 很大时极度保守——例如 \(m = 1000\) 时，\(\alpha_{\text{adjusted}} = 0.00005\)
案例：A 股强势股票筛选
import pandas as pd # 导入pandas用于数据操作
import os # 导入os用于路径处理
DATA_DIR = 'C:/qiufei/data' if os.name == 'nt' else '/home/ubuntu/r2_data_mount/qiufei/data' # 跨平台数据路径
# 读取前复权股价数据
stock_price_data = pd.read_hdf(os.path.join(DATA_DIR, 'stock/stock_price_pre_adjusted.h5')).reset_index() # 读取并重置索引
stock_price_2020 = stock_price_data[ # 筛选2020年度数据
(stock_price_data['date'] >= '2020-01-01') & (stock_price_data['date'] < '2021-01-01') # 限定时间范围
] # 用于强势股票筛选分析
# 计算日收益率
stock_price_2020 = stock_price_2020.copy() # 复制以避免SettingWithCopyWarning
stock_price_2020['daily_return'] = stock_price_2020.groupby('order_book_id')['close'].pct_change() # 按股票计算日收益率
# 筛选交易日数 >= 100 的股票
valid_stock_counts = stock_price_2020.groupby('order_book_id')['daily_return'].count() # 各股票的有效交易日数
valid_stock_list = valid_stock_counts[valid_stock_counts >= 100].index.tolist() # 筛选数据充足的股票
print(f'2020年有效股票数量: {len(valid_stock_list)}') # 输出有效股票数
未校正 vs Bonferroni 校正效果对比
Code
from scipy import stats # 导入统计库
from statsmodels.stats.multitest import multipletests # 导入多重检验校正工具
# 对每只股票执行单侧 t 检验 (H1: μ > 0)
p_value_list = [] # 存储原始p值
t_stat_list = [] # 存储t统计量
for stock_code in valid_stock_list[:500]: # 取前500只股票（控制计算时间）
stock_returns = stock_price_2020[stock_price_2020['order_book_id'] == stock_code]['daily_return'].dropna() # 提取收益率
t_statistic, two_sided_p = stats.ttest_1samp(stock_returns, 0) # 单样本t检验
one_sided_p = two_sided_p / 2 if t_statistic > 0 else 1 - two_sided_p / 2 # 转为单侧p值
p_value_list.append(one_sided_p) # 记录p值
t_stat_list.append(t_statistic) # 记录t统计量
p_values_array = np.array(p_value_list) # 转为numpy数组
# Bonferroni 校正
bonferroni_reject, bonferroni_p_adjusted, _, _ = multipletests(p_values_array, alpha=0.05, method='bonferroni') # Bonferroni校正
fig, axes = plt.subplots(1, 2, figsize=(13, 5)) # 创建1行2列子图
# 左图：未校正
is_significant_raw = p_values_array < 0.05 # 未校正下的显著性判断
axes[0].scatter(range(len(p_values_array)), -np.log10(p_values_array), # 绘制-log10(p)散点
c=['#E3120B' if s else 'gray' for s in is_significant_raw], s=8, alpha=0.6) # 显著=红色
axes[0].axhline(y=-np.log10(0.05), color='blue', linestyle='--', label='α=0.05') # 原始阈值线
axes[0].set_title(f'未校正（显著: {sum(is_significant_raw)} 只）') # 标题
axes[0].set_xlabel('股票编号') # x轴
axes[0].set_ylabel('$-\\log_{10}(p)$') # y轴
axes[0].legend() # 图例
# 右图：Bonferroni 校正后
axes[1].scatter(range(len(p_values_array)), -np.log10(p_values_array), # 同样的-log10(p)
c=['#E3120B' if s else 'gray' for s in bonferroni_reject], s=8, alpha=0.6) # Bonferroni显著=红色
bonferroni_threshold = 0.05 / len(p_values_array) # 计算Bonferroni校正阈值
axes[1].axhline(y=-np.log10(bonferroni_threshold), color='blue', linestyle='--', label=f'α/m={bonferroni_threshold:.2e}') # 校正阈值线
axes[1].set_title(f'Bonferroni 校正（显著: {sum(bonferroni_reject)} 只）') # 标题
axes[1].set_xlabel('股票编号') # x轴
axes[1].set_ylabel('$-\\log_{10}(p)$') # y轴
axes[1].legend() # 图例
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
为什么 BH 一定要先给 p 值排序
BH 方法的阈值不是看“这个假设是谁”，而是看“这个 p 值排第几名”。
先排序的原因：
- 最小的 p 值最像真信号
- 中间的 p 值真假难分
- 最大的 p 值更像噪声
BH 的做法是：
- 按从小到大排序
- 给第 1 名最严格但仍可接受的阈值
- 随着名次往后，阈值线性放宽
- 找到“最后一个还能过线的位置”
所以，BH 不是逐个独立判定，而是在整批结果里寻找一条整体截断线。
BH 实战：注入 Alpha 信号的检测
# 使用2021年数据
stock_price_2021 = stock_price_data[ # 筛选2021年数据
(stock_price_data['date'] >= '2021-01-01') & (stock_price_data['date'] < '2022-01-01') # 限定时间范围
].copy() # 创建副本
stock_price_2021['daily_return'] = stock_price_2021.groupby('order_book_id')['close'].pct_change() # 计算日收益率
# 选取数据充足的股票
valid_2021_counts = stock_price_2021.groupby('order_book_id')['daily_return'].count() # 各股票交易日数
valid_2021_stocks = valid_2021_counts[valid_2021_counts >= 100].index.tolist() # 筛选有效股票
# 随机注入100只股票的正Alpha（+0.5%/天）
np.random.seed(2024) # 设置随机种子
injected_alpha_stocks = np.random.choice(valid_2021_stocks, size=100, replace=False) # 随机选100只注入Alpha
injected_stock_set = set(injected_alpha_stocks) # 转为集合便于查询
# 为注入Alpha的股票修改收益率
stock_price_2021_modified = stock_price_2021.copy() # 复制数据
alpha_mask = stock_price_2021_modified['order_book_id'].isin(injected_stock_set) # 标记被注入的股票
stock_price_2021_modified.loc[alpha_mask, 'daily_return'] += 0.005 # 注入每日+0.5%的正Alpha
print(f'2021年有效股票总数: {len(valid_2021_stocks)}') # 输出总股票数
print(f'注入 Alpha 的股票数: {len(injected_alpha_stocks)}') # 输出注入Alpha的股票数
2021年有效股票总数: 4428
注入 Alpha 的股票数: 100
三种方法对比：BH 的优势一目了然
# 计算各方法的关键指标
def calculate_performance_metrics(reject_decisions, true_labels): # 定义性能计算函数
true_positive = np.sum(reject_decisions & true_labels) # 真阳性（正确拒绝）
false_positive = np.sum(reject_decisions & ~true_labels) # 假阳性（错误拒绝）
total_rejections = np.sum(reject_decisions) # 总拒绝数
fdr = false_positive / total_rejections if total_rejections > 0 else 0 # 错误发现率
tpr = true_positive / np.sum(true_labels) if np.sum(true_labels) > 0 else 0 # 真阳性率（检验效力）
return {'TP': true_positive, 'FP': false_positive, 'R': total_rejections, 'FDR': f'{fdr:.2%}', 'TPR': f'{tpr:.2%}'} # 返回指标字典
results_comparison = pd.DataFrame({ # 创建对比表
'未校正 (α=0.05)': calculate_performance_metrics(raw_reject, truth_arr), # 未校正结果
'Bonferroni': calculate_performance_metrics(bonf_reject, truth_arr), # Bonferroni结果
'BH (q=0.05)': calculate_performance_metrics(bh_reject, truth_arr) # BH结果
}) # 三列三行的性能对比表
results_comparison # 显示对比表格
可视化：BH 自适应阈值线与火山图
Code
fig, axes = plt.subplots(1, 2, figsize=(13, 5)) # 创建1行2列子图
# 左图：BH 自适应阈值线
sorted_indices = np.argsort(p_arr) # 按p值从小到大排序的索引
sorted_p = p_arr[sorted_indices] # 排序后的p值
m_total = len(p_arr) # 总检验次数
bh_threshold_line = np.arange(1, m_total + 1) / m_total * 0.05 # BH阈值线: k/(m*q)，线性递增
axes[0].scatter(range(m_total), sorted_p, s=3, c='gray', alpha=0.5, label='排序后的p值') # 排序p值散点图
axes[0].plot(range(m_total), bh_threshold_line, color='#008080', linewidth=2, label='BH 阈值线') # BH阈值线
axes[0].set_xlabel('排序位置 k') # x轴
axes[0].set_ylabel('p 值') # y轴
axes[0].set_title('BH 自适应阈值线') # 标题
axes[0].set_ylim(0, 0.1) # 限制y轴范围便于观察
axes[0].legend() # 图例
# 右图：火山图
for i in range(len(t_arr)): # 遍历每个检验
color = 'gray' # 默认灰色（不显著）
if bh_reject[i] and t_arr[i] > 0: # BH显著且t为正
color = '#E3120B' # 红色（显著正效应）
elif bh_reject[i] and t_arr[i] < 0: # BH显著且t为负
color = '#008080' # 青色（显著负效应）
axes[1].scatter(t_arr[i], -np.log10(p_arr[i]), s=6, c=color, alpha=0.5) # 绘制散点
axes[1].set_xlabel('t 统计量') # x轴
axes[1].set_ylabel('$-\\log_{10}(p)$') # y轴
axes[1].set_title('火山图（BH 校正）') # 标题
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
校正后的 p 值应该怎样解释
原始 p 值回答的是：单个检验在原假设下有多极端。
校正后的 p 值回答的是：当你把这组检验放在一起看时，它是否仍然足够显著。
阅读规则：
- 校正后 p 值通常会比原始 p 值更大
- 如果校正后仍小于阈值，说明它在“整轮筛选”里仍站得住
- 如果原始显著但校正后不显著，说明它更像“筛选过程里的幸运结果”
报告结果时，最好同时给出原始 p 值与校正后 p 值，而不是只汇报最漂亮的那个版本。
P-Hacking 与数据偷窥的危害
P-Hacking（p值黑客）是量化金融研究中的严重问题：
- 在大量回测中”挖掘”出看起来显著的策略
- Harvey, Liu & Zhu (2016) 发现：超过一半的金融”异象因子”可能是假的
- 他们建议将 t 统计量的显著性阈值提高到 t > 3.0
常见的数据偷窥手法：
- 在同一数据集上测试大量策略，只报告显著的
- 不断调整回测参数直到出现好结果
- 对同一假设尝试多种统计检验方法
防范措施：
- 预先注册研究计划和假设
- 样本外验证：将数据分为训练集和测试集
- 报告所有结果，包括不显著的检验
- 使用 BH 或 Bonferroni 进行系统性的多重检验校正
样本外验证为什么不能替代多重检验校正
很多同学会问：‘我已经做了样本外验证，还需要多重检验校正吗？’
答案通常是：仍然需要。
| 多重检验校正 |
同一轮筛选里测试太多，假阳性累积 |
| 样本外验证 |
历史规律能否跨时间继续成立 |
这两个工具不能互相替代，因为：
- 你可能在样本内先试了很多策略，已经发生了选择偏差
- 即使最后只拿 1 个策略去做样本外，前面的“筛选运气”仍可能残留
- 最稳妥的做法是：先做多重检验控制，再做样本外验证
案例：海康威视技术指标的多重检验
Code
# 获取海康威视(002415.XSHE)的股价数据
hikvision_data = stock_price_data[ # 筛选海康威视
(stock_price_data['order_book_id'] == '002415.XSHE') & # 股票代码
(stock_price_data['date'] >= '2020-01-01') # 2020年以来
].copy().sort_values('date') # 按日期排序
hikvision_data['daily_return'] = hikvision_data['close'].pct_change() # 计算日收益率
# 测试 50 条均线策略 (MA5 到 MA54)
ma_periods = range(5, 55) # MA周期从5到54
strategy_p_values = [] # 存储各策略的p值
for ma_period in ma_periods: # 遍历每个MA周期
hikvision_data[f'MA_{ma_period}'] = hikvision_data['close'].rolling(ma_period).mean() # 计算移动平均
signal_mask = hikvision_data['close'] > hikvision_data[f'MA_{ma_period}'] # 价格在均线上方时持仓
strategy_returns = hikvision_data.loc[signal_mask, 'daily_return'].dropna() # 持仓期间的收益率
if len(strategy_returns) > 30: # 样本量充足时
t_stat, p_val = stats.ttest_1samp(strategy_returns, 0) # t检验：均值是否>0
strategy_p_values.append(p_val / 2 if t_stat > 0 else 1) # 转为单侧p值
else: # 样本量不足
strategy_p_values.append(1.0) # 设为不显著
strategy_p_arr = np.array(strategy_p_values) # 转numpy数组
bonf_reject_ma, _, _, _ = multipletests(strategy_p_arr, alpha=0.05, method='bonferroni') # Bonferroni校正
bh_reject_ma, _, _, _ = multipletests(strategy_p_arr, alpha=0.05, method='fdr_bh') # BH校正
fig, ax = plt.subplots(figsize=(11, 5)) # 创建图表
x_pos = list(ma_periods) # x轴位置
ax.bar(x_pos, -np.log10(strategy_p_arr), color=[ # 绘制-log10(p)柱状图
'#E3120B' if strategy_p_arr[i] < 0.05 else 'gray' for i in range(len(strategy_p_arr)) # 未校正显著=红色
], alpha=0.7) # 透明度0.7
ax.axhline(y=-np.log10(0.05), color='blue', linestyle='--', label='未校正 α=0.05') # 未校正阈值
ax.axhline(y=-np.log10(0.05/len(strategy_p_arr)), color='#E3120B', linestyle='-', label='Bonferroni') # Bonferroni阈值
ax.set_xlabel('移动平均周期 (MA)') # x轴
ax.set_ylabel('$-\\log_{10}(p)$') # y轴
ax.set_title(f'海康威视50条MA策略: 未校正显著 {sum(strategy_p_arr < 0.05)} 条 → 校正后 {sum(bonf_reject_ma)} 条') # 标题
ax.legend() # 图例
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
关键发现：大量未校正下”显著”的策略，在 Bonferroni 或 BH 校正后全部变为不显著 → 典型的数据偷窥结果！
案例：行业 Alpha 的多重检验
Code
# 三种校正
industry_raw_sig = industry_p_arr < 0.05 # 未校正显著判断
industry_bonf_rej, _, _, _ = multipletests(industry_p_arr, alpha=0.05, method='bonferroni') # Bonferroni
industry_bh_rej, _, _, _ = multipletests(industry_p_arr, alpha=0.05, method='fdr_bh') # BH
fig, axes = plt.subplots(1, 3, figsize=(15, 5)) # 1行3列子图
method_configs = [ # 三种方法的配置
('未校正 (α=0.05)', industry_raw_sig, int(sum(industry_raw_sig))), # 未校正配置
('Bonferroni', industry_bonf_rej, int(sum(industry_bonf_rej))), # Bonferroni配置
('BH (q=0.05)', industry_bh_rej, int(sum(industry_bh_rej))) # BH配置
] # 三种方法
for ax_idx, (method_name, reject_mask, count) in enumerate(method_configs): # 遍历三种方法
colors = ['#E3120B' if r else 'gray' for r in reject_mask] # 显著=红色，不显著=灰色
axes[ax_idx].barh(range(len(industry_p_arr)), -np.log10(industry_p_arr), color=colors, alpha=0.7) # 水平柱状图
axes[ax_idx].set_title(f'{method_name}（显著: {count}）') # 标题含显著数量
axes[ax_idx].set_xlabel('$-\\log_{10}(p)$') # x轴标签
plt.tight_layout() # 自动调整间距
plt.show() # 显示图表
方法选择决策框架
多重检验在不同金融场景的应用
场景1：因子发现
- 测试上百个潜在Alpha因子
- 推荐：BH方法，控制FDR在5%-10%
- 目标：在控制假发现率的前提下，尽量多发现有效因子
场景2：风险监控
- 同时监控多个交易策略的风险指标
- 推荐：Bonferroni，控制FWER
- 理由：任何一个假警报都可能导致不必要的平仓
场景3：学术研究
- 资产定价实证中的多重假设检验
- 推荐：Holm方法或BH方法
- 需要在论文中明确报告校正方法
核心原则：测试的假设越多，校正越重要
