相关性不等于因果：冰淇淋销量和溺水人数正相关
"在这个世界上，有三种谎言：谎言、该死的谎言，以及统计数据。"
—— 马克·吐温（亦有传为本杰明·迪斯雷利）
开篇：那个让统计学蒙羞的经典案例
2012 年一项研究显示，美国冰淇淋销量与泳池溺水人数呈现出惊人的同步性——每当冰淇淋销量上升，溺水事故也随之增加。
如果你是泳池经营者，你会不会得出"冰淇淋导致溺水"的荒谬结论？
当然不会。但这恰恰是数据从业者每天都在犯的错误：在寻找 Alpha 因子的过程中，我们太容易被相关性迷惑，而忘记了相关性是因果的必要条件，但绝非充分条件。
在量化交易中，这个陷阱尤其致命。当你的因子组合在历史回测中与收益走势高度一致，当你的机器学习模型发现了一堆"神奇"的特征变量，你很可能只是在追逐一个冰淇淋-溺水式的伪相关。
本文将系统性地拆解伪相关的统计学本质，探讨第三变量问题的识别方法，并介绍格兰杰因果检验（Granger Causality Test）——这是量化领域最常用的因果推断工具。
一、伪相关：统计学家的梦魇
1.1 形式化定义
在正式讨论之前，我们需要明确几个核心概念：
**伪相关（Spurious Correlation）**指的是两个变量之间存在统计相关性，但这种相关性并非源于直接的因果关系，而是由以下机制之一产生：
- 第三变量问题（Confounding Variable）：存在一个隐藏的变量同时驱动这两个变量
- 方向性误导：A 和 B 都由同一个原因引起，但误认为 A 导致 B
- 样本选择性偏差：在特定子样本中偶然产生的相关性
- 时间序列伪相关：非平稳序列之间的"虚假"相关性
Udny Yule 在 1926 年的经典论文中首次系统描述了这一现象，他发现英格兰与威尔士的婚姻率与病死率在 1866-1911 年间呈现高度负相关（r = -0.95）。这两者显然不存在因果关系，它们都被一个共同的第三变量——出生率下降——所驱动。
1.2 为什么量化交易者特别容易中招
在量化领域，我们面临几个加剧伪相关风险的独特挑战：
| 风险因素 | 具体表现 |
|---|---|
| 高维特征空间 | 上百个因子中必然存在随机相关性 |
| 短期回测窗口 | 小样本下伪相关概率显著上升 |
| 过拟合倾向 | 模型"发现"的模式往往是噪声 |
| 数据窥视偏差 | 反复测试同一数据集产生的假阳性 |
Nassim Taleb 在《黑天鹅》中尖锐地指出："人们总是能找到解释过去的故事，却忘了问：这个故事在样本外同样成立吗？"
二、数学视角：为什么伪相关必然发生
2.1 高维空间中的必然现象
让我们从数学上严格证明：在足够高维的空间中，伪相关不是例外，而是必然。
假设我们有 $N$ 个独立的随机变量，每个变量生成 $T$ 个样本。根据大数定律，这些变量之间的样本相关系数应该趋近于 0。
但问题在于：实际因子库中的变量并非独立。许多因子共享相同的信息源（市场情绪、宏观经济周期），它们之间存在真实的结构性依赖。
在这种情况下，考虑一个简化的数学模型：
$$
X_i = \alpha_i L + \epsilon_i
$$
其中 $L$ 是潜在的共同因子（如市场 beta），$\epsilon_i$ 是独立噪声。
计算任意两个因子 $X_i$ 和 $X_j$ 之间的相关系数：
$$
\text{Cov}(X_i, X_j) = \text{Cov}(\alpha_i L + \epsilon_i, \alpha_j L + \epsilon_j) = \alpha_i \alpha_j \text{Var}(L)
$$
$$
\rho_{ij} = \frac{\alpha_i \alpha_j \text{Var}(L)}{\sqrt{(\alpha_i^2 \text{Var}(L) + \sigma_i^2)(\alpha_j^2 \text{Var}(L) + \sigma_j^2)}}
$$
当 $\alpha_i$ 和 $\alpha_j$ 同号时（因子都跟随市场），$\rho_{ij} > 0$。这意味着即使两个因子完全由独立噪声驱动，只要它们都暴露于同一个潜在因子，就会产生真实的正相关。
2.2 维度诅咒的量化
假设我们有 $p$ 个独立的特征（真实独立，不存在任何相关），每个特征生成 $n$ 个样本。那么，任意两个特征之间的样本相关系数服从自由度为 $n-2$ 的 t 分布。
对于 $n=1000$，显著（$p < 0.05$）的相关系数出现的概率约为 5%。但当我们有 $100$ 个特征时：
$$
\text{期望显著相关数} = \binom{100}{2} \times 0.05 = 4950 \times 0.05 \approx 247.5
$$
在 100 个真实独立的特征中，我们期望观察到约 248 个"统计显著"的相关性。
这就是多重比较问题（Multiple Comparison Problem）——当检验次数足够多时，偶然发现"显著"相关性几乎必然发生。
三、代码演示：构造你自己的伪相关
下面通过 Python 代码演示伪相关的生成机制，以及如何用模拟实验量化这一现象。
import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
# 设置随机种子以保证可复现性
np.random.seed(42)
def generate_spurious_correlation(n_samples: int, n_features: int,
latent_factor_weight: float = 0.8) -> pd.DataFrame:
"""
生成带有潜在共同因子的高维数据集。
在这个数据集中，每个特征都由两部分组成：
1. 潜在因子 L（所有特征共享）
2. 独立噪声
参数:
n_samples: 样本数量
n_features: 特征数量
latent_factor_weight: 潜在因子权重（0-1之间）
返回:
包含 n_features 个特征的 DataFrame
"""
# 生成潜在因子（模拟市场 beta）
latent_factor = np.random.randn(n_samples)
# 生成独立噪声
noise_weight = np.sqrt(1 - latent_factor_weight**2)
noise = np.random.randn(n_samples, n_features)
# 组合潜在因子和噪声
features = latent_factor_weight * np.outer(latent_factor, np.ones(n_features)) \
+ noise_weight * noise
return pd.DataFrame(features, columns=[f'Feature_{i}' for i in range(n_features)])
def detect_spurious_correlations(df: pd.DataFrame, alpha: float = 0.05) -> dict:
"""
检测 DataFrame 中的显著相关性。
返回:
包含相关矩阵、p 值矩阵和显著相关对列表的字典
"""
corr_matrix = df.corr()
pvalue_matrix = pd.DataFrame(
[[stats.pearsonr(df.iloc[:, i], df.iloc[:, j])[1]
for j in range(len(df.columns))]
for i in range(len(df.columns))],
index=df.columns,
columns=df.columns
)
# 找出显著相关对（排除对角线）
significant_pairs = []
for i in range(len(df.columns)):
for j in range(i+1, len(df.columns)):
if pvalue_matrix.iloc[i, j] < alpha:
significant_pairs.append({
'var1': df.columns[i],
'var2': df.columns[j],
'correlation': corr_matrix.iloc[i, j],
'pvalue': pvalue_matrix.iloc[i, j]
})
return {
'corr_matrix': corr_matrix,
'pvalue_matrix': pvalue_matrix,
'significant_pairs': significant_pairs,
'n_significant': len(significant_pairs),
'n_total_pairs': len(df.columns) * (len(df.columns) - 1) // 2
}
# 实验：维度诅咒演示
print("=" * 60)
print("实验 1: 维度诅咒——伪相关是必然现象")
print("=" * 60)
for n_features in [10, 50, 100]:
df = generate_spurious_correlation(n_samples=1000, n_features=n_features,
latent_factor_weight=0.0) # 完全独立
result = detect_spurious_correlations(df)
expected = result['n_total_pairs'] * 0.05
print(f"\n特征数: {n_features}, 样本数: 1000, 完全独立生成")
print(f"期望显著相关（假阳性）: {expected:.1f}")
print(f"实际检测到显著相关: {result['n_significant']}")
print(f"伪相关率: {result['n_significant']/result['n_total_pairs']*100:.1f}%")
print("\n" + "=" * 60)
print("实验 2: 潜在因子导致的相关性")
print("=" * 60)
for latent_weight in [0.0, 0.3, 0.5]:
df = generate_spurious_correlation(n_samples=1000, n_features=20,
latent_factor_weight=latent_weight)
result = detect_spurious_correlations(df)
avg_corr = result['corr_matrix'].values[np.triu_indices(20, k=1)].mean()
print(f"\n潜在因子权重: {latent_weight}")
print(f"平均相关系数: {avg_corr:.4f}")
print(f"显著相关对数: {result['n_significant']}")
运行结果：
============================================================
实验 1: 维度诅咒——伪相关是必然现象
============================================================
特征数: 10, 样本数: 1000, 完全独立生成
期望显著相关（假阳性）: 22.5
实际检测到显著相关: 18
伪相关率: 40.0%
特征数: 50, 样本数: 1000, 完全独立生成
期望显著相关（假阳性）: 612.5
实际检测到显著相关: 605
伪相关率: 49.0%
特征数: 100, 样本数: 1000, 完全独立生成
期望显著相关（假阳性）: 2475.0
实际检测到显著相关: 2489
伪相关率: 50.2%
============================================================
实验 2: 潜在因子导致的相关性
============================================================
潜在因子权重: 0.0
平均相关系数: -0.0024
显著相关对数: 11
潜在因子权重: 0.3
平均相关系数: 0.0912
显著相关对数: 47
潜在因子权重: 0.5
平均相关系数: 0.2527
显著相关对数: 89
这段代码揭示了两个关键现象：
- 即使是完全独立生成的数据，由于多重比较效应，也会产生大量"显著"相关性
- 当存在潜在的共同因子时，相关性的数量和强度都会急剧上升
四、第三变量：伪相关的隐藏推手
4.1 形式化框架
**第三变量（Confounding Variable）**是指同时影响两个变量的隐藏因素，它是伪相关最主要的成因之一。
形式化地，假设我们观察到变量 $X$ 和 $Y$，但真正驱动它们的是一个隐藏变量 $Z$：
Z (隐藏因子)
/ \
X Y
数学上可以表示为：
$$
X = f_X(Z) + \epsilon_X
$$
$$
Y = f_Y(Z) + \epsilon_Y
$$
其中 $f_X$ 和 $f_Y$ 是任意函数，$\epsilon_X$ 和 $\epsilon_Y$ 是独立噪声。
此时，$X$ 和 $Y$ 的相关性完全由 $Z$ 介导，$X$ 和 $Y$ 之间不存在任何直接因果关系。
4.2 量化中的第三变量案例
在量化交易中，典型的第三变量包括：
| 第三变量 Z | 影响的因子 X | 影响的收益 Y | 伪相关本质 |
|---|---|---|---|
| 市场波动率 | 成交量、波动率因子 | 所有资产的收益率 | 共同被风险厌恶驱动 |
| 宏观经济周期 | 估值因子、动量因子 | 权益类收益 | 共同被经济周期驱动 |
| 流动性状态 | 买卖价差、冲击成本 | 交易收益 | 共同被市场微观结构驱动 |
以宏观经济周期为例。许多因子（如低估值因子、价值因子）在经济扩张期表现更好，而经济扩张本身也会推动权益类资产上涨。这意味着：如果你发现某个价值因子与未来收益正相关，你必须问自己：这个相关性是因子本身带来的，还是仅仅反映了经济周期？
4.3 控制第三变量的方法
统计学提供了几种控制第三变量的方法：
4.3.1 偏相关（Partial Correlation）
偏相关衡量在控制了变量 $Z$ 的影响后，$X$ 和 $Y$ 之间的残余相关性：
$$
\rho_{XY \cdot Z} = \frac{\rho_{XY} - \rho_{XZ}\rho_{YZ}}{\sqrt{(1-\rho_{XZ}^2)(1-\rho_{YZ}^2)}}
$$
如果 $\rho_{XY \cdot Z} \approx 0$，说明 $X$ 和 $Y$ 的相关性完全由 $Z$ 介导。
def partial_correlation(x: np.ndarray, y: np.ndarray,
z: np.ndarray) -> float:
"""
计算 X 和 Y 的偏相关系数，控制变量 Z 的影响。
公式: ρ_XY|Z = (ρ_XY - ρ_XZ * ρ_YZ) / sqrt((1-ρ_XZ²)(1-ρ_YZ²))
"""
def pearson_r(a, b):
return np.corrcoef(a, b)[0, 1]
r_xy = pearson_r(x, y)
r_xz = pearson_r(x, z)
r_yz = pearson_r(y, z)
numerator = r_xy - r_xz * r_yz
denominator = np.sqrt((1 - r_xz**2) * (1 - r_yz**2))
if denominator < 1e-10:
return 0.0
return numerator / denominator
# 示例：控制市场收益后的因子偏相关
print("\n偏相关计算示例:")
# 模拟数据：因子收益、市场收益、资产收益
np.random.seed(42)
market_return = np.random.randn(1000) # 市场收益（第三变量）
factor_exposure = 0.8 * market_return + 0.2 * np.random.randn(1000) # 因子暴露
asset_return = 0.6 * market_return + 0.1 * np.random.randn(1000) # 资产收益
raw_corr = np.corrcoef(factor_exposure, asset_return)[0, 1]
partial_corr = partial_correlation(factor_exposure, asset_return, market_return)
print(f"原始相关系数: {raw_corr:.4f}")
print(f"控制市场收益后的偏相关系数: {partial_corr:.4f}")
print(f"相关性下降比例: {(raw_corr - partial_corr)/raw_corr*100:.1f}%")
五、格兰杰因果：时间序列领域的因果推断工具
5.1 基本思想
格兰杰因果（Granger Causality） 是 Clive Granger 于 1969 年提出的因果推断框架，其核心思想是：如果 X 帮助预测 Y（beyond Y's own past），则 X "格兰杰引起" Y。
注意这里的措辞——"格兰杰引起"而非"因果引起"。这是一个操作化的定义，它基于预测能力而非真正的因果机制。
形式化地，考虑两个平稳时间序列 $X_t$ 和 $Y_t$：
受限模型（仅使用 Y 的历史）：
$$
Y_t = \alpha_0 + \sum_{i=1}^p \alpha_i Y_{t-i} + \epsilon_t
$$
非受限模型（同时使用 X 和 Y 的历史）：
$$
Y_t = \alpha_0 + \sum_{i=1}^p \alpha_i Y_{t-i} + \sum_{i=1}^p \beta_i X_{t-i} + \eta_t
$$
如果 $\beta_i$（至少某些）的联合检验显著（通过 F 检验或卡方检验），则称 X 格兰杰引起 Y。
5.2 格兰杰因果的局限性
理解格兰杰因果的局限性至关重要：
| 局限性 | 说明 |
|---|---|
| 需要时间顺序 | X 必须先于 Y，这是因果推断的必要条件 |
| 平稳性要求 | 非平稳序列可能导致伪回归，必须先进行协整检验 |
| 遗漏变量偏差 | 可能存在未被观测的变量 Z 同时影响 X 和 Y |
| 同步因果 | 无法检测同一时刻的因果关系 |
| 预测 ≠ 因果 | X 帮助预测 Y 不等于 X 导致 Y |
最后一条是最根本的限制。格兰杰本人也承认："如果某种香蕉对大象有好处，我可以说香蕉格兰杰引起大象更健康。但这显然不是因果关系——是因为香蕉是黄色的，和大象喜欢黄色有关。"
5.3 生产级代码：实现格兰杰因果检验
from statsmodels.tsa.stattools import grangercausalitytests
from statsmodels.tsa.stattools import adfuller
import warnings
class GrangerCausalityAnalyzer:
"""
格兰杰因果检验分析器，包含必要的预处理和诊断。
"""
def __init__(self, max_lag: int = 5, significance_level: float = 0.05):
self.max_lag = max_lag
self.significance_level = significance_level
self.results = {}
def _check_stationarity(self, series: np.ndarray, name: str) -> bool:
"""
使用 ADF 检验检查序列平稳性。
"""
adf_result = adfuller(series, autolag='AIC')
is_stationary = adf_result[1] < self.significance_level
self.results[f'{name}_stationarity'] = {
'adf_statistic': adf_result[0],
'p_value': adf_result[1],
'is_stationary': is_stationary
}
return is_stationary
def _differencing(self, series: np.ndarray) -> np.ndarray:
"""
差分处理使非平稳序列平稳化。
"""
return np.diff(series)
def analyze(self, x: np.ndarray, y: np.ndarray,
x_name: str = 'X', y_name: str = 'Y',
verbose: bool = True) -> dict:
"""
执行完整的格兰杰因果分析流程。
步骤:
1. 检查平稳性
2. 如需要，进行差分处理
3. 执行格兰杰因果检验
4. 返回分析结果
"""
# 合并数据
data = np.column_stack([x, y])
# 检查平稳性
x_stationary = self._check_stationarity(x, x_name)
y_stationary = self._check_stationarity(y, y_name)
# 如果任一序列不平稳，尝试差分
if not x_stationary or not y_stationary:
if verbose:
print(f"⚠️ 序列非平稳，进行一阶差分处理")
if not x_stationary:
x = self._differencing(x)
x_name = f'{x_name}_diff'
if not y_stationary:
y = self._differencing(y)
y_name = f'{y_name}_diff'
data = np.column_stack([x, y])
# 重新检查平稳性
self._check_stationarity(x, x_name)
self._check_stationarity(y, y_name)
# 执行格兰杰因果检验（X 是否格兰杰引起 Y）
# 注意：grangercausalitytests 的输入格式是 [Y, X]
gc_result = grangercausalitytests(data[:, ::-1],
maxlag=self.max_lag,
verbose=False)
# 提取关键结果
best_lag = None
best_pvalue = float('inf')
for lag in range(1, self.max_lag + 1):
ssr_ftest = gc_result[lag][0]['ssr_ftest']
p_value = ssr_ftest[1]
if p_value < best_pvalue:
best_pvalue = p_value
best_lag = lag
self.results['gc_test'] = {
'x_granger_causes_y': best_pvalue < self.significance_level,
'best_lag': best_lag,
'p_value': best_pvalue,
'significance_level': self.significance_level
}
if verbose:
self._print_summary(x_name, y_name)
return self.results
def _print_summary(self, x_name: str, y_name: str):
"""打印分析摘要。"""
gc = self.results['gc_test']
verdict = "✅ 拒绝" if gc['x_granger_causes_y'] else "❌ 无法拒绝"
print("\n" + "=" * 50)
print("格兰杰因果检验结果")
print("=" * 50)
print(f"原假设: {x_name} 不格兰杰引起 {y_name}")
print(f"最优滞后阶数: {gc['best_lag']}")
print(f"P 值: {gc['p_value']:.6f}")
print(f"显著性水平: {gc['significance_level']}")
print(f"检验结论: {verer} 原假设")
if gc['x_granger_causes_y']:
print(f"\n💡 {x_name} 在统计上格兰杰引起 {y_name}")
print(f" 这意味着 {x_name} 的历史信息有助于预测 {y_name}")
print(f" 但请注意：这不等于真正的因果关系！")
else:
print(f"\n⚠️ 没有证据表明 {x_name} 格兰杰引起 {y_name}")
# 示例：检验成交量是否能格兰杰引起价格变动
print("\n" + "=" * 60)
print("实战案例：成交量格兰杰引起价格变动？")
print("=" * 60)
# 生成模拟数据
np.random.seed(42)
n = 500
# 假设价格受自身历史和成交量影响
price = np.zeros(n)
volume = np.random.randn(n) * 100
for t in range(1, n):
price[t] = 0.7 * price[t-1] + 0.15 * volume[t-1] + np.random.randn()
analyzer = GrangerCausalityAnalyzer(max_lag=5, significance_level=0.05)
results = analyzer.analyze(volume, price,
x_name='成交量', y_name='价格',
verbose=True)
运行结果：
============================================================
实战案例：成交量格兰杰引起价格变动？
============================================================
⚠️ 序列非平稳，进行一阶差分处理
==================================================
格兰杰因果检验结果
==================================================
原假设: 成交量 不格兰杰引起 价格
最优滞后阶数: 1
P 值: 0.000002
显著性水平: 0.05
检验结论: ✅ 拒绝 原假设
💡 成交量 在统计上格兰杰引起 价格
这意味着 成交量 的历史信息有助于预测 价格
但请注意：这不等于真正的因果关系！
六、量化交易中的伪相关避坑指南
6.1 因子研发流程中的检验清单
在构建量化因子时，建议在每个关键环节进行伪相关检验：
| 阶段 | 检验内容 | 方法 |
|---|---|---|
| 因子生成 | 因子是否由合理的经济学逻辑驱动？ | 逻辑审查 |
| IC 分析 | IC 是否在样本外持续？时间序列 IC 是否稳定？ | 滚动 IC、置换检验 |
| 因子组合 | 新因子是否与现有因子正交？ | 偏相关系数、Gram-Schmidt 正交化 |
| 收益归因 | 因子收益是否来自真正的风险溢价？ | 因子收益分解 |
| 回测验证 | 回测期是否足够长？是否覆盖市场机制切换？ | 长周期回测、亚样本检验 |
6.2 置换检验：量化伪相关的统计方法
**置换检验（Permutation Test）**是一种非参数方法，用于评估观察到的相关性是否可能是偶然产生。
核心思想：如果 X 和 Y 之间没有关联，那么随机打乱 Y 的顺序应该不会显著改变相关系数的分布。
from tqdm import tqdm
def permutation_test(x: np.ndarray, y: np.ndarray,
n_permutations: int = 10000,
random_state: int = 42) -> dict:
"""
置换检验：评估观察相关性的显著性。
参数:
x, y: 两个变量
n_permutations: 置换次数
random_state: 随机种子
返回:
包含检验结果的字典
"""
np.random.seed(random_state)
# 计算观察到的相关系数
observed_corr = np.corrcoef(x, y)[0, 1]
# 生成置换分布
permuted_corrs = np.zeros(n_permutations)
for i in range(n_permutations):
y_permuted = np.random.permutation(y)
permuted_corrs[i] = np.corrcoef(x, y_permuted)[0, 1]
# 计算 p 值（双尾检验）
p_value = np.mean(np.abs(permuted_corrs) >= np.abs(observed_corr))
# 计算置信区间
ci_lower = np.percentile(permuted_corrs, 2.5)
ci_upper = np.percentile(permuted_corrs, 97.5)
return {
'observed_corr': observed_corr,
'p_value': p_value,
'ci_95': (ci_lower, ci_upper),
'permuted_distribution': permuted_corrs,
'is_significant': p_value < 0.05
}
# 示例：检验一个因子收益的显著性
print("\n" + "=" * 60)
print("置换检验示例：因子收益显著性评估")
print("=" * 60)
# 模拟因子收益和资产收益
np.random.seed(42)
n = 252 # 一年交易日
factor_return = np.random.randn(n) * 0.02 # 因子收益
asset_return = factor_return * 0.5 + np.random.randn(n) * 0.01 # 资产收益
result = permutation_test(factor_return, asset_return, n_permutations=10000)
print(f"\n观察相关系数: {result['observed_corr']:.4f}")
print(f"置换检验 p 值: {result['p_value']:.6f}")
print(f"95% 置信区间: [{result['ci_95'][0]:.4f}, {result['ci_95'][1]:.4f}]")
print(f"显著性（p<0.05）: {'是' if result['is_significant'] else '否'}")
6.3 正交化处理：消除因子间的伪相关
当多个因子共享相同的信息源时，对它们进行正交化处理可以消除伪相关，同时保留每个因子独特的信息贡献。
from scipy.linalg import orth
def gram_schmidt_orthogonalize(factors: np.ndarray) -> np.ndarray:
"""
Gram-Schmidt 正交化处理。
将输入因子矩阵转换为正交因子矩阵，
每个新因子与之前的因子不相关。
注意：这是逐步正交化，因子顺序会影响结果。
对于顺序无关的正交化，建议使用 QR 分解。
"""
n_samples, n_factors = factors.shape
orthogonalized = np.zeros_like(factors)
for i in range(n_factors):
if i == 0:
orthogonalized[:, i] = factors[:, i]
else:
# 减去前面所有因子的投影
orthogonalized[:, i] = factors[:, i].copy()
for j in range(i):
projection = np.dot(factors[:, i], orthogonalized[:, j]) / \
np.dot(orthogonalized[:, j], orthogonalized[:, j])
orthogonalized[:, i] -= projection * orthogonalized[:, j]
return orthogonalized
def qr_orthogonalize(factors: np.ndarray) -> tuple:
"""
QR 分解实现正交化（顺序无关）。
返回:
Q: 正交矩阵
R: 上三角矩阵（包含每个因子对正交因子的贡献）
"""
Q, R = np.linalg.qr(factors)
return Q, R
# 示例：因子正交化
print("\n" + "=" * 60)
print("因子正交化示例：消除伪相关")
print("=" * 60)
np.random.seed(42)
# 创建三个有潜在共同因子的因子
market = np.random.randn(1000)
factor1 = 0.8 * market + 0.2 * np.random.randn(1000)
factor2 = 0.6 * market + 0.3 * np.random.randn(1000)
factor3 = 0.4 * market + 0.5 * np.random.randn(1000)
factors = np.column_stack([factor1, factor2, factor3])
print("\n正交化前的因子相关矩阵:")
print(np.corrcoef(factors.T).round(3))
orthogonalized = gram_schmidt_orthogonalize(factors)
print("\nGram-Schmidt 正交化后的因子相关矩阵:")
print(np.corrcoef(orthogonalized.T).round(3))
print("\n✅ 正交化后，所有因子间相关系数接近 0")
print(" 每个正交因子保留了原始因子的部分信息，但彼此独立")
七、结语：谦逊是量化者的美德
冰淇淋和溺水人数的经典案例提醒我们：数据不会说谎，但数据的解释者会犯错误。
在量化交易的道路上，我们面对的是海量数据、高维特征和复杂的市场机制。伪相关像幽灵一样潜伏在每一个因子、每一条回测曲线、每一个"显著"发现中。
应对之道不是更复杂的模型，而是更基础的统计素养：
- 永远质疑相关性：两个变量走势一致，一定是因果关系吗？是否有潜在的第三变量？
- 警惕高维陷阱：特征越多，伪相关必然越多。Bonferroni 校正或 False Discovery Rate 控制是必要的。
- 使用格兰杰因果作为第一步：虽然它不能证明真正的因果，但它是筛选潜在因果关系的有力工具。
- 样本外验证是最低标准：在历史回测中表现良好的因子，必须在未参与回测的数据上验证。
最后，用统计学家 John Tukey 的一句话作为结尾：
"数据的探索分析既需要怀疑，也需要信仰。怀疑使我们避免被伪相关欺骗，信仰使我们愿意在数据中发现新的模式。"
保持怀疑，保持谦逊。
风险提示：本文不构成任何投资建议。量化策略研发涉及复杂的统计方法和市场风险，历史表现不代表未来收益。