12.4 经典案例与Python实现

从理论到实践：完整复现California Proposition 99

---

本节目标

• 完整理解California Proposition 99的研究背景
• 掌握SCM的Python手动实现（从零开始）
• 学会使用SparseSC库
• 解读权重表和协变量平衡表
• 绘制Gap Plot和安慰剂图

---

经典案例：California Proposition 99

背景

Abadie, Diamond & Hainmueller (2010)

政策: 1988年11月,加州通过Proposition 99:

• 香烟税从10美分/包 → 35美分/包（提高250%）
• 全美最高烟草税
• 目标:降低吸烟率

研究问题: Prop 99对加州香烟消费量的因果效应?

挑战:

• 只有1个处理单位（加州）
• 加州在政策前就与其他州不同（更自由、更健康意识）
• 全国趋势:吸烟率普遍下降

数据

结果变量: 人均香烟消费量（包/人/年）

时间: 1970-2000 (政策时点:1988)

单位: 39个州（排除:有大规模烟草控制政策的州）

协变量():

1. 政策前香烟消费量（1975, 1980, 1988）
2. 啤酒消费量（1984-1988平均）
3. 人均收入（1980, 1988）
4. 15-24岁人口占比（1980, 1988）
5. 香烟零售价（1990）

合成加州的构建

最优权重:

州	权重
犹他	0.276
内华达	0.234
科罗拉多	0.200
蒙大拿	0.154
新墨西哥	0.136
其他34个州	0.000

解释: 只有5个州获得正权重！

为什么这5个州？

特征	真实加州	合成加州	差异
人均收入（1988）	$21.5k	$21.2k	1.4%
啤酒消费量	24.3	24.1	0.8%
15-24岁占比（1988）	17.6%	17.3%	1.7%
香烟消费（1980）	120.2	120.5	0.2%
香烟消费（1988）	90.1	91.0	1.0%

惊人的拟合！

结果：因果效应

图示: 真实加州 vs 合成加州

香烟消费量（包/人/年）

140|
120|     ●────●────●────●●●●●●●  真实加州
100|     ◯────◯────◯────◯╲      合成加州
 80|                      ╲◯────◯────◯
 60|                        ╲________↓
 40|                         政策效应
 20|
  0|_____|_____|_____|_____|_____|_____|_____
   1970  1975  1980  1985  1990  1995  2000
                           ↑
                        Prop 99

定量结果:

1988-2000年平均效应:

相对效应:

解释: Prop 99使加州香烟消费量下降21.5%！

---

Python实现：完整代码

手动实现：从零开始

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize, LinearConstraint
from cvxopt import matrix, solvers

# 中文设置
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams['axes.unicode_minus'] = False
sns.set_style("whitegrid")

# =================================================================
# 第1步：权重优化（二次规划）
# =================================================================

def optimize_weights(X1, X0, V):
    """
    求解最优权重W

    min (X1 - X0 @ W)' @ V @ (X1 - X0 @ W)
    s.t. W >= 0, sum(W) = 1

    参数:
    ------
    X1: (K,) 处理单位特征
    X0: (K, J) 对照单位特征矩阵
    V: (K, K) 权重矩阵

    返回:
    ------
    W: (J,) 最优权重
    """
    K, J = X0.shape

    # 使用cvxopt求解二次规划
    # 构造QP: min 0.5 * w' @ P @ w + q' @ w
    X0_weighted = V @ X0
    P = 2 * (X0.T @ X0_weighted)  # (J, J)
    q = -2 * (X1.T @ X0_weighted)  # (J,)

    # 约束: sum(W) = 1, W >= 0
    A = np.ones((1, J))
    b = np.ones(1)
    G = -np.eye(J)
    h = np.zeros(J)

    # 转换为cvxopt格式
    P = matrix(P)
    q = matrix(q)
    G = matrix(G)
    h = matrix(h)
    A = matrix(A)
    b = matrix(b)

    # 求解
    solvers.options['show_progress'] = False
    sol = solvers.qp(P, q, G, h, A, b)

    W = np.array(sol['x']).flatten()
    return W

# =================================================================
# 第2步：优化V矩阵
# =================================================================

def optimize_V(Y1, Y0, X1, X0, T0):
    """
    优化V使政策前MSPE最小

    参数:
    ------
    Y1: (T,) 处理单位结果
    Y0: (T, J) 对照单位结果矩阵
    X1, X0: 特征矩阵
    T0: 政策时点

    返回:
    ------
    V: (K, K) 最优V矩阵
    """
    K = X1.shape[0]

    def loss(v):
        """MSPE损失函数"""
        # 确保V是对角矩阵且非负
        v = np.abs(v)
        V = np.diag(v / v.sum())  # 归一化

        # 求最优W
        W = optimize_weights(X1, X0, V)

        # 计算政策前预测误差
        Y_synthetic = Y0[:T0] @ W
        mspe = np.mean((Y1[:T0] - Y_synthetic)**2)

        return mspe

    # 初始V
    v_init = np.ones(K)

    # 优化（使用Nelder-Mead）
    result = minimize(loss, v_init, method='Nelder-Mead')

    v_opt = np.abs(result.x)
    V_opt = np.diag(v_opt / v_opt.sum())

    return V_opt

# =================================================================
# 第3步：合成控制主函数
# =================================================================

def synthetic_control(Y1, Y0, X1, X0, T0, optimize_v=True):
    """
    合成控制法主函数

    参数:
    ------
    Y1: (T,) 处理单位结果时间序列
    Y0: (T, J) 对照单位结果矩阵
    X1: (K,) 处理单位特征
    X0: (K, J) 对照单位特征矩阵
    T0: 政策时点索引
    optimize_v: 是否优化V（否则用单位矩阵）

    返回:
    ------
    results: dict, 包含权重、合成结果、效应等
    """
    T = len(Y1)
    K, J = X0.shape

    # 优化V
    if optimize_v:
        print("优化V矩阵...")
        V = optimize_V(Y1, Y0, X1, X0, T0)
    else:
        V = np.eye(K)

    # 计算最优W
    print("计算最优权重...")
    W = optimize_weights(X1, X0, V)

    # 构造合成控制
    Y_synthetic = Y0 @ W

    # 计算效应
    effects = Y1 - Y_synthetic

    # MSPE
    mspe_pre = np.mean((Y1[:T0] - Y_synthetic[:T0])**2)
    mspe_post = np.mean((Y1[T0:] - Y_synthetic[T0:])**2)

    # 平均效应
    att = np.mean(effects[T0:])

    results = {
        'weights': W,
        'V': V,
        'Y_synthetic': Y_synthetic,
        'effects': effects,
        'mspe_pre': mspe_pre,
        'mspe_post': mspe_post,
        'ATT': att
    }

    return results

# =================================================================
# 第4步：安慰剂检验（Placebo Test）
# =================================================================

def placebo_test(Y, X, T0, treated_idx=0):
    """
    安慰剂检验：对每个对照单位做"假处理"

    参数:
    ------
    Y: (T, J+1) 所有单位结果矩阵（第0列是真实处理单位）
    X: (K, J+1) 所有单位特征矩阵
    T0: 政策时点
    treated_idx: 真实处理单位索引

    返回:
    ------
    placebo_effects: (T, J+1) 所有单位的"效应"
    """
    T, J_plus_1 = Y.shape
    placebo_effects = np.zeros((T, J_plus_1))
    mspe_ratios = []

    for j in range(J_plus_1):
        # 将单位j作为"处理单位"
        Y1_fake = Y[:, j]
        Y0_fake = np.delete(Y, j, axis=1)

        X1_fake = X[:, j]
        X0_fake = np.delete(X, j, axis=1)

        # 运行SCM
        try:
            results = synthetic_control(Y1_fake, Y0_fake, X1_fake, X0_fake, T0, optimize_v=False)
            placebo_effects[:, j] = results['effects']

            # MSPE比率
            ratio = results['mspe_post'] / results['mspe_pre']
            mspe_ratios.append(ratio)
        except:
            placebo_effects[:, j] = np.nan
            mspe_ratios.append(np.nan)

    return placebo_effects, np.array(mspe_ratios)

# =================================================================
# 第5步：可视化
# =================================================================

def plot_synthetic_control(years, Y1, Y_synthetic, T0, title="合成控制法结果"):
    """绘制合成控制结果"""
    fig, axes = plt.subplots(2, 1, figsize=(14, 10), sharex=True)

    # 子图1：趋势对比
    ax1 = axes[0]
    ax1.plot(years, Y1, 'o-', color='black', linewidth=2,
             label='真实处理单位', markersize=6)
    ax1.plot(years, Y_synthetic, 's--', color='blue', linewidth=2,
             label='合成对照', markersize=5, alpha=0.7)
    ax1.axvline(x=years[T0], color='red', linestyle='--', linewidth=2, alpha=0.7)
    ax1.text(years[T0]+0.5, ax1.get_ylim()[1]*0.9, '政策实施',
             fontsize=12, color='red', fontweight='bold')
    ax1.set_ylabel('结果变量', fontsize=12, fontweight='bold')
    ax1.set_title(title, fontsize=14, fontweight='bold')
    ax1.legend(loc='best', fontsize=11)
    ax1.grid(True, alpha=0.3)

    # 子图2：处理效应（Gap Plot）
    ax2 = axes[1]
    effects = Y1 - Y_synthetic
    ax2.plot(years, effects, 'o-', color='darkgreen', linewidth=2, markersize=6)
    ax2.axhline(y=0, color='black', linestyle='-', linewidth=1, alpha=0.5)
    ax2.axvline(x=years[T0], color='red', linestyle='--', linewidth=2, alpha=0.7)
    ax2.fill_between(years[T0:], 0, effects[T0:], alpha=0.3, color='green')
    ax2.set_xlabel('年份', fontsize=12, fontweight='bold')
    ax2.set_ylabel('处理效应', fontsize=12, fontweight='bold')
    ax2.set_title('估计的因果效应（Gap Plot）', fontsize=13, fontweight='bold')
    ax2.grid(True, alpha=0.3)

    plt.tight_layout()
    return fig

def plot_placebo_test(years, placebo_effects, T0, treated_idx=0):
    """绘制安慰剂检验"""
    fig, ax = plt.subplots(figsize=(14, 8))

    T, J = placebo_effects.shape

    # 绘制所有安慰剂单位（灰色）
    for j in range(J):
        if j != treated_idx:
            ax.plot(years, placebo_effects[:, j],
                   color='gray', alpha=0.3, linewidth=1)

    # 绘制真实处理单位（粗黑线）
    ax.plot(years, placebo_effects[:, treated_idx],
           color='black', linewidth=3, label='真实处理单位')

    ax.axhline(y=0, color='black', linestyle='-', linewidth=1)
    ax.axvline(x=years[T0], color='red', linestyle='--', linewidth=2, alpha=0.7)
    ax.text(years[T0]+0.5, ax.get_ylim()[1]*0.9, '政策实施',
           fontsize=12, color='red', fontweight='bold')

    ax.set_xlabel('年份', fontsize=12, fontweight='bold')
    ax.set_ylabel('处理效应', fontsize=12, fontweight='bold')
    ax.set_title('安慰剂检验：真实处理单位 vs 假处理单位', fontsize=14, fontweight='bold')
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)

    plt.tight_layout()
    return fig

# =================================================================
# 应用示例：模拟Prop 99
# =================================================================

# 设置随机种子
np.random.seed(42)

# 模拟数据
T = 31  # 1970-2000
T0 = 19  # 1988
J = 38  # 38个对照州

years = np.arange(1970, 2001)

# 生成协变量（K=5个特征）
X1 = np.array([120, 90, 24, 21.5, 0.176])  # 加州特征
X0 = np.random.randn(5, J) * 10 + X1.reshape(-1, 1) + np.random.randn(5, J) * 5

# 生成结果变量
# 基础趋势：线性下降
baseline = 130 - 2 * np.arange(T)

# 加州：政策后额外下降
Y1 = baseline + np.random.randn(T) * 3
Y1[T0:] -= 20  # 政策效应

# 对照州：只有基础趋势
Y0 = np.zeros((T, J))
for j in range(J):
    state_effect = np.random.randn() * 5
    Y0[:, j] = baseline + state_effect + np.random.randn(T) * 3

# 运行SCM
print("="*60)
print("合成控制法分析")
print("="*60)

results = synthetic_control(Y1, Y0, X1, X0, T0, optimize_v=True)

# 输出结果
print(f"\n政策前MSPE: {results['mspe_pre']:.3f}")
print(f"政策后MSPE: {results['mspe_post']:.3f}")
print(f"平均处理效应(ATT): {results['ATT']:.3f}")

# 显示权重
weights_df = pd.DataFrame({
    '州': [f'州{j+1}' for j in range(J)],
    '权重': results['weights']
})
weights_df = weights_df[weights_df['权重'] > 0.01].sort_values('权重', ascending=False)
print(f"\n主要权重(>1%):")
print(weights_df.to_string(index=False))

# 可视化
fig1 = plot_synthetic_control(years, Y1, results['Y_synthetic'], T0,
                               "加州Proposition 99的效应：合成控制法")
plt.show()

# 安慰剂检验
print("\n执行安慰剂检验...")
Y_all = np.column_stack([Y1, Y0])
X_all = np.column_stack([X1.reshape(-1, 1), X0])
placebo_effects, mspe_ratios = placebo_test(Y_all, X_all, T0, treated_idx=0)

fig2 = plot_placebo_test(years, placebo_effects, T0, treated_idx=0)
plt.show()

# 计算p值
true_mspe_ratio = mspe_ratios[0]
p_value = np.mean(mspe_ratios >= true_mspe_ratio)
print(f"\nMSPE比率 p值: {p_value:.3f}")
print(f"结论: {'拒绝原假设(有显著效应)' if p_value < 0.05 else '不能拒绝原假设'}")

---

使用SparseSC库

# 安装: pip install SparseSC

from SparseSC import fit

# 准备数据
outcomes = Y_all  # (T, J+1)
treated_units = [0]  # 加州是第0列
T0_sc = T0

# 拟合
sc_result = fit(
    outcomes=outcomes,
    treated_units=treated_units,
    T0=T0_sc,
    verbose=1
)

# 提取结果
Y_synthetic_sparse = sc_result.synthetic_outcomes[:, 0]
effects_sparse = sc_result.effects[:, 0]

print(f"\n[SparseSC] 平均处理效应: {np.mean(effects_sparse[T0:]):.3f}")

---

Python工具

• 手动实现: cvxopt, scipy.optimize
• 现成库: SparseSC, augsynth
• 新方法: synthdid

---

本节小结

核心要点

1. Prop 99案例：SCM成功识别了加州烟草税政策降低香烟消费21.5%的因果效应
2. 手动实现：权重优化的核心是二次规划（QP），可用cvxopt求解
3. 安慰剂检验：对所有对照单位做"假处理"，构建RMSPE比率的参考分布
4. 可视化：趋势对比图和Gap Plot是SCM论文的标准图表
5. SparseSC：现成Python库，适合快速原型实现

---

<< 上一节：12.3 统计推断与稳健性 | 下一节：12.5 本章小结 >>