簡單介紹一下實證論文中雙重差分法（DID）的安慰劑檢驗（Placebo Test）在Stata中如何操作。
（本文首發(fā)于個人微信公眾號DMETP，是往期兩篇推文的合輯幅聘，歡迎關(guān)注Ｆ计簟）
下面的內(nèi)容根據(jù)實際使用的數(shù)據(jù)集分為兩個部分。
一是以一個截面數(shù)據(jù)集為例蚤假，介紹一下安慰劑檢驗的整個思路與流程。這里使用的是系統(tǒng)數(shù)據(jù)集
auto.dta吧兔，由于是簡單介紹思路磷仰，因此該部分并沒有第二部分面板數(shù)據(jù)那么復(fù)雜，且模型中不包括DID的交互項境蔼，僅僅是對一個核心變量rep78進(jìn)行1,000次隨機(jī)抽樣灶平；二是以一個面板數(shù)據(jù)集為例，介紹一下面板數(shù)據(jù)DID中安慰劑檢驗的整個流程欧穴。這里使用的數(shù)據(jù)集是石大千等（2018）公布在《中國工業(yè)經(jīng)濟(jì)》網(wǎng)站上的附件內(nèi)容民逼。論文使用的模型是普通DID模型，也即政策發(fā)生時點（2012年）固定涮帘，處理組與控制組的設(shè)置也固定拼苍。
參考文獻(xiàn)：
石大千, 丁海, 衛(wèi)平, 劉建江. 智慧城市建設(shè)能否降低環(huán)境污染[J]. 中國工業(yè)經(jīng)濟(jì), 2018(06): 117-135.
一、安慰劑檢驗方法選擇
通常而言，DID中的安慰劑檢驗方法包括兩種疮鲫。
一是改變政策發(fā)生時點吆你，具體又包括前置處理組的政策發(fā)生時點，此時安慰劑檢驗的作用與平行趨勢檢驗相同俊犯，都是考察政策發(fā)生前基礎(chǔ)回歸中時間虛擬變量與處理組交互項系數(shù)（
F(-1)妇多、F(-2)、F(-3)燕侠、......）的顯著性者祖，如果不顯著說明檢驗通過；還包括將政策發(fā)生時點隨機(jī)化绢彤，也即將時點前置或后置七问，這是一種更一般化的做法；二是將處理組隨機(jī)化茫舶，對處理組變量進(jìn)行一定次數(shù)的隨機(jī)抽樣械巡，然后再觀測隨機(jī)化后的DID項系數(shù)或觀測值的核密度圖是否集中分布于0附近，以及是否顯著偏離其真實值饶氏；
第二種方法更為常見讥耗，第一種方法的不足在于：如果樣本期間較短，導(dǎo)致隨機(jī)抽樣的時段區(qū)間過短疹启，得出的結(jié)論不一定真實古程，雖然抽樣次數(shù)可以很多，但抽樣空間過短將影響結(jié)論的穩(wěn)健性皮仁。綜合來說籍琳，對處理組進(jìn)行隨機(jī)化處理是一種更為合適的做法菲宴。當(dāng)然贷祈，具體論文要具體分析。
二喝峦、截面數(shù)據(jù)集的安慰劑檢驗
這部分代碼使用的是Stata系統(tǒng)自帶的數(shù)據(jù)集auto.dta势誊，該數(shù)據(jù)集是截面數(shù)據(jù)且不包含DID項，在實際使用中谣蠢，可以將reg改為面板數(shù)據(jù)回歸命令（如xtreg粟耻、reghdfe等），同時將這里的核心解釋變量rep78改為論文中需要進(jìn)行隨機(jī)化處理的關(guān)鍵變量眉踱。
需要提醒的是挤忙，陸菁等（2021）在隨機(jī)化處理時，提到“DID模型估計要求政策實施年份前后至少有一年數(shù)據(jù)”谈喳，因此在時間窗口不長且需要進(jìn)行安慰劑檢驗的情況下册烈，需要特別注意這一點。
參考文獻(xiàn)：
陸菁, 鄢云, 王韜璇. 綠色信貸政策的微觀效應(yīng)研究——基于技術(shù)創(chuàng)新與資源再配置的視角[J]. 中國工業(yè)經(jīng)濟(jì), 2021(01): 174-192.
此外婿禽，代碼部分還參考了李青原和章尹賽楠（2021）發(fā)表在中國工業(yè)經(jīng)濟(jì)上的一篇文章的附件赏僧，同時還參考了簡書的一篇文章大猛。
參考文獻(xiàn)：
李青原, 章尹賽楠. 金融開放與資源配置效率——來自外資銀行進(jìn)入中國的證據(jù)[J]. 中國工業(yè)經(jīng)濟(jì), 2021(05): 95-113.
2.1 整體思路
第一步：在原始數(shù)據(jù)集
auto.dta中單獨(dú)剔除核心變量rep78的樣本數(shù)據(jù)；第二步：將剔除出來的
rep78隨機(jī)打亂順序淀零，再將隨機(jī)化的rep78合并至已被處理過的原始數(shù)據(jù)集中挽绩；第三步：將隨機(jī)化的
rep78放入回歸方程中進(jìn)行回歸；第四步：以上操作步驟重復(fù)1,000次驾中；
第五步：單獨(dú)提取出1,000次回歸結(jié)果中
rep78的系數(shù)與標(biāo)準(zhǔn)誤唉堪，最后分別繪制系數(shù)和t值的核密度分布圖以及P值 - 系數(shù)散點圖。
2.2 代碼實現(xiàn)
*==============================================================================*
* 雙重差分法 | 安慰劑檢驗 *
*==============================================================================*
** Stata Version: 16 | 17
** 【文章首發(fā)】這篇文章首發(fā)于本人微信公眾號『DMETP』,是兩篇推文的合輯肩民，歡迎關(guān)注巨坊！
cd "C:\Users\KEMOSABE\Desktop\placebo_test"
**# 一、截面數(shù)據(jù)的安慰劑檢驗
**# 1.1 根據(jù)原始樣本進(jìn)行基礎(chǔ)回歸
sysuse auto.dta, clear
global ctrlvar1 "mpg headroom trunk weight length"
reg price rep78 $ctrlvar , r
*- 基礎(chǔ)回歸中核心變量rep78的真實系數(shù)與真實t值分別為：
*- 真實系數(shù) = 889.6715
*- 真實t值 = 3.2206（可手動計算此改，也即889.6715 / 276.2438）
clear all
**# 1.2 對rep78進(jìn)行1,000次隨機(jī)抽樣趾撵、回歸并繪制核密度圖
*- 思路：
*- a. 在原始數(shù)據(jù)集auto.dta中單獨(dú)剔除核心變量rep78的樣本數(shù)據(jù)
*- b. 將剔除出來的rep78隨機(jī)打亂順序，再將隨機(jī)化的rep78合并至已被處理過的原始數(shù)據(jù)集中
*- c. 將隨機(jī)化的rep78放入回歸方程中進(jìn)行回歸
*- d. 以上操作步驟重復(fù)1,000次
*- e. 單獨(dú)提取出1,000次回歸結(jié)果中rep78的系數(shù)與標(biāo)準(zhǔn)誤共啃，最后分別繪制系數(shù)和t值的核密度估計圖以及P值與系數(shù)的散點圖
set seed 13579 // 設(shè)置隨機(jī)種子數(shù)
forvalue i = 1/1000 {
sysuse auto.dta, clear
preserve
gen randomvar = runiform()
sort randomvar
gen id = _n
label var id "用于匹配的樣本序號"
keep rep78 id
save rep78_random.dta, replace // 該數(shù)據(jù)集中僅保留rep78和隨機(jī)化的id
// id用于將該數(shù)據(jù)集中的rep78橫向合并（merge）至原數(shù)據(jù)集
restore
gen id = _n
label var id "用于匹配的樣本序號"
drop rep78
save rep78_dropped.dta, replace // 原數(shù)據(jù)集中已剔除rep78字段
use rep78_dropped.dta, clear
merge 1:1 id using rep78_random.dta, keepusing(rep78) // 將隨機(jī)化排序的rep78合并至原數(shù)據(jù)集
qui reg price rep78 $ctrlvar1 , r
gen b_rep78 = _b[rep78] // 提取回歸后rep78的系數(shù)
gen se_rep78 = _se[rep78] // 提取回歸后rep78的標(biāo)準(zhǔn)誤
keep b_rep78 se_rep78
duplicates drop b_rep78, force
save placebo_`i'.dta, replace
}
erase rep78_random.dta
erase rep78_dropped.dta
use placebo_1.dta, clear
forvalue k = 2/1000 {
append using placebo_`k'.dta // 將1,000次回歸中rep78的系數(shù)和標(biāo)準(zhǔn)誤縱向合并（append）至單獨(dú)的數(shù)據(jù)集（placebo_1.dta）
erase placebo_`k'.dta
}
graph set window fontface "Times New Roman"
graph set window fontfacesans "宋體" // 設(shè)置圖形輸出的字體
gen tvalue = b_rep78 / se_rep78
gen pvalue = 2 * ttail(e(df_r), abs(tvalue)) // 計算t值和P值
**# 1.2.1 系數(shù)
sum b_rep78, detail
twoway(kdensity b_rep78, ///
xline(`r(mean)', lpattern(dash) lcolor(black)) ///
xline(889.6715 , lpattern(solid) lcolor(black)) ///
scheme(qleanmono) ///
xtitle("{stSans:系數(shù)}" , 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的系數(shù)的核密度圖
graph export "placebo_test_Coefficient1.png", replace // 導(dǎo)出為矢量圖占调，方便論文中的圖形展示；medlarge可以改為large以將標(biāo)題文字加大
**# 1.2.2 P值
sum b_rep78, detail
twoway(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:系數(shù)}" , 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, detail
twoway(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", replace
erase placebo_1.dta
discard
clear all
cls
2.3 運(yùn)行結(jié)果與解讀
以上這段代碼的運(yùn)行結(jié)果是3張圖移剪，如下所示究珊。其中圖 1是系數(shù)的核密度估計圖；圖 2是P值 - 系數(shù)散點圖纵苛；圖 3是t值的核密度估計圖剿涮。
針對圖 1至圖 3的解讀如下：
隨機(jī)化核心解釋變量后系數(shù)與t值的核密度估計值的均值都接近于0（分別為1.2233和0.0029）；
隨機(jī)化后系數(shù)與t值的核密度估計值的均值都大大偏離其真實值（真實值分別為889.6715和3.2206）攻人；
隨機(jī)化后多數(shù)系數(shù)的P值位于
P value = 0.1線以上取试，說明多數(shù)系數(shù)至少在10%的水平下不顯著；以上三點均說明
rep78對price的影響不是由其他不可觀測因素（或遺漏變量）推動的怀吻；設(shè)置隨機(jī)種子數(shù)為13,579時瞬浓，可重復(fù)以上結(jié)果并得出一致結(jié)論；
-
從P值的散點圖可以得到以下兩點信息：
- 第一蓬坡，更多的散點集中分布于0附近猿棉，而位于真實值垂直線上的散點只有幾個，這說明在隨機(jī)化后真實值是一個異常點屑咳；
- 第二萨赁，雖然多數(shù)散點集中于0附近，但這些散點所代表的系數(shù)至少在10%的水平上是不顯著的兆龙。
三杖爽、面板數(shù)據(jù)集的安慰劑檢驗
前面一部分介紹了安慰劑檢驗的具體操作，但都是以一個截面數(shù)據(jù)集（auto.dta）作為示例的，且模型中沒有加入DID的交互項掂林，因此嚴(yán)格來說這個例子還不太恰當(dāng)臣缀。這里用一個具體例子介紹面板數(shù)據(jù)雙重差分模型中的安慰劑檢驗，這個例子是一個普通DID模型泻帮，政策發(fā)生時點固定精置，處理組和控制組也是固定的，相對而言模型設(shè)置比較簡單锣杂，但也可以延伸至相對復(fù)雜的DID模型中（如多期DID脂倦、連續(xù)DID和廣義DID等），所需的可能僅僅是要發(fā)揮更多的想象力元莫。
原始數(shù)據(jù)來源于石大千(2018)赖阻，但這篇文章和中國工業(yè)經(jīng)濟(jì)官網(wǎng)放出來的附件都沒有詳細(xì)解釋處理組和控制組的具體設(shè)置，因此雖然可以用gen dt = (year >= 2012)生成政策發(fā)生時間虛擬變量（dt）踱蠢，但處理組虛擬變量(du)無法生成火欧，因此這里使用的是公眾號『功夫計量經(jīng)濟(jì)學(xué)』提供的數(shù)據(jù)，數(shù)據(jù)有效性本人無法保證茎截，這里只作為一個參考示例苇侵。
關(guān)于論文的具體內(nèi)容詳見參考文獻(xiàn)，這里不做介紹企锌，也可以快速瀏覽『功夫計量經(jīng)濟(jì)學(xué)』的相關(guān)推文榆浓。值得一提的是，『功夫計量經(jīng)濟(jì)學(xué)』給出了另外一種隨機(jī)抽樣的方法撕攒，可以與本推送給出的方案對照閱讀陡鹃。
這里設(shè)置了一個隨機(jī)種子（seed），方便復(fù)現(xiàn)結(jié)果與推送內(nèi)容保持一致抖坪，隨機(jī)種子數(shù)是223萍鲸，至于為什么是這個數(shù)，純粹是試錯試出來的柳击，因為設(shè)置成這個數(shù)畫出來的圖最好看猿推。
3.1 整體思路
第一步：在原始數(shù)據(jù)集
smart_city2018.dta中單獨(dú)剔除變量id的樣本數(shù)據(jù)片习；第二步：將剔除出來的
id隨機(jī)打亂順序捌肴，再將隨機(jī)化的id合并至已被處理過的原始數(shù)據(jù)集中；第三步：將隨機(jī)化的
treat與dt的交互項（did）放入回歸方程中進(jìn)行回歸藕咏；第四步：以上操作步驟重復(fù)1,000次状知；
第五步：單獨(dú)提取出1,000次回歸結(jié)果后
did的系數(shù)與標(biāo)準(zhǔn)誤，最后分別繪制系數(shù)和t值的核密度估計圖以及P值與系數(shù)的散點圖孽查。
3.2 代碼實現(xiàn)
*==============================================================================*
* 雙重差分法 | 安慰劑檢驗 *
*==============================================================================*
** Stata Version: 16 | 17
** 【文章首發(fā)】這篇文章首發(fā)于本人微信公眾號『DMETP』,是兩篇推文的合輯饥悴，歡迎關(guān)注！
**# 二、面板數(shù)據(jù)的安慰劑檢驗
*- 【原始數(shù)據(jù)來源】石大千, 丁海, 衛(wèi)平, 劉建江. 智慧城市建設(shè)能否降低環(huán)境污染[J]. 中國工業(yè)經(jīng)濟(jì), 2018(06): 117-135.
* 《中國工業(yè)經(jīng)濟(jì)》附件下載地址：http://ciejournal.ajcass.org/Magazine/show/?id=54281
*- 【加工數(shù)據(jù)來源】微信公眾號『功夫計量經(jīng)濟(jì)學(xué)』2021.4.20推送文章“雙重差分法（DID）安慰劑檢驗的做法：隨機(jī)抽取500次”
* 微信公眾號推文地址：https://mp.weixin.qq.com/s/06v6s90G1pp-yLju_yAy1Q
cd "C:\Users\KEMOSABE\Desktop\placebo_test\example"
**# 2.1 基礎(chǔ)回歸
use smart_city2018.dta, clear
global ctrlvars2 "lnrgdp lninno lnurb lnopen lnss"
reghdfe lnrso dudt $ctrlvars2 , absorb(id year) cluster(id)
*- 基礎(chǔ)回歸中核心變量dudt的真實系數(shù)與真實t值分別為：
*- 真實系數(shù) = -0.1706
*- 真實t值 = -2.1245（可手動計算西设，即-0.1706 / 0.0803）
discard
**# 2.2 安慰劑檢驗
*- 思路：
*- a. 在原始數(shù)據(jù)集smart_city2018.dta中單獨(dú)剔除變量id的樣本數(shù)據(jù)
*- b. 將剔除出來的id隨機(jī)打亂順序瓣铣，再將隨機(jī)化的id合并至已被處理過的原始數(shù)據(jù)集中
*- c. 將隨機(jī)化的treat與dt的交互項（did）放入回歸方程中進(jìn)行回歸
*- d. 以上操作步驟重復(fù)1,000次
*- e. 單獨(dú)提取出1,000次回歸結(jié)果后did的系數(shù)與標(biāo)準(zhǔn)誤，最后分別繪制did系數(shù)和t值的核密度估計圖以及P值與系數(shù)的散點圖
set seed 223 // 設(shè)置隨機(jī)種子數(shù)
forvalue i = 1/1000 {
use smart_city2018.dta, clear
preserve
keep if year == 2005
gen randomvar = runiform()
sort randomvar
keep if id in 1/32
keep id
save id_random.dta, replace // 數(shù)據(jù)集僅保留隨機(jī)化后的id
restore
merge m:1 id using id_random.dta // 將隨機(jī)化id合并至原數(shù)據(jù)集
gen treat = (_merge == 3)
gen did = treat * dt
qui reghdfe lnrso did $ctrlvars2 , absorb(id year) cluster(id)
gen b_did = _b[did] // 提取單次回歸后did的系數(shù)
gen se_did = _se[did] // 提取單次回歸后did的標(biāo)準(zhǔn)誤
keep b_did se_did
duplicates drop b_did, force
save placebo_`i'.dta, replace
}
erase id_random.dta
use placebo_1.dta, clear
forvalue k = 2/1000 {
append using placebo_`k'.dta // 將1,000次回歸后did的系數(shù)和標(biāo)準(zhǔn)誤縱向合并（append）至單獨(dú)的數(shù)據(jù)集（placebo_1.dta）
erase placebo_`k'.dta
}
graph set window fontface "Times New Roman"
graph set window fontfacesans "宋體" // 設(shè)置圖形輸出的字體
gen tvalue = b_did / se_did
gen pvalue = 2 * ttail(e(df_r), abs(tvalue)) // 計算t值和P值
**# 2.2.1 系數(shù)
sum b_did, detail
twoway(kdensity b_did, ///
xline(`r(mean)', lpattern(dash) lcolor(black)) ///
xline(-0.1706 , lpattern(solid) lcolor(black)) ///
scheme(qleanmono) ///
xtitle("{stSans:系數(shù)}" , 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的系數(shù)的核密度圖
graph export "placebo_test_Coefficient2.png", replace // 導(dǎo)出為矢量圖贷揽，方便論文中的圖形展示棠笑；可改medlarge為large以加大字體
**# 2.2.2 P值
sum b_did, detail
twoway(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:系數(shù)}" , 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, detail
twoway(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", replace
erase placebo_1.dta
3.3 運(yùn)行結(jié)果與解讀
同樣，以上代碼的運(yùn)行結(jié)果是3張圖禽绪，如下圖 4至圖 6蓖救。
針對以上3張圖，有如下幾點解讀印屁。
第一循捺，
圖 4是隨機(jī)化處理組后did項回歸系數(shù)的核密度估計圖，其中實線是基礎(chǔ)回歸估計出來的真實系數(shù)雄人，虛線是1,000個“虛擬”系數(shù)的均值从橘；第二，
圖 5是t值的核密度估計圖础钠，其中實線是真實t值洋满，虛線是均值，兩根短虛線分別代表t = -1.65和t = 1.65（也即大樣本下10%的顯著性水平所對應(yīng)的t值珍坊，t值的絕對值小于該數(shù)說明至少在10%的水平下不顯著）牺勾；第三，
圖 6是P值的散點圖阵漏，其中水平短虛線是P = 0.1驻民，散點位于該虛線以下說明系數(shù)至少在10%的水平下顯著，反之則不顯著履怯；第四回还，
圖 4和圖 5都說明了一個基本事實，那就是絕大多數(shù)系數(shù)和t值均集中分布在0附近叹洲，均值與真實值的距離較遠(yuǎn)柠硕，且絕大多數(shù)估計系數(shù)并不顯著，這意味著智慧城市試點對環(huán)境污染治理的政策效應(yīng)沒有受到其他未被觀測因素的影響运提。這個基本事實其實完全可以從P值的散點圖（圖 6）中得知蝗柔，如散點集中分布在0附近，且遠(yuǎn)離其真實值民泵，多數(shù)散點都位于虛線以上癣丧，同時說明在10%的水平下不顯著，也就是說栈妆，P值散點圖包含的信息其實更多更凝練胁编。
