8.5 经典案例与 Python 实现
"The best way to learn econometrics is to do econometrics.""学习计量经济学最好的方式就是动手做计量经济学。"— Clive Granger, Nobel Laureate
从经典论文到可复现的 Python 代码
📌 本节概要
在本节中，我们将：
- 案例 1：复现 Angrist & Krueger (1991) — 出生季度 → 教育 → 工资
- 案例 2：复现 Card (1995) — 距离学院 → 教育 → 工资
- 使用模拟数据还原论文核心结构
- 完整的 Python 实现：从数据生成到结果解读
- 诊断检验和稳健性分析
📚 案例 1：Angrist & Krueger (1991)
📖 论文背景
论文："Does Compulsory School Attendance Affect Schooling and Earnings?"发表：Quarterly Journal of Economics, 1991
核心问题：义务教育法是否真的增加了教育年限，并进而提高了工资？
内生性问题：
- 教育与工资的 OLS 回归受到"能力偏误"的困扰
- 高能力者同时倾向于接受更多教育和获得更高工资
工具变量：出生季度（Quarter of Birth, QOB）
制度逻辑：
义务教育法规定：年满 16 岁才能辍学
入学规则：在当年某日期前满 6 岁 → 进入一年级
结果：
- 1 月出生 → 满 6 岁时进入一年级 → 16 岁时在 10 年级
- 12 月出生 → 满 6 岁时进入一年级 → 16 岁时在 9 年级末
→ 第一季度出生的人，平均受教育年限略低于第四季度出生的人🔢 模拟数据生成
我们创建一个结构上模拟 Angrist & Krueger (1991) 核心发现的数据集。
python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.formula.api as smf
from linearmodels.iv import IV2SLS
from scipy import stats
# 设置
plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
sns.set_style("whitegrid")
np.random.seed(1991) # 致敬论文年份
n = 50000 # 大样本（模仿 Census 数据）
# --- 数据生成过程 ---
# 出生季度（1, 2, 3, 4）
quarter_of_birth = np.random.choice([1, 2, 3, 4], size=n)
# 出生年份（1930-1949，模仿原论文的队列）
birth_year = np.random.choice(range(1930, 1950), size=n)
# 不可观测的能力
ability = np.random.normal(0, 1, n)
# 教育年限：受出生季度和能力影响
# 核心机制：Q1 出生的人受教育年限略低
qob_effect = np.where(quarter_of_birth == 1, -0.15,
np.where(quarter_of_birth == 2, -0.05,
np.where(quarter_of_birth == 3, 0.05, 0.15)))
education = 12 + qob_effect + 1.5 * ability + np.random.normal(0, 1.5, n)
education = np.clip(education, 6, 20) # 限制在合理范围
# 工资（对数周工资）：受教育和能力影响
# 真实教育回报 = 0.08（每增加一年教育，工资增加约 8%）
beta_true = 0.08
log_wage = 5.5 + beta_true * education + 0.15 * ability + \
0.005 * (birth_year - 1940) + np.random.normal(0, 0.4, n)
# 创建虚拟变量
df_ak = pd.DataFrame({
'log_wage': log_wage,
'education': education,
'quarter_of_birth': quarter_of_birth,
'birth_year': birth_year,
'ability': ability,
'Q1': (quarter_of_birth == 1).astype(int),
'Q2': (quarter_of_birth == 2).astype(int),
'Q3': (quarter_of_birth == 3).astype(int),
'Q4': (quarter_of_birth == 4).astype(int),
'age': 1991 - birth_year
})
# 出生年份虚拟变量
for year in range(1930, 1950):
df_ak[f'yr_{year}'] = (df_ak['birth_year'] == year).astype(int)
print("=" * 70)
print("案例 1：Angrist & Krueger (1991) — 出生季度 IV")
print("=" * 70)
print(f"样本量: {n:,}")
print(f"真实教育回报: {beta_true} (每年教育增加 {beta_true*100:.0f}% 工资)")
print(f"\n按出生季度的平均教育年限:")
print(df_ak.groupby('quarter_of_birth')['education'].mean().round(3))
print(f"\n按出生季度的平均对数工资:")
print(df_ak.groupby('quarter_of_birth')['log_wage'].mean().round(4))📊 关键图表：出生季度与教育年限
python
# 图 1：Angrist & Krueger 的经典图——按出生年份-季度的平均教育年限
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
# 按出生季度的平均教育年限
qob_educ = df_ak.groupby('quarter_of_birth')['education'].mean()
colors = ['#e74c3c', '#3498db', '#2ecc71', '#f39c12']
axes[0].bar(qob_educ.index, qob_educ.values, color=colors,
edgecolor='black', width=0.6)
axes[0].set_xlabel('出生季度', fontsize=13)
axes[0].set_ylabel('平均教育年限', fontsize=13)
axes[0].set_title('出生季度与教育年限\n（Q1 出生者教育年限最低）', fontsize=14, fontweight='bold')
axes[0].set_xticks([1, 2, 3, 4])
axes[0].set_xticklabels(['Q1\n(1-3月)', 'Q2\n(4-6月)', 'Q3\n(7-9月)', 'Q4\n(10-12月)'])
# 标注差异
q1_mean = qob_educ[1]
q4_mean = qob_educ[4]
axes[0].annotate(f'Q4-Q1 差异\n= {q4_mean - q1_mean:.3f} 年',
xy=(2.5, (q1_mean + q4_mean) / 2),
fontsize=12, fontweight='bold', color='purple',
bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.3))
# 按出生季度的平均对数工资
qob_wage = df_ak.groupby('quarter_of_birth')['log_wage'].mean()
axes[1].bar(qob_wage.index, qob_wage.values, color=colors,
edgecolor='black', width=0.6)
axes[1].set_xlabel('出生季度', fontsize=13)
axes[1].set_ylabel('平均对数周工资', fontsize=13)
axes[1].set_title('出生季度与对数工资\n（Q1 出生者工资最低）', fontsize=14, fontweight='bold')
axes[1].set_xticks([1, 2, 3, 4])
axes[1].set_xticklabels(['Q1\n(1-3月)', 'Q2\n(4-6月)', 'Q3\n(7-9月)', 'Q4\n(10-12月)'])
plt.suptitle('Angrist & Krueger (1991) 核心发现的再现', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('ak1991_qob_patterns.png', dpi=300, bbox_inches='tight')
plt.show()📈 OLS 估计
python
# --- OLS 回归 ---
# 模型 1：简单回归
ols1 = smf.ols('log_wage ~ education', data=df_ak).fit()
# 模型 2：加入出生年份控制
year_controls = ' + '.join([f'yr_{y}' for y in range(1931, 1950)])
ols2 = smf.ols(f'log_wage ~ education + {year_controls}', data=df_ak).fit()
print("\n" + "=" * 70)
print("OLS 估计结果")
print("=" * 70)
print(f"\n模型 1（简单 OLS）：")
print(f" 教育回报: {ols1.params['education']:.5f} ({ols1.params['education']*100:.2f}%)")
print(f" 标准误: {ols1.bse['education']:.5f}")
print(f" 95% CI: [{ols1.conf_int().loc['education', 0]:.5f}, {ols1.conf_int().loc['education', 1]:.5f}]")
print(f"\n模型 2（控制出生年份）：")
print(f" 教育回报: {ols2.params['education']:.5f} ({ols2.params['education']*100:.2f}%)")
print(f" 标准误: {ols2.bse['education']:.5f}")
print(f"\n真实效应: {beta_true} ({beta_true*100:.0f}%)")
print(f"OLS 偏差: {ols1.params['education'] - beta_true:.5f}（因为能力混淆，OLS 高估）")🔧 IV/2SLS 估计
python
# --- 第一阶段：education ~ Q1 + Q2 + Q3 (Q4 为参照组) ---
first_stage = smf.ols('education ~ Q1 + Q2 + Q3', data=df_ak).fit()
print("\n" + "=" * 70)
print("第一阶段回归")
print("=" * 70)
print(f"{'变量':<8} {'系数':<10} {'标准误':<10} {'t 统计量':<10} {'p 值':<10}")
print("-" * 48)
for var in ['Q1', 'Q2', 'Q3']:
print(f"{var:<8} {first_stage.params[var]:<10.4f} {first_stage.bse[var]:<10.4f} "
f"{first_stage.tvalues[var]:<10.2f} {first_stage.pvalues[var]:<10.6f}")
print(f"\n第一阶段 F 统计量: {first_stage.fvalue:.2f}")
print(f"第一阶段 R²: {first_stage.rsquared:.6f}")
if first_stage.fvalue > 10:
print(f"✅ F = {first_stage.fvalue:.2f} > 10，工具变量足够强")
else:
print(f"⚠️ F = {first_stage.fvalue:.2f} < 10，弱工具变量警告！")
# --- IV/2SLS 估计 ---
# 方法 A：使用 Q1 作为单一 IV
iv_single = IV2SLS.from_formula('log_wage ~ 1 [education ~ Q1]', data=df_ak).fit(cov_type='robust')
# 方法 B：使用 Q1, Q2, Q3 作为多个 IV
iv_multi = IV2SLS.from_formula('log_wage ~ 1 [education ~ Q1 + Q2 + Q3]', data=df_ak).fit(cov_type='robust')
# 方法 C：控制出生年份
iv_controls = IV2SLS.from_formula(
f'log_wage ~ 1 + {year_controls} [education ~ Q1 + Q2 + Q3]', data=df_ak
).fit(cov_type='robust')
print("\n" + "=" * 70)
print("IV/2SLS 估计结果")
print("=" * 70)
print(f"\n{'模型':<30} {'教育回报':<12} {'标准误':<12} {'95% CI':<25}")
print("-" * 79)
print(f"{'真实值':<30} {beta_true:<12.5f} {'—':<12} {'—':<25}")
print(f"{'OLS（有偏）':<28} {ols1.params['education']:<12.5f} {ols1.bse['education']:<12.5f} "
f"[{ols1.conf_int().loc['education', 0]:.5f}, {ols1.conf_int().loc['education', 1]:.5f}]")
print(f"{'IV (Q1 单一IV)':<28} {iv_single.params['education']:<12.5f} {iv_single.std_errors['education']:<12.5f} "
f"[{iv_single.conf_int().loc['education', 'lower']:.5f}, {iv_single.conf_int().loc['education', 'upper']:.5f}]")
print(f"{'IV (Q1+Q2+Q3 多IV)':<27} {iv_multi.params['education']:<12.5f} {iv_multi.std_errors['education']:<12.5f} "
f"[{iv_multi.conf_int().loc['education', 'lower']:.5f}, {iv_multi.conf_int().loc['education', 'upper']:.5f}]")
print(f"{'IV (多IV+年份控制)':<27} {iv_controls.params['education']:<12.5f} {iv_controls.std_errors['education']:<12.5f} "
f"[{iv_controls.conf_int().loc['education', 'lower']:.5f}, {iv_controls.conf_int().loc['education', 'upper']:.5f}]")📋 诊断检验
python
print("\n" + "=" * 70)
print("诊断检验")
print("=" * 70)
# 第一阶段诊断
print("\n📊 第一阶段诊断：")
print(iv_multi.first_stage.diagnostics)
# 过度识别检验（3 个 IV - 1 个内生变量 = 2 个过度识别约束）
print("\n📏 过度识别检验 (Wooldridge)：")
try:
overid = iv_multi.wooldridge_overid
print(f" 统计量: {overid.stat:.4f}")
print(f" p 值: {overid.pval:.4f}")
if overid.pval > 0.05:
print(f" ✅ 不拒绝 H₀：所有 IV 可能有效")
else:
print(f" ❌ 拒绝 H₀：至少一个 IV 可能无效")
except Exception as e:
print(f" {e}")
# 内生性检验
print("\n🔍 内生性检验 (Wu-Hausman)：")
try:
wh = iv_multi.wu_hausman()
print(f" 统计量: {wh.stat:.4f}")
print(f" p 值: {wh.pval:.6f}")
if wh.pval < 0.05:
print(f" ❌ 拒绝 H₀：教育是内生的，需要 IV")
else:
print(f" ✅ 不拒绝 H₀：教育可能是外生的")
except Exception as e:
print(f" {e}")💡 结果讨论
python
print("\n" + "=" * 70)
print("结果讨论")
print("=" * 70)
print("""
📌 主要发现：
1. OLS 估计的教育回报高于真实值
- 原因：能力偏误（高能力 → 高教育 + 高工资）
- OLS 把"能力的效应"也算进了"教育的效应"
2. IV 估计更接近真实值
- 出生季度是一个合理的 IV：它通过义务教育法影响教育
- 但出生季度本身不直接影响工资
3. IV 的标准误比 OLS 大得多
- 这是 IV 方法的代价：用效率换一致性
- 出生季度对教育的影响很微弱 → 标准误较大
4. 弱工具变量问题
- 出生季度对教育的影响非常小（F 统计量可能不够大）
- Bound, Jaeger & Baker (1995) 指出了这个问题
- 这是 IV 方法中经典的"弱 IV 争论"
5. LATE 解释
- IV 估计的是"Compliers"的效应
- 即那些因为出生季度不同而改变教育年限的人
- 这些人是"处于义务教育边缘的人"——他们可能恰好是教育回报较高的群体
""")📚 案例 2：Card (1995)
📖 论文背景
论文："Using Geographic Variation in College Proximity to Estimate the Return to Schooling"发表：Aspects of Labour Market Behaviour: Essays in Honour of John Vanderkamp, 1995
核心问题：教育的因果回报率是多少？
工具变量：是否在四年制学院附近长大（Proximity to College）
制度逻辑：
住在学院附近 → 上大学的成本更低（通勤、信息）
→ 更可能上大学
→ 更高的教育年限
但住在学院附近本身不直接影响能力或工资
（需要控制地区特征）🔢 模拟数据生成
python
np.random.seed(1995) # 致敬论文年份
n = 30000
# --- 个体特征 ---
ability = np.random.normal(0, 1, n)
family_income = np.random.normal(50, 15, n) # 家庭收入（千美元）
family_income = np.clip(family_income, 10, 120)
black = np.random.binomial(1, 0.15, n)
south = np.random.binomial(1, 0.35, n)
urban = np.random.binomial(1, 0.65, n)
# --- 工具变量：是否在学院附近 ---
# 城市地区更可能在学院附近
near_college = np.random.binomial(1, 0.4 + 0.2 * urban)
# --- 教育年限 ---
# 受能力、家庭收入、种族、是否在学院附近影响
education = (12
+ 0.6 * near_college # IV 的效应
+ 1.2 * ability # 能力的效应
+ 0.02 * family_income # 家庭收入
- 0.5 * black # 种族差异
- 0.3 * south # 地区差异
+ np.random.normal(0, 1.5, n))
education = np.clip(education, 6, 20)
education = np.round(education, 1)
# --- 工资 ---
# 真实教育回报 = 0.10（10%）
beta_true_card = 0.10
log_wage = (5.0
+ beta_true_card * education # 教育的因果效应
+ 0.20 * ability # 能力的效应（混淆！）
+ 0.003 * family_income # 家庭收入
- 0.15 * black # 种族工资差距
+ 0.10 * urban # 城市工资溢价
- 0.05 * south # 南方工资折扣
+ np.random.normal(0, 0.35, n))
df_card = pd.DataFrame({
'log_wage': log_wage,
'education': education,
'near_college': near_college,
'ability': ability,
'family_income': family_income,
'black': black,
'south': south,
'urban': urban
})
print("=" * 70)
print("案例 2：Card (1995) — 学院距离 IV")
print("=" * 70)
print(f"样本量: {n:,}")
print(f"真实教育回报: {beta_true_card} ({beta_true_card*100:.0f}%)")
print(f"\n数据概览:")
print(df_card.describe().round(3))
print(f"\n学院附近 vs 不在学院附近:")
print(f" 在学院附近的比例: {near_college.mean():.3f}")
print(f" 在学院附近：平均教育 = {df_card.loc[df_card['near_college']==1, 'education'].mean():.3f}")
print(f" 不在学院附近：平均教育 = {df_card.loc[df_card['near_college']==0, 'education'].mean():.3f}")
print(f" 差异: {df_card.loc[df_card['near_college']==1, 'education'].mean() - df_card.loc[df_card['near_college']==0, 'education'].mean():.3f}")📊 可视化
python
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# 图 1：学院距离与教育年限
group_educ = df_card.groupby('near_college')['education'].mean()
axes[0].bar(['远离学院\n(near=0)', '学院附近\n(near=1)'],
group_educ.values, color=['steelblue', 'coral'],
edgecolor='black', width=0.5)
axes[0].set_ylabel('平均教育年限', fontsize=12)
axes[0].set_title('第一阶段：学院距离 → 教育', fontsize=14, fontweight='bold')
diff_educ = group_educ[1] - group_educ[0]
axes[0].annotate(f'差异 = {diff_educ:.3f} 年',
xy=(0.5, group_educ.mean()),
fontsize=12, fontweight='bold', color='purple',
ha='center',
bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.3))
# 图 2：学院距离与工资
group_wage = df_card.groupby('near_college')['log_wage'].mean()
axes[1].bar(['远离学院\n(near=0)', '学院附近\n(near=1)'],
group_wage.values, color=['steelblue', 'coral'],
edgecolor='black', width=0.5)
axes[1].set_ylabel('平均对数工资', fontsize=12)
axes[1].set_title('简化形式：学院距离 → 工资', fontsize=14, fontweight='bold')
diff_wage = group_wage[1] - group_wage[0]
axes[1].annotate(f'差异 = {diff_wage:.4f}',
xy=(0.5, group_wage.mean()),
fontsize=12, fontweight='bold', color='purple',
ha='center',
bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.3))
# 图 3：教育年限的分布（按是否在学院附近）
df_near = df_card[df_card['near_college'] == 1]['education']
df_far = df_card[df_card['near_college'] == 0]['education']
axes[2].hist(df_far, bins=30, alpha=0.6, color='steelblue',
label='远离学院', density=True)
axes[2].hist(df_near, bins=30, alpha=0.6, color='coral',
label='学院附近', density=True)
axes[2].set_xlabel('教育年限', fontsize=12)
axes[2].set_ylabel('密度', fontsize=12)
axes[2].set_title('教育年限分布', fontsize=14, fontweight='bold')
axes[2].legend(fontsize=11)
plt.suptitle('Card (1995)：学院距离作为教育的工具变量', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('card1995_illustration.png', dpi=300, bbox_inches='tight')
plt.show()📈 逐步估计
python
# --- Step 1：朴素 OLS ---
ols_naive = smf.ols('log_wage ~ education', data=df_card).fit()
# --- Step 2：带控制变量的 OLS ---
ols_controls = smf.ols('log_wage ~ education + black + south + urban + family_income',
data=df_card).fit()
# --- Step 3：第一阶段回归 ---
first_stage_card = smf.ols(
'education ~ near_college + black + south + urban + family_income',
data=df_card
).fit()
# --- Step 4：IV/2SLS（不带控制变量）---
iv_simple = IV2SLS.from_formula(
'log_wage ~ 1 [education ~ near_college]',
data=df_card
).fit(cov_type='robust')
# --- Step 5：IV/2SLS（带控制变量）---
iv_controls_card = IV2SLS.from_formula(
'log_wage ~ 1 + black + south + urban + family_income [education ~ near_college]',
data=df_card
).fit(cov_type='robust')
# 结果汇总
print("\n" + "=" * 70)
print("Card (1995) 回归结果汇总")
print("=" * 70)
print(f"\n{'模型':<35} {'教育系数':<12} {'标准误':<12} {'备注':<20}")
print("-" * 79)
print(f"{'真实值':<35} {beta_true_card:<12.5f} {'—':<12} {'—':<20}")
print(f"{'(1) OLS（简单）':<33} {ols_naive.params['education']:<12.5f} {ols_naive.bse['education']:<12.5f} {'有偏（高估）':<20}")
print(f"{'(2) OLS（+控制变量）':<31} {ols_controls.params['education']:<12.5f} {ols_controls.bse['education']:<12.5f} {'仍有偏':<20}")
print(f"{'(3) IV（简单）':<33} {iv_simple.params['education']:<12.5f} {iv_simple.std_errors['education']:<12.5f} {'可能有偏':<20}")
print(f"{'(4) IV（+控制变量）':<31} {iv_controls_card.params['education']:<12.5f} {iv_controls_card.std_errors['education']:<12.5f} {'推荐':<20}")
# 第一阶段诊断
print(f"\n--- 第一阶段诊断 ---")
print(f"near_college 系数: {first_stage_card.params['near_college']:.4f}")
print(f"t 统计量: {first_stage_card.tvalues['near_college']:.2f}")
F_card = first_stage_card.tvalues['near_college'] ** 2
print(f"F 统计量 (≈ t²): {F_card:.2f}")
if F_card > 10:
print(f"✅ F = {F_card:.2f} > 10，工具变量足够强")
else:
print(f"⚠️ F = {F_card:.2f} < 10，弱工具变量警告！")📋 完整诊断
python
print("\n" + "=" * 70)
print("完整诊断报告")
print("=" * 70)
# Wald 估计量验证
wald_card = (df_card.loc[df_card['near_college']==1, 'log_wage'].mean() -
df_card.loc[df_card['near_college']==0, 'log_wage'].mean()) / \
(df_card.loc[df_card['near_college']==1, 'education'].mean() -
df_card.loc[df_card['near_college']==0, 'education'].mean())
print(f"\nWald 估计量（不含控制变量）: {wald_card:.5f}")
print(f"IV 估计量（不含控制变量）: {iv_simple.params['education']:.5f}")
print(f"（两者应相等）")
# 内生性检验
print(f"\n🔍 内生性检验 (Wu-Hausman):")
try:
wh_card = iv_controls_card.wu_hausman()
print(f" 统计量: {wh_card.stat:.4f}")
print(f" p 值: {wh_card.pval:.6f}")
if wh_card.pval < 0.05:
print(f" ❌ 拒绝 H₀：教育是内生的")
else:
print(f" ✅ 不拒绝 H₀：教育可能是外生的")
except Exception as e:
print(f" {e}")💡 Card (1995) 的关键讨论
python
print("\n" + "=" * 70)
print("Card (1995) 关键讨论")
print("=" * 70)
print("""
📌 Card (1995) 的主要发现和争论：
1. IV 估计值 vs OLS 估计值
- 在原论文中，IV 估计（~13%）竟然高于 OLS（~7%）
- 这与"能力偏误导致 OLS 高估"的直觉矛盾！
2. 可能的解释：
a) LATE ≠ ATE
- IV 估计的是 Compliers 的效应
- Compliers = 因为住在学院附近才上大学的人
- 这些人可能是"受教育回报最高"的群体（高边际回报）
- 对他们来说，教育的回报确实更高
b) 测量误差
- 教育年限的自报数据存在测量误差
- 经典测量误差 → 衰减偏差（attenuation bias）
- OLS 因此低估了真实效应
- IV 可以修正这个衰减偏差
c) 排他性约束的争议
- 住在学院附近 ≈ 住在城市
- 城市本身可能有更高的工资
- 如果排他性被违反，IV 估计是有偏的
3. 教训：
- IV 的有效性高度依赖排他性约束的可信性
- 需要仔细论证 + 控制尽可能多的混淆变量
- LATE 的政策含义需要谨慎解读
""")🔧 完整 Python 工作流模板
一站式 IV 分析模板
python
def complete_iv_analysis(df, dep_var, endog_var, instruments, controls=None,
title="IV/2SLS 分析"):
"""
一站式 IV 分析函数
Parameters
----------
df : DataFrame, 数据
dep_var : str, 因变量
endog_var : str, 内生变量
instruments : list, 工具变量列表
controls : list, 控制变量列表（可选）
title : str, 分析标题
"""
from linearmodels.iv import IV2SLS
import statsmodels.formula.api as smf
print("=" * 70)
print(f"📊 {title}")
print("=" * 70)
# --- 构建公式 ---
ctrl_str = ' + '.join(controls) if controls else ''
iv_str = ' + '.join(instruments)
if controls:
ols_formula = f'{dep_var} ~ {endog_var} + {ctrl_str}'
fs_formula = f'{endog_var} ~ {iv_str} + {ctrl_str}'
iv_formula = f'{dep_var} ~ 1 + {ctrl_str} [{endog_var} ~ {iv_str}]'
else:
ols_formula = f'{dep_var} ~ {endog_var}'
fs_formula = f'{endog_var} ~ {iv_str}'
iv_formula = f'{dep_var} ~ 1 [{endog_var} ~ {iv_str}]'
# --- 1. OLS ---
print("\n📈 [1/5] OLS 估计")
ols = smf.ols(ols_formula, data=df).fit()
print(f" {endog_var}: {ols.params[endog_var]:.5f} (SE={ols.bse[endog_var]:.5f})")
# --- 2. 第一阶段 ---
print(f"\n🔧 [2/5] 第一阶段回归")
fs = smf.ols(fs_formula, data=df).fit()
for iv in instruments:
print(f" {iv}: coef={fs.params[iv]:.4f}, t={fs.tvalues[iv]:.2f}, p={fs.pvalues[iv]:.6f}")
if len(instruments) == 1:
F = fs.tvalues[instruments[0]] ** 2
else:
# 部分 F 统计量
fs_restricted = smf.ols(
f'{endog_var} ~ {ctrl_str}' if controls else f'{endog_var} ~ 1',
data=df
).fit()
q = len(instruments)
F = ((fs_restricted.ssr - fs.ssr) / q) / (fs.ssr / fs.df_resid)
print(f"\n 第一阶段 F 统计量: {F:.2f}")
print(f" {'✅ F > 10' if F > 10 else '⚠️ F < 10（弱 IV 警告）'}")
# --- 3. IV/2SLS ---
print(f"\n📐 [3/5] IV/2SLS 估计")
iv_result = IV2SLS.from_formula(iv_formula, data=df).fit(cov_type='robust')
ci = iv_result.conf_int()
print(f" {endog_var}: {iv_result.params[endog_var]:.5f} (SE={iv_result.std_errors[endog_var]:.5f})")
print(f" 95% CI: [{ci.loc[endog_var, 'lower']:.5f}, {ci.loc[endog_var, 'upper']:.5f}]")
# --- 4. 内生性检验 ---
print(f"\n🔍 [4/5] 内生性检验")
try:
wh = iv_result.wu_hausman()
print(f" Wu-Hausman: stat={wh.stat:.4f}, p={wh.pval:.6f}")
print(f" {'❌ D 是内生的' if wh.pval < 0.05 else '✅ D 可能是外生的'}")
except:
print(" 无法计算")
# --- 5. 过度识别检验 ---
print(f"\n📏 [5/5] 过度识别检验")
if len(instruments) > 1:
try:
overid = iv_result.wooldridge_overid
print(f" J 统计量: {overid.stat:.4f}, p={overid.pval:.4f}")
print(f" {'✅ IV 可能有效' if overid.pval > 0.05 else '❌ 至少一个 IV 可能无效'}")
except:
print(" 无法计算")
else:
print(" （恰好识别，无法进行 J 检验）")
# --- 汇总表 ---
print(f"\n{'='*50}")
print(f"{'方法':<15} {'系数':<12} {'标准误':<12}")
print(f"{'-'*39}")
print(f"{'OLS':<15} {ols.params[endog_var]:<12.5f} {ols.bse[endog_var]:<12.5f}")
print(f"{'IV/2SLS':<15} {iv_result.params[endog_var]:<12.5f} {iv_result.std_errors[endog_var]:<12.5f}")
print(f"{'差异':<15} {iv_result.params[endog_var]-ols.params[endog_var]:<12.5f}")
print(f"{'='*50}")
return iv_result, ols, fs
# --- 使用模板 ---
print("\n\n" + "#" * 70)
print("使用模板分析 Card (1995) 数据")
print("#" * 70)
iv_r, ols_r, fs_r = complete_iv_analysis(
df=df_card,
dep_var='log_wage',
endog_var='education',
instruments=['near_college'],
controls=['black', 'south', 'urban', 'family_income'],
title="Card (1995): 教育的因果回报"
)📊 结果可视化总结
python
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# --- 图 1：两个案例的 OLS vs IV ---
cases = ['AK (1991)', 'Card (1995)']
ols_vals = [ols1.params['education'], ols_controls.params['education']]
iv_vals = [iv_multi.params['education'], iv_controls_card.params['education']]
true_vals = [beta_true, beta_true_card]
x = np.arange(len(cases))
width = 0.25
bars1 = axes[0].bar(x - width, ols_vals, width, label='OLS', color='salmon',
edgecolor='black', alpha=0.8)
bars2 = axes[0].bar(x, iv_vals, width, label='IV/2SLS', color='mediumseagreen',
edgecolor='black', alpha=0.8)
bars3 = axes[0].bar(x + width, true_vals, width, label='真实值', color='gold',
edgecolor='black', alpha=0.8)
axes[0].set_xticks(x)
axes[0].set_xticklabels(cases, fontsize=12)
axes[0].set_ylabel('教育回报率', fontsize=12)
axes[0].set_title('OLS vs IV/2SLS vs 真实值', fontsize=14, fontweight='bold')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# --- 图 2：IV 估计的置信区间 ---
methods = ['OLS\n(AK)', 'IV\n(AK)', 'OLS\n(Card)', 'IV\n(Card)']
estimates = [ols1.params['education'],
iv_multi.params['education'],
ols_controls.params['education'],
iv_controls_card.params['education']]
ses = [ols1.bse['education'],
iv_multi.std_errors['education'],
ols_controls.bse['education'],
iv_controls_card.std_errors['education']]
colors_ci = ['salmon', 'mediumseagreen', 'salmon', 'mediumseagreen']
for i, (est, se, color) in enumerate(zip(estimates, ses, colors_ci)):
axes[1].errorbar(i, est, yerr=1.96*se, fmt='o', color=color,
markersize=10, capsize=5, linewidth=2, capthick=2)
axes[1].axhline(y=beta_true, color='blue', linestyle='--', alpha=0.5, label=f'AK 真实值 ({beta_true})')
axes[1].axhline(y=beta_true_card, color='red', linestyle='--', alpha=0.5, label=f'Card 真实值 ({beta_true_card})')
axes[1].set_xticks(range(len(methods)))
axes[1].set_xticklabels(methods, fontsize=11)
axes[1].set_ylabel('教育回报率 (95% CI)', fontsize=12)
axes[1].set_title('估计值与 95% 置信区间', fontsize=14, fontweight='bold')
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
plt.suptitle('经典 IV 研究结果比较', fontsize=16, fontweight='bold')
plt.tight_layout()
plt.savefig('iv_classic_comparison.png', dpi=300, bbox_inches='tight')
plt.show()📝 本节要点
Angrist & Krueger (1991)：出生季度作为教育的 IV
- 通过义务教育法的制度性变异提供外生变异
- 弱 IV 问题是该研究最大的争议点
- 开创了现代 IV 研究的范式
Card (1995)：学院距离作为教育的 IV
- 地理变异提供了"准实验"
- IV 估计高于 OLS → LATE vs ATE、测量误差修正
- 排他性约束需要仔细论证
实践经验：
- 始终报告第一阶段 F 统计量
- 始终比较 OLS 和 IV 估计值
- 讨论 LATE 的含义和外部有效性
- 做好全套诊断检验
下一节：8.6 本章小结 — 核心概念回顾、方法比较和推荐文献。