双重差分法之安慰剂检验
简单介绍一下实证论文中双重差分法(DID)的安慰剂检验(Placebo Test)在Stata
中如何操作。
(本文首发于个人微信公众号DMETP
,是往期两篇推文的合辑,欢迎关注!)
下面的内容根据实际使用的数据集分为两个部分。
一是以一个截面数据集为例,介绍一下安慰剂检验的整个思路与流程。这里使用的是系统数据集
auto.dta
,由于是简单介绍思路,因此该部分并没有第二部分面板数据那么复杂,且模型中不包括DID的交互项,仅仅是对一个核心变量rep78
进行1,000次随机抽样;二是以一个面板数据集为例,介绍一下面板数据DID中安慰剂检验的整个流程。这里使用的数据集是石大千等(2018)公布在《中国工业经济》网站上的附件内容。论文使用的模型是普通DID模型,也即政策发生时点(2012年)固定,处理组与控制组的设置也固定。
参考文献:
石大千, 丁海, 卫平, 刘建江. 智慧城市建设能否降低环境污染[J]. 中国工业经济, 2018(06): 117-135.
一、安慰剂检验方法选择
通常而言,DID中的安慰剂检验方法包括两种。
一是改变政策发生时点,具体又包括前置处理组的政策发生时点,此时安慰剂检验的作用与平行趋势检验相同,都是考察政策发生前基础回归中时间虚拟变量与处理组交互项系数(
F(-1)
、F(-2)
、F(-3)
、…)的显著性,如果不显著说明检验通过;还包括将政策发生时点随机化,也即将时点前置或后置,这是一种更一般化的做法;二是将处理组随机化,对处理组变量进行一定次数的随机抽样,然后再观测随机化后的DID项系数或观测值的核密度图是否集中分布于0附近,以及是否显著偏离其真实值;
第二种方法更为常见,第一种方法的不足在于:如果样本期间较短,导致随机抽样的时段区间过短,得出的结论不一定真实,虽然抽样次数可以很多,但抽样空间过短将影响结论的稳健性。综合来说,对处理组进行随机化处理是一种更为合适的做法。当然,具体论文要具体分析。
二、截面数据集的安慰剂检验
这部分代码使用的是Stata
系统自带的数据集auto.dta
,该数据集是截面数据且不包含DID项,在实际使用中,可以将reg
改为面板数据回归命令(如xtreg
、reghdfe
等),同时将这里的核心解释变量rep78
改为论文中需要进行随机化处理的关键变量。
需要提醒的是,陆菁等(2021)在随机化处理时,提到“DID模型估计要求政策实施年份前后至少有一年数据”,因此在时间窗口不长且需要进行安慰剂检验的情况下,需要特别注意这一点。
参考文献:
陆菁, 鄢云, 王韬璇. 绿色信贷政策的微观效应研究——基于技术创新与资源再配置的视角[J]. 中国工业经济, 2021(01): 174-192.
此外,代码部分还参考了李青原和章尹赛楠(2021)发表在中国工业经济上的一篇文章的附件,同时还参考了简书的一篇文章。
参考文献:
李青原, 章尹赛楠. 金融开放与资源配置效率——来自外资银行进入中国的证据[J]. 中国工业经济, 2021(05): 95-113.
2.1 整体思路
第一步:在原始数据集
auto.dta
中单独剔除核心变量rep78
的样本数据;第二步:将剔除出来的
rep78
随机打乱顺序,再将随机化的rep78
合并至已被处理过的原始数据集中;第三步:将随机化的
rep78
放入回归方程中进行回归;第四步:以上操作步骤重复1,000次;
第五步:单独提取出1,000次回归结果中
rep78
的系数与标准误,最后分别绘制系数和t值的核密度分布图以及P值 - 系数散点图。
2.2 代码实现
*==============================================================================*
* 双重差分法 | 安慰剂检验 *
*==============================================================================*** Stata Version: 16 | 17** 【文章首发】这篇文章首发于本人微信公众号『DMETP』,是两篇推文的合辑,欢迎关注!cd "C:\Users\KEMOSABE\Desktop\placebo_test"**# 一、截面数据的安慰剂检验**# 1.1 根据原始样本进行基础回归sysuse auto.dta, clearglobal ctrlvar1 "mpg headroom trunk weight length"reg price rep78 $ctrlvar , r*- 基础回归中核心变量rep78的真实系数与真实t值分别为:
*- 真实系数 = 889.6715
*- 真实t值 = 3.2206(可手动计算,也即889.6715 / 276.2438)clear all**# 1.2 对rep78进行1,000次随机抽样、回归并绘制核密度图*- 思路:
*- a. 在原始数据集auto.dta中单独剔除核心变量rep78的样本数据
*- b. 将剔除出来的rep78随机打乱顺序,再将随机化的rep78合并至已被处理过的原始数据集中
*- c. 将随机化的rep78放入回归方程中进行回归
*- d. 以上操作步骤重复1,000次
*- e. 单独提取出1,000次回归结果中rep78的系数与标准误,最后分别绘制系数和t值的核密度估计图以及P值与系数的散点图set seed 13579 // 设置随机种子数forvalue i = 1/1000 {sysuse auto.dta, clearpreservegen randomvar = runiform()sort randomvargen id = _nlabel var id "用于匹配的样本序号"keep rep78 idsave rep78_random.dta, replace // 该数据集中仅保留rep78和随机化的id// id用于将该数据集中的rep78横向合并(merge)至原数据集restoregen id = _nlabel var id "用于匹配的样本序号"drop rep78save rep78_dropped.dta, replace // 原数据集中已剔除rep78字段use rep78_dropped.dta, clearmerge 1:1 id using rep78_random.dta, keepusing(rep78) // 将随机化排序的rep78合并至原数据集qui reg price rep78 $ctrlvar1 , rgen b_rep78 = _b[rep78] // 提取回归后rep78的系数gen se_rep78 = _se[rep78] // 提取回归后rep78的标准误keep b_rep78 se_rep78duplicates drop b_rep78, forcesave placebo_`i'.dta, replace}erase rep78_random.dtaerase rep78_dropped.dtause placebo_1.dta, clearforvalue k = 2/1000 {append using placebo_`k'.dta // 将1,000次回归中rep78的系数和标准误纵向合并(append)至单独的数据集(placebo_1.dta)erase placebo_`k'.dta}graph set window fontface "Times New Roman" graph set window fontfacesans "宋体" // 设置图形输出的字体gen tvalue = b_rep78 / se_rep78gen pvalue = 2 * ttail(e(df_r), abs(tvalue)) // 计算t值和P值**# 1.2.1 系数sum b_rep78, detailtwoway(kdensity b_rep78, ///xline(`r(mean)', lpattern(dash) lcolor(black)) ///xline(889.6715 , lpattern(solid) lcolor(black)) ///scheme(qleanmono) ///xtitle("{stSans:系数}" , size(medlarge)) ///ytitle("{stSans:核}""{stSans:密}""{stSans:度}", size(medlarge) orientation(h)) ///saving(placebo_test_Coefficient1, replace)), ///xlabel(, labsize(medlarge)) ///ylabel(, labsize(medlarge) format(%05.4f)) // 绘制1,000次回归rep78的系数的核密度图graph export "placebo_test_Coefficient1.png", replace // 导出为矢量图,方便论文中的图形展示;medlarge可以改为large以将标题文字加大**# 1.2.2 P值sum b_rep78, detailtwoway(scatter pvalue b_rep78, ///msy(oh) mcolor(black) ///xline(`r(mean)', lpattern(dash) lcolor(black)) ///xline(889.6715 , lpattern(solid) lcolor(black)) ///yline(0.1 , lpattern(shortdash) lcolor(black)) ///scheme(qleanmono) ///xtitle("{stSans:系数}" , size(medlarge)) ///ytitle("{stSans:P}""{stSans:值}" , size(medlarge) orientation(h)) ///saving(placebo_test_Pvalue1, replace)), ///xlabel( , labsize(medlarge)) ///ylabel(0(0.25)1, labsize(medlarge) format(%03.2f))graph export "placebo_test_Pvalue1.png", replace**# 1.2.3 t值sum tvalue, detailtwoway(kdensity tvalue, ///xline(`r(mean)', lpattern(dash) lcolor(black)) ///xline(3.2206 , lpattern(solid) lcolor(black)) ///xline(-1.65 , lpattern(shortdash) lcolor(black)) ///xline( 1.65 , lpattern(shortdash) lcolor(black)) ///scheme(qleanmono) ///xtitle("{stSans:t值}" , size(medlarge)) ///ytitle("{stSans:核}""{stSans:密}""{stSans:度}", size(medlarge) orientation(h)) ///saving(placebo_test_Tvalue1, replace)), ///xlabel(, labsize(medlarge)) ///ylabel(, labsize(medlarge) format(%02.1f)) // 绘制1,000次回归rep78的t值的核密度图graph export "placebo_test_Tvalue1.png", replaceerase placebo_1.dtadiscardclear allcls
2.3 运行结果与解读
以上这段代码的运行结果是3张图,如下所示。其中图 1
是系数的核密度估计图;图 2
是P值 - 系数散点图;图 3
是t值的核密度估计图。
针对图 1
至图 3
的解读如下:
随机化核心解释变量后系数与t值的核密度估计值的均值都接近于0(分别为1.2233和0.0029);
随机化后系数与t值的核密度估计值的均值都大大偏离其真实值(真实值分别为889.6715和3.2206);
随机化后多数系数的P值位于
P value = 0.1
线以上,说明多数系数至少在10%的水平下不显著;以上三点均说明
rep78
对price
的影响不是由其他不可观测因素(或遗漏变量)推动的;设置随机种子数为13,579时,可重复以上结果并得出一致结论;
从P值的散点图可以得到以下两点信息:
- 第一,更多的散点集中分布于0附近,而位于真实值垂直线上的散点只有几个,这说明在随机化后真实值是一个异常点;
- 第二,虽然多数散点集中于0附近,但这些散点所代表的系数至少在10%的水平上是不显著的。
三、面板数据集的安慰剂检验
前面一部分介绍了安慰剂检验的具体操作,但都是以一个截面数据集(auto.dta
)作为示例的,且模型中没有加入DID的交互项,因此严格来说这个例子还不太恰当。这里用一个具体例子介绍面板数据双重差分模型中的安慰剂检验,这个例子是一个普通DID模型,政策发生时点固定,处理组和控制组也是固定的,相对而言模型设置比较简单,但也可以延伸至相对复杂的DID模型中(如多期DID、连续DID和广义DID等),所需的可能仅仅是要发挥更多的想象力。
原始数据来源于石大千(2018),但这篇文章和中国工业经济官网放出来的附件都没有详细解释处理组和控制组的具体设置,因此虽然可以用gen dt = (year >= 2012)
生成政策发生时间虚拟变量(dt
),但处理组虚拟变量(du
)无法生成,因此这里使用的是公众号『功夫计量经济学』提供的数据,数据有效性本人无法保证,这里只作为一个参考示例。
关于论文的具体内容详见参考文献,这里不做介绍,也可以快速浏览『功夫计量经济学』的相关推文。值得一提的是,『功夫计量经济学』给出了另外一种随机抽样的方法,可以与本推送给出的方案对照阅读。
这里设置了一个随机种子(seed
),方便复现结果与推送内容保持一致,随机种子数是223,至于为什么是这个数,纯粹是试错试出来的,因为设置成这个数画出来的图最好看。
3.1 整体思路
第一步:在原始数据集
smart_city2018.dta
中单独剔除变量id
的样本数据;第二步:将剔除出来的
id
随机打乱顺序,再将随机化的id
合并至已被处理过的原始数据集中;第三步:将随机化的
treat
与dt
的交互项(did
)放入回归方程中进行回归;第四步:以上操作步骤重复1,000次;
第五步:单独提取出1,000次回归结果后
did
的系数与标准误,最后分别绘制系数和t值的核密度估计图以及P值与系数的散点图。
3.2 代码实现
*==============================================================================*
* 双重差分法 | 安慰剂检验 *
*==============================================================================*** Stata Version: 16 | 17** 【文章首发】这篇文章首发于本人微信公众号『DMETP』,是两篇推文的合辑,欢迎关注!**# 二、面板数据的安慰剂检验*- 【原始数据来源】石大千, 丁海, 卫平, 刘建江. 智慧城市建设能否降低环境污染[J]. 中国工业经济, 2018(06): 117-135.
* 《中国工业经济》附件下载地址:http://ciejournal.ajcass.org/Magazine/show/?id=54281*- 【加工数据来源】微信公众号『功夫计量经济学』2021.4.20推送文章“双重差分法(DID)安慰剂检验的做法:随机抽取500次”
* 微信公众号推文地址:https://mp.weixin.qq.com/s/06v6s90G1pp-yLju_yAy1Qcd "C:\Users\KEMOSABE\Desktop\placebo_test\example"**# 2.1 基础回归use smart_city2018.dta, clearglobal ctrlvars2 "lnrgdp lninno lnurb lnopen lnss"reghdfe lnrso dudt $ctrlvars2 , absorb(id year) cluster(id)*- 基础回归中核心变量dudt的真实系数与真实t值分别为:
*- 真实系数 = -0.1706
*- 真实t值 = -2.1245(可手动计算,即-0.1706 / 0.0803)discard**# 2.2 安慰剂检验*- 思路:
*- a. 在原始数据集smart_city2018.dta中单独剔除变量id的样本数据
*- b. 将剔除出来的id随机打乱顺序,再将随机化的id合并至已被处理过的原始数据集中
*- c. 将随机化的treat与dt的交互项(did)放入回归方程中进行回归
*- d. 以上操作步骤重复1,000次
*- e. 单独提取出1,000次回归结果后did的系数与标准误,最后分别绘制did系数和t值的核密度估计图以及P值与系数的散点图set seed 223 // 设置随机种子数forvalue i = 1/1000 {use smart_city2018.dta, clearpreservekeep if year == 2005gen randomvar = runiform()sort randomvarkeep if id in 1/32keep idsave id_random.dta, replace // 数据集仅保留随机化后的idrestoremerge m:1 id using id_random.dta // 将随机化id合并至原数据集gen treat = (_merge == 3)gen did = treat * dtqui reghdfe lnrso did $ctrlvars2 , absorb(id year) cluster(id)gen b_did = _b[did] // 提取单次回归后did的系数gen se_did = _se[did] // 提取单次回归后did的标准误keep b_did se_didduplicates drop b_did, forcesave placebo_`i'.dta, replace}erase id_random.dtause placebo_1.dta, clearforvalue k = 2/1000 {append using placebo_`k'.dta // 将1,000次回归后did的系数和标准误纵向合并(append)至单独的数据集(placebo_1.dta)erase placebo_`k'.dta}graph set window fontface "Times New Roman"graph set window fontfacesans "宋体" // 设置图形输出的字体gen tvalue = b_did / se_didgen pvalue = 2 * ttail(e(df_r), abs(tvalue)) // 计算t值和P值**# 2.2.1 系数sum b_did, detailtwoway(kdensity b_did, ///xline(`r(mean)', lpattern(dash) lcolor(black)) ///xline(-0.1706 , lpattern(solid) lcolor(black)) ///scheme(qleanmono) ///xtitle("{stSans:系数}" , size(medlarge)) ///ytitle("{stSans:核}""{stSans:密}""{stSans:度}", size(medlarge) orientation(h)) ///saving(placebo_test_Coefficient2, replace)), ///xlabel(, labsize(medlarge) format(%02.1f)) ///ylabel(, labsize(medlarge) format(%02.1f)) // 绘制1,000次回归did的系数的核密度图graph export "placebo_test_Coefficient2.png", replace // 导出为矢量图,方便论文中的图形展示;可改medlarge为large以加大字体**# 2.2.2 P值sum b_did, detailtwoway(scatter pvalue b_did, ///msy(oh) mcolor(black) ///xline(`r(mean)', lpattern(dash) lcolor(black)) ///xline(-0.1706 , lpattern(solid) lcolor(black)) ///yline( 0.1 , lpattern(shortdash) lcolor(black)) ///scheme(qleanmono) ///xtitle("{stSans:系数}" , size(medlarge)) ///ytitle("{stSans:P}""{stSans:值}" , size(medlarge) orientation(h)) ///saving(placebo_test_Pvalue2, replace)), ///xlabel( , labsize(medlarge) format(%02.1f)) ///ylabel(0(0.25)1, labsize(medlarge) format(%03.2f))graph export "placebo_test_Pvalue2.png", replace**# 2.2.3 t值sum tvalue, detailtwoway(kdensity tvalue, ///xline(`r(mean)', lpattern(dash) lcolor(black)) ///xline(-2.1245 , lpattern(solid) lcolor(black)) ///xline(-1.65 , lpattern(shortdash) lcolor(black)) ///xline( 1.65 , lpattern(shortdash) lcolor(black)) ///scheme(qleanmono) ///xtitle("{stSans:t值}" , size(medlarge)) ///ytitle("{stSans:核}""{stSans:密}""{stSans:度}", size(medlarge) orientation(h)) ///saving(placebo_test_Tvalue2, replace)), ///xlabel(, labsize(medlarge)) ///ylabel(, labsize(medlarge) format(%02.1f)) // 绘制1,000次回归did的t值的核密度图graph export "placebo_test_Tvalue2.png", replaceerase placebo_1.dta
3.3 运行结果与解读
同样,以上代码的运行结果是3张图,如下图 4
至图 6
。
针对以上3张图,有如下几点解读。
第一,
图 4
是随机化处理组后did
项回归系数的核密度估计图,其中实线是基础回归估计出来的真实系数,虚线是1,000个“虚拟”系数的均值;第二,
图 5
是t值的核密度估计图,其中实线是真实t值,虚线是均值,两根短虚线分别代表t = -1.65
和t = 1.65
(也即大样本下10%的显著性水平所对应的t值,t值的绝对值小于该数说明至少在10%的水平下不显著);第三,
图 6
是P值的散点图,其中水平短虚线是P = 0.1
,散点位于该虚线以下说明系数至少在10%的水平下显著,反之则不显著;第四,
图 4
和图 5
都说明了一个基本事实,那就是绝大多数系数和t值均集中分布在0附近,均值与真实值的距离较远,且绝大多数估计系数并不显著,这意味着智慧城市试点对环境污染治理的政策效应没有受到其他未被观测因素的影响。这个基本事实其实完全可以从P值的散点图(图 6
)中得知,如散点集中分布在0附近,且远离其真实值,多数散点都位于虚线以上,同时说明在10%的水平下不显著,也就是说,P值散点图包含的信息其实更多更凝练。
双重差分法之安慰剂检验相关推荐
- 双重差分法|DID|PSM|平行趋势检验|安慰剂检验|Stata代码
双重差分与稳健性检验原理 双重差分法 以两期面板为例 数据采用cardkrueger1994.dta,可在陈强老师个人主页<高级计量经济学及Stata应用>, 第2版下载数据集 基本模型 ...
- 双重差分模型能做固定效应吗_数据分析之道 | 双重差分法(DID)
Picture from Internet DID是什么? 双重差分法(DID)又被称为"倍差法",小名"差中差",是种专门用于分析政策效果的计量方法. 我国最 ...
- python包怎么做双重差分did分析_【研究方法】双重差分法(DID)介绍
双重差分法,英文名Differences-in-Differences,别名"倍差法",小名"差中差".作为政策效应评估方法中的一大利器,双重差分法受到越来越多 ...
- 因果推断 | 双重差分法笔记补充
换了新的环境后,一直在适应(其实是一直被推着走),所以停更了笔记好久啦.这一周周末终于有点得空,当然也是因为疫情,哪里都不能去,哈哈,所以来冒个泡~ 整理了最近pre的作业,分享一下双重差分法的一些笔 ...
- 双重差分法(DID):标准化流程和stata代码实现
文章目录 标准化流程 平行假设检验 效果评估 安慰剂检验 标准化流程 此前的文章介绍了双重差分法(difference-in-differences,DID)的原理,并说明了其是算法策略效果评估的有效 ...
- python 双重差分_双重差分法(DID)介绍
双重差分法,英文名Differences-in-Differences,别名"倍差法",小名"差中差".作为政策效应评估方法中的一大利器,双重差分法受到越来越多 ...
- did双重差分法_互助问答第252期:双重差分平行趋势检验等问题
老师您好,关于双重差分法的适用性问题:基于现有样本数据,做出平行趋势图如下,设想将A组作为控制组,B组作为实验组,采用双重差分法检验政策出台对B组有抑制效应,请问老师基于这张图是否可行? 2. 关于三 ...
- did双重差分法_Stata中双重差分操流程及代码
01 简介 现代计量经济学和统计学的发展为我们的研究提供了可行的工具.倍差法来源于计量经济学的综列数据模型,是政策分析和工程评估中广为使用的一种计量经济方法.主要是应用于在混合截面数据集中,评价某一事 ...
- did双重差分法_互助问答第47期:政策时点不一致DID的问题
问题: 各位老师好,在政策效果评估时经常会用到双重差分法,我想请问一个关于政策时点不一致DID的问题.我知道DID需要满足两个前提,一是政策外生性,二是平行趋势假设.但是我国绝大多数政策其实是采取&q ...
最新文章
- 我用 YOLOv5 做情感识别!
- 一步一步学习hadoop(七)
- @这位没带口罩的朋友,你让我感染新冠的风险升高百倍!马普所建模计算结果,认真的...
- 【错误记录】Tinker 热修复示例运行报错 ( Execution failed for task ‘:app:tinkerProcessD‘ . tinkerId is not set!!! )
- “面试不败计划”:多线程
- php如何做浏览量,php+ajax实现的点击浏览量加1
- leetcode237 删除链表中的节点(你意想不到的做法,注意细节)
- 第一阶段 03Java的基本数据类型
- java第七章jdbc课后简答题_java学习路线流程
- 中英文搜索引擎收录口整理
- 未来的5年内,我为什么不看好“AI+教育”
- dev-mysql_GitHub - intergrate-dev/mysql-elasticsearch
- 微软 HoloLens 2 的幕后故事
- mysql获取当前时间,前一天,后一天
- Hue添加Spark notebook
- 嵌入式相关开源项目、库、资料
- juk互粉攻略set结构体
- phalapi可以依赖注入么_3.2 PhalApi 配置
- vue面试常见问题小结
- CocosCreator开发笔记(4)-Windows搭建幼麟麻将运行环境
热门文章
- AJAX实验(添加+模糊查询 图书)
- vue2中h(“router-view“) vue3如何写?
- 【Golang 中的 type A = XXX 与 type A XXXX的区别】
- ubuntu将cuda卸载干净
- 线性代数基础10--特征值与特征向量,行列式的空间关系
- android环信群聊显名称,Android环信群聊插入头像和昵称
- 重载和重写的区别???
- 如何防止破解?MCU加密技术揭秘
- 前端报表导出成word文档(含echarts图表)
- Java Pair的使用