合成控制法 Synthetic Control Method：为单个处理对象构造"合成对照组"

* ============================================================
* 合成控制法 SCM 完整演示
* 场景：某省（id=1）在 2010 年实施自贸区政策，其余 20 省未实施
* 目标：估计自贸区对 lnGDP 的因果效应
* ============================================================
clear all
set seed 20240101
set obs 2100                              // 21 地区 × 100 年
gen id = ceil(_n/100)                      // 地区编号 1..21
bysort id: gen year = 1950 + _n            // 年份 1950..2049

* --- 1. 模拟结果变量：地区固定效应 + 共同趋势 + 政策效应 ---
gen u_i = rnormal(0, 0.5) in 1/21          // 地区固定效应（每地区抽一次）
bysort id: replace u_i = u_i[1]
gen t = year - 1950
gen y = 2 + 0.05*t + u_i + rnormal(0, 0.2) // 无政策时的 lnGDP

* 处理单位 id=1 在 2010 年起获得 +0.3 的政策效应
replace y = y + 0.3 if id==1 & year>=2010

* --- 2. 预测变量（政策前 2000-2009 年的均值）---
* 为简化，这里用几个常数型变量 + 政策前结果变量均值
gen x1 = rnormal(0, 1) in 1/21
bysort id: replace x1 = x1[1]
gen x2 = rnormal(0, 1) in 1/21
bysort id: replace x2 = x2[1]

* --- 3. 安装外部命令 ---
* ssc install synth, replace
* ssc install synth_runner, replace

* --- 4. 跑 synth（经典版）---
* synth 语法：synth 结果变量 预测变量(英文变量列表), trunit(处理单位) trperiod(政策期)
* 这里我们用 x1 x2 + 政策前若干年的结果变量作为预测变量
preserve
    keep if year <= 2009                  // 只用政策前数据估计权重
    * 构造政策前结果变量的若干年份均值作为预测变量
    * synth 要求数据是 wide 形式：每地区一行，列是各年 y
    collapse (mean) y, by(id year)
    reshape wide y, i(id) j(year)
    * 用 2000-2009 年的 y 作为预测变量（典型 SCM 做法）
    synth y2009 y2008 y2007 y2006 y2005 x1 x2, ///
        trunit(1) trperiod(2010) keep(scm_result)
restore

* --- 5. 用 synth_runner 跑完整流程 + 安慰剂检验 ---
* synth_runner 是更现代的封装，自动画 gap plot、算 RMSPE ratio
preserve
    keep if year <= 2019                  // 政策后 10 年窗口
    reshape wide y x1 x2, i(id) j(year)
    * synth_runner 语法（需先 install）
    * synth_runner y, trunit(1) trperiod(2010) ///
    *     pre(1980/2009) post(2010/2019) ///
    *    covariates(x1 x2) ///
    *    mspeperiod(1980/2009) ///
    *    figure_name(gap_plot)
    * 输出：
    *   - _data.dta     : 处理 vs 合成的逐年差距
    *   - placebo_results: 所有单位的伪效应
    *   - gap plot      : 真实效应是否跳出信封
restore

* --- 6. 手动解读关键输出 ---
* synth 输出表 1：权重 W*
*   关注：哪几个 donor 权重 > 0.01？加起来是否 = 1？
* synth 输出表 2：预测变量均值对比
*   关注：处理单位 vs 合成单位在每个预测变量上的差距是否 < 0.1 SD
* synth 输出表 3：RMSPE
*   pre-RMSPE 应 < post-RMSPE / 5

# ============================================================
# 合成控制法 Python 实现（pysyncon）
# ============================================================
import numpy as np
import pandas as pd
from pysyncon import Dataprep, Synth

rng = np.random.default_rng(20240101)

# --- 1. 构造面板：21 地区 × 70 年 ---
n_units, n_years = 21, 70
years = np.arange(1980, 2050)
ids = np.arange(1, n_units + 1)

data = []
for i in ids:
    u_i = rng.normal(0, 0.5)
    x1 = rng.normal(0, 1)
    x2 = rng.normal(0, 1)
    for t, yr in enumerate(years):
        y = 2 + 0.05 * t + u_i + rng.normal(0, 0.2)
        if i == 1 and yr >= 2010:
            y += 0.3                       # 政策效应
        data.append({"id": i, "year": yr, "y": y, "x1": x1, "x2": x2})
df = pd.DataFrame(data)

# --- 2. Dataprep：指定处理单位、政策期、预测变量 ---
dataprep = Dataprep(
    foo=df,
    predictors=["x1", "x2"],
    predictors_op="mean",                  # 预测变量在 pre 期取均值
    time_predictors_prior=np.arange(2000, 2010),
    dependent="y",
    unit_variable="id",
    time_variable="year",
    treatment_identifier=1,                # 处理单位 id=1
    controls_identifier=list(range(2, 22)),  # donor pool
    time_optimize_ssr=np.arange(1990, 2010),  # 用 1990-2009 拟合
)

# --- 3. 跑 SCM ---
synth = Synth()
synth.fit(dataprep, optim_method="Nelder-Mead")
print(synth.weights())                    # 各 donor 的权重 W*
print(synth.summary())                    # 预测变量均值对比 + RMSPE

# --- 4. 画处理 vs 合成轨迹图 ---
synth.plot_treatment_outcome()

# --- 5. 安慰剂检验（placebo in-space）---
from pysyncon import Placebo
placebo = Placebo(dataprep, synth)
placebo.plot_placebo(figsize=(10, 6))     # gap plot
print(placebo.p_value(period_post=np.arange(2010, 2020)))