logo
搜索
我的图书馆

发文章

发文工具

其他工具

留言交流
搜索
分享

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

 张春强2022 2023-05-05 发布于安徽

目录

  • 1. 改变政策时间

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

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

  • 2. 虚构处理组

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

  • 4. 相关推文

安慰剂检验核心思想在于虚构处理组或者虚构政策时间进行估计，如果不同虚构方式下的估计量的回归结果依然显著，那么就说明说明原来的估计结果很有可能出现了偏误。现有的关于 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./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) 在研究国家级新区设立对经济增长的影响时，采用以下方法随机提前政策时点：

假定设立国家级新区的城市不变，如果现实中城市 在 年设立了国家级新区，那么从 [2003，t-1] 的时间范围内随机抽取任意 1 年作为城市 设立国家级新区的时间，据此利用新的样本重新估计基准模型便可以得到变量 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) 的做法：将处理组中设立国家级新区的城市视为新的控制组；在控制组中选取与处理组样本数相同的城市作为新的处理组，并且保持新处理组样本中城市的政策实施时间、批次与原处理组一致。

即，如果在 年有 个城市设立了国家级新区，则新控制组中应当有 个城市在 年实施政策，在此基础上利用新的样本重新估计基准模型，并重复上述操作 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 附近，绝大多数 值大于 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模型的求解与模拟
  • 专题：论文重现
    • 论文复现：基于组级纵向数据评估政策的试验模拟方法
图片

New！ Stata 搜索神器：lianxh 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

图片

    本站是提供个人知识管理的网络存储空间，所有内容均由用户发布，不代表本站观点。请注意甄别内容中的联系方式、诱导购买等信息，谨防诈骗。如发现有害或侵权内容，请点击一键举报。
    转藏 分享 献花（0）

    来自： 张春强2022 > 《共同趋势检验》

    举报/认领

    0条评论

    发表

    请遵守用户 评论公约

    类似文章