计量：事件研究法 2

clear

*-生成双重差分面板随机数据

set seed 2024 
set obs 100 
gen i = _n 
gen treated = i > 50 
gen E_i = 11 if treated == 1 
expand 20 
sort i 
qui by i: gen t = _n 
xtset i t 


*-设置处理效应

gen D = (t == E_i) 
gen etime = (t - E_i) 
label var i "个体编号"
label var t "日期"
label var D "did变量=time*post"
label var E_i "事件冲击时间"


*-画出实验组和处理组的分布

replace E_i=16 if treated==0
replace E_i=5 if treated==1 & i>=51 & i<=61
replace E_i=6 if treated==1 & i>=62 & i<=68
replace E_i=7 if treated==1 & i>=69 & i<=70
replace E_i=10 if treated==1 & i>=71 & i<=75
replace E_i=11 if treated==1 & i>=76 & i<=90
replace E_i=12 if treated==1 & i>=91


#delimit ; 
/*
分界符设置，修改分号；为命令分行符号。
如此以来，除非末尾加上分号，即便分行写也看作一排
画图命令这样使用会方面很多
*/
twoway  (histogram E_i if treated==1, discrete frequency start(4) barwidth(0.8) xscale(range(0 17)) yscale(range(0 1200))   color(cranberry%100)) 
(histogram E_i if treated==0, discrete frequency start(4) barwidth(0.8) xscale(range(0 17)) yscale(range(0 1200))   color(navy%40)),  
title("实验组和对照组事件冲击分布图", 
color(black) size(v.small)) legend(off)
xlabel(  5 "5" 6 "6" 7 "7" 8 "8" 9 "9" 10 "10" 11 "11" 12 "12" 16 "NA")
text(1100 15.5 "从来没接受干预", size(small))
name(histogram4 , replace);

#delimit cr

*-画出实验组和处理组的分布
preserve //保存当前数据集

#delimit ;
cdfplot E_i,  
    title("事件时间的累计分布图",
color(black) size(v.small)) legend(off)
ytitle(百分比)
name(CDF4, replace)
note(Hybrid事件数据类型，实验组事件时间不同，同时对照组从未受到事件冲击)
;
#delimit cr
restore //重新加载保存数据集


*-实验组和处理组的分布
replace treated=1
replace E_i=5 if treated==1 & i>=1 & i<=20
replace E_i=6 if treated==1 & i>=21 & i<=32
replace E_i=7 if treated==1 & i>=33 & i<=37
replace E_i=10 if treated==1 & i>=38 & i<=58
replace E_i=11 if treated==1 & i>=59 & i<=64
replace E_i=12 if treated==1 & i>=65 & i<=90
replace E_i=15 if treated==1 & i>90
#delimit ;
twoway  (histogram E_i if treated==1, discrete frequency start(4) barwidth(0.8) xscale(range(0 17)) yscale(range(0 800))   color(cranberry%100)),  
    title("实验组和对照组事件冲击分布图", 
    color(black) size(v.small)) legend(off)
ylabel(0(100)800)
name(histogram5 , replace)
;
#delimit cr

*-画出实验组和处理组的累计分布
#delimit ;
cdfplot E_i,   
    title("事件时间的累计分布图",
color(black) size(v.small)) legend(off)
ytitle(百分比)
note(Timing-based事件数据类型，由于实验组事件时间不同，内部也能形成比较)
name(CDF5, replace)
    ;

#delimit cr
*-合并展示以上图
graph combine histogram4 CDF4 histogram5 CDF5 
iscale 0.25 //图片比例缩小为四分之一

#delimit ;

cap log close ;

set trace off ;
set more off ;
clear all ;
set seed 101 ;

cap prog drop runme ; 
prog def runme ;
args list ctrl trends ;

qui gendata "`list'" `ctrl' `trends' ;

qui forvalues i = 0/20 { ;
gen D_p`i' = (etime == `i') ;  // 生成处理组正向虚拟变量
gen D_m`i' = (etime == -1 * `i') ;  // 生成处理组负向虚拟变量
summ D_p`i' ;  // 汇总处理组正向虚拟变量
local dropme = (r(mean) == 0) ;
if `dropme' == 1 { ;
drop D_p`i'  ;  // 如果均值为0，则删除处理组正向虚拟变量
} ;
summ D_m`i' ;  // 汇总处理组负向虚拟变量
local dropme = (r(mean) == 0) ;
if `dropme' == 1 { ;
drop D_m`i'  ;  // 如果均值为0，则删除处理组负向虚拟变量
} ;
} ;

drop D_m0  ;  // 删除处理组负向虚拟变量的基准情况

egen unittype = group(E_i) , miss ;  // 根据 E_i 生成单位类型变量，处理缺失值
qui tab unittype , gen(udum) ;  // 汇总单位类型变量，生成对应的虚拟变量
qui tab t , gen(tdum) ;  // 汇总时间变量，生成对应的虚拟变量

unab vars : D_* udum* tdum* ;  // 生成所有与 D_*、udum* 和 tdum* 相关的变量列表
mata: count() ;  // 调用 Mata 函数计算变量的统计信息
local c1 = count[1,1] ;  // 获取统计结果中的第一列信息
local c2 = count[1,2] ;  // 获取统计结果中的第二列信息
local c3 = count[1,3] ;  // 获取统计结果中的第三列信息

di "`list'" _column(20) `c1' _column(25) `c2'  
_column(35) `c3' ;  // 显示变量列表及其统计信息

end ;

mata: ;

void count() 
{  ;
v = st_local("vars") ;  // 获取 Stata 定义的局部宏变量 "vars"
X = st_data(.,v) ;  // 从当前数据集中获取变量数据

mycols = cols(X) ;  // 计算数据矩阵的列数
myrank = rank(X) ;  // 计算数据矩阵每列的非缺失值数量
gap = mycols - myrank ;  // 计算数据矩阵每列的缺失值数量
output = (mycols,myrank,gap) ;  // 将统计结果组合成元组
st_matrix("count",output) ;  // 将统计结果存储到名为 "count" 的矩阵中
} ;
end ;

cap prog drop gendata ;
program define gendata ;
args list ctrl trends ;

drop _all ;

local numobs = 48 ;  // 定义观察值数量
local treatedunittypes = wordcount("`list'") ;  // 计算处理单元类型数量

/* "ctrl" is 0/1 whether we have control (untreated) units */
/* if we have them, put half of data into control units */
local numtreatedobs = round(`numobs' / (1 + `ctrl')) ;  // 计算处理单元观察值数量
local numuntreatedobs = `numobs' - `numtreatedobs' ;  // 计算对照单元观察值数量

local obspertype = floor(`numtreatedobs' /
`treatedunittypes') ;  // 计算每种处理单元类型的观察值数量

forvalues i = 1/`treatedunittypes' { ;
local eventdate_`i' = word("`list'",`i') ;  // 获取每种处理单元类型的事件日期
} ;

set obs `numobs' ;/* 观察值数量 */
gen i = _n ;

gen treated = i > `numuntreatedobs' ;/* 处理组单元 */

gen E_i = . ;
forvalues i = 1/`treatedunittypes' { ;
local start = (`i'-1)*`obspertype' + `numuntreatedobs' + 1 ;
local stop = (`i')*`obspertype' + `numuntreatedobs' ;
replace E_i = `eventdate_`i'' if i >= `start' & i <= `stop' ;  // 根据每种处理单元类型设置事件日期
} ;
replace E_i = `eventdate_`treatedunittypes'' if i > `stop' ;
tab E_i , miss ;  // 汇总事件日期变量，显示缺失值情况

expand 20 ;/* 时间期限 */

sort i ;
qui by i: gen t = _n ;

xtset i t ;  // 定义面板数据的单位和时间变量

/* make variables that determine the DGP */
gen etime = (t - E_i) ;/* 事件时间 */

/* gen TE = 1 * (etime >= 0) ;/* 步骤功能处理效应 */
*/
gen TE = (etime >= 0) * (etime+1) ;/* 无尽的斜坡功能处理效应 */
replace TE = 0 if E_i == . ;

gen Y0_pure = 0 ;/* 最简单的反事实情况 */
*gen Y0_pure = 4 * treated +  0.3 * treated * t ;/* 处理组有前趋势 ... */
*gen Y0_pure = 4 * treated +  0.1 * treated * (t-10) * (E_i - 9);/* 基于 E_i 的前趋势 ... */

gen eps = sqrt(0.2) * rnormal() ;
gen actual = Y0_pure + TE * treated ;
gen y = actual + eps ;/* 观察到的 Y */


end ;

cap prog drop bigone ;
prog def bigone ;

runme "10 11" 0 0 ;

forvalues i = 4/16 { ;
runme "4 16 `i' " 0 0 ;
local rest_m3_`i' = count[1,3] - 2 ;
runme "4 16 `i' 10" 0 0 ;
local rest_m4_`i' = count[1,3] - 2 ;
} ;

drop _all ;
set obs 7 ;
gen E_i = _n + 3 ;
qui gen units3 = . ;
qui gen units4 = . ;
qui forvalues i = 4/10 { ;
replace units3 = `rest_m3_`i'' if E_i == `i' ;
replace units4 = `rest_m4_`i'' if E_i == `i' ;
} ;

gen E4_i = E_i + 0.15 ;
graph twoway (scatter units3 E_i , msymbol(o) msize(medsmall)) 
(scatter units4 E4_i , msymbol(triangle) msize(medsmall)) 
,
ti("需要额外参数限制") 
xtitle("E_3 的位置") 
note("# DiD 限制之外需要的限制。无对照组。T=20。E_1 = 4, E_2 = 16。"
"圆圈显示 3 种单位类型，作为第三种事件日期的函数。三角形显示 4 种单位类型，其中 E_4 = 10。") 
;
graph export figures/apdx_fig_A03.png , replace ;


end ;



bigone ;

#delimit ;

cap log close ;

set trace off ;
set more off ;
clear ;
set seed 2024202431 ;

cap prog drop runme ; 
prog def runme ;


gendata ;
/* 创建估计变量 */
/* 事件时间虚拟变量 */

forvalues i = 0/20 { ;
gen D_p`i' = (etime == `i') ;
gen D_m`i' = (etime == -1 * `i') ;
summ D_p`i' ;
local dropme = (r(mean) == 0) ;
if `dropme' == 1 { ;
drop D_p`i'  ;
} ;
summ D_m`i' ;
local dropme = (r(mean) == 0) ;
if `dropme' == 1 { ;
drop D_m`i'  ;
} ;
} ;

drop D_m0  ;

gen byte group2 = E_i == 11 ;
gen trend_group2 = t * group2 ;

summ ; 


/* 1 规范化事件事件-1的虚拟变量为0 */
constraint define 1 D_m1 = 0 ;

/* 2, 所有组的固定效应加总为0 */
constraint define 2 1.i + 2.i + 3.i + 4.i + 5.i + 
6.i + 7.i + 8.i + 9.i + 10.i = 0 ;

/* 7 处理前的事件时间虚拟变量的平均值设为零*/
constraint define 7 D_m10 + D_m9 + D_m8 + D_m7 + D_m6 + D_m5 + D_m4 + 
D_m3 + D_m2 + D_m1 = 0 ; 

cnsreg y D_m* D_p* 
ibn.t ibn.i , cluster(i) nocons constraints(2 1 )  collinear ;
matrix myb_ES1 = e(b) ;
matrix myV_ES1 = e(V) ;

cnsreg y D_m* D_p* 
ibn.t ibn.i , cluster(i) nocons constraints(2 7 )  collinear ;
matrix myb_ES2 = e(b) ;
matrix myV_ES2 = e(V) ;


tempfile main ;
save `main' ;

/* 画图 */

*开始调试;
tempfile pooled ;
forvalues i = 0/10 { ;
drop _all ;
set obs 2 ;
gen label = "m`i'" in 1 ;
replace label = "p`i'" in 2 ;

gen etime = -1 * `i' in 1 ;
replace etime = `i' in 2 ;

foreach j in 1 2 { ;
gen cf_m`j' = . ;

/* 设置治疗效应的估计方程 */
gen truth_m`j' = (0) in 1 ;
replace truth_m`j' = (1) in 2 ;

local myb = 0 ;
cap local myb = myb_ES`j'[1,"D_m`i'"] ;
gen ES_b_m`j' = `myb' in 1 ;
local myv = . ;
cap local myv = myV_ES`j'["D_m`i'","D_m`i'"] ;
local myse = sqrt(`myv') ;
gen ES_se_m`j' = `myse' in 1 ;

local myb = 0 ;
cap local myb = myb_ES`j'[1,"D_p`i'"] ;
replace ES_b_m`j' = `myb' in 2 ;
local myv = . ;
cap local myv = myV_ES`j'["D_p`i'","D_p`i'"] ;
local myse = sqrt(`myv') ;
replace ES_se_m`j' = `myse' in 2 ;
} ;


capture append using `pooled' ;
save `pooled' , replace ;
} ;




drop if label == "m0" ;
drop if label == "p10" ;
sort etime ;
save `pooled' , replace ;


local note1 "gam(-1)=0" ;
local note2 "Avg(gam(-))=0." ;

forvalues i = 1/2 { ;

gen top_`i' = ES_b_m`i' + 1.96 * ES_se_m`i' ; 
gen bot_`i' = ES_b_m`i' - 1.96 * ES_se_m`i' ; 

graph twoway (connected ES_b_m`i' truth_m`i' etime ,  msize(medsmall medium) 
msymbol(o oh) lpattern(solid dot) )
(line cf_m`i' etime , lpattern(dash) )
(line top_`i' bot_`i' etime, lpattern(dash dash) lcolor(brown brown))
, legend(off) xline(-0.5) 
note("`i': `note`i''") 
name(g`i' , replace) yscale(range(-3,3.5)) yline(0 , lpattern(dash) lcolor(black)) ;
} ;


graph combine g1 g2 , ti("不同等标准化估计会影响结果") 
note("All models: Avg(alpha)=0; DiD data structure. T=20, E1=11.") ;
graph export figures/apdx_fig_A04.png , replace ;


end ;

cap prog drop gendata ;
program define gendata ;

/*****************
主要数据生成过程选项

/* 可能的治疗效应类型 */
1 - 零治疗效应
2 - 阶梯函数治疗效应
3 - 永久爬坡
4 - 爬坡到平台
5 - AR(1)类型

/* E_i的分布 */
T = 19; E_i ~ U(6,14)

/* 是否存在从未被处理的单位？ */
是; 否

/* Y0（潜在结果的"信号"）动态 */
1 - 基本: Y0 = 0
2 - 水平变化: Y0_i,t = E_i
3 - 趋势变化: Y0_i,t = t * E_i

主要选项结束
******************/

drop _all ;

/*
/* 经典 ES 数据结构 */
set obs 9 ;/* 组数 */
gen i = _n ;
gen E_i = 5 + i ;
gen treated = 1 ;/* 全部实验组 */
expand 20 ;/* 时间段数 */
*/


set obs 10 ;/* 组数 */
gen i = _n ;

/* NxT DiD 数据结构*/
gen treated = i > 5 ;/* 一半实验组 */
gen E_i = 11 if treated == 1 ;


expand 20 ;/* 时间段数 */

sort i ;
qui by i: gen t = _n ;

xtset i t ;

/* 创建决定 DGP 的变量 */
gen D = (t == E_i) ; /* 事件"脉冲" */
gen etime = (t - E_i) ;/* 事件虚拟时间 */

gen TE = 1 * (etime >= 0) ;/* 实验组判断 */
*gen TE = (etime >= 0) * (etime+1) ;/* endless ramp function for treatment effect */
replace TE = 0 if E_i == . ;

gen treated_post = etime >= 0 ;

gen Y0_pure = 0 + 1 * treated ;/* 偏移量 counterfactual */
*gen Y0_pure = 4 * treated +  0.3 * treated * t ;/* treated have a pre-trend ... */
*gen Y0_pure = 4 * treated +  0.01 * treated * (t-10) * (E_i - 9);/* pre-trend based on E_i ... */

gen eps = sqrt(0.3) * rnormal() ;
gen actual = Y0_pure + TE * treated ;
gen y = actual + eps ;/* observed Y */


end ;


runme ;

/* 创建反事实预测值 cf_ES2 */
gen cf_ES2 = y ;  // 将 cf_ES2 初始化为观测变量 y

/* 按时间范围循环处理 */
qui summ etime ;  // 求 etime 的汇总统计
local mymin = r(min) ;  // 获取最小值
local mymax = r(max) ;  // 获取最大值

forvalues i = `mymin'/`mymax' {  // 循环遍历时间范围
    if `i' < 0 {  // 如果时间小于0
        local j = abs(`i') ;  // 取绝对值
        replace cf_ES2 = cf_ES2 - _b[D_m`j'] * D_m`j' ;  // 减去处理效应系数乘以对应的 dummy 变量
    }
    else {  // 如果时间大于等于0
        replace cf_ES2 = cf_ES2 - _b[D_p`i'] * D_p`i' ;  // 减去处理效应系数乘以对应的 dummy 变量
    }
}

/* 按单元类型和时间求反事实预测均值 */
egen meancf_ES2 = mean(cf_ES2), by(t unittype) ;  // 按 t 和 unittype 求 cf_ES2 的均值

#delimit; 
/*两期一组*/
constraint define 9 D_m10=D_m9;
constraint define 10 D_m8=D_m7;
constraint define 11 D_m6=D_m5;
constraint define 12 D_m4=D_m3;
constraint define 13 D_m2=D_m1;
constraint define 14 D_p0=D_p1;
constraint define 15 D_p2=D_p3;
constraint define 16 D_p4=D_p5;
constraint define 17 D_p6=D_p7;
constraint define 18 D_p8=D_p9;

/*三期一组 */
#delimit;
constraint define 21 D_m1=D_m2;
constraint define 22 D_m1=D_m3;

constraint define 23 D_m4=D_m5;
constraint define 24 D_m4=D_m6;

constraint define 25 D_m7=D_m8;
constraint define 26 D_m7=D_m9;

constraint define 27 D_p0=D_p1;
constraint define 28 D_p2=D_p1;

constraint define 29 D_p3=D_p4;
constraint define 30 D_p3=D_p5;

constraint define 31 D_p6=D_p7;
constraint define 32 D_p6=D_p8;