因果推断
预计学习时间:2-3 小时
难度:⭐⭐⭐⭐
核心问题:“相关不等于因果”——但量化中我们经常需要回答因果问题。如何从观测数据中识别因果效应?
概述
你看到”冰激凌销量高的时候,溺水事故多”。你会说”吃冰激凌导致溺水”吗?当然不会。真正的原因是:夏天到了,天热 → 冰激凌销量增加,同时游泳的人增加 → 溺水事故增加。
这就是**混杂因素(Confounder)**的问题。在很多场景中,混杂因素不像”天气”这么明显,你需要严谨的方法来识别真正的因果关系。
量化中,因果推断的典型问题包括:
- 价值因子的高收益,是因为”承担了风险”,还是”行为偏差”?
- 交易所降费后,流动性改善是因为降费本身,还是市场环境变化?
- 被纳入指数后,股票流动性提升是因为纳入指数,还是因为纳入时市场刚好在涨?
一、相关 vs 因果
1.1 经典例子
| 现象 | 相关 | 因果? | 真正原因 |
|---|---|---|---|
| 冰激凌销量与溺水事故 | 高 | 不是 | 天气(混杂因素) |
| 小学成绩和鞋子尺码 | 正相关 | 不是 | 年龄(混杂因素) |
| 医院越多,死亡率越高 | 正相关 | 不是 | 严重疾病集中在大医院 |
| 教育水平和收入 | 正相关 | 大概率是 | 但能力可能也是混杂因素 |
1.2 量化中的相关 vs 因果
| 量化现象 | 你看到的相关 | 可能的混杂因素 |
|---|---|---|
| 高 BM 股票收益高 | 价值因子有效 | 杠杆约束、机构偏好 |
| 低波动股票长期收益不差 | 低波动异象 | 杠杆限制、做空限制 |
| 加入指数后价格上涨 | 指数效应 | 信息不对称降低、流动性提升 |
| 分析师上调评级后价格上涨 | 评级效应 | 信息泄露、选择性覆盖 |
核心教训:在量化研究中,如果你只报告”因子 X 和收益 Y 相关”,而不考虑混杂因素,你的结论可能是错的。
二、潜在结果框架(Rubin Causal Model)
2.1 核心思想
白话解释:对于每个个体,存在两个”潜在结果”——接受处理时的结果和不接受处理时的结果。因果效应就是两者的差。
对于个体 :
- :如果接受处理的潜在结果
- :如果不接受处理的潜在结果
个体因果效应:
问题:现实中,你只能观测到一个结果(个体要么接受了处理,要么没有)。你永远无法同时观测到 和 。这就是因果推断的根本难题。
2.2 平均处理效应(ATE)
既然个体效应不可观测,我们退而求其次,估计群体的平均处理效应:
如果处理是随机分配的(如随机对照试验 RCT),ATE 可以直接计算:
但现实中的大多数数据不是随机分配的。这就是因果推断方法要解决的问题。
2.3 选择偏差
白话解释:选择偏差说的是”接受处理的个体和未接受处理的个体本身就有差异”。
例子:你想研究”上大学对收入的影响”。但上大学的和不上大学的人,本身在能力、家庭背景上就不同。你看到的收入差异,可能不是大学教育的效果,而是本来就有的能力差异。
2.4 Python 示例:模拟选择偏差
import numpy as np
import pandas as pd
np.random.seed(42)
n = 1000
# ============================================================
# 模拟数据:教育的"因果效应"
# ============================================================
# 能力(混杂因素,影响是否接受教育和收入)
ability = np.random.normal(100, 15, n)
# 是否上大学(和能力相关 → 选择偏差!)
# 能力高的人更容易上大学
prob_college = 1 / (1 + np.exp(-(ability - 100) / 10))
college = np.random.binomial(1, prob_college)
# 收入 = 能力效应 + 教育效应 + 噪声
# 真实的教育因果效应 = 5000 元
true_effect = 5000
income = (200 * ability + true_effect * college +
np.random.normal(0, 10000, n))
# ============================================================
# 朴素估计(忽略选择偏差)
# ============================================================
naive_ate = np.mean(income[college == 1]) - np.mean(income[college == 0])
# ============================================================
# 控制混杂因素后的估计
# ============================================================
import statsmodels.api as sm
X = sm.add_constant(np.column_stack([college, ability]))
model = sm.OLS(income, X).fit()
controlled_ate = model.params[1]
print("=" * 50)
print("选择偏差演示")
print("=" * 50)
print(f" 真实因果效应: {true_effect}")
print(f" 朴素估计(有偏): {naive_ate:.0f}")
print(f" 控制能力后: {controlled_ate:.0f}")
print(f"\n 选择偏差 = {naive_ate - true_effect:.0f}")
print(f" 能力高的人更可能上大学,也更容易高收入")
print(f" 所以朴素估计把'能力效应'也算进了'教育效应'")三、断点回归(RDD)
3.1 核心思想
白话解释:在某个阈值附近,处理分配”就像随机分组”一样。
例子:假设政策规定”市值 < 100 亿的股票可以享受税收优惠”。你比较市值 99 亿和 101 亿的股票。这两组股票在各方面都差不多(除了一个刚好在阈值下面,一个刚好在上面),唯一的区别就是税收待遇。
其中 是断点(阈值), 是驱动变量(running variable)。
为什么有效? 在阈值附近,处理组和非处理组的个体几乎一样——差异只来源于处理本身。
3.2 Python 示例
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
np.random.seed(42)
n = 1000
# ============================================================
# 模拟 RDD 数据
# ============================================================
# 驱动变量:股票市值(10 亿为单位),断点在 100 亿
market_cap = np.random.uniform(80, 120, n)
# 处理变量:市值 < 100 亿 → 税收优惠
treatment = (market_cap < 100).astype(int)
# 结果变量:流动性(买卖价差的倒数)
# 真实因果效应:税收优惠提升流动性 0.5 个单位
base_liquidity = 5 + 0.02 * market_cap + np.random.normal(0, 0.5, n)
true_rdd_effect = 0.5
liquidity = base_liquidity + true_rdd_effect * treatment
df = pd.DataFrame({
'market_cap': market_cap,
'treatment': treatment,
'liquidity': liquidity
})
# ============================================================
# RDD 估计:只使用断点附近的观测值
# ============================================================
bandwidth = 5 # 只看断点 ±5 的范围
df_narrow = df[(df['market_cap'] >= 95) & (df['market_cap'] <= 105)]
# 方法 1:简单均值差(在带宽内)
simple_rdd = (df_narrow.loc[df_narrow['treatment'] == 1, 'liquidity'].mean() -
df_narrow.loc[df_narrow['treatment'] == 0, 'liquidity'].mean())
# 方法 2:局部线性回归(更精确)
# 分别在断点左右做回归
left = df_narrow[df_narrow['market_cap'] < 100]
right = df_narrow[df_narrow['market_cap'] >= 100]
# 左侧:Y = a + b*(X - 100)
X_left = sm.add_constant(left['market_cap'] - 100)
model_left = sm.OLS(left['liquidity'], X_left).fit()
# 右侧:Y = c + d*(X - 100)
X_right = sm.add_constant(right['market_cap'] - 100)
model_right = sm.OLS(right['liquidity'], X_right).fit()
# RDD 效应 = 左侧在断点处的拟合值 - 右侧在断点处的拟合值
local_rdd = model_left.params[0] - model_right.params[0]
print("=" * 50)
print("断点回归(RDD)估计")
print("=" * 50)
print(f" 断点: 市值 = 100 亿")
print(f" 带宽: ±{bandwidth} 亿")
print(f" 带宽内样本量: {len(df_narrow)}")
print(f"\n 真实因果效应: {true_rdd_effect}")
print(f" 简单均值差: {simple_rdd:.4f}")
print(f" 局部线性回归: {local_rdd:.4f}")
print(f"\n 结论: 税收优惠显著提升了小盘股流动性")
# ============================================================
# 可视化
# ============================================================
fig, ax = plt.subplots(figsize=(12, 6))
# 散点图
ax.scatter(df['market_cap'], df['liquidity'], alpha=0.3, s=10,
c=df['treatment'], cmap='coolwarm')
# 拟合线
x_left_fit = np.linspace(80, 100, 50)
x_right_fit = np.linspace(100, 120, 50)
y_left_fit = model_left.params[0] + model_left.params[1] * (x_left_fit - 100)
y_right_fit = model_right.params[0] + model_right.params[1] * (x_right_fit - 100)
ax.plot(x_left_fit, y_left_fit, 'r-', linewidth=2, label='断点左侧拟合')
ax.plot(x_right_fit, y_right_fit, 'b-', linewidth=2, label='断点右侧拟合')
ax.axvline(x=100, color='gray', linestyle='--', label='断点')
ax.set_xlabel('市值(亿元)')
ax.set_ylabel('流动性')
ax.set_title('断点回归(RDD)')
ax.legend()
plt.tight_layout()
plt.show()四、合成控制法
4.1 核心思想
白话解释:为处理组”合成”一个对照组——用多个未处理的个体加权组合,模拟处理组”如果没被处理”会怎样。
适用场景:只有一个(或很少)处理个体,无法用 DID。
例子:2015 年上海证券交易所实施了某项创新政策。你想评估这个政策的效果。但只有一个上交所,没有”另一个上交所”作为对照组。
解决方法:用深交所、中金所等”类似但未处理”的市场,按某种权重组合,“合成”一个”如果没有政策的话上交所应该长什么样”。
4.2 Python 示例
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import minimize
np.random.seed(42)
# ============================================================
# 模拟数据:一个"处理单位"和多个"候选对照单位"
# ============================================================
T_pre = 12 # 政策前 12 个月
T_post = 12 # 政策后 12 个月
T = T_pre + T_post
# 处理单位(比如上交所)
# 政策前和对照组有相似的轨迹
common_trend = np.cumsum(np.random.normal(0, 1, T))
treated = 2 * common_trend + np.random.normal(0, 0.5, T)
# 政策后有一个正向冲击
treated[T_pre:] += 5 # 政策效果
# 5 个候选对照单位
n_controls = 5
controls = np.zeros((n_controls, T))
for i in range(n_controls):
weight_true = np.random.uniform(-0.5, 0.5)
controls[i] = weight_true * common_trend + np.random.normal(0, 0.5 + 0.2*i, T)
# ============================================================
# 合成控制法:找到最优权重
# ============================================================
# 目标:在政策前,让加权和尽可能接近处理单位
def mse_loss(weights):
"""政策前的均方误差"""
synthetic = weights @ controls[:, :T_pre]
return np.mean((treated[:T_pre] - synthetic)**2)
# 约束:权重非负,权重之和为 1
constraints = {'type': 'eq', 'fun': lambda w: np.sum(w) - 1}
bounds = [(0, 1)] * n_controls
initial_weights = np.array([1/n_controls] * n_controls)
result = minimize(mse_loss, initial_weights,
method='SLSQP',
bounds=bounds,
constraints=constraints)
optimal_weights = result.x
synthetic = optimal_weights @ controls
# ============================================================
# 结果
# ============================================================
print("=" * 50)
print("合成控制法")
print("=" * 50)
print(f"\n最优权重:")
for i, w in enumerate(optimal_weights):
if w > 0.01:
print(f" 对照单位 {i+1}: {w:.4f}")
# 政策前的拟合优度
pre_mse = np.mean((treated[:T_pre] - synthetic[:T_pre])**2)
print(f"\n政策前 MSE: {pre_mse:.4f}")
# 政策效果估计
effect = np.mean(treated[T_pre:] - synthetic[T_pre:])
print(f"政策后平均效果: {effect:.4f}")
# ============================================================
# 可视化
# ============================================================
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左图:处理组 vs 合成对照组
months = np.arange(T)
axes[0].plot(months, treated, 'b-o', label='处理单位', markersize=3)
axes[0].plot(months, synthetic, 'r--', label='合成对照组', linewidth=2)
axes[0].axvline(x=T_pre - 0.5, color='gray', linestyle=':', label='政策实施')
axes[0].set_xlabel('时间')
axes[0].set_ylabel('结果变量')
axes[0].set_title('处理组 vs 合成对照组')
axes[0].legend()
# 右图:处理效应(差异)
diff = treated - synthetic
axes[1].plot(months, diff, 'g-o', markersize=4)
axes[1].axvline(x=T_pre - 0.5, color='gray', linestyle=':', label='政策实施')
axes[1].axhline(y=0, color='black', linestyle='-', linewidth=0.5)
axes[1].set_xlabel('时间')
axes[1].set_ylabel('处理效应(差异)')
axes[1].set_title('政策效应')
axes[1].legend()
plt.tight_layout()
plt.show()五、事件研究法(Event Study)
5.1 核心思想
白话解释:围绕一个事件(如财报发布、政策变化、评级调整),看股票的”异常收益”是否显著偏离零。
关键概念:
- 正常收益:根据市场模型预测的收益(“如果没有事件,股票应该赚多少”)
- 异常收益(AR):实际收益 - 正常收益
- 累积异常收益(CAR):事件窗口内异常收益的累加
如果 CAR 显著不为零,说明事件有因果效应。
5.2 事件窗口选择
估计窗口 事件窗口
|<---------->| |<------>|<-------->|
[-250, -30] [-10, 0] [1, 10]
事件前 事件后
- 估计窗口:用来估计”正常收益”模型(通常是市场模型)
- 事件窗口:用来计算异常收益
- 避免在估计窗口中包含事件本身
5.3 量化中的应用场景
| 事件类型 | 典型事件 | 因果问题 |
|---|---|---|
| 公司事件 | 财报、分红、回购、并购 | 事件是否创造了价值? |
| 政策事件 | 监管政策、税收变化、准入限制 | 政策是否影响了市场? |
| 评级事件 | 分析师评级、信用评级变化 | 评级是否包含了新信息? |
| 指数事件 | 纳入/剔除指数 | 指数效应是否存在? |
5.4 Python 完整代码
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
np.random.seed(42)
# ============================================================
# 模拟事件研究数据
# ============================================================
n_stocks = 50 # 50 只发生了事件的股票
T_est = 200 # 估计窗口 200 天
T_event = 21 # 事件窗口 [-10, +10]
T_total = T_est + T_event
results = []
for i in range(n_stocks):
# 市场收益
market = np.random.normal(0.001, 0.015, T_total)
# 个股 beta(对市场的敏感度)
beta_i = np.random.uniform(0.5, 1.5)
alpha_i = np.random.normal(0, 0.001)
# 个股正常收益 = alpha + beta * market + noise
normal_return = alpha_i + beta_i * market + np.random.normal(0, 0.02, T_total)
# 在事件日添加异常收益(模拟事件效应)
event_day = T_est # 事件发生在估计窗口之后
# 事件效应:从事件前一天到事件后 3 天,累计异常收益 = 3%
event_effect = np.zeros(T_total)
event_effect[event_day - 1] = 0.01
event_effect[event_day] = 0.015
event_effect[event_day + 1] = 0.005
actual_return = normal_return + event_effect
# 估计窗口:估计市场模型
est_X = sm.add_constant(market[:T_est])
est_model = sm.OLS(actual_return[:T_est], est_X).fit()
# 事件窗口:计算异常收益
event_X = sm.add_constant(market[T_est:])
event_predicted = est_model.predict(event_X)
abnormal_return = actual_return[T_est:] - event_predicted
# 累积异常收益
cumulative_return = np.cumsum(abnormal_return)
results.append({
'stock': f'S{i:03d}',
'alpha': est_model.params[0],
'beta': est_model.params[1],
'abnormal_returns': abnormal_return,
'cumulative_returns': cumulative_return,
})
# ============================================================
# 汇总分析
# ============================================================
# 计算平均异常收益(AAR)和平均累积异常收益(CAAR)
ar_matrix = np.array([r['abnormal_returns'] for r in results]) # (n_stocks, T_event)
car_matrix = np.array([r['cumulative_returns'] for r in results])
aar = np.mean(ar_matrix, axis=0) # 每天的平均异常收益
caar = np.mean(car_matrix, axis=0) # 每天的平均累积异常收益
# t 检验
ar_std = np.std(ar_matrix, axis=0, ddof=1)
ar_t = aar / (ar_std / np.sqrt(n_stocks))
car_std = np.std(car_matrix, axis=0, ddof=1)
car_t = caar / (car_std / np.sqrt(n_stocks))
event_days = np.arange(-10, 11)
print("=" * 60)
print("事件研究结果")
print("=" * 60)
print(f"\n平均异常收益(AAR)和 t 检验:")
print(f"{'事件日':>6} {'AAR':>10} {'t值':>8} {'显著':>6}")
print("-" * 35)
for d in [0, -1, 1, 2, 3, -2, 5]:
idx = d + 10
sig = "***" if abs(ar_t[idx]) > 2.58 else "**" if abs(ar_t[idx]) > 1.96 else "*"
print(f"{d:>6} {aar[idx]:>10.4f} {ar_t[idx]:>8.2f} {sig:>6}")
# 事件窗口 CAAR
print(f"\n累积异常收益(CAAR):")
print(f"{'窗口':>12} {'CAAR':>10} {'t值':>8}")
print("-" * 35)
for window_name, start, end in [('[0,0]', 10, 10),
('[-1,0]', 9, 10),
('[-1,+1]', 9, 11),
('[0,+3]', 10, 13),
('[-5,+5]', 5, 15)]:
window_caar = np.mean(car_matrix[:, start:end+1], axis=1)[:, -1]
window_t = np.mean(window_caar) / (np.std(window_caar, ddof=1) / np.sqrt(n_stocks))
print(f"{window_name:>12} {np.mean(window_caar):>10.4f} {window_t:>8.2f}")
# ============================================================
# 可视化
# ============================================================
fig, axes = plt.subplots(2, 1, figsize=(12, 8))
# 上图:CAAR
axes[0].plot(event_days, caar, 'b-o', markersize=4, linewidth=1.5)
axes[0].axhline(y=0, color='black', linestyle='-', linewidth=0.5)
axes[0].axvline(x=0, color='red', linestyle='--', alpha=0.5, label='事件日')
# 置信区间
axes[0].fill_between(event_days,
caar - 1.96 * car_std / np.sqrt(n_stocks),
caar + 1.96 * car_std / np.sqrt(n_stocks),
alpha=0.2, color='blue')
axes[0].set_xlabel('事件日')
axes[0].set_ylabel('CAAR')
axes[0].set_title('平均累积异常收益(CAAR)')
axes[0].legend()
# 下图:AAR 柱状图
colors = ['red' if abs(t) > 1.96 else 'gray' for t in ar_t]
axes[1].bar(event_days, aar, color=colors, alpha=0.7)
axes[1].axhline(y=0, color='black', linestyle='-', linewidth=0.5)
axes[1].axvline(x=0, color='red', linestyle='--', alpha=0.5)
axes[1].set_xlabel('事件日')
axes[1].set_ylabel('AAR')
axes[1].set_title('平均异常收益(AAR),红色柱表示 t > 1.96')
plt.tight_layout()
plt.show()5.5 事件研究的注意事项
- 估计窗口要足够长:至少 100-250 个交易日,确保市场模型估计稳定
- 事件窗口不宜过长:过长会混入其他事件的影响
- 注意事件重叠:如果多只股票的事件日相近,它们的异常收益可能相关
- 考虑事件泄露:如果事件信息在正式公布前就泄露了,CAR 会在事件日前就开始累积
- 调整市场模型:可以用多因子模型代替简单的市场模型
六、Granger 因果 vs 真实因果
6.1 关键区别
Granger 因果检验是时间序列分析中常用的工具,但它不等于真正的因果关系。
Granger 因果的定义:如果用 X 的过去值能帮助更好地预测 Y(相对于只用 Y 的过去值),就说 “X Granger causes Y”。
数学上:
如果 联合显著(F 检验),则 X Granger causes Y。
注意:
- “Granger causes” 只是”在统计上有预测力”
- 可能存在第三变量 Z 同时影响 X 和 Y
- 可能是 Y → X 而不是 X → Y(方向可能搞反)
- Granger 因果只检验时间上的先后关系,不是真正的因果
6.2 Python 示例
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import grangercausalitytests
np.random.seed(42)
n = 500
# ============================================================
# 情况 1:X 真正 Granger causes Y
# ============================================================
x1 = np.random.normal(0, 1, n)
y1 = np.zeros(n)
for t in range(2, n):
y1[t] = 0.5 * y1[t-1] + 0.3 * x1[t-1] + np.random.normal(0, 1)
data1 = pd.DataFrame({'y': y1, 'x': x1})
print("=" * 50)
print("Granger 因果检验")
print("=" * 50)
print(f"\n情况 1:X 真正 Granger causes Y")
print(f"检验: X → Y")
grangercausalitytests(data1[['y', 'x']], maxlag=4)
print(f"\n检验: Y → X(应该不显著)")
grangercausalitytests(data1[['x', 'y']], maxlag=4)
# ============================================================
# 情况 2:X 和 Y 有共同驱动因子 Z(虚假 Granger 因果)
# ============================================================
z = np.random.normal(0, 1, n)
x2 = 0.5 * z + np.random.normal(0, 0.5, n)
y2 = 0.8 * z + np.random.normal(0, 0.5, n) # Y 不依赖 X 的过去
data2 = pd.DataFrame({'y': y2, 'x': x2})
print(f"\n情况 2:X 和 Y 有共同驱动因子 Z(虚假 Granger 因果)")
print(f"检验: X → Y(可能'显著',但这是虚假的!)")
grangercausalitytests(data2[['y', 'x']], maxlag=4)6.3 量化中的误用
Granger 因果检验在量化中经常被误用:
-
误用:发现利率 Granger causes 股市收益,就认为”利率决定股市”
- 可能是宏观因素同时驱动两者
-
误用:用 Granger 因果来判断因子是否有”因果效应”
- 因子只是预测变量,不是因果变量
-
误用:忽略滞后阶数的选择
- 滞后阶数不同,结论可能完全不同
正确理解:Granger 因果只说明”X 的过去信息对预测 Y 有帮助”,不要过度解读为因果关系。
七、因果推断在量化中的应用场景
7.1 因子有效性检验
问题:价值因子的高收益,是”风险补偿”还是”数据挖掘”?
方法:
- 用 Fama-MacBeth + 多种稳健性检验
- 控制已知因子后看增量解释力
- 样本外验证 + 子样本分析
7.2 策略归因
问题:策略的超额收益来自 Alpha 还是 Beta?
方法:
- 多因子回归归因
- 控制市场、风格因子后看残差 Alpha
- 事件研究法评估特定事件贡献
7.3 政策评估
问题:某监管政策对市场质量的影响?
方法:
- DID:对比受影响和不受影响的资产
- RDD:利用政策阈值附近
- 合成控制法:构造反事实
7.4 执行质量评估
问题:算法交易是否真的改善了执行质量?
方法:
- DID:比较使用算法前后的执行成本
- 控制订单特征和市场条件
本文件要点回顾
| 概念 | 一句话总结 |
|---|---|
| 相关 vs 因果 | 相关不等于因果——混杂因素可能同时影响两者 |
| 潜在结果框架 | 因果 = 接受处理的结果 - 不接受处理的结果(但永远看不到后者) |
| 选择偏差 | 接受处理的个体本身就和未接受的不同 |
| RDD | 在阈值附近”就像随机分组” |
| 合成控制法 | 为处理组构造一个”合成对照组” |
| 事件研究法 | 围绕事件看累积异常收益 |
| Granger 因果 | ”X 的过去能帮助预测 Y”——不等于真正的因果 |
| CAR | 累积异常收益,事件研究的核心指标 |
整个计量经济学模块总结
你现在掌握了传统量化研究所需的核心计量工具:
| 文件 | 核心能力 |
|---|---|
| 经典回归 | 知道回归结果何时可信、何时不可信 |
| 时间序列 | 知道如何正确分析时间序列、避免虚假回归 |
| 面板数据 | 知道如何利用截面-时间二维结构 |
| 因果推断 | 知道”相关不等于因果”,掌握识别因果的方法 |
最终的检验标准:你是否能在面对一个研究问题时,自动判断:
- 该用什么方法(截面、时间序列、面板)
- 可能有哪些统计陷阱(异方差、自相关、虚假回归、选择偏差)
- 该怎么修正(HAC、Fama-MacBeth、IV、DID)
- 结果是否稳健(不同设定下是否一致)
这就是计量经济学给你的核心能力——不是跑更多的回归,而是知道怎样跑才是对的。