hahahoby 发表于 2020-7-1 07:55 
赵老师 您好，最近在用合成控制法做研究 有一个不太会，在安慰剂检验中不少人用等间距方法进行随机抽样作 ...
画安慰剂图，都需要通过编程的方式实现，下面我提供一个我写的程序供大家参考
*======================================
*合成控制法假设检验(Placebo test and draw graph)*
*======================================
set more off
use smoking,clear
tsset state year
* 获取一些后文要用的参数
qui tab state
local n = r(r) // 州数
qui tab year
local n_year = r(r) // 年份数
*======================================
*根据个人研究需要，调整这些参数
*======================================
local date_t = "1989" // 干预时间点
local m = 2 // 限制MSPE为干预州MSPE的m倍，m=0表示无限制
*local slow = "nested" // 取消*使用nested选项，计算量大，拟合更好
local id_t=3 // 干预州的id或行号
local treat_name ="California" // 图中显示的干预组名称
local ctrl_name="Control States" // 图中显示的控制组名称
local xtitle "year" // 横轴变量名称
local ytitle "gap in per-capita cigarette sales (in packs)" //纵轴变量名称
local saving "syn_plot" //保存安慰剂检验图
*======================================
tempname resmat
forvalues i=1/`n' {
synth cigsale beer lnincome retprice age15to24 cigsale(1988) cigsale(1980) cigsale(1975) , ///
trunit(`i') trperiod(`date_t') xperiod(1980(1)1988) `slow' keep(tmp`i', replace)
//上述循环命令分别对所有州作为干预组进行合成, tmp`I'.dta保存合成结果
local rmspe = e(RMSPE)[1,1] //取RMSPE
use tmp`i',clear
keep _Y_treated _Y_synthetic _time
gen te = _Y_treated- _Y_synthetic
gen id = `i'
keep in 1/`n_year' //1970-2000, there are 31 years, which is keep in the first 31 obs.
gen te2 = te*te // use it to calculate MSPE
local n_before = `date_t' - _time[1] //取干预期之前对应位置或序号
local n_after = `n_before' + 1 //干预期起点
qui sum te2 in 1/`n_before' // MSPE
local mspe_pre = r(mean) // 干预前的MSPE
qui sum te2 in `n_after'/`n_year'
local mspe_post = r(mean) // 干预后的MSPE
local r = `mspe_post'/`mspe_pre' //计算Abadie-R统计量
matrix `resmat' = nullmat(`resmat')\(`rmspe', `mspe_pre', `mspe_post', `r') //resmat saves the RMSPE for each model
local names `"`names'`"`i'"'"' // names of each
save tmp`i', replace
use smoking,clear
tsset state year
}
mat colnames `resmat' = "RMSPE" "MSPE_pre" "MSPE_post" "Abadie_R"
mat rownames `resmat' = `names'
matlist `resmat', row("Treated Unit")
*Placebo Graphs - Draw Figure 3
*Get the RMSPE of the treated unit
local RMSPE_t=`resmat'[`id_t',1]
use tmp1, clear
local num = 0 // # of units includes in the graph
forvalues i=2/`n' {
if `m'==0 {
append using tmp`i'
local num = `num' + 1
}
else if `resmat'[`i',1]^2<=`m'*`RMSPE_t'^2 { // MSPE comparation
append using tmp`i'
local num = `num' + 1
}
}
*======================================
*画安慰剂图1
local s="" // string to store the graph command
local controls = "" //string to store the id of control units used
local num_t = `num'+1 // # postion to identify the treated unit
levelsof id, local(levels)
foreach l of local levels {
if `l'!=`id_t' {
local s = "`s'"+"(line te _time if id==`l', lc(gs13))"
local controls = "`controls'"+" "+"`l'"
}
}
local date_before = `date_t'-1
two `s'(line te _time if id==`id_t', lc(black)), ///
legend(order(`num_t' "`treat_name'" `num' "`ctrl_name'") cols(1) pos(11) ring(0)) xline(`date_before', lp(dot) lc(black)) yline(0, lp(dash) lc(black)) ///
xlabel(1970(5)2000) xtitle("`xtitle'") ytitle("`ytitle'") saving(`saving'_`m', replace)
di "# of controls after limit `m' times of RMSPE of treated unit: " `num' //显示保留的控制组数量
di "ID of controls:" "`controls'" //显示保留的控制组id或序号
*======================================
*画出Abadie-R统计量分布图，Abadie et al. (2010)
*======================================
clear
svmat `resmat', names(col)
save tmp_R, replace //unstar this line if you want to save the file
histogram Abadie_R, freq width(1) text(1 77 "California {&rarr}", placement(s)) xtitle("post/pre-Proposition 99 mean squared prediction error")
*======================================
*删除所有临时文件
!del tmp*
set more on
exit
发帖
雷达卡
扫码
京公网安备 11010802022788号
论坛法律顾问：王进律师
知识产权保护声明
免责及隐私声明
