显著性操作手法库与学术规范判别 Significance Toolkit & Research-Integrity Discipline

* A1 样本筛选：合法(事前制度窗口) vs 红线(事后挑子样本)
clear all
set seed 101
set obs 1000
gen id = _n
gen year = 2005 + mod(_n,15)          // 2005-2019
gen soe = runiform()>0.5              // 国企虚拟
gen x  = rnormal()
gen y  = 0.2*x + rnormal(0,1.2)       // 真实效应偏弱

* (合法)事前由政策时点限定：2011年起
reg y x if year>=2011
* (合法)事前理论关心国企子样本，并报告交互项组间差异
reg y x if soe==1
gen x_soe = x*soe
reg y x soe x_soe
test x_soe                            // 必须用交互项验组间差异

* (红线演示，勿用于正式研究) 事后从2005..2019逐年试窗口
* forvalues s=2005/2019 {
*   quietly reg y x if year>=`s'
*   display "起始年=`s'  系数=" %6.3f _b[x] "  p=" %6.3f (2*ttail(e(df_r),abs(_b[x]/_se[x])))
* }
* 上面这段"扫窗口"在论文里只能作为敏感性曲线(specification curve)整体报告，
* 绝不能只挑其中一个显著的窗口写进正文。

* A2 样本量：同一小效应，N=200 不显著，N=20000 显著
clear all
set seed 202
forvalues n=200(4950)10000 {
    preserve
    quietly set obs `n'
    gen x = rnormal()
    gen y = 0.03*x + rnormal(0,1)     // 真实斜率极小(0.03)
    quietly reg y x
    display "N=" %6.0f `n' "  b=" %6.4f _b[x] "  t=" %6.2f (_b[x]/_se[x])
    restore
}
* 观察：N 越大 t 越大，但效应量始终是0.03
* 结论：大样本下的"显著"必须配经济显著性换算，否则是误导

* A3 删高缺失控制变量导致 N 变化
clear all
set seed 303
set obs 8000
gen x = rnormal()
gen y = 0.3*x + rnormal()
gen ctrl_missing = .
replace ctrl_missing = rnormal() in 1/7500      // 仅前7500个有值，500个缺失
gen ctrl_other = rnormal()

* 含高缺失变量：N 只剩 7500
reg y x ctrl_missing ctrl_other
* 删掉高缺失变量：N 回到 8000，样本变大、SE 变小
reg y x ctrl_other
* 诚实做法：报告两套 N，并说明删除理由
* 不诚实做法：发现删完后 x 变显著，就悄悄把 ctrl_missing 从表里删掉

* A4 winsor2 压回 vs trim 删除；可 by(Year)/by(Ind) 分组
* ssc install winsor2
clear all
set seed 404
set obs 2000
gen year = 2010 + mod(_n,10)
gen ind  = ceil(runiform()*20)
gen x = rnormal()
replace x = 50 in 1/20                 // 20个极端值
gen y = 0.3*x + rnormal(0,1.5)

* 压回(winsor)：保留观测数，把尾部压到分位
winsor2 x, cuts(1 99) suffix(_w)
* 按年份分组缩尾（消除跨年量纲漂移）
winsor2 x, cuts(1 99) suffix(_wy) by(year)
* 按行业分组缩尾
winsor2 x, cuts(1 99) suffix(_wi) by(ind)

reg y x      ; est store raw
reg y x_w    ; est store w1
reg y x_wy   ; est store w_year
esttab raw w1 w_year, keep(x x_w x_wy) b(%9.3f) se star(*.1 **.05 ***.01)
* 红线：不要从 cuts(1 99)/(2 98)/(5 95) 里挑一个显著的写正文

* A5 先 log 还是先 winsor：固定 pipeline，不要试
clear all
set seed 505
set obs 3000
gen rev = exp(rnormal(3,1.5))          // 右偏的营收
gen x   = rnormal()
gen y   = 0.2*x + 0.1*ln(rev+1) + rnormal()

* 顺序1：先缩尾后 log
winsor2 rev, cuts(1 99) suffix(_w)
gen ln_rev1 = ln(rev_w + 1)
* 顺序2：先 log 后缩尾
gen ln_rev0 = ln(rev + 1)
winsor2 ln_rev0, cuts(1 99) suffix(_w)

reg y x ln_rev1   ; est store order1
reg y x ln_rev0_w ; est store order2
esttab order1 order2, b(%9.3f) se
* 把你采用的那套顺序写进数据清洗章节，另一套放附录即可

* A6 Cook's D / DFBETA：极端值≠有影响力
clear all
set seed 606
set obs 1200
gen id = _n
gen x = rnormal()
replace x = 12 in 1/5                  // 5个高杠杆点
gen y = 0.3*x + rnormal(0,1.5)

reg y x
predict cook, cooksd                   // Cook's D > 4/N 视为高影响
predict lev, leverage
dfbeta                                 // 逐点删去后 x 系数变化
list id x y cook if cook > 4/_N

* 合法：经查是错误观测，有依据地剔除
reg y x if cook <= 4/_N
* 红线方向：把"按 b_x 排序删观测"做成循环（见 F1）——切勿使用

* B1 控制变量逐级加入，esttab 输出 Model1-3
clear all
set seed 707
set obs 1500
gen x = rnormal()
gen c1 = rnormal()                      // 企业控制
gen c2 = rnormal()                      // 行业控制(制造缺失)
replace c2 = . in 1/200
gen y = 0.25*x + 0.3*c1 + 0.5*c2 + rnormal()

reg y x                              ; est store m1
reg y x c1                           ; est store m2
reg y x c1 c2                        ; est store m3
esttab m1 m2 m3, keep(x c1 c2) b(%9.3f) se star(*.1 **.05 ***.01) ///
    mtitles("Model1 无控制" "Model2 企业控制" "Model3 全控制")
* 注意 m3 的 N 比 m2 小200（因 c2 缺失）——正文必须写明
* 红线：逐格试 c1..c10 哪个加进去 x 就显著，只留那个

* B2 核心变量多种操作化，事前锁定主定义
clear all
set seed 808
set obs 1200
gen id = ceil(_n/6)
gen year = 2010 + mod(_n,6)
xtset id year
gen policy = rnormal()
gen y = 0.2*l.policy + rnormal()

gen policy_l1 = l.policy               // 滞后一期
gen policy_d  = d.policy               // 一阶差分
gen policy_dum = policy>0              // 0-1 虚拟
gen policy_stage = cond(year<=2012,0,cond(year<=2015,1,2))  // 阶段

reg y policy      ; est store def_level
reg y policy_l1  ; est store def_lag
reg y policy_d    ; est store def_diff
reg y i.policy_stage#c.policy        ; est store def_stage
esttab def_level def_lag def_diff, b(%9.3f) se star(*.1 **.05 ***.01)

* B3 改善分布(log/标准化) vs 改变构念(连续切分类)
clear all
set seed 909
set obs 1500
gen x_raw = exp(rnormal(0,1))          // 右偏
gen y = 0.3*ln(x_raw) + rnormal()

* 合法：log 改善分布
gen ln_x = ln(x_raw)
reg y ln_x
* 合法：标准化
sum ln_x
gen z_x = (ln_x-r(mean))/r(sd)
reg y z_x
* 红线方向：在10/25/50/75分位数上逐个试切点
* xtile x_bin = x_raw, nq(2)  // 切点必须事前定，不能事后挑

* B4 逆向指标正向化：事前定义一次，不再改
clear all
set seed 111
set obs 1000
gen cost = rnormal(0,1)                 // 融资成本(越高越差)
gen y = -0.3*cost + rnormal()           // 成本越高，y越低

* 事前定义：融资便利 = -成本（正向指标）
gen ease = -cost
reg y ease                              // 预期为正，与故事一致
* 红线：为了"系数显著为负"，有时报 cost 有时报 -cost

* B5 企业年龄控制：ln(1+age)
clear all
set seed 222
set obs 2000
gen listyear = 1990 + floor(runiform()*30)
gen year = 2020
gen age = year - listyear               // 上市年限
gen lnage = ln(1+age)                   // 取对数+1 缓解偏态
gen x = rnormal()
gen y = 0.2*x - 0.01*age + rnormal()
reg y x age
reg y x lnage                           // 推荐：对数化后更稳健

* C1 按因变量性质选模型，不要试
clear all
set seed 333
set obs 1500
gen id = ceil(_n/5)
gen x = rnormal()
gen y_latent = 0.5*x + rnormal()
gen y_bin = (y_latent>0)                // 二值因变量
gen y_count = poisson(0.6+0.4*x)        // 计数因变量

reg y_bin x        ; est store ols_lpm      // LPM(可作稳健性)
logit y_bin x      ; est store logit
probit y_bin x     ; est store probit
xtset id
xtlogit y_bin x, fe                         // FE logit 只用组内
poisson y_count x  ; est store poisson
nbreg y_count x    ; est store negbin
esttab ols_lpm logit probit, b(%9.3f) se star(*.1 **.05 ***.01)

* C2 固定效应组合：reghdfe 吸收多维 FE
clear all
set seed 444
set obs 2000
gen firm = ceil(_n/10)
gen year = 2010 + mod(_n,10)
gen ind  = ceil(runiform()*15)
gen x = rnormal() + 0.3*rnormal(ind)    // 行业层面冲击
gen y = 0.3*x + rnormal(firm) + rnormal(year)

reghdfe y x, absorb(firm) vce(cluster firm)            ; est store fe_firm
reghdfe y x, absorb(firm year) vce(cluster firm)       ; est store fe_2way
reghdfe y x, absorb(ind#year) vce(cluster ind)         ; est store fe_indyear
esttab fe_firm fe_2way fe_indyear, keep(x) b(%9.3f) se star(*.1 **.05 ***.01)

* C3 行业代码粒度 Ind1/Ind2/Ind3
clear all
set seed 555
set obs 2000
gen ind1 = ceil(runiform()*20)          // 一位
gen ind2 = ind1*10 + ceil(runiform()*5)// 两位(100个)
gen ind3 = ind2*10 + ceil(runiform()*3)// 三位(300个)
gen x = rnormal()
gen y = 0.3*x + rnormal(ind2) + rnormal()

reghdfe y x, absorb(i.ind1) vce(cluster ind1) ; est store g1
reghdfe y x, absorb(i.ind2) vce(cluster ind2) ; est store g2
reghdfe y x, absorb(i.ind3) vce(cluster ind3) ; est store g3
esttab g1 g2 g3, keep(x) b(%9.3f) se star(*.1 **.05 ***.01)

* C4 二次项识别 U 型/倒 U；交互识别调节
clear all
set seed 666
set obs 1500
gen x = rnormal(-2,2)
gen y = 1 - 0.4*x^2 + rnormal()        // 真实倒U

reg y x                                 // 线性误设：平均斜率≈0
gen x2 = x^2
reg y c.x##c.x                          // = reg y x c.x#c.x
test x2                                 // 二次项显著=非线性
* 调节/阶段：c.x##i.Stage
gen stage = cond(x<0,0,1)
reg y c.x##i.stage
margins, dydx(x) at(stage=(0 1))       // 两阶段各自边际效应

* C5 xtreg vs reghdfe；R² 分项；singleton 样本
clear all
set seed 777
set obs 3000
gen firm = ceil(_n/20)
gen year = 2010 + mod(_n,10)
gen x = rnormal()
gen y = 0.3*x + rnormal(firm) + rnormal(year)

xtset firm year
xtreg y x i.year, fe                    ; est store xtreg_fe
reghdfe y x, absorb(firm year)          ; est store reghdfe_r2
esttab xtreg_fe reghdfe_r2, keep(x) b(%9.3f) se
* reghdfe 输出会显示 absorbed FE 个数与 R² 分解
* 检查 singleton：被 absorb 掉的单观测组
display "实际估计 N = " e(N)

* C6 稀有事件因变量：firthlogit / relogit
* ssc install firthlogit
clear all
set seed 888
set obs 3000
gen x = rnormal()
gen y_latent = -3 + 0.5*x + rnormal()
gen rare = (y_latent>0)                 // 事件仅约5%

sum rare                                 // 确认稀有程度
logit rare x           ; est store plain_logit
firthlogit rare x, iter(100)            ; est store firth
esttab plain_logit firth, b(%9.3f) se star(*.1 **.05 ***.01)
* 稀有事件下普通 logit 系数通常偏大，firth 更可靠

* D1 同一系数，四种标准误，p 值天差地别
clear all
set seed 999
set obs 2000
gen firm = ceil(_n/5)
gen ind  = ceil(runiform()*20)
gen x = rnormal()
gen y = 0.3*x + rnormal(ind) + rnormal()   // 同行业残差相关

reg y x, vce(ols)                    ; est store s_ols
reg y x, vce(robust)                 ; est store s_rob
reg y x, vce(cluster firm)           ; est store s_firm
reghdfe y x, absorb(ind) vce(cluster ind) ; est store s_ind
* 双向聚类
reghdfe y x, absorb(ind) vce(cluster ind firm) ; est store s_two
esttab s_ols s_rob s_firm s_ind s_two, keep(x) b(%9.3f) se star(*.1 **.05 ***.01)
* 红线：哪一列最窄就只报哪一列；合法：论证后报告最保守且可辩护的那列

* D2 加权：必须有抽样设计/理论依据
clear all
set seed 1010
set obs 1500
gen country = ceil(runiform()*10)
gen pop = runiform()*1000              // 各国人口权重
gen x = rnormal()
gen y = 0.2*x + rnormal(country)

* 不加权：每个国家等权
reg y x
* 人口加权：估计"人口意义上的"平均效应
reg y x [pw=pop]
* 合法：因为目标总体是全球人口加权平均，理论上本就该加权重
* 红线：发现不加权不显著、加权后才显著，就只报加权版

* E1 滞后结构：理论定阶，不逐期扫显著
clear all
set seed 1212
set obs 1200
gen id = ceil(_n/6)
gen year = 2010 + mod(_n,6)
xtset id year
gen x = rnormal()
gen y = 0.3*l2.x + rnormal()          // 真实效应滞后两期

reg y x        ; est store lag0        // 当期：弱
reg y l1.x     ; est store lag1
reg y l2.x     ; est store lag2        // 理论预期的传导期
reg y l3.x     ; est store lag3
esttab lag0 lag1 lag2 lag3, b(%9.3f) se star(*.1 **.05 ***.01)
* 合法：理论说"政策两期后起效"，主报 L2，其余稳健性
* 红线：L1..L5 全扫一遍，哪个星多报哪个

* E2 平行趋势：完整事件研究，不缝合
clear all
set seed 1313
set obs 1500
gen firm = ceil(_n/5)
gen year = 2005 + mod(_n,10)
gen treat = (firm>100)
gen post  = year>=2012
gen x = rnormal()
gen y = 0.4*treat*post + 0.2*treat*(year-2012) + rnormal()  // 事前就有趋势

* 事件研究：相对政策年的虚拟变量
gen rel = year - 2012
forvalues k=-4/3 {
    gen e`k' = treat*(rel==`k')
}
reghdfe y e-4(1/0) e1 e2 e3, absorb(firm year) vce(cluster firm)
* 看 e-4..e-1 是否显著（若显著=平行趋势不通过）
* 合法：报告全部事前系数；红线：只画 e0 之后的部分"显得平行"

* E3 逐步回归：红线演示，仅用于理解其风险
clear all
set seed 1414
set obs 1200
gen y = rnormal()
* 制造20个纯噪声控制变量
forvalues j=1/20 {
    gen c`j' = rnormal()
}
gen x = rnormal() + 0.2*c1
replace y = y + 0.1*x

* (红线)按 p 值逐步加入控制：会把噪声变量也选进来
stepwise , pe(0.05): reg y x c1-c20
* 观察：被选入的控制变量里混着纯噪声，核心 x 的 p 值被玩弄
* 诚实替代：理论事前定控制集，如 reg y x c1 c3 c7
reg y x c1 c3 c7

* G1 事件研究图：eventdd 一键 vs reghdfe 手写 leads/lags
* ssc install eventdd, replace
* ssc install coefplot,  replace
clear all
set seed 2301
set obs 600
gen id = ceil(_n/6)
gen year = 2005 + mod(_n,10)
gen g = cond(id<=60, 2012, 0)            // 前60个2012年处理，其余never-treated
xtset id year
gen treat = (g!=0 & year>=g)
gen rel = year - g if g!=0
gen y = 0.8*treat + rnormal(id) + rnormal(year) + rnormal()

* 方法1：eventdd 一键出图（含95% CI）
*   lead(4) lag(3) ：政策前4期、政策后3期
*   baseline(-1)  ：以 k=-1 为基准
*   accum(.)      ：pre期系数自动汇总
eventdd y i.year, timevar(rel) lead(4) lag(3) baseline(-1) ///
    controls(i.id) cluster(id) graph_op( ///
    yline(0, lcolor(red)) xline(0.5, lpattern(dash)) ///
    title("事件研究图：leads/lags 与95%置信带") ///
    ytitle("估计系数") xtitle("相对政策时点 k"))

* 方法2：reghdfe 手写 leads/lags（更透明，推荐主回归用）
forvalues k=-4/3 {
    local kk = `k' + 100
    gen e`kk' = (rel==`k') if g!=0
    replace e`kk' = 0 if g==0
}
* 省略 e099 (即 k=-1) 作为基准
reghdfe y e096 e097 e098 e100 e101 e102 e103, ///
    absorb(id year) vce(cluster id)
* e096=k=-4, e097=k=-3, e098=k=-2, e100=k=0...
coefplot, keep(e096 e097 e098 e100 e101 e102 e103) vertical ///
    yline(0) xline(4.5, lp(dash)) ///
    order(e096 e097 e098 e100 e101 e102 e103) ///
    title("手写 leads/lags 事件研究图")
* 看图：pre期 e096/e097/e098 应贴0；post期 e100..e103 应显著为正

* G2 pre-trend 联合检验 + Roth (2022) pretrends 功效分析
* ssc install pretrends, replace
* 沿用 G1 模拟数据：reghdfe 后已有 e096(k=-4) e097(k=-3) e098(k=-2)

* 1) 朴素做法：pre 期 leads 联合 F 检验
reghdfe y e096 e097 e098 e100 e101 e102 e103, ///
    absorb(id year) vce(cluster id)
test e096 e097 e098          // H0: pre期系数全=0
* 若 p>0.05，只能说"不能拒绝平行趋势"，不能说"平行趋势成立"

* 2) Roth (2022) pretrends：自动评估检验功效
*    pretrends 会:
*      - 自动识别 pre 期虚拟变量
*      - 估计"如果真实存在线性 pre-trend，本检验有多大概率抓住"
*      - 报告 " detectable bias " —— 在当前功效下，
*        多大的 pre-trend 漂移才会被检测到
*   语法（在 event study 回归之后跑）：
* pretrends, leads(4) lag(0) post(3)
*   输出关键量:
*     Power of pre-trend test (5% level) ≈ 0.3x
*     Partial R^2 of pre-trend ...
* 解读:
*   - 若 power 只有 0.35：意味着真实 pre-trend 即使把你的
*     ATT 偏误 20%，你也有 65% 概率检验不出来。
*   - 论文写法:
*     "We fail to reject parallel trends (joint F p=0.21),
*      but the pre-trend test has limited power (~35% for a
*      trend that would bias ATT by 20%), so we additionally
*      report Rambachan-Roth sensitivity bounds."

* 3) 观测量合并检验：当某事件期观测太少（如远端 leads 只有
*    20 个处理组），单独系数 SE 巨大，应合并相邻期：
*    例如把 k=-4 与 k=-3 合并成 "leads early"
gen e_pre_early = (rel<=-3) if g!=0
replace e_pre_early = 0 if g==0
reghdfe y e_pre_early e098 e100 e101 e102 e103, ///
    absorb(id year) vce(cluster id)
* 注意：合并必须在看结果前由"该期观测太少"这一客观理由决定，
*      不能发现合并后"恰好平行"才合并。

* G3 安慰剂检验：伪时点 + 伪处理组随机化
* 沿用 G1 模拟数据
clear all
set seed 2303
set obs 600
gen id = ceil(_n/6)
gen year = 2005 + mod(_n,10)
gen g = cond(id<=60, 2012, 0)
xtset id year
gen treat = (g!=0 & year>=g)
gen y = 0.8*treat + rnormal(id) + rnormal(year) + rnormal()

* (a) 虚构处理时间：把政策时点提前到 2008（伪政策）
*     若 2008 版 DID 显著，说明 pre-trend 或隐藏趋势存在
gen treat_fake = (g!=0 & year>=2008)
reghdfe y treat_fake, absorb(id year) vce(cluster id)
* 期望：treat_fake 系数不显著（因为真实政策是 2012）
* 若显著 → 平行趋势可疑

* (b) 随机分配伪处理组：500 次蒙特卡洛
permute treat, reps(500) seed(2303) ///
    saving("placebo_dist.dta", replace): ///
    reghdfe y treat, absorb(id year) vce(cluster id)
* 真实系数 0.8 应落在 placebo 分布的 95% 分位之外
use "placebo_dist.dta", clear
sum _b_treat, detail
* 查看 p10/p50/p90 分位；真实 0.8 若在 p99 右侧 → 安慰剂通过

* G4 post 期动态效应：联合显著 + 趋势
* 沿用 G1 回归：e100=k=0, e101=k=1, e102=k=2, e103=k=3
reghdfe y e096 e097 e098 e100 e101 e102 e103, ///
    absorb(id year) vce(cluster id)

* (i) post 期联合显著性
test e100 e101 e102 e103
* (ii) post 期线性趋势：e101+2*e102+3*e103 是否单调
test (e101 - e100) = (e102 - e101) = (e103 - e102)
* 若拒绝 → post 期不是线性单调，需解释曲线形状

* 诚实写法：把 e100..e103 全部列在主回归表里，
* 不要只挑 e101 一个显著系数说"政策第2年起作用"。

* G5 交叠 DID 的 pre-trend 诊断：csdid 事件期 leads 联合检验
* ssc install csdid, replace
* ssc install drdid,  replace
clear all
set seed 2305
set obs 200
gen id = _n
* 三队列：g=2008 (50个), g=2012 (50个), g=0 never-treated (100个)
gen g = cond(_n<=50, 2008, cond(_n<=100, 2012, 0))
expand 11
bysort id: gen year = 1999 + _n
xtset id year
gen treat = (g!=0 & year>=g)
gen y = 0.8*treat + rnormal(id) + rnormal(year) + rnormal()

* csdid 主估计
csdid y, ivar(id) time(year) gvar(g) method(dripw)
* 事件时间聚合（含 leads/lags）
estat event
* csdid_plot 画干净的事件研究图（不含负权重污染）
csdid_plot, title("csdid 事件研究图（交叠 DID 专用）") ///
    legend(off) yline(0)

* 关键：在 csdid 输出上做 pre 期 leads 联合检验
*   estat event 返回的系数按事件期排列；
*   pre 期（k<0）系数联合为 0 的检验
* 实务写法：
*   "We report the pre-trend joint test from csdid:
*    the leads k=-4..-1 are jointly insignificant
*    (Wald p=0.31), and pretrends (Roth 2022) indicates
*    power of approximately 45% to detect a trend that
*    would bias ATT by 20%."
*
* 对照：传统 TWFE 事件研究图在此处会被污染，
*       不能作为交叠 DID 的 pre-trend 证据。

* G6 平行趋势不通过时的合法处理（红线示例注释在内）
* 沿用 G5 模拟数据

* 合法方案1：局部窗口 —— 仅用政策前后±2年
*   窗口必须在看结果前由制度依据决定
preserve
    keep if year>=2006 & year<=2014   // 政策(2008/2012)前后各±4年
    csdid y, ivar(id) time(year) gvar(g) method(dripw)
    estat event
restore

* 合法方案2：csdid 加协变量做条件平行趋势
*   csdid 会按 X 做 IPW 平衡处理组/对照组的可观测特征
csdid y alpha_i gamma_t, ivar(id) time(year) gvar(g) ///
    method(dripw) covariates(alpha_i gamma_t)
estat event

* 合法方案3：Bacon 分解找出污染队列
* bacondecomp y treat, ddetail
* 若"Later T vs Earlier C"权重>0.3，说明某批队列在污染

* 合法方案4：Rambachan-Roth 敏感性界（Stata 中需配合
*   R 包 Honor & Roth, 或手动给出不同 Mbar 下的 CI）
*   Mbar=0 表示严格平行趋势；Mbar=1 允许 pre 期有线性偏离。
*   报告 "Even when Mbar=1, the lower bound of the ATT
*   CI remains positive at 0.21, so our conclusion is
*   robust to modest pre-trend violations."

* ⛔ 红线（严禁）：
*   - 先跑全样本 pre-trend 不通过，回头把窗口改成±2年"凑通过"
*   - 把漂移严重的队列从样本里删掉，理由是"该组异常"
*   - 只在正文报"pre-trend 通过"，附录才出现不通过的版本
*   审稿人只要你交出全部 event-time 系数 + do 文件，立刻识破。

* H1 坏控制识别：处理后变量不能进主回归
clear all
set seed 2311
set obs 2000
gen X = rnormal()
gen U = rnormal()                       // 不可观测混淆
gen M = 0.5*X + U + rnormal()          // 中介：X→M（处理后）
gen Y = 0.4*X + 0.3*M + U + rnormal()   // Y = X + M + U

* 主回归（合法）：只控制前定变量，不控制 M
reg Y X U
* 真实总效应 X→Y ≈ 0.4 + 0.3*0.5 = 0.55

* 坏控制（红线）：把中介 M 当控制
reg Y X M U
* 系数会从 0.55 掉到 ≈0.4 —— 你"控制掉"了 X→M→Y 这条通路
* 这不是"更严格"，是机制偏误

* 合法机制检验：两步法
reg M X U           ; est store stage1   // X→M
reg Y X U           ; est store total    // X→Y 总效应
* 中介效应 = 间接 = total - direct；不要在 Y 回归里放 M

* H2 事前交互项 vs 事后挑子样本
clear all
set seed 2312
set obs 2000
gen X = rnormal()
gen H = runiform()>0.5                  // 异质性维度(事前定义)
gen Y = 0.3*X + 0.1*X*H + rnormal()     // 真实异质性很小

* 合法：全样本交互项 + test 组间差异
reg Y c.X##i.H
test 1.H#c.X                            // 组间差异是否显著
* 关键：即使 test 不显著，也要如实报

* 红线：分样本各跑一次，挑带星的那张讲故事
* reg Y X if H==1   // 国企子样本
* reg Y X if H==0   // 民企子样本
* 两组系数 t=1.65 vs t=2.10，就写"异质性显著"
* 错误：两个独立回归的系数不能直接比；必须用交互项 test。

* H3 多重检验校正：Bonferroni / Holm / BH-FDR
* ssc install multproc, replace
clear all
set seed 2313
set obs 1500
gen X = rnormal()
* 构造 6 个"因变量"，其中 2 个真有信号、4 个纯噪声
forvalues j=1/6 {
    gen Y`j' = cond(`j'<=2, 0.25*X, 0) + rnormal()
}

* 未校正：6 个回归里挑 p<0.05 的报
forvalues j=1/6 {
    reg Y`j' X
    display "Y`j': p=" %6.4n (2*ttail(e(df_r),_b[X]/_se[X]))
}

* 把 6 个 p 值存成矩阵，做三种校正
* multproc 用法示例（按实际输出填 p 向量）：
* matrix pvec = (p1,p2,p3,p4,p5,p6)
* multproc, p(pvec) methods(bonferroni holm benjamini)
* 输出：每个原始 p 对应的 Bonferroni/Holm/BH 调整后 p

* 实务建议：
*   - 主假设(1个)：不校正
*   - 次要/探索性(6-20个)：报告 BH-FDR 校正后 p
*   - Bonferroni 太保守，只在 p 本来就极显著(如 p<.001)时才安全

* H4 Specification Curve：Simonsohn et al. (2020)
* 把所有合理规格的系数按大小排成曲线
clear all
set seed 2314
set obs 2000
gen x = rnormal()
gen c1 = rnormal()
gen c2 = rnormal()
gen c3 = rnormal()
gen y = 0.2*x + 0.3*c1 + rnormal()

* 定义规格：是否加 c1/c2/c3、缩尾 1%/5%、聚类/不聚类
* （简化演示：跑 8 个规格）
matrix sc = J(1, 8, .)
local j = 0
foreach cset in " " "c1" "c1 c2" "c1 c2 c3" {
    foreach clu in "robust" "cluster id" {
        local ++j
        quietly reg y x `cset', vce(`clu')
        matrix sc[1, `j'] = _b[x]
    }
}
* 把矩阵列拉成变量，画曲线
matrix list sc
* 横轴 = 规格编号；纵轴 = 系数；一条水平 y=0 线
* 诚实结论：
*   "The specification curve across 8 reasonable
*    specifications shows the median coefficient is 0.21,
*    and 7/8 specifications have the expected sign."
* 红线：只挑 8 个里 1 个显著的写正文，其余不提。

* H5 功效分析：事前 + 事后
* Stata 13+ 内置 power 命令

* 事前：双样本均值差异，要检出 δ=0.3 SD，power=0.8，α=0.05
power twomeans 0, diff(0.3) power(0.8) alpha(0.05)

* 事后：N=200 时，能以 80% 功效检出多大的处理效应？
*   输出 MDE (minimum detectable effect)
power twomeans 0, n(200) power(0.8) alpha(0.05)

* DID 语境下更实用的事后功效：
*   假设 SE(β)=0.15，要 80% 功效、双侧 5%，
*   临界 t≈1.96+0.84=2.80，MDE = 2.80*0.15 = 0.42
*   含义：只有真实 |ATT|>0.42，你才有 80% 概率检出。
*   若你的点估计是 0.15，诚实写法:
*   "We cannot rule out an effect as large as 0.42
*    with 80% power; our null result should not be
*    interpreted as evidence of no effect."

* H6 TOST 等价检验：把"不显著"升级为"统计上等价于0"
* ssc install ttost, replace
clear all
set seed 2316
set obs 2000
gen x = rnormal()
gen y = 0.05*x + rnormal(0,1)         // 真实效应很小(0.05)

* 常规回归：p 不显著
reg y x
* 看到 p=0.18，错误写法："政策无影响"

* TOST：设定等价界值 δ=0.1（即 |β|<0.1 视为无经济意义）
*   ttost 会做两个单侧检验：
*     H0a: β >= 0.1    vs H1a: β < 0.1
*     H0b: β <= -0.1   vs H1b: β > -0.1
*   两个都拒绝 → β 在 [-0.1, 0.1] 内 → 等价于0
* ttost coef=x, 0.1     // 或用回归后的 e() 直接做

* 关键：等价界值 0.1 必须在事前由经济意义定，
*      不能跑完 TOST 发现不通过，再换 δ=0.2 凑通过。
* 诚实写法:
*   "The 95% CI for β is [-0.03, 0.13], which lies
*    within our pre-specified equivalence bounds of
*    ±0.10; we therefore conclude the policy has at
*    most a negligible effect."

*=================================================================
* F1 【红线手法 · 严禁用于正式研究】
* 留一法回归 + 按系数排序删样本，直到核心变量显著
* 收录目的：让你认识它、在审稿/复现时能识破它
*=================================================================
ssc install rangestat
gen _i = 1
gen _t = 1
* 对每个观测做"剔除该观测后"的回归，得到留一法系数 b_x
* (by(_i) 分组、i(_t 0 0) 在单点上做窗口；excludeself 排除自身)
rangestat (reg) y x controls, by(_i) i(_t 0 0) excludeself

* --- 若要把系数"调显著为负"：按系数升序排列，删系数最小(最拖后腿)的那批 ---
gsort b_x
drop if _n<=20          // 一次删20个；不显著就继续删，直到显著
save data1.dta, replace

* --- 若要把系数"调显著为正"：按系数降序排列 ---
gsort -b_x
drop if _n<=120
save data2.dta, replace

* 识别特征(审稿人视角)：
*   - 删的数量(20/120)与"经济理由"毫无关系，纯由显著性决定
*   - 全样本回归不显著，删完才显著；删除模式单调对齐系数
*   - 复现代码里会出现 gsort b_x + drop if _n<=N 这种语句
* 诚实替代：Cook's D/DFBETA 只删"已确认录入错误"的观测，
*           并在附录列出删除清单与理由；样本定义事前预注册。