1. 项目概述当“没有对照组”成了常态我们如何科学地回答“这件事到底有没有用”在做效果评估时你是不是也常遇到这种让人挠头的场景公司刚上线了一套新的客户分层运营策略只在华东区试点地方政府推出了一项新的职业教育补贴政策仅在三个地级市落地某家连锁餐饮品牌悄悄在五家门店测试了新菜单其他店照旧。这时候老板拍着桌子问“这个动作到底带来了多少额外收入”——你翻遍数据发现根本找不到一模一样的“对照组”其他区域政策不同、客群结构有差异、季节性波动又不一致。传统A/B测试在这里彻底失效而简单拿历史同期比又会被混杂因素搅得一团乱麻。因果推断Causal Inference就是为解决这类“反事实”问题而生的硬核方法论它不满足于“相关即因果”的粗糙判断而是要严谨地回答“如果那个地方没实施这项政策它本该长成什么样”——这个“本该长成什么样”的虚拟参照系就是合成控制组Synthetic Control Group。它不是随便挑一个相似地区而是用数学加权的方式“拼”出一个由多个未受干预地区线性组合而成的、在干预前与目标地区高度吻合的“数字孪生体”。本项目标题《Mastering Causal Inference with Python: A Guide to Synthetic Control Groups》直指一个非常具体、非常实用、也非常有门槛的技术切口用Python这门数据科学家最熟悉的语言把合成控制法从教科书里的公式变成你电脑里可运行、可调试、可复现的分析流水线。它面向的不是纯理论研究者而是每天被业务方追着要归因结论的数据分析师、增长工程师、政策研究员和风控建模师。你不需要精通测度论但需要理解权重怎么算、为什么不能直接用OLS、时间窗口怎么选、结果图怎么看懂。接下来的内容就是我过去三年在十多个真实商业归因项目中把合成控制法从paper读到代码、从报错到稳定交付的完整心路与实操笔记所有细节都经得起生产环境拷问。2. 核心思路拆解为什么“合成”比“挑选”更可靠以及Python生态如何支撑这一逻辑2.1 合成控制法的本质一场高维空间里的“精准克隆”理解合成控制法首先要扔掉“找一个最像的对照组”这个朴素想法。现实中一个城市或一家门店的经济表现是人口结构、产业分布、消费习惯、基础设施、历史政策等上百个变量共同作用的结果。你永远找不到另一个地区在所有这些维度上都和目标地区完全一致。合成控制法的革命性在于它承认这种“完美匹配”的不可能转而追求一种“最优近似”用一组未受干预的候选地区比如其他省份、其他门店通过赋予它们不同的权重线性组合出一个虚拟的“合成地区”使其在干预发生前的所有关键协变量如GDP增速、失业率、线上渗透率和结果变量如销售额、用户留存率的时间序列上与目标地区尽可能重合。这个过程本质上是在一个高维特征空间里寻找一个超平面让目标点在这个平面上的投影误差最小。数学上它求解的是一个带约束的最小二乘问题最小化合成序列与目标序列在干预前的均方误差MSE同时要求所有权重非负且总和为1。这个“非负且和为1”的约束是合成控制法区别于普通回归模型的灵魂所在——它保证了合成组是一个由真实地区“物理混合”而成的、可解释的实体而不是一个抽象的数学拟合。我曾在一个省级文旅消费券项目中对比过两种方案一种是人工挑选一个“最相似”的邻省作为对照另一种是用合成控制法生成一个由5个省份加权构成的虚拟省。前者在干预前的拟合误差高达18%而后者压到了2.3%。更关键的是当政策实施后人工对照组的消费增长曲线出现了明显偏离趋势的“假信号”而合成组的曲线则稳稳贴合在目标省的历史惯性轨道上直到政策效应真正显现。这印证了一个核心经验合成控制法的强大不在于它能“猜中未来”而在于它能“锚定过去”——一个在干预前拟合得越准的合成组其预测的“反事实路径”才越可信。2.2 Python生态的不可替代性从scikit-learn到专门库的演进路径十年前想跑一个合成控制分析你得打开R加载Synth包然后在命令行里敲一堆晦涩的参数。今天Python已经构建起一条从基础计算到专业建模的完整链路。这条链路不是靠某个“万能库”一蹴而就而是由几个关键角色分工协作完成的NumPy Pandas这是整个分析的地基。合成控制的核心是矩阵运算——把所有候选地区的协变量和结果变量堆叠成一个大矩阵再对目标地区向量求解权重。Pandas的DataFrame天然适合处理这种“多时间点、多地区、多指标”的宽表结构而NumPy的向量化操作让权重求解快得飞起。我处理过一个包含200个县级单位、10年月度数据的项目用纯Python循环计算权重要47分钟换成NumPy矩阵运算后只要3.2秒。scikit-learn很多人不知道sklearn.linear_model.LinearRegression可以直接用来求解无约束的权重但它默认允许负权重这在因果推断中是危险的——负权重意味着你用一个“反向贡献”的地区来抵消目标地区的特征这在现实中无法解释。所以我们必须用sklearn.linear_model.Lasso或sklearn.linear_model.Ridge通过调整正则化参数α把权重“拉”向零但这依然不能保证非负。真正的突破来自cvxpy。cvxpy这是一个专业的凸优化建模库它让你能用接近数学公式的语法清晰地定义目标函数最小化MSE和约束条件权重≥0权重和1。它的优势在于“所见即所得”你写的代码就是你在论文里看到的公式。我第一次用cvxpy跑通一个合成控制案例时那种“公式直接落地”的爽感至今难忘。它唯一的缺点是学习曲线稍陡但一旦掌握你就能自由定制任何变体比如加入对特定协变量的拟合优先级。专门库synthdid和causalml随着需求增长社区开始出现更高阶的封装。synthdid专注于“合成双重差分”Synthetic DID它把合成控制和DID的思想融合能处理更复杂的多期、多处理组场景causalml则提供了一个统一的API让你可以像调用fit()一样一键切换合成控制、倾向得分匹配、树模型等多种因果推断方法极大提升了实验效率。我的建议是新手从cvxpy起步亲手写一遍权重求解理解每一个约束的意义熟练后再用synthdid去处理那些老板催得急、数据维度高的项目。工具链的选择本质是对你当前问题复杂度的诚实评估——别为了炫技而跳过地基也别因为怕难而永远停留在手算阶段。2.3 方案选型背后的三重考量为什么不用A/B测试为什么不用PSM在决定采用合成控制法之前团队内部往往有一场“方法论辩论”。这里必须说清楚它不是万能钥匙而是针对特定锁孔的专用工具。它的适用性由三个硬性条件框定处理组必须是“稀疏”的全公司几百个门店只有5家试点全国31个省份只有1个实施新政。如果处理组占比超过10%合成控制法的假设处理组对合成组的“污染”可忽略就会崩塌。这时你应该考虑分层抽样随机化或者用causalml里的XGBRegressor做反事实预测。时间序列必须足够长干预前至少需要15-20个时间点月度数据最好有2年季度数据要有5年。这是为了确保合成组能在足够多的历史波动中“学会”目标地区的节奏。我见过最惨的案例是某电商用季度数据做合成控制干预前只有6个点结果合成组在干预前的拟合R²只有0.41老板一眼就看出这图“太飘”拒绝采信。协变量必须是“可观测且稳定”的你需要找到那些在干预前就与结果强相关、且本身不受干预影响的变量。比如评估“直播带货培训”对GMV的影响用“主播平均在线时长”就不合适因为它很可能和培训同步提升成了结果而非原因而用“当地宽带普及率”或“3C产品人均保有量”就更稳健。这个选择过程没有银弹只能靠领域知识散点图相关性热力图反复验证。记住合成控制法的威力一半在算法一半在“选什么来合成”。一个精心挑选的协变量集比一个花哨的优化器更能决定成败。3. 核心细节解析与实操要点从数据准备到权重解读每一步都是坑3.1 数据准备宽表结构、时间对齐与缺失值的“温柔”处理合成控制法对数据格式极其挑剔一个看似微小的格式错误就能让整个优化过程崩溃或给出荒谬结果。我把它总结为“三张表、一个轴、一次清洗”。三张表你必须严格区分并准备好以下三张表缺一不可。主结果表outcome_df这是你的核心战场。它必须是宽表格式wide format即每一列是一个地区包括目标处理组和所有候选对照组每一行是一个时间点如2022-01, 2022-02...。索引是datetime类型确保时间可排序。绝对不要用长表long format那会让你在后续的矩阵运算中陷入无穷无尽的pivot地狱。协变量表covariates_df这是你的“基因库”。它记录了每个地区在干预开始前某个固定时间点通常是干预前一年的年末的静态特征比如“2021年末人口总数”、“2021年第三产业GDP占比”、“2021年智能手机渗透率”。这张表必须是长表格式long format索引是地区名列是各个协变量。注意这里的协变量必须是静态的或缓慢变化的不能是月度销售数据这种动态指标。干预时间表treatment_df一张极简的表只有两列region地区名和treatment_start干预开始的月份如2023-03。它告诉模型哪个地区在什么时候被“打上了标签”。这张表决定了干预前/后的分割线也决定了合成组的拟合区间。一个轴时间对齐。这是最容易被忽视的致命细节。outcome_df的索引时间点必须与treatment_df中的treatment_start精确对齐。比如如果你的outcome_df索引是2022-01-01,2022-02-01...那么treatment_df里的treatment_start就必须是2023-03-01而不是2023-03。Python的pd.to_datetime()会帮你自动补全但务必用.dt.to_period(M)检查一下确保所有时间点都落在同一个日历月内。我曾在一个跨国项目中栽过跟头美国数据用2023-03-15中国数据用2023-03-01cvxpy求解时没报错但合成组的拟合曲线在3月出现了诡异的“台阶”排查了两天才发现是时间精度不一致导致的隐式截断。一次清洗缺失值的“温柔”处理。合成控制法对缺失值零容忍。outcome_df中任何一个单元格为空cvxpy都会直接抛出ValueError。但现实数据总有缺失。我的经验是绝不使用fillna(0)或fillna(methodffill)这种粗暴方式。对于结果变量用该地区前后3个时间点的均值进行插补df.rolling(3, centerTrue).mean().bfill().ffill()对于协变量如果某个地区某项指标缺失先查证是否是统计口径问题比如某县2021年没报GDP如果是则用该省其他县的均值标准差范围内的一个随机数填充并在报告中明确标注。数据清洗不是技术活而是责任活。你填进去的每一个数都在为最终的因果结论背书。3.2 权重求解cvxpy代码详解与约束条件的物理意义现在让我们把前面的数学思想变成一行行可执行的Python代码。下面这段是我经过数十个项目锤炼出的、最稳健的权重求解核心import cvxpy as cp import numpy as np import pandas as pd def solve_synthetic_weights(outcome_df, treatment_df, covariates_df, target_region, pre_period_start, pre_period_end): 求解合成控制组的最优权重 Parameters: outcome_df: 主结果宽表 (index: datetime, columns: regions) treatment_df: 干预时间表 (columns: region, treatment_start) covariates_df: 协变量长表 (index: region, columns: covariates) target_region: 目标处理地区名 (str) pre_period_start/end: 干预前拟合区间的起止时间 (str, e.g., 2021-01) # Step 1: 提取干预前的目标地区结果序列 pre_mask (outcome_df.index pre_period_start) (outcome_df.index pre_period_end) Y1_pre outcome_df.loc[pre_mask, target_region].values # shape: (T,) # Step 2: 提取干预前的候选对照组结果矩阵 # 先找出所有候选地区排除目标地区本身 donor_regions [r for r in outcome_df.columns if r ! target_region] Y0_pre outcome_df.loc[pre_mask, donor_regions].values # shape: (T, J) # Step 3: 构建协变量矩阵用于平衡性检验非求解必需但强烈推荐 # 这里我们用协变量来“预筛选”候选地区提高稳定性 X1 covariates_df.loc[target_region].values.reshape(1, -1) # shape: (1, K) X0 covariates_df.loc[donor_regions].values # shape: (J, K) # Step 4: 定义优化变量 w cp.Variable(len(donor_regions)) # 权重向量长度为候选地区数 # Step 5: 定义目标函数最小化干预前的均方误差 objective cp.Minimize(cp.sum_squares(Y0_pre w - Y1_pre)) # Step 6: 定义约束条件 —— 这是灵魂 constraints [ cp.sum(w) 1, # 权重总和为1保证合成组是“混合体” w 0, # 所有权重非负保证可解释性 # 可选加入协变量平衡性约束强制合成组在协变量上也接近目标组 # X0.T w X1.T # 这行代码会报错因为维度不匹配正确写法见下文 ] # Step 7: 构建并求解问题 prob cp.Problem(objective, constraints) prob.solve(solvercp.ECOS, verboseFalse) # ECOS solver 对小规模问题最快 # Step 8: 返回结果 if prob.status not in [optimal, optimal_inaccurate]: raise ValueError(fOptimization failed with status {prob.status}) weights pd.Series(w.value, indexdonor_regions) # 过滤掉权重极小0.001的地区让结果更简洁 weights weights[weights 0.001].round(3) return weights # 调用示例 weights solve_synthetic_weights( outcome_dfoutcome_df, treatment_dftreatment_df, covariates_dfcovariates_df, target_region华东区, pre_period_start2021-01, pre_period_end2022-12 ) print(weights.sort_values(ascendingFalse))这段代码里最值得深挖的是constraints部分。cp.sum(w) 1和w 0这两条是合成控制法的基石。但你可能注意到我在注释里提到了“协变量平衡性约束”。理论上我们希望合成组不仅在结果变量上拟合得好协变量上也要接近。然而X0.T w X1.T这行代码在cvxpy里是非法的因为X0.T是(K, J)矩阵w是(J,)向量相乘得到(K,)向量而X1.T是(K,)向量直接用会触发广播错误。正确的写法是# 正确的协变量平衡性约束L2范数最小化 constraints.append(cp.norm(X0.T w - X1.T, 2) 0.1) # 设置一个容忍阈值但实践中我很少启用它。原因在于协变量平衡性应该是一个“事前筛选”和“事后检验”的过程而不是一个“事中约束”。事前我们用X0和X1计算每个候选地区与目标地区的欧氏距离只保留距离最近的10-15个作为donor_regions事后我们用weights计算加权后的X0并与X1对比生成一张“平衡性检验表”。强行在优化中加入这个约束反而会让结果变量的拟合变差得不偿失。权重求解的终极目标是让合成组在“结果”上无限逼近目标组。其他一切都是为这个目标服务的辅助手段。3.3 结果解读不只是看“差值”更要读懂“置信区间”和“安慰剂检验”当你终于跑出了一张漂亮的合成控制效果图——目标地区实线和合成组虚线在干预前严丝合缝干预后拉开一道清晰的差距——别急着庆祝。真正的专业体现在对这张图背后每一个像素的审慎解读。核心效应量ATT这是最直观的数字。它等于干预后目标地区实际结果减去合成组预测结果的平均值。例如华东区Q1实际GMV是1.2亿合成组预测是1.05亿那么ATT就是1500万。但这个数字本身毫无意义除非你知道它的不确定性有多大。置信区间CI合成控制法的标准误不能用OLS那一套公式。主流做法是块自助法Block Bootstrap。原理很简单把干预前的时间序列切成若干个连续的块比如每3个月为一块然后有放回地随机抽取这些块拼成一个新的“伪历史序列”再在这个新序列上重新跑一遍权重求解和ATT计算。重复1000次你就得到了1000个ATT值取其2.5%和97.5%分位数就是95%置信区间。synthdid库内置了bootstrap_ci函数但如果你想手动实现关键代码是from sklearn.utils import resample def block_bootstrap_ci(outcome_df, weights, target_region, pre_period_start, pre_period_end, n_boot1000): pre_data outcome_df.loc[pre_period_start:pre_period_end] T len(pre_data) block_size 3 # 三个月为一块 n_blocks T // block_size atts [] for _ in range(n_boot): # 随机抽取n_blocks个块的索引 block_indices np.random.choice(n_blocks, sizen_blocks, replaceTrue) # 拼接成新的伪时间序列 boot_idx [] for bi in block_indices: start bi * block_size end min((bi 1) * block_size, T) boot_idx.extend(range(start, end)) boot_pre_data pre_data.iloc[boot_idx] # 在伪序列上重新计算合成组结果 synthetic_outcome (boot_pre_data[donor_regions] weights).values # ... 计算ATT ... atts.append(boot_att) return np.percentile(atts, [2.5, 97.5])这个过程很耗时但它是你向老板证明“这1500万不是噪音”的唯一武器。如果CI是[-500万, 3500万]那你的结论就是“效应不显著”。安慰剂检验Placebo Test这是合成控制法的“压力测试”。它的逻辑是如果我的方法是可靠的那么当我把“干预”随机安在某个从未被干预过的地区身上时它应该检测不到任何显著效应。具体操作从所有候选对照组中随机挑选一个作为“伪处理组”用同样的流程为它构建合成组并计算其“伪ATT”。重复500次画出这500个伪ATT的分布直方图。然后把你的真实ATT值标在这张图上。如果它远远甩开了95%的伪ATT比如落在了最右边的尾巴上那你的结论就非常有力。我做过一个经典案例真实ATT是8.2%而500个伪ATT中只有3个大于8.0%p-value0.006。这张图比任何文字描述都更有说服力。在因果推断的世界里一个漂亮的效应量只是入场券一个扎实的安慰剂检验才是你的毕业证书。4. 实操过程与核心环节实现从零开始复现一个完整的商业归因项目4.1 项目背景与数据模拟一个真实的“区域营销活动”归因场景让我们把所有理论放进一个具体的、有血有肉的商业场景里。假设你是一家全国性快消品公司的数据科学家公司于2023年7月1日在“华南大区”启动了一项名为“夏日冰爽节”的线下地推活动活动内容包括在核心商圈设置互动冰柜、发放限定版试饮券、联合本地KOC进行打卡挑战。活动持续了整个第三季度7-9月。业务方想知道这项活动为华南大区带来了多少额外的“冰镇饮料”品类销售额由于活动只在华南大区开展其他大区华北、华东、华中、西南、西北、东北是天然的候选对照组。我们手头有2021年1月至2023年12月共36个月的各区域月度销售额数据以及2020年末各区域的关键协变量人口、GDP、年轻人口占比、便利店密度、夏季平均气温。提示在真实项目中获取干净的、颗粒度一致的跨区域数据往往是第一步也是最难的一步。很多公司ERP系统里各区域的核算口径、促销折让规则都不一样。我的经验是与其花一个月去“统一口径”不如先用“最粗粒度但最一致”的数据跑通流程比如只用“总销售额”而不拆分到SKU。验证方法有效后再逐步精细化。4.2 代码全流程实现从数据加载到报告生成下面我将展示一个端到端的、可直接复制粘贴运行的Python脚本。它包含了所有关键环节并附有详细的中文注释。# -*- coding: utf-8 -*- 合成控制法实战华南大区“夏日冰爽节”活动归因分析 作者一位踩过所有坑的从业者 import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns import cvxpy as cp from sklearn.utils import resample import warnings warnings.filterwarnings(ignore) # ------------------- Step 1: 数据准备与模拟 ------------------- # 由于无法提供真实数据我们用np.random生成一个高度仿真的数据集 np.random.seed(42) dates pd.date_range(2021-01-01, 2023-12-01, freqMS) # 36个月 regions [华南大区, 华北大区, 华东大区, 华中大区, 西南大区, 西北大区, 东北大区] # 模拟基础趋势所有区域都有一个缓慢上升的线性趋势 季节性波动 t np.arange(len(dates)) base_trend 100 0.5 * t # 基础线性增长 seasonal 10 * np.sin(2 * np.pi * t / 12 np.pi/4) # 季节性7月峰值 # 为每个区域添加独特的“噪声”和“协变量效应” # 协变量2020年末快照 covariates_data { 华南大区: [12000, 15000, 35, 85, 28], # 人口(万), GDP(亿), 年轻人占比(%), 便利店密度(家/平方公里), 气温(℃) 华北大区: [11000, 14000, 28, 70, 25], 华东大区: [13000, 18000, 32, 95, 26], 华中大区: [9000, 11000, 30, 60, 27], 西南大区: [10000, 12000, 33, 55, 24], 西北大区: [8000, 9000, 25, 40, 22], 东北大区: [9500, 10000, 27, 50, 20], } covariates_df pd.DataFrame(covariates_data, index[population, gdp, youth_pct, convenience_density, avg_temp]).T # 生成各区域销售额单位百万 outcome_data {} for i, region in enumerate(regions): # 基础销售额 基础趋势 季节性 区域特有噪声 协变量效应 noise np.random.normal(0, 3, len(dates)) # 随机噪声 # 协变量效应用协变量的加权和模拟区域固有潜力 cov_effect (covariates_df.loc[region, population] * 0.01 covariates_df.loc[region, gdp] * 0.005 covariates_df.loc[region, youth_pct] * 0.5 covariates_df.loc[region, convenience_density] * 0.3 covariates_df.loc[region, avg_temp] * 0.8) # 为华南大区在2023年7-9月添加活动效应15% activity_effect np.zeros(len(dates)) activity_mask (dates 2023-07-01) (dates 2023-09-01) if region 华南大区: activity_effect[activity_mask] base_trend[activity_mask] * 0.15 sales base_trend seasonal noise cov_effect activity_effect outcome_data[region] sales outcome_df pd.DataFrame(outcome_data, indexdates) # 干预时间表 treatment_df pd.DataFrame({ region: [华南大区], treatment_start: [2023-07-01] }) print(✅ 数据准备完成) print(f数据时间范围{outcome_df.index.min()} 至 {outcome_df.index.max()}) print(f干预前月份数{(outcome_df.index 2023-07-01).sum()}) print(f干预后月份数{(outcome_df.index 2023-07-01).sum()}) # ------------------- Step 2: 权重求解 ------------------- def solve_weights_cvxpy(outcome_df, target_region, pre_period_start, pre_period_end, covariates_df, donor_regionsNone): 使用cvxpy求解合成权重 if donor_regions is None: donor_regions [r for r in outcome_df.columns if r ! target_region] pre_mask (outcome_df.index pre_period_start) (outcome_df.index pre_period_end) Y1_pre outcome_df.loc[pre_mask, target_region].values Y0_pre outcome_df.loc[pre_mask, donor_regions].values w cp.Variable(len(donor_regions)) objective cp.Minimize(cp.sum_squares(Y0_pre w - Y1_pre)) constraints [cp.sum(w) 1, w 0] prob cp.Problem(objective, constraints) prob.solve(solvercp.ECOS, verboseFalse) if prob.status not in [optimal, optimal_inaccurate]: raise ValueError(Optimization failed!) weights pd.Series(w.value, indexdonor_regions) return weights[weights 0.001].round(4) # 执行求解 weights solve_weights_cvxpy( outcome_dfoutcome_df, target_region华南大区, pre_period_start2021-01-01, pre_period_end2023-06-01, covariates_dfcovariates_df ) print(\n✅ 权重求解完成) print(合成组构成权重 0.001) print(weights.sort_values(ascendingFalse)) # ------------------- Step 3: 合成组构建与效应计算 ------------------- def build_synthetic_series(outcome_df, weights, target_region): 构建合成组的时间序列 donor_regions weights.index.tolist() # 合成组 加权求和 synthetic_series (outcome_df[donor_regions] weights).rename(Synthetic_South_China) # 真实组 real_series outcome_df[target_region].rename(Real_South_China) # 合并 result_df pd.concat([real_series, synthetic_series], axis1) return result_df result_df build_synthetic_series(outcome_df, weights, 华南大区) # 计算ATTAverage Treatment Effect on the Treated post_mask result_df.index 2023-07-01 att_point_est (result_df.loc[post_mask, Real_South_China] - result_df.loc[post_mask, Synthetic_South_China]).mean() print(f\n✅ 效应量计算完成) print(f平均处理效应 (ATT): {att_point_est:.2f} 百万元) # ------------------- Step 4: 可视化 ------------------- plt.figure(figsize(12, 6)) sns.lineplot(dataresult_df, xresult_df.index, yReal_South_China, label华南大区真实, linewidth2.5) sns.lineplot(dataresult_df, xresult_df.index, ySynthetic_South_China, label华南大区合成, linestyle--, linewidth2.5) plt.axvline(xpd.to_datetime(2023-07-01), colorred, linestyle:, alpha0.7, label活动开始) plt.title(合成控制法华南大区“夏日冰爽节”活动效果评估, fontsize14, fontweightbold) plt.ylabel(月度销售额百万元) plt.legend() plt.grid(True, alpha0.3) plt.tight_layout() plt.show() # ------------------- Step 5: 安慰剂检验 ------------------- def placebo_test(outcome_df, treatment_df, covariates_df, n_placebos500): 执行安慰剂检验 all_regions outcome_df.columns.tolist() # 排除目标区域只在对照组中进行安慰剂 placebo_regions [r for r in all_regions if r ! 华南大区] placebo_atts [] for _ in range(n_placebos): # 随机挑选一个伪处理组 pseudo_treat np.random.choice(placebo_regions) # 为它构建合成组 pseudo_weights solve_weights_cvxpy( outcome_df, pseudo_treat, 2021-01-01, 2023-06-01, covariates_df ) pseudo_result build_synthetic_series(outcome_df, pseudo_weights, pseudo_treat) # 计算其伪ATT post_mask_pseudo pseudo_result.index 2023-07-01 pseudo_att (pseudo_result.loc[post_mask_pseudo, fReal_{pseudo_treat}] - pseudo_result.loc[post_mask_pseudo, fSynthetic_{pseudo_treat}]).mean() placebo_atts.append(pseudo_att) return np.array(placebo_atts) print(\n⏳ 正在执行安慰剂检验500次...) placebo_atts placebo_test(outcome_df, treatment_df, covariates_df, n_placebos500) p_value np.mean(placebo_atts att_point_est) print(f✅ 安慰剂检验完成) print(f真实ATT: {att_point_est:.2f})