EDA 探索性分析与样本构造 Exploratory Data Analysis

* 三件套：散点 + lowess + binscatter
twoway (scatter y x) (lfit y x), name(sc1, replace)
twoway (lowess y x, bw(0.3)), name(sc2, replace)
binscatter y x, nquantiles(20) name(sc3, replace)

* 并排放三张图
graph combine sc1 sc2 sc3, rows(1) ycommon
graph export "output/figures/bivariate_eda.png", replace

* 按 (year, treat) 计算 Y 的均值
collapse (mean) y, by(year treat)
twoway (line y year if treat==1) ///
       (line y year if treat==0), ///
       legend(label(1 "处理组") label(2 "对照组")) ///
       xline(2012, lpattern(dash) lcolor(red)) ///  // 政策实施年份
       title("政策前后两组均值趋势")
* 如果政策前两组基本平行，DID 平行趋势假设看起来合理
* 如果政策前两组已经在发散，DID 就有问题

* 假设 kept==1 是保留下来的样本
gen byte kept = !missing(Y, X, controls)

* 对比被删组 vs 保留组在 X、Y、控制变量上的均值
foreach v in X Y age size leverage {
    quietly ttest `v', by(kept)
    display "变量 `v': 被删组均值=" r(mu_1) " 保留组均值=" r(mu_2) ///
            "  p值=" r(p)
}
* 如果 p 值很小，说明"能否进样本"和这些变量系统相关
* 你的 Sample 层有选择偏误，需要在论文里讨论

*==============================================================
* 14_eda_full.do：模拟数据的全套 EDA
* 目标：在跑任何回归之前，先用眼睛看一遍数据
*==============================================================
clear all
set more off
set seed 20260911

* --- 1. 生成模拟面板数据 ---
* 1000 家企业 × 10 年 = 10000 观测
set obs 1000
gen firm_id = _n
gen treat = (runiform() < 0.5)            // 处理组 50%
expand 10
bysort firm_id: gen year = 2010 + _n - 1

* 处理变量：2015 年及以后处理组受到政策冲击
gen post = (year >= 2015)
gen did = treat * post

* 协变量：企业规模（对数正态）、年龄
gen size = exp(rnormal(2, 1))
gen age  = floor(runiform()*20 + 1)

* 结果变量：企业利润（对数）
* 真实模型：lprofit = 0.5*did + 0.3*ln(size) - 0.01*age + 噪声
gen ln_size = log(size)
gen lprofit = 0.5*did + 0.3*ln_size - 0.01*age + rnormal(0, 1)

* 故意加几个异常值（1% 的极端大值）
replace lprofit = lprofit + 8 if runiform() < 0.01

* 故意让部分企业在 2010 年缺失利润（模拟选择偏误）
replace lprofit = . if year==2010 & treat==1 & runiform()<0.3

*==============================================================
* 2. 单变量分布 EDA
*==============================================================
* 2.1 直方图 + 正态叠加
histogram lprofit, normal name(h_lp, replace) ///
    title("利润对数分布")
graph export "eda_1_hist_lprofit.png", replace

* 2.2 核密度（处理组 vs 对照组）
twoway (kdensity lprofit if treat==1) ///
       (kdensity lprofit if treat==0), ///
       legend(label(1 "处理组") label(2 "对照组")) ///
       name(kd_grp, replace)
graph export "eda_2_kdensity_group.png", replace

* 2.3 箱线图：按处理组
graph box lprofit, over(treat) name(box_grp, replace)
graph export "eda_3_boxplot.png", replace

*==============================================================
* 3. 双变量关系 EDA
*==============================================================
* 3.1 binscatter：ln(size) vs lprofit
binscatter lprofit ln_size, nquantiles(30) ///
    name(bin_size, replace) ///
    title("企业规模 vs 利润")
graph export "eda_4_binscatter_size.png", replace

* 3.2 lowess：年龄 vs lprofit（看是否非线性）
lowess lprofit age, bw(0.3) name(lowess_age, replace)
graph export "eda_5_lowess_age.png", replace

*==============================================================
* 4. 时间趋势 EDA（DID 平行趋势预检）
*==============================================================
preserve
    collapse (mean) lprofit, by(year treat)
    twoway (line lprofit year if treat==1) ///
           (line lprofit year if treat==0), ///
           legend(label(1 "处理组") label(2 "对照组")) ///
           xline(2015, lpattern(dash) lcolor(red)) ///
           name(trend, replace) ///
           title("政策前后利润趋势（2015 为政策年）")
    graph export "eda_6_trend.png", replace
restore

*==============================================================
* 5. 分组均值 + 置信区间
*==============================================================
* 2015 年前 vs 后，处理组 vs 对照组
preserve
    keep if !missing(lprofit)
    ci means lprofit, by(treat post)
restore

*==============================================================
* 6. 样本构造 flow chart
*==============================================================
count
scalar N0 = r(N)
display "Step0 原始观测: " N0

* 只保留 2010 年以后（其实数据就是 2010 起）
keep if year >= 2010
count
scalar N1 = r(N)
display "Step1 保留2010年后: " N1

* 删除 lprofit 缺失
drop if missing(lprofit)
count
scalar N2 = r(N)
display "Step2 删除lprofit缺失: " N2

* 缩尾 1%/99%（注意：缩尾不改变 N）
capture which winsor2
if _rc ssc install winsor2, replace
winsor2 lprofit ln_size, cuts(1 99) replace

count
scalar N3 = r(N)
display "Step3 缩尾后: " N3 "（N 不变）"

display _n "===== Flow Chart ====="
display "N0 = " N0
display "N1 = " N1 " (删了 " N0-N1 ")"
display "N2 = " N2 " (删了 " N1-N2 ")"
display "N3 = " N3 " (最终回归 N)"

*==============================================================
* 7. 样本选择诊断：被删的是谁？
*==============================================================
* 重新读一遍原始数据，标记哪些被删了
preserve
    * 这里略：实际项目中用 kept 哑变量 + ttest 对比
    display "请参考本页第 07 节的 ttest 模板"
restore

*==============================================================
* 8. 最后才跑回归
*==============================================================
reg lprofit did ln_size i.age, vce(cluster firm_id)
estimates table, b(%9.3f) se(%9.3f) stats(N r2)

display _n "===== EDA 全部完成 ====="
display "先看图，再看回归表。两者必须对得上。"