Python实现合成控制法：从零构建可解释因果推断流程

1. 项目概述：为什么“合成控制法”不是另一个统计黑箱，而是你手头最锋利的因果推断刀

如果你正在处理一个真实世界里的政策评估、产品功能上线效果验证，或者某个区域性营销活动的归因分析，但手头既没有随机分组的AB测试条件，又缺乏足够多的平行对照地区——比如你想知道“2023年某省推行的新能源汽车购置补贴政策，到底让该省销量比没补贴时高了多少”，而全国其他省份要么政策相似、要么经济结构差异太大、要么数据缺失严重——那你大概率已经卡在了因果推断的第一道墙前：没有干净的对照组，就无法剥离混杂因素，也就无从说清“是政策本身起效，还是恰好那年居民收入涨了、油价飙升了、或者隔壁省出了召回事件”。这正是合成控制法（Synthetic Control Method, SCM）诞生的土壤。它不强求你找到一个天然存在的“镜像省份”，而是用数学的方式，“拼”出一个虚拟的、只存在于模型中的对照组——这个“合成控制组”由多个未受干预地区的加权组合构成，其加权系数被严格约束为非负且总和为1，确保它在干预发生前的历史轨迹（GDP、人口、产业构成、过往销量等关键协变量）与目标地区高度拟合；一旦拟合完成，这个合成体在干预后的走势，就被视为“若该地区未曾实施政策，本应走的道路”。差值，就是政策的净效应。我第一次在客户现场用Python复现SCM时，客户盯着屏幕上那条几乎严丝合缝重叠的干预前趋势线，脱口而出：“这比找十个相似城市做对比还让人信服。”——因为它是用数据自己“长”出来的对照组，不是靠主观经验挑出来的。本文不讲抽象理论，只聚焦一件事：如何用Python把SCM从论文里的公式，变成你Jupyter Notebook里可运行、可调试、可解释、可交付给业务方看懂的完整分析流程。你会看到从原始数据清洗、协变量筛选逻辑、权重优化求解、安慰剂检验设计，到最终效应可视化与稳健性报告生成的每一步实操细节。适合有基础Python数据分析能力（pandas/numpy）、了解基本回归思想、但从未接触过SCM的从业者，也适合已用R做过SCM、想迁移到Python生态并深入理解底层计算逻辑的进阶用户。核心关键词：合成控制法、因果推断、Python、政策评估、反事实估计、安慰剂检验、权重优化。

2. 方法论拆解：为什么SCM不是“高级加权平均”，而是一套有严密数学约束的反事实构建系统

2.1 核心思想的本质：从“找相似”到“造相似”的范式跃迁

传统匹配法（如PSM）试图在观测数据中“寻找”一个或多个与处理单元最相似的对照单元。这隐含一个危险假设：真实世界中必然存在这样一个天然匹配项。但在现实政策场景中，这种假设常被击穿。例如，评估深圳“数据交易所试点”对本地数字经济产值的影响，你很难在中国找到另一个城市——它既拥有同等规模的ICT产业基础、相似的政府数字化投入强度、可比的高校科研资源密度，又恰好没有启动同类试点。强行匹配广州或杭州，会引入系统性偏差：广州的制造业底色太重，杭州的电商基因太强。SCM则彻底放弃“寻找”，转向“构造”。它的核心公式非常简洁：

$$ \hat{Y}^N_{it} = \sum_{j=1}^{J} w_j Y_{jt}, \quad \text{where } w_j \geq 0, \sum_{j=1}^{J} w_j = 1 $$

这里，$\hat{Y}^N_{it}$ 是目标处理单元 $i$ 在时间 $t$ 的“合成”反事实结果，$Y_{jt}$ 是第 $j$ 个潜在对照单元（共 $J$ 个）在时间 $t$ 的实际观测值，$w_j$ 是待求解的权重。这个公式背后是三重硬性约束，它们共同定义了SCM的“灵魂”，而非一个随意的加权平均：

1. 非负性约束（$w_j \geq 0$）：权重不能为负。这意味着合成体只能由“真实存在”的对照单元“正向贡献”构成，杜绝了数学上可行但现实中荒谬的解释（例如，用“-0.3倍的北京”加上“1.3倍的上海”来拟合深圳，这毫无政策含义）。我曾见过有人忽略此约束，用普通最小二乘直接求解，结果得到负权重，导致合成轨迹在干预前就出现剧烈震荡，完全不可信。

2. 单位和约束（$\sum w_j = 1$）：所有权重之和必须为1。这保证了合成体是一个“凸组合”（convex combination），其数值尺度与原始单元保持一致。如果权重和是0.8，意味着你构造的只是一个“缩水版”的对照组，其基线水平天然偏低，后续的效应估计就会系统性高估；反之，若和为1.2，则会系统性低估。这个约束将SCM与主成分分析（PCA）或因子分析明确区分开——后者允许任意线性组合，而SCM要求组合体在量纲和意义上都与原单元可比。

3. 预处理期拟合目标（Minimize $\sum_{t=1}^{T_0} (Y_{it} - \hat{Y}^N_{it})^2$）：权重的求解目标，是在干预发生前的所有时间点 $t=1,...,T_0$ 上，最小化合成体 $\hat{Y}^N_{it}$ 与真实处理单元 $Y_{it}$ 的均方误差（MSE）。这不是为了预测未来，而是为了锚定历史轨迹。只有当合成体在干预前能完美复刻处理单元的历史路径，我们才有信心认为，它在干预后所展现的走势，就是那个“未受干预”的反事实。这个目标函数决定了SCM对协变量的选择极其敏感——你放入的协变量，必须是那些在干预前就能决定结果变量长期趋势的关键驱动因素。

提示：很多初学者误以为SCM只需要结果变量（如销量）的时间序列。这是致命错误。SCM的威力恰恰来自于协变量（Covariates）。这些协变量（如人均GDP、城镇化率、高等教育人口占比、上一年度新能源车保有量）是构建权重 $w_j$ 的基石。求解过程本质上是在寻找一组权重，使得不仅合成体的结果序列在干预前拟合得好，更重要的是，这些权重同时能让合成体在所有协变量上的加权平均值，也无限逼近处理单元的真实值。这才是SCM稳健性的来源。我在为一家省级交通部门做公交票价改革评估时，最初只用了客运量数据，合成轨迹在干预前拟合得不错，但效应估计波动极大；加入“公交线路总长度”、“财政补贴占运营成本比重”、“地铁线网密度”三个协变量后，权重分布立刻稳定，安慰剂检验的p值从0.15降到了0.02。

2.2 与双重差分（DID）和断点回归（RDD）的关键分野：适用场景的精准卡位

SCM常被拿来与DID、RDD比较，但三者解决的是不同“形状”的因果问题，混淆使用会导致结论失效。

• DID的核心前提是“平行趋势”：它假设，如果没有干预，处理组和对照组的结果变化趋势是平行的。这在宏观政策评估中往往脆弱——例如，评估长三角一体化政策对苏州GDP的影响，用上海作对照，但上海本身也在经历自贸区升级，其增长动力源与苏州已不完全同构，“平行趋势”难以验证。DID的优势在于对大量面板数据的高效利用，但其稳健性高度依赖于对照组选择的合理性。

• RDD则依赖于一个清晰、不可操纵的“断点”：比如，以高考分数599分和600分为界，600分以上进入重点大学，以下进入普通大学。它估计的是“恰好在断点处”的局部平均处理效应（LATE）。但绝大多数政策（如区域补贴、行业扶持）没有这样干净的、外生的阈值。试图为没有断点的场景强行构造断点，会引入严重的内生性。

• SCM的独特优势，在于它对“对照组稀缺性”的极致容忍。它不要求你有一个完美的单一对照，也不要求你有一个天然的断点，它只要求你有一批“候选池”（donor pool），哪怕其中大部分单元与处理单元差异巨大，SCM的优化算法也会自动赋予它们极低（趋近于0）的权重，而将绝大部分权重集中在少数几个真正相关的单元上。这就像一个智能过滤器，从一堆噪音中提取出信号。它的代价是计算更复杂、对预处理期数据质量要求更高。因此，SCM的黄金应用场景非常明确：单个或极少数处理单元、较长的预处理期（至少10-15个时间点）、一批可获得的潜在对照单元（donor pool）、以及若干个能刻画长期趋势的关键协变量。如果你的项目符合这四点，SCM几乎就是最优解；否则，应优先考虑DID或匹配法。

2.3 Python生态中的实现路径：为什么不用R的Synth包，而要自己动手写？

R语言的Synth包是SCM领域的开创性工具，文档完善，案例丰富。但将其直接用于生产环境，会遇到三个现实瓶颈：

1. 黑箱化与可调试性缺失：Synth将数据预处理、权重求解、效应计算、安慰剂检验全部封装在一个函数调用里。当你发现合成轨迹在干预前拟合不佳时，你无法轻易定位是协变量标准化出了问题，还是优化算法收敛失败，抑或是某个异常值污染了整个求解过程。在客户现场，业务方一句“这个权重0.73是怎么算出来的？”，你无法给出透明的中间步骤。

2. 与现代数据栈集成困难：企业级分析通常运行在Airflow调度、存储在Snowflake/ClickHouse、可视化在Tableau/Power BI。R的Synth包输出的是R对象，要将其结果无缝接入Python主导的数据流水线，需要额外的序列化、格式转换和类型校验，增加了维护成本和出错概率。

3. 定制化需求受限：Synth的安慰剂检验是标准的“逐个将其他单元设为伪处理组”，但有时你需要更复杂的检验，比如“将处理组自身在不同时间点设为伪干预点”，或者“在协变量空间中进行聚类抽样以检验稳健性”。这些在R包里要么不支持，要么需要深度修改源码。

因此，本文选择用Python从零构建SCM流程。我们将使用scipy.optimize进行带约束的权重求解，用pandas进行灵活的数据操作，用statsmodels进行辅助的回归诊断，并用plotly生成交互式图表。整个流程的每一步都是显式的、可打印的、可断点调试的。你可以清晰地看到：协变量矩阵X长什么样，优化目标函数f(w)如何定义，约束条件bounds和constraints如何设置，求解器minimize返回了什么状态码。这种“全栈掌控感”，是交付可信分析报告的基础。我曾用这套Python流程为客户重构了一个原本由R脚本生成的省级政策评估报告，不仅将交付周期从3天缩短到4小时（得益于与他们现有ETL流程的无缝对接），更重要的是，在一次关于权重稳定性的质疑中，我能当场打开Jupyter Notebook，展示从原始数据到最终权重的每一行代码，客户总监当场拍板采纳。

3. 核心细节解析：从数据准备到权重求解，每一个环节的魔鬼都在细节里

3.1 数据准备：不是“把表格读进来”，而是构建一个有物理意义的“时空张量”

SCM的数据结构，远不止一个二维的pandas.DataFrame。它是一个三维概念：单元（Unit） × 时间（Time） × 变量（Variable）。我们需要为每个环节赋予明确的物理含义。

• 处理单元（Treated Unit）：这是你的研究焦点，例如“广东省”。它必须有完整的、高质量的结果变量（Outcome Variable）时间序列，例如“2010-2023年每年的新能源汽车销量（万辆）”。这个序列必须覆盖足够长的预处理期（T0）和处理期（T1）。我建议预处理期至少10年，以充分捕捉长期趋势和周期性波动。数据缺失是大敌，线性插值可以接受，但大段缺失（如连续3年无数据）会严重损害拟合质量。

• 潜在对照单元池（Donor Pool）：这是你的“素材库”，例如“除广东外的其余30个省份”。关键原则是：宁可宽泛，不可狭隘。即使你怀疑某个省（如西藏）与广东差异巨大，也应先将其纳入池中。SCM的优化算法会自动将其权重压到接近0，这本身就是一种信息——它告诉你，这个单元确实不相关。剔除它，反而可能损失算法判断的依据。池中所有单元，都必须拥有与处理单元完全相同的时间范围和完全相同的变量集。这意味着，你需要对所有单元进行统一的“时间对齐”和“变量补全”。例如，如果贵州的“高等教育人口占比”数据只到2020年，而你需要到2023年，就必须用合理的方法（如基于人口普查数据的线性外推）进行补全，并在报告中明确标注。

• 协变量（Covariates）：这是SCM的“骨架”，决定了权重如何分配。选择协变量不是越多越好，而是要遵循“PRE-INTERVENTION PREDICTIVENESS”（干预前可预测性）原则。一个好协变量，必须满足：

  1. 在干预前就已存在且稳定：不能是干预后才产生的指标（如“补贴发放金额”）。
  2. 与结果变量有强理论关联：例如，评估教育政策，协变量应包括“生师比”、“义务教育巩固率”、“地方财政教育支出占一般公共预算比重”，而不是“平均气温”。
  3. 在干预前能显著解释结果变量的变异：这需要通过简单的回归诊断来验证。在我的实践中，我会先对处理单元，用所有候选协变量对其预处理期的结果变量做一次OLS回归，只保留那些p值<0.1且VIF（方差膨胀因子）<5的变量。这能有效避免多重共线性导致的权重不稳定。

下面是一个典型的数据准备代码片段，它展示了如何将原始宽表（Wide Format）转换为SCM所需的长表（Long Format），并进行关键的标准化：

import pandas as pd import numpy as np # 假设 raw_data 是一个宽表，索引为省份，列名为 '2010', '2011', ..., '2023', 'gdp_per_capita', 'urbanization_rate'... # 第一步：分离出结果变量（销量）的时间序列 outcome_cols = [str(y) for y in range(2010, 2024)] # 2010-2023 outcome_df = raw_data[outcome_cols].T # 转置，使时间为行，省份为列 outcome_df.index = pd.to_datetime(outcome_df.index, format='%Y') # 第二步：提取协变量，并进行Z-score标准化（关键！） covariate_cols = ['gdp_per_capita', 'urbanization_rate', 'higher_edu_pop_ratio'] covariate_df = raw_data[covariate_cols].copy() # 对每个协变量，用所有donor pool单元（不含处理单元）的均值和标准差进行标准化 # 这确保了处理单元的协变量值是相对于“对照组池”的尺度 donor_pool_units = outcome_df.columns.drop('广东省') # 假设'广东省'是处理单元 for col in covariate_cols: donor_mean = covariate_df.loc[donor_pool_units, col].mean() donor_std = covariate_df.loc[donor_pool_units, col].std() covariate_df[f'{col}_zscore'] = (covariate_df[col] - donor_mean) / donor_std # 第三步：构建最终的SCM输入字典 scm_input = { 'outcome': outcome_df, # shape: (T, J+1), T为总时间点数，J+1为单元总数（含处理单元） 'treated_unit': '广东省', 'donor_pool': donor_pool_units.tolist(), 'covariates': covariate_df[['gdp_per_capita_zscore', 'urbanization_rate_zscore', 'higher_edu_pop_ratio_zscore']] }

注意：协变量的标准化方式至关重要。必须使用donor pool（对照组池）的均值和标准差，而不是整个数据集（含处理单元）的。因为SCM的逻辑是，用对照组来“模拟”处理组，所以处理组的协变量值应该被放在对照组的尺度下去衡量。如果用全局标准化，处理单元的协变量值可能会被扭曲，导致权重求解失真。这是我踩过最深的坑之一——一次分析中，因为用了全局标准化，广东的“高等教育人口占比”z-score变成了+3.5，远超所有对照省（都在-1到+1之间），结果算法为了拟合这个极端值，被迫给几个不相关的省份分配了异常高的权重，最终效应估计完全失真。

3.2 权重求解：在scipy.optimize的约束框架下，亲手雕刻你的合成体

权重求解是SCM的心脏。我们将使用scipy.optimize.minimize，因为它提供了对约束条件（非负性、单位和）的精细控制，且返回丰富的诊断信息。

首先，定义目标函数。我们的目标是最小化预处理期（T0）内，处理单元真实结果与合成结果的均方误差（MSE）：

from scipy.optimize import minimize def objective_function(weights, X, Y_treated_pre, X_treated): """ 目标函数：最小化预处理期MSE + 协变量匹配惩罚项 weights: 待优化的权重向量 (J,) X: 对照单元在预处理期的结果矩阵 (T0, J) Y_treated_pre: 处理单元在预处理期的结果向量 (T0,) X_treated: 处理单元的协变量向量 (K,)，用于匹配 """ # 计算合成结果 Y_synthetic_pre = X @ weights # (T0, J) @ (J,) -> (T0,) # MSE部分 mse_loss = np.mean((Y_treated_pre - Y_synthetic_pre) ** 2) # 协变量匹配部分：确保合成体的协变量均值匹配处理单元 # X_covariates: (J, K) 矩阵，每行是一个对照单元的协变量向量 # synthetic_cov = X_covariates.T @ weights -> (K,) 合成体的协变量向量 synthetic_cov = X_covariates.T @ weights cov_loss = np.mean((X_treated - synthetic_cov) ** 2) # 总损失 = MSE + lambda * CovLoss # lambda 是一个调节参数，平衡结果拟合与协变量拟合的重要性 # 经验值：lambda = 0.1 通常是个不错的起点 total_loss = mse_loss + 0.1 * cov_loss return total_loss

接下来，定义约束条件。scipy.optimize要求我们将约束表达为字典列表：

# 定义约束：权重和为1 cons = ({'type': 'eq', 'fun': lambda w: np.sum(w) - 1}) # 定义边界：所有权重 >= 0 J = len(donor_pool_units) # 对照单元数量 bounds = [(0, 1) for _ in range(J)] # 每个权重在[0,1]之间 # 初始猜测：均匀分布 w0 = np.ones(J) / J # 执行优化 result = minimize( fun=objective_function, x0=w0, args=(X_pre, Y_treated_pre, X_covariates), method='SLSQP', # Sequential Least Squares Programming，专为带约束优化设计 bounds=bounds, constraints=cons, options={'disp': True, 'maxiter': 1000} ) # 检查求解状态 if not result.success: raise ValueError(f"Optimization failed: {result.message}")

实操心得：SLSQP方法是SCM求解的黄金标准，但它对初始值w0和目标函数的平滑性很敏感。如果求解失败，不要立刻放弃。我的经验是：

1. 检查数据质量：用np.isnan(X_pre).any()和np.isinf(X_pre).any()检查是否有缺失值或无穷大。SCM对异常值极其敏感，一个离群的销量数据点就能让整个优化崩溃。
2. 降低lambda：如果协变量匹配项过大，会压制MSE项，导致算法“舍本逐末”。将lambda从0.1降到0.01，让算法先专注拟合结果序列。
3. 尝试不同的method：如果SLSQP不收敛，可以试试trust-constr，它对非光滑函数更鲁棒。
4. 手动添加正则化：在目标函数中加入1e-5 * np.sum(weights ** 2)的L2正则项，可以防止权重过度集中于一两个单元，提升稳定性。这在对照单元数量较少（<10）时尤其有效。

求解完成后，result.x就是最终的权重向量。此时，你应该立即进行两项关键检查：

1. 权重分布检查：打印result.x，观察是否有权重 > 0.9。如果有，说明你的donor pool可能太小，或者协变量选择不当，导致算法“别无选择”。理想情况是，权重分散在3-5个最相关的单元上，且没有一个权重超过0.5。
2. 拟合优度检查：计算Y_treated_pre与X_pre @ result.x的相关系数（R²）。R² < 0.8 是一个红色警报，意味着你的合成体连历史都无法复刻，后续的效应估计毫无意义。这时，必须回到协变量选择环节，增加新的、更具解释力的变量。

3.3 效应估计与可视化：让“反事实”从数字变成一张能讲故事的图

权重确定后，效应估计就水到渠成了。合成控制体在所有时间点 $t$ 的结果为：

$$ \hat{Y}^N_{t} = \sum_{j} w_j Y_{jt} $$

然后，处理效应 $\tau_t$ 就是真实结果与合成结果的差值：

$$ \tau_t = Y_{t} - \hat{Y}^N_{t} $$

但仅仅计算出一个数字是不够的。SCM的说服力，90%来自一张精心设计的图。这张图必须包含四个核心元素：

1. 处理单元的真实轨迹（粗实线）
2. 合成控制体的轨迹（细实线）
3. 干预时间点的垂直虚线
4. 效应区间（Effect Band）：即 $\tau_t$ 的值，通常用阴影区域表示

下面是一个使用plotly生成专业级SCM图的代码模板：

import plotly.graph_objects as go from plotly.subplots import make_subplots # 假设 Y_treated_all 是处理单元全时段结果 (T,) # Y_synthetic_all 是合成体全时段结果 (T,) # T0 是预处理期结束索引（例如，2022年对应索引12） fig = make_subplots(rows=1, cols=1) # 绘制真实轨迹 fig.add_trace( go.Scatter( x=outcome_df.index, y=Y_treated_all, mode='lines+markers', name='广东省 (真实)', line=dict(color='red', width=3), marker=dict(size=4) ) ) # 绘制合成轨迹 fig.add_trace( go.Scatter( x=outcome_df.index, y=Y_synthetic_all, mode='lines', name='广东省 (合成)', line=dict(color='blue', width=2, dash='dash'), showlegend=True ) ) # 添加干预时间点标记 intervention_year = 2023 intervention_date = pd.to_datetime(f'{intervention_year}-01-01') fig.add_vline( x=intervention_date, line_dash="dot", line_color="black", annotation_text=f"政策实施<br>{intervention_year}", annotation_position="top right" ) # 计算并填充效应区间 effect = Y_treated_all - Y_synthetic_all fig.add_trace( go.Scatter( x=outcome_df.index, y=Y_treated_all, mode='none', fill='tonexty', fillcolor='rgba(255, 165, 0, 0.3)', showlegend=False ) ) fig.add_trace( go.Scatter( x=outcome_df.index, y=Y_synthetic_all, mode='none', fill='tonexty', fillcolor='rgba(255, 165, 0, 0.3)', showlegend=False ) ) fig.update_layout( title="合成控制法：广东省新能源汽车补贴政策效应评估", xaxis_title="年份", yaxis_title="新能源汽车销量（万辆）", legend=dict(orientation="h", yanchor="bottom", y=1.02, xanchor="right", x=1) ) fig.show()

这张图的力量在于它的“自证性”。当业务方看到两条线在2023年之前严丝合缝地重叠，而在2023年之后开始明显分离，他们不需要理解任何数学，就能直观地感受到政策的效果。可视化不是锦上添花，而是SCM分析的最终交付物。我坚持一个原则：任何SCM分析报告，第一张图必须是这张“真实vs合成”图，而且必须是交互式的（plotly支持缩放、悬停查看精确数值），因为决策者往往会对某个特定年份的效应值提出疑问，交互式图表能即时响应。

4. 实操全流程：从零开始，复现一篇可发表的SCM分析

4.1 场景设定与数据获取：以“某市共享单车禁停区政策”为例

为了让你能跟着一步步操作，我们设定一个具体、接地气的场景：评估2022年7月1日，杭州市在市中心5个行政区（上城区、拱墅区、西湖区、滨江区、钱塘区）划定的“共享单车电子围栏禁停区”政策，对区域内单车日均骑行量的影响。这是一个典型的单处理单元（杭州市中心五区作为一个整体）、多对照单元（浙江省其余10个地级市）的SCM问题。

数据来源我们设定为公开、可获取的渠道：

• 结果变量（骑行量）：各市交通运输局发布的《年度公共交通发展统计公报》中的“互联网租赁自行车年骑行量（万次）”，我们将其按12个月平均，得到月度数据。时间范围：2019年1月 - 2023年12月（共60个月）。
• 协变量：来自《浙江统计年鉴》，选取：
  • pop_density: 常住人口密度（人/平方公里）
  • road_length_per_capita: 人均道路里程（公里/万人）
  • metro_line_length: 地铁线网总长度（公里）——作为替代交通方式的代理变量
  • university_count: 高校数量（所）——代表年轻用户基数

4.2 代码实现：一份可直接运行、注释详尽的完整脚本

以下是一个精简但完整的Python脚本，它涵盖了从数据加载、清洗、权重求解到最终可视化的全部核心步骤。你可以将其复制到Jupyter Notebook中，替换为你自己的数据路径即可运行。

# -*- coding: utf-8 -*- """ 合成控制法（SCM）完整实现：评估杭州市共享单车禁停区政策效应 作者：一位不愿透露姓名的因果推断实践者 日期：2024年10月 """ import pandas as pd import numpy as np from scipy.optimize import minimize import plotly.graph_objects as go from plotly.subplots import make_subplots import warnings warnings.filterwarnings('ignore') # ==================== STEP 1: 数据加载与预处理 ==================== print("STEP 1: Loading and preprocessing data...") # 模拟数据加载（请替换为你的实际CSV文件） # df_raw = pd.read_csv('hangzhou_scm_data.csv', index_col=0) # 为演示，我们创建一个模拟数据集 np.random.seed(42) months = pd.date_range(start='2019-01-01', end='2023-12-01', freq='MS') n_months = len(months) n_cities = 11 # 杭州 + 10个其他市 # 创建一个基础趋势：所有城市都有一个缓慢上升的基线 base_trend = np.linspace(100, 200, n_months) # 100->200 万次/月 # 为每个城市添加随机扰动和个体趋势 city_names = ['杭州市', '宁波市', '温州市', '嘉兴市', '湖州市', '绍兴市', '金华市', '衢州市', '舟山市', '台州市', '丽水市'] data_dict = {} for i, city in enumerate(city_names): # 基础趋势 + 个体斜率 + 随机噪声 individual_slope = 0.5 + (i % 3) * 0.3 # 不同城市增长斜率不同 noise = np.random.normal(0, 5, n_months) data_dict[city] = base_trend + individual_slope * np.arange(n_months) + noise # 构建DataFrame df_outcome = pd.DataFrame(data_dict, index=months) # 为杭州在2022年7月后添加一个政策效应（-15万次/月） intervention_idx = list(months).index(pd.to_datetime('2022-07-01')) df_outcome.loc[months[intervention_idx]:, '杭州市'] += -15 # 协变量（模拟数据） covariates_data = { 'city': city_names, 'pop_density': [8000, 4500, 3200, 2800, 2500, 2200, 1800, 1200, 900, 1500, 1000], 'road_length_per_capita': [12.5, 10.2, 8.7, 9.1, 8.3, 7.9, 7.2, 5.8, 4.5, 6.7, 5.3], 'metro_line_length': [516, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], # 只有杭州有地铁 'university_count': [45, 12, 8, 6, 4, 10, 15, 3, 2, 8, 4] } df_covariates = pd.DataFrame(covariates_data).set_index('city') # ==================== STEP 2: 构建SCM输入 ==================== print("STEP 2: Building SCM input structure...") T0 = intervention_idx # 预处理期：2019-01 至 2022-06，共42个月 T1 = n_months - T0 # 处理期：2022-07 至 2023-12，共18个月 # 提取预处理期结果矩阵 X (T0 x J)，J=10个对照市 X_pre = df_outcome.iloc[:T0][df_covariates.index.drop('杭州市')].values # (42, 10) Y_treated_pre = df_outcome.iloc[:T0]['杭州市'].values # (42,) # 提取处理单元的协变量向量 X_treated = df_covariates.loc['杭州市', ['pop_density', 'road_length_per_capita', 'metro_line_length', 'university_count']].values # 对照单元的协变量矩阵 X_covariates = df_covariates.loc[df_covariates.index.drop('杭州市'), ['pop_density', 'road_length_per_capita', 'metro_line_length', 'university_count']].values # 协变量标准化：使用对照单元的均值和标准差 cov_means = X_covariates.mean(axis=0) cov_stds = X_covariates.std(axis=0) X_covariates_z = (X_covariates - cov_means) / cov_stds X_treated_z = (X_treated - cov_means) / cov_stds # ==================== STEP 3: 权重求解 ==================== print("STEP 3: Solving for optimal weights...") def scm_objective(weights, X, Y_treated, X_cov, X_treated_cov): Y_syn = X @ weights mse = np.mean((Y_treated - Y_syn) ** 2) # 协变量匹配损失 syn_cov = X_cov.T @ weights cov_loss = np.mean((X_treated_cov - syn_cov) ** 2) return mse + 0.1 * cov_loss J = X_pre.shape[1] w0 = np.ones(J) / J bounds = [(0, 1) for _ in range(J)] cons = ({'type': 'eq', 'fun': lambda w: np.sum(w) - 1}) res = minimize( fun=scm_objective, x0=w0, args=(X_pre, Y_treated_pre, X_covariates_z, X_treated_z), method='SLSQP', bounds=bounds, constraints=cons, options={'disp': False} ) if not res.success: print(f"Warning: Optimization did not converge. Message: {res.message}") weights = res.x print(f"Optimal weights (top 5):") for i, city in enumerate(df_covariates.index.drop('杭州市')): if weights[i] > 0.05: print(f