Stata实战：双重差分模型(DID)的完整检验流程与可视化呈现
1. 双重差分模型(DID)入门：从原理到应用场景
我第一次接触双重差分模型是在分析一项教育政策效果时，当时需要评估某省推行的"重点高中扩招"政策对学生升学率的影响。传统方法难以区分政策效果和其他因素，而DID完美解决了这个问题。简单来说，DID就像在实验室里做对照实验——我们比较政策实施前后，实验组和对照组的变化差异。
核心逻辑可以打个比方：假设你想知道新肥料对玉米产量的影响，可以选两块地（A用新肥料，B用旧肥料），比较它们产量变化的差值。在Stata中实现这个思路需要三个关键变量：
- 处理组虚拟变量（treated=1/0）
- 时间虚拟变量（post=1/0）
- 两者的交互项（did = treated × post）
最基础的DID模型公式如下：
reg outcome did treated post, vce(robust)这个简单模型却藏着几个容易踩坑的细节：
- 必须满足平行趋势假设——政策前实验组和对照组的变化趋势应该一致
- 要警惕溢出效应——对照组成员可能间接受到政策影响
- 时间虚拟变量设置要避开政策实施当年（避免混淆短期冲击和长期效果）
2. 数据预处理：为DID分析打好基础
2.1 数据清洗实战技巧
拿到原始数据时，我通常会先运行browse快速浏览数据结构。最近处理的一个县域经济数据集就发现日期格式混乱的问题——有的记录是"2015-12-31"，有的是"2015年12月"。这时需要用字符串函数统一处理：
replace date = subinstr(date,"2015-2015-12-31","2015",.) destring date, replace force异常值处理我偏好使用截尾法(winsor2)而非简单删除，特别是当某些极端值可能反映真实情况时。比如县域GDP数据中，某个县因新发现矿产资源导致GDP暴增，这本身就是研究价值所在：
ssc install winsor2 winsor2 gdp, cuts(1 99) replace // 缩尾处理 winsor2 gdp, cuts(1 99) trim replace // 截尾处理2.2 面板数据声明与检查
DID分析必须正确定义面板结构。有次我忘记声明面板数据，结果标准误计算完全错误。正确的操作流程应该是：
xtset id year // 声明面板 xtdescribe // 检查平衡性如果发现非平衡面板，可以用tsfill补全时间维度，或者用reghdfe直接处理非平衡数据。我常用的固定效应回归命令如下：
reghdfe y did, absorb(id year) vce(cluster id)3. 核心模型估计与结果解读
3.1 基础模型构建
在评估环保政策对企业排放的影响时，我构建了这样一个模型：
reghdfe pollution did size profit, absorb(province#year industry#year) vce(cluster firm)这里有几个关键点：
- 通过
province#year和industry#year控制了省份-年份和行业-年份的双重固定效应 vce(cluster firm)考虑了企业层面的聚类标准误- 加入企业规模(size)和利润(profit)作为控制变量
结果解读要特别注意系数大小和经济意义。有次我发现did系数显著为负，高兴地得出政策有效的结论，后来才发现是因为因变量单位是"万吨污染物"，应该取对数才对！
3.2 动态效应检验
静态DID只能看平均效应，而政策效果往往是渐进的。我推荐使用事件研究法：
forvalues i = 2/5 { gen event_`i' = (year == policy_year + `i') * treated } reghdfe y event_* treated, absorb(id year) vce(cluster id)然后用coefplot可视化动态路径：
coefplot, keep(event_*) vertical /// yline(0) xline(3) /// ytitle("处理效应") xtitle("相对政策年份") /// ciopts(recast(rcap))4. 稳健性检验：让结果经得起推敲
4.1 安慰剂检验的三种玩法
时间安慰剂是我最常用的方法——假装政策提前实施：
gen fake_year = year - 3 // 假设政策提前3年 gen fake_post = (fake_year >= policy_year) gen fake_did = treated * fake_post reghdfe y fake_did, absorb(id year) vce(cluster id)空间安慰剂也很有意思。有次我做省级政策评估，随机抽取部分县作为假定的"处理组"：
preserve keep if year == 2010 // 选政策前一年 sample 50, count // 随机抽50个县作为假处理组 save fake_treated.dta, replace restore merge m:1 county using fake_treated.dta, gen(fake_treated) reghdfe y fake_treated##post, absorb(id year) vce(cluster id)4.2 随机抽样检验
这个检验能告诉我们结果是否可能出于偶然。我通常做500次抽样：
permute did beta = _b[did] se = _se[did], /// reps(500) saving(simulations.dta): /// reghdfe y did, absorb(id year) vce(cluster id)然后看真实估计值在模拟分布中的位置：
use simulations.dta, clear sum beta local real_beta = 0.12 // 替换为你的真实估计值 count if abs(beta) > abs(`real_beta')5. 结果呈现：让图表自己讲故事
5.1 专业回归表格输出
我习惯用esttab输出三种格式的表格：
esttab using results.rtf, b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) /// scalars(N r2) mtitles("模型1" "模型2" "模型3") replace esttab using results.csv, csv b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) /// scalars(N r2) mtitles("模型1" "模型2" "模型3") replace esttab using results.tex, tex b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) /// scalars(N r2) mtitles("模型1" "模型2" "模型3") replace5.2 高级可视化技巧
系数图可以这样美化：
coefplot, baselevels keep(event_*) /// coeflabels(event_2 = "-2年" event_3 = "-1年" /// event_4 = "基准年" event_5 = "+1年") /// yline(0, lpattern(dash)) xline(3, lpattern(dash)) /// ytitle("政策效应") xtitle("相对时间") /// ciopts(color(black)) msymbol(D) mcolor(blue) /// graphregion(color(white)) plotregion(color(white))平行趋势检验图建议这样做：
egen time_mean = mean(y), by(year treated) twoway (connected time_mean year if treated==1) /// (connected time_mean year if treated==0), /// xline(2016) legend(label(1 "处理组") label(2 "对照组")) /// ytitle("结果变量均值") xtitle("年份")6. 实战中的疑难解答
6.1 处理效应异质性分析
有时候政策对不同类型的个体效果不同。我常用分样本回归：
reghdfe y did if size > median_size, absorb(id year) vce(cluster id) estimates store large_firm reghdfe y did if size <= median_size, absorb(id year) vce(cluster id) estimates store small_firm esttab large_firm small_firm, b(3) se(3) star(* 0.1 ** 0.05 *** 0.01)更优雅的做法是加入三重交互项：
gen large = (size > median_size) reghdfe y did##large, absorb(id year) vce(cluster id)6.2 多期DID处理
当政策在不同时间分批实施时，传统DID需要调整。我推荐使用Cengiz等(2019)的方法：
gen first_treat = year if treated == 1 & L.treated == 0 bysort id (year): replace first_treat = first_treat[1] if missing(first_treat) gen rel_year = year - first_treat然后生成事件虚拟变量：
forvalues k = -5/5 { if `k' >= -1 { gen event_`k' = (rel_year == `k') } } drop event_-1 // 以-2期为基准 reghdfe y event_*, absorb(id year) vce(cluster id)7. 完整案例：最低工资政策就业效应评估
最近我用DID分析了某省最低工资调整对服务业就业的影响，完整流程如下：
- 数据准备
import excel "data.xlsx", sheet("就业数据") firstrow case(lower) destring emp_rate min_wage, replace force winsor2 emp_rate, cuts(1 99) replace- 模型设定
gen post = (year >= 2017) // 政策在2017年实施 gen did = treated * post- 基准回归
reghdfe emp_rate did, absorb(city year) vce(cluster city)- 动态效应
forvalues y = 2015/2020 { gen y`y' = (year == `y') * treated } reghdfe emp_rate y2016 y2018 y2019 y2020, absorb(city year) vce(cluster city)- 稳健性检验
// 安慰剂检验 gen fake_year = year - 2 gen fake_post = (fake_year >= 2017) gen fake_did = treated * fake_post reghdfe emp_rate fake_did, absorb(city year) vce(cluster city) // 更换控制组 reghdfe emp_rate did if industry != "餐饮业", absorb(city year) vce(cluster city)- 结果输出
esttab using minwage.rtf, b(3) se(3) star(* 0.1 ** 0.05 *** 0.01) /// scalars(N r2) mtitles("基准" "动态" "安慰剂") replace这个案例中，我发现政策实施后第二年就业率显著下降约3.2%，但到第三年影响就消失了，说明企业可能通过调整用工结构适应了工资成本上涨。