多时点DID保姆级教程(下)-安慰剂检验

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
  . ssc install lianxh
详情参见帮助文件 (有惊喜)：
  . help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg
 

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

• 全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
• 计量专题 | 论文写作 | 数据分享 | 专题课程
• 结果输出 | Stata绘图 | 数据处理 | Stata程序
• 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
• 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
• Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
• 文本分析-爬虫 | Python-R-Matlab | 机器学习
• Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：李闯 (新疆大学)
邮箱：chuangli0808@163.com

目录

安慰剂检验核心思想在于虚构处理组或者虚构政策时间进行估计，如果不同虚构方式下的估计量的回归结果依然显著，那么就说明说明原来的估计结果很有可能出现了偏误。现有的关于 DID 的安慰剂检验主要有以下几种方式：

1. 改变政策时间

改变政策发生时点，具体包括政策时点提前、滞后两种情况，还包括将政策时间随机化这种更一般化的处理方式 (多用于政策时点随机提前)。

政策时点提前，若回归系数不显著，则说明是该政策起作用，反向验证了结论的稳健性。政策时点滞后，系数依旧显著，一般而言绝对值增大，表明政策效果明显，并且具有一定的持续性。

1.1 政策时点提前、滞后一定期数

数据来自石大千等 (2018)《智慧城市建设能否降低环境污染》这篇文章，严格意义来说本文只使用了 2012 年的一期智慧城市试点，不属于多期DID，但是方法相同。

use smart_city, clear
merge m:1 id using smart_pt //匹配政策变量，生成政策时间变量
drop _merge
replace pt=0 if pt==.
g did_1=0
replace did_1=1 if pt==2012&year>=2011 //政策时点提前1期
g did_2=0
replace did_2=1 if pt==2012&year>=2010 //政策时点提前2期
g did1=0
replace did1=1 if pt==2012&year>=2013 //政策时点滞后1期
g did2=0
replace did2=1 if pt==2012&year>=2014 //政策时点滞后2期
xtreg lnrso did_1 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
xtreg lnrso did_2 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
xtreg lnrso did1 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
xtreg lnrso did2 lnrgdp lntgdp lninno lnurb lnopen lnss,fe //以人均废气排放量为例

数据表明在假定政策时点提前 1 期、2 期后，智慧城市政策对人均废气排放量均无显著抑制作用，从侧面证实了智慧试点政策的环境友好效应。同时在政策分别滞后 1 期、2 期时，系数依旧显著为负，且绝对值增大，表明政策的促进环境友好发展作用在加强。

. lxhuse smart_city, clear
. merge m:1 id using https://file.lianxh.cn/data/s/smart_pt.dta // 匹配政策变量，生成政策时间变量
. drop _merge
. replace pt=0 if pt==.
. g did_1=0
. replace did_1=1 if pt==2012&year>=2011 // 政策时点提前1期
. g did_2=0
. replace did_2=1 if pt==2012&year>=2010 // 政策时点提前2期
. g did1=0
. replace did1=1 if pt==2012&year>=2013  // 政策时点滞后1期
. g did2=0
. replace did2=1 if pt==2012&year>=2014  // 政策时点滞后2期

. xtset id year
. g lntgdp=lnrgdp*lnrgdp
. gen lnww=log(ww)
. gen lnso=log(so)
. g lnrww=log(ww/pop)
. xtreg lnrso did_1 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
. est  store m1
. xtreg lnrso did_2 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
. est  store m2
. xtreg lnrso did1 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
. est  store m3
. xtreg lnrso did2 lnrgdp lntgdp lninno lnurb lnopen lnss,fe
. est  store m4

. esttab m1 m2 m3 m4, star(* 0.1 ** 0.05 *** 0.01) b(%6.3f) ///
>     order(did_1 did_2 did1 did2) compress nogap

--------------------------------------------------------------
                 (1)          (2)          (3)          (4)   
               lnrso        lnrso        lnrso        lnrso   
--------------------------------------------------------------
did_1         -0.051                                          
             (-0.91)                                          
did_2                      -0.091                             
                          (-1.63)                             
did1                                    -0.114*               
                                       (-1.86)                
did2                                                 -0.128*  
                                                    (-1.84)   
lnrgdp         2.695***     2.621***     2.600***     2.619***
              (5.79)       (5.66)       (5.63)       (5.69)   
lntgdp        -0.137***    -0.133***    -0.132***    -0.133***
             (-5.96)      (-5.82)      (-5.79)      (-5.86)   
lninno        -0.037*      -0.037*      -0.037*      -0.038** 
             (-1.90)      (-1.90)      (-1.89)      (-1.97)   
lnurb          0.172***     0.172***     0.174***     0.174***
              (5.06)       (5.06)       (5.11)       (5.11)   
lnopen         0.030        0.030        0.030        0.031   
              (1.23)       (1.25)       (1.27)       (1.28)   
lnss           0.033        0.033        0.026        0.027   
              (0.24)       (0.25)       (0.19)       (0.20)   
_cons        -14.771***   -14.430***   -14.326***   -14.415***
             (-6.14)      (-6.03)      (-6.00)      (-6.06)   
--------------------------------------------------------------
N               2167         2167         2167         2167   
--------------------------------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01

1.2 处理组政策时点随机提前

曹清峰 (2020) 在研究国家级新区设立对经济增长的影响时，采用以下方法随机提前政策时点：

假定设立国家级新区的城市不变，如果现实中城市 $i$ 在 $t$ 年设立了国家级新区，那么从 [2003，t-1] 的时间范围内随机抽取任意 1 年作为城市 $i$ 设立国家级新区的时间，据此利用新的样本重新估计基准模型便可以得到变量 did 的估计系数。将上述过程重复 1000 次，得到其概率分布图如下：

mat b = J(1000,1,0)  // 系数矩阵
mat se = J(1000,1,0) // 标准误矩阵 
mat p = J(1000,1,0)  // P值矩阵

* 循环1000次
lxhget cao_treat.dta, replace 
lxhget cao_control.dta, replace
forvalues i=1/1000{
    use cao_treat.dta, clear
    qui xtset id year  
    gen DID1=0
    replace DID1=1 if pt==2005&year<2005
    replace DID1=1 if pt==2009&year<2009
    replace DID1=1 if pt==2012&year<2012
    replace DID1=1 if pt==2013&year<2013
    replace DID1=1 if pt==2014&year<2014
    replace DID1=1 if pt==2015&year<2015
    replace DID1=1 if pt==2016&year<2016   // 选取政策前年份
    keep if DID1==1
    sample 1, count by(id) // 根据id分组，每组随机抽取一个年份 
    keep id year           // 得到所抽取样本的id编号
    rename year policy_year
    save match_id_year.dta, replace   // 另存id编号数据
    merge 1:m id using cao_treat.dta  // 与原数据匹配
    qui xtset id year
    gen treatment = (_merge == 3)      // 生成政策分组虚拟变量
    gen period = (year >= policy_year) //生成政策时间虚拟变量
    gen did0 = treatment*period
    append using cao_control.dta
    xtreg gdpr did0 invest consume export gov second agg innov i.year, fe cluster(id)
    * 将回归结果赋值到对应矩阵的对应位置
    mat b[`i',1] = _b[did]   // 系数矩阵
    mat se[`i',1] = _se[did] // 标准误矩阵
    mat p[`i',1] = 2*ttail(e(df_r), abs(_b[did]/_se[did])) // 计算P值并赋值于矩阵
}
* 矩阵转化为向量
svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)
* 删除空值并添加标签
drop if pvalue1 == .
label var pvalue1 p值
label var coef1 估计系数
keep coef1 se1 pvalue1  
gen tvalue = coef1/se1     // 计算t值
save placebo.dta, replace   //关于p值，估计系数的文件，要用作画图
kdensity coef1, normal scheme(qleanmono) //绘制核密度分布图    

由图可知，系数基本满足正态分布，而且平均值远小于真实系数估计值，即随机提前国家级新区的设立时间会导致国家新区对所在城市经济增长的带动效应出现明显下降， 这也从反事实角度证实了国家级新区设立后确实提高了所在城市的经济增长率。

2. 虚构处理组

同样，借鉴曹清峰 (2020) 的做法：将处理组中设立国家级新区的城市视为新的控制组；在控制组中选取与处理组样本数相同的城市作为新的处理组，并且保持新处理组样本中城市的政策实施时间、批次与原处理组一致。

即，如果在 $t$ 年有 $n$ 个城市设立了国家级新区，则新控制组中应当有 $n$ 个城市在 $t$ 年实施政策，在此基础上利用新的样本重新估计基准模型，并重复上述操作 1000 次，便可以完成 1 次安慰剂检验。

mat b = J(1000,1,0)  // 系数矩阵
mat se = J(1000,1,0) // 标准误矩阵 
mat p = J(1000,1,0)  // P值矩阵
lxhget new_area.dta, replace
use new_area, clear 
xtset id year
forvalues i=1/100 {
    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 1, count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2005  // 将新随机样本政策时间设为2005年
    save random_control_2005,replace

    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 1,count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2009  // 将新随机样本政策时间设为2009年
    save random_control_2009,replace

    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 2,count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2012  // 将新随机样本政策时间设为2012年
    save random_control_2012,replace

    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 2,count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2013  // 将新随机样本政策时间设为2013年
    save random_control_2013,replace

    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 4,count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2014  // 将新随机样本政策时间设为2014年
    save random_control_2014,replace

    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 5,count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2015  // 将新随机样本政策时间设为2015年
    save random_control_2015,replace

    use new_area, clear 
    drop if treat==1 // 删除原来的处理组
    drop treat
    gen random=1
    sample 2,count 
    gen treat=1 // 将新随机样本treat设为1
    gen t=2016  // 将新随机样本政策时间设为2016年
    save random_control_2016,replace

    use random_control_2005,clear
    append using random_control_2009
    append using random_control_2012
    append using random_control_2013
    append using random_control_2014
    append using random_control_2015
    append using random_control_2016
    duplicates drop id year,force
    save random_control_all.dta,replace
    use new_area, clear 
    drop treat
    gen random=0
    gen start=0
    append using random_control_all
    by id year,sort: gen test=_N
    drop if test==2 & random==0
    by id,sort: egen t2=max(t)
    by id,sort: egen treat2=max(treat)
    gen DID0=1 if year>=t2&treat2==1
    replace DID0=0 if DID0==.
    xtreg gdpr DID0 invest consume export gov second agg innov i.year, fe cluster(id)
    * 将回归结果赋值到对应矩阵的对应位置
    mat b[`i',1] = _b[DID0]   // 系数矩阵
    mat se[`i',1] = _se[DID0] // 标准误矩阵
    * 计算P值并赋值于矩阵
    mat p[`i',1] = 2*ttail(e(df_r), abs(_b[DID0]/_se[DID0]))
}
* 矩阵转化为向量
svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)
* 删除空值并添加标签
drop if pvalue1 == .
label var pvalue1 p值
label var coef1 估计系数
keep coef1 se1 pvalue1  
gen tvalue = coef1/se1     // 计算t值
save placebo.dta,replace   // 关于p值，估计系数的文件，要用作画图
sum coef1  // 计算系数均值
kdensity coef1,normal scheme(qleanmono)

估计结果显示，此次变量 did 系数的均值为 -0.34，远小于基准回归结果 1.51，这表明国家级新区的政策效应表现出了明显的区位导向性，对设立国家级新区城市经济增长的带动效应最显著。

3. 同时随机化政策时间与处理组样本

以白俊红等 (2022) 为例，从所有样本中随机选取与原处理组相同数量的样本，并随机生成政策实施时间，构建城市与政策时间均随机的新处理组。基于此，重新估计基准回归模型，并随机重复 500 次实验，操作 Stata 代码如下：

clear
set matsize 5000
mat b = J(500,1,0)
mat se = J(500,1,0)
mat p = J(500,1,0)
lxhget inno_policy.dta, replace
forvalues i=1/500{
	use "inno_policy.dta" , clear 
	xtset city year
	keep if year==2004
	sample 71, count
	keep city
	save "atchcity.dta",replace
	merge 1:m city using "inno_policy.dta"
	gen groupnew=(_merge==3) // 生成伪处理组的虚拟变量
	save "matchcity`i'.dta",replace
		
	* 伪政策虚拟变量
	use "inno_policy.dta",clear 
	bsample 1, strata(city)
	keep year
	save "matchyear.dta", replace
	mkmat year, matrix(sampleyear)
	use "matchcity`i'.dta",replace
	xtset city year
	gen time = 0
	foreach j of numlist 1/280 {
		replace time = 1 if (city == `j' & year >= sampleyear[`j',1])
	}	
	gen  did=time*groupnew
	global xlist  "lnagdp indust_stru finance ainternet market "	
	xtreg entre_activation did  $xlist  i.year, fe robust
	
	mat b[`i',1] = _b[did]
	mat se[`i',1] = _se[did]
	scalar df_r = e(N) - e(df_m) -1
	mat p[`i',1] = 2*ttail(df_r,abs(_b[did]/_se[did]))
}
svmat b, names(coef)
svmat se, names(se)
svmat p, names(pvalue)

drop if pvalue1 == .
label var pvalue1 p值
label var coef1 估计系数
save placebo.dta,replace  // 保存数据

* 删除临时文件
forvalue i=1/500{
    erase  "matchcity`i'.dta"
}
graph set window fontface "Times New Roman"
graph set window fontfacesans "宋体"  // 设置图形输出的字体
woway (scatter pvalue1 coef1, xlabel(-0.2(0.05)0.4, format(%03.2f) angle(0)) ///
    yline(0.1,lp(shortdash)) xline(0.2997,lp(shortdash))                     ///
    xtitle("{stSans:估计系数}") ytitle("{stSans:P}" "{stSans:值}",           ///
    orientation(h)) msymbol(smcircle_hollow) mcolor(grey) legend(off))       ///
    (kdensity coef1, title("{stSans:安慰剂检验}"))

据下图可知，随机化后的 DID 项系数集中分布于 0 附近，绝大多数 $p$ 值大于 0.1，并且随机系数基本位于真实值 0.2997 左侧，表明双重随机处理后，政策效果在显著性与作用强度方面均有大幅削弱，间接证实了原文结论的稳健性。

4. 相关推文

Note：产生如下推文列表的 Stata 命令为：
  lianxh 安慰剂 模拟, m
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：Stata教程

  • Stata博文系列：边模拟边学计量经济学-Econometrics-By-Simulation
• 专题：数据处理

  • Stata：蒙特卡洛模拟新命令-simulate2
• 专题：Stata程序

  • Stata：蒙特卡洛模拟A-(Monte-Carlo-Simulation)没那么神秘
  • Stata程序：Monte-Carlo-模拟之产生符合特定分布的随机数
  • Stata：蒙特卡洛模拟分析 (Monte Carlo Simulation)
• 专题：回归分析

  • Stata：一行代码实现安慰剂检验-permute
  • aoeplacebo：地理安慰剂检验
  • Stata：在线可视化模拟-OLS-的性质
• 专题：倍分法DID

  • 多期DID之安慰剂检验、平行趋势检验
• 专题：内生性-因果推断

  • 安慰剂检验！安慰剂检验！
  • B-Stata模拟：控制变量！控制变量！Good-Controls-Bad-Controls
• 专题：时间序列

  • Stata: VAR - 模拟、估计和推断
• 专题：Python-R-Matlab

  • MATLAB：Ramsey模型的求解与模拟
• 专题：其它

  • 世代交叠OLG模型的求解与模拟
• 专题：论文重现

  • 论文复现：基于组级纵向数据评估政策的试验模拟方法

相关课程

免费公开课

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版
• Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

最新课程-直播课

• Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

⛳ 课程主页

⛳ 课程主页

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 连享会-主页 和 知乎专栏，400+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会-常见问题解答：
✨ https://gitee.com/lianxh/Course/wikis

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh