*-------------8.2个体安慰剂检验：随机化政策时间和处理组样本---------*/       
*生成备用矩阵
clear
mat b=J(500,1,0)   
mat se=J(500,1,0)
mat p=J(500,1,0)
forvalues i = 1/500{
*第二步，设置为伪处理组
use input_stata2.dta,clear
xtset id date
keep if date==13  
*13随便选一年   
sample 12,count
*12为实验组ID数量      
keep id
save matchcity.dta,replace
merge 1:m id using "input_stata2.dta"
gen treat1=(_merge==3)
save matchcity`i'.dta,replace
*第三步，设置伪虚拟变量
use input_stata2.dta,clear
bsample 1,strata(id)
keep date
save matchyear.dta,replace
mkmat date,matrix(sampleyear)
*第四步，设置为政策的伪虚拟变量
use matchcity`i'.dta,replace
xtset id date
gen time1=0
foreach j of numlist 1/18{
*18为所有ID数量
replace time1=1 if (id==`j' & date>=sampleyear[`j',1])
}
gen DID1=time1*treat1
*第五步，回归
global X1 " lntrade lniifa lninstruction lncapacity cpi rpi"
qui xtreg lnfreight DID1 $X1 i.date,fe r
mat b[`i',1]= _b[DID1]
mat se[`i',1]= _se[DID1]
scalar df_r=e(N)-e(df_m)-1
mat p[`i',1]=2*ttail(df_r,abs(_b[DID1]/_se[DID1]))
}
*绘图
svmat b,names(coef)
svmat se,names(se)
svmat p,names(pvalue)

drop if pvalue==.
label var pvalue p值
label var coef1 估计系数

twoway(scatter pvalue1 coef1,xlabel(-0.15(0.05)0.15,grid) yline(0.1,lp(shortdash)) xline(-0.0193,lp(shortdash)) xtitle(估计系数) ytitle(p值) msymbol(smcircle_hollow) mcolor(blue) legend(off))(kdensity coef1,title(安慰剂检验))
*-0.0193为DID回归的结果
forvalue i=1/500{
erase matchcity`i'.dta
}
*安慰剂结果解释见https://www.bilibili.com/video/av292757788
* 竖的虚线不在范围之内；大多数在横虚线之上；关于0对称