关于 Full-empirical-analysis-skill-Stata
面向应用经济学研究者的 Stata 技能，提供涵盖经济学、流行病学和因果机器学习的八步实证分析流程。重点是经过验证的可复现 .do 文件，以及面向论文的表格和图表。
主要功能
- 01八步实证分析流程
- 02经济学、流行病学和因果机器学习模式
- 03可复现的 Stata .do 文件
- 04诊断与稳健性检验
- 05导出表格和图表
适用场景
- 01开展 DID、IV、RD 或事件研究
- 02准备流行病学或公共卫生分析
- 03估计 DML 或因果森林模型
- 04为论文生成表格和图表
安装
复制下面的命令，在终端中运行。
$
npx skills add https:/安装后在 Claude Code 中输入/full-empirical-analysis-skill-stata确认已加载。
包含的文件
除 SKILL.md 外，安装时一起下载的参考资料和脚本。
README.mdreferences/04-statistical-tests.md08-tables-plots.md01-data-cleaning.md02-data-transformation.md07-further-analysis.md06-robustness.md03-descriptive-stats.md05-modeling.md详细介绍
完整实证分析——传统 Stata 工作流
这项技能是应用经济学家处理每篇实证论文时采用的规范性 8 步流程，基于传统 Stata 生态系统编写——原生 Stata，以及 20 多个已成为事实标准的社区命令（reghdfe、ivreg2、csdid、did_imputation、eventstudyinteract、sdid、rdrobust、rddensity、synth、synth_runner、psmatch2、teffects、ebalance、coefplot、esttab、outreg2、boottest、ritest、rwolf、bacondecomp、honestdid、binscatter）。
配套技能：如果用户希望使用 Python 实现相同流程，请转到 00-StatsPAI_skill（原生支持 Agent 的 DSL）或 00.1-Full-empirical-analysis-skill（明确的 Python 技术栈）。本技能是对应的 Stata 版本——每一步都会生成一个 .do 文件，可交给期刊的复现办公室，或交给不愿离开 Stata 的合作者。
理念
- 采用 Stata 惯用写法，而非翻译成 Python。 使用
reghdfe，而不是“reghdfe的 statsmodels 类比项”。使用esttab，而不是“Stata 的 stargazer”。 - 确保
.do文件可复现。 以下每个代码块都可以在运行use data.dta, clear后执行。没有 Jupyter，也没有 notebook——只有 do 文件和日志文件。 - 覆盖完整流程，而不只是回归。 Stata 用户历来过度投入第 5 步（建模），而对第 1–4 步和第 6–8 步投入不足。本技能将这些步骤同等对待。
- 提供丰富的输出。 每一步至少生成一个表格（
.tex/.rtf）或图形（.pdf/.png）——绝不只在 Results 窗口中打印一个系数，然后将其遗忘。 - 渐进式披露。
SKILL.md提供每一步的规范命令；references/则包含针对不同变体的深入内容（数十种检验、针对特定估计量的诊断、图形配方）。
SkillOpt 风格的执行关卡
将这份详细指南用作种子技能，而不是必须逐项穷尽执行的脚本。遵循 SkillOpt 原则：将每项局部 .do 文件改动视为候选补丁，只有通过筛选检查和留出检查后，才能将其转为可复用的样板代码。编写或修订 .do 文件前，先将用户的请求压缩成一张任务专属的 best_skill 卡片：
best_skill: <mode + design + artifact target>
train_signal: <current failure, user goal, or missing evidence>
selection_split: <focal dataset/spec/output used to judge the candidate>
heldout_gate: <checks the patch must pass beyond the focal example>
accepted_patterns: <rules to reuse after validation>
rejected_patterns: <failed shortcuts not to retry without new evidence>
patch_scope: <one estimator/sample/export/robustness change>
reject_if: <conditions that force rollback to the last passing spec>
- 路由卡片：记录模式（
econ、epi或ml-causal）、估计目标、识别设计、核心结果变量/处理变量、Stata 命令系列和所需产物。 - 限制改动范围：每次只更改一个决策（样本规则、估计量命令、固定效应、聚类、导出格式或稳健性检验）。优先采用能够通过验证的最小补丁。
- 遵守筛选拆分原则：将用户眼下遇到的故障或所需产物视为筛选集。至少留出一个备选结果变量、样本窗口、估计量命令系列或导出目标，作为留出关卡。
- 留出关卡：运行代码前定义检查项：
isid/重复记录状态、合并断言、处理变量支持情况、缺失值阈值、预期的表格/图形文件，以及一项不能被此次改动破坏的非核心稳健性检验/设定。 - 拒绝缓冲区：如果候选设定未通过关卡，就在
analysis_log.md或 Stata 日志中记录失败情况、.do差异和关卡输出；恢复到上一个通过的设定，不要重复尝试相同的未检查模式。 - 缓慢/元层面更新：任务结束时，根据执行过程记录
accepted_patterns和rejected_patterns。不要仅凭一次通过的运行就扩展规范.do模板。 - 仅在验证后推广：只有在当前数据以及至少一个备选结果变量/样本/设定上通过后，才能将一次性修复转为可复用的
.do样板代码。
三种领域模式（默认模式 = AER 经济学；备选模式 = 流行病学与 ML 因果推断）
上面的默认流程是 AER 风格的应用计量经济学——遵循 AEA 惯例：写出完整的估计方程、识别假设、设计方案对比，以及全面的稳健性检验。该技能还提供了另外两种因果推断传统对应的并行子流程，它们都沿用相同的步骤 1–4（清理 / 构造 / 表 1 / 诊断）和步骤 8（表格 / 图表）——只有步骤 5（估计量）以及步骤 6/7 会替换命令：
| 模式 | 读者惯例 | 步骤 5 估计量组合 | 报告组合 | 跳转到 |
|---|---|---|---|---|
| 默认模式 — 应用经济学（AER / QJE / AEJ） | “展示方程 + 识别假设 + 设计方案对比；控制变量清晰可见；聚类标准误” | DID / IV / RD / SCM / 匹配 / reghdfe HDFE | AER 风格的多列 esttab / outreg2 / coefplot + 8 部分论文布局 | 步骤 1 → 8（下方完整流程） |
| 模式 A — 流行病学 / 公共卫生 | “STROBE / TRIPOD-AI；目标试验方案；双重稳健估计目标；绝对与相对风险；KM 生存分析” | 目标试验模拟 · IPTW（teffects ipw）· IPWRA / AIPW（teffects ipwra / teffects aipw）· g-公式（gformula）· TMLE（eltmle）· 孟德尔随机化（mrrobust IVW / mregger / mrpresso）· KM/Cox/AFT（sts/stcox/streg/strmst2） | 相同的 esttab + 风险差 / 风险比 / E 值行 | §A. 流行病学流程 |
| 模式 B — ML 因果推断 | “DML / 元学习器 / 因果森林 / DR 学习器；CATE 分布；策略价值” | DML（ddml / pdslasso）· S/T/X/R/DR 学习器（ddml interactive）· GRF 因果森林（crforest / cforest）· BART / BCF（通过 python_user 调用外部 Python）· 矩阵补全（外部工具） | esttab ML 方案对比 + crforest CATE 图 + 策略价值表 | §B. ML 因果流程 |
如何调用非默认模式（Claude / Agent 会根据用户的措辞选择模式）：
| 用户说…… | 技能切换到的模式 |
|---|---|
| “运行 DID / IV / RD / 事件研究”、“AER 表格”、“应用微观经济学” | 默认模式（AER 经济学）— 步骤 1 → 8 |
| “目标试验模拟”、“g-公式”、“IPTW”、“TMLE”、“孟德尔随机化”、“STROBE / TRIPOD”、“公共健康 / 流行病学”、“epi 流程”、“RWE 研究”、“队列研究”、“病例对照” | 模式 A（流行病学）— §A |
| “DML”、“双重机器学习”、“ddml”、“因果森林”、“crforest”、“元学习器”、“CATE”、“策略学习”、“ML 因果推断”、“因果机器学习” | 模式 B（ML 因果推断）— §B |
| “混合”（例如“估计 DID，然后对异质性进行 ML CATE 分析”） | 依次运行默认模式 + 模式 B — 每个估计量都通过 eststo 存储结果，再将结果全部放入同一个 esttab 中，作为方案对比列 |
这三种模式共用 相同的步骤 1–4 清理 / 表 1 / 诊断框架、相同的步骤 8 导出组合，以及相同的 DAG 优先识别逻辑——切换模式只会改变步骤 5 所用的命令系列。如果你只需要描述性统计 / 表 1 / 平衡性检验，步骤 3 中的 AER tabstat / balancetable / asdoc 命令在三种模式下都同样适用。
模式 B 的 Stata 特别说明：Stata 原生的 ML 因果推断支持不如 Python/R 完善。对于 Dragonnet / TARNet / CEVAE / cfcausal / 公平性审计，可通过 Stata 18 的
python:/python script代码块调用 Python（或通过 shell 调用配套的.py文件），再通过frame或import delimited读回结果。该技能优先使用原生 Stata 命令（ddml、pdslasso、crforest），并在 §B 中明确标出 Python 调用。
默认输出规范 — 经济学实证论文
该技能默认遵循 应用经济学论文惯例。除非用户明确要求只给出单个点估计，否则每次运行 .do 文件都会生成以下完整的、可供发表的输出集。将其视为步骤 8 的约定——必须执行，而非可选项。
必须生成的表格（始终生成）
| # | 表格 | Stata 来源 | 保存到 |
|---|---|---|---|
| T1 | 汇总统计与平衡性（处理组 vs 对照组，含 SMD / p 值） | balancetable + asdoc sum（第 3 步） | tables/table1_balance.{tex,rtf,xlsx,docx} |
| T2 ★ | 主要结果 — 多列回归 M1→M6（逐步加入控制变量 + FE） | eststo 6 个设定 → esttab（第 5–6 步） | tables/table2_main.{tex,rtf,xlsx,docx} |
| T3 | 机制 / 结果阶梯 — 相同处理变量，3 个或更多结果并排展示 | 遍历各结果执行 eststo: reghdfe → esttab（第 7 步） | tables/table3_mechanism.{tex,rtf,xlsx,docx} |
| T4 | 异质性 — 子组 × 主要系数（性别、年龄、地区等） | 子组 eststo + suest Wald → esttab（第 7 步） | tables/table4_heterogeneity.{tex,rtf,xlsx,docx} |
| T5 | 稳健性检验组 — 替代 SE / 聚类 / 样本 / 安慰剂，放在同一张表中 | eststo × 各变体 → esttab（第 6 步） | tables/table5_robustness.{tex,rtf,xlsx,docx} |
★ 表 2 是每篇经济学论文的核心。 它是一张多列回归表，展示读者如何从原始相关性（M1）逐步走向完整设定的设计（M6：双向 FE + 交互 FE + 聚类稳健 SE）。不要将其合并为单列。不要只报告主要系数。这个逐步推进过程本身就是可信度论证：如果 M1→M6 的结果单调且稳定，说明设计很可能能够识别；如果加入 FE 后结果消失，那本身就是研究结果。
标准的 6 列，顺序如下：
- M1 原始双变量回归（
reg y treat）- M2 + 人口统计变量（
+ age + edu）- M3 + 行业控制变量（
+ tenure / firm_size）- M4 + 个体 FE（
reghdfe ..., absorb(unit)）- M5 + 双向 FE（
absorb(unit year)）- M6 + 交互 FE（
absorb(unit year i.industry#i.year)），使用vce(cluster unit)
必须生成的图表（始终生成）
| # | 图表 | Stata 来源 | 保存到 |
|---|---|---|---|
| F1 | 趋势 / 研究动机 — 处理组与对照组随时间变化的趋势，含政策线 | collapse (mean) y, by(year treat) → twoway line（第 3 步） | figures/fig1_trend.pdf（及 .png） |
| F2 | 事件研究系数及 95% CI，基期为 –1 | eventstudyinteract / csdid / coefplot, keep(*.rel)（第 5 步） | figures/fig2_event_study.pdf |
| F3 | M1→M6 各设定的系数图 | coefplot m1 m2 m3 m4 m5 m6, keep(treat) vertical（第 8 步） | figures/fig3_coefplot.pdf |
| F4 | 稳健性 / 敏感性曲线 — bacondecomp 图、honestdid 图或聚类比较森林图 | 依具体情形而定（第 6 步） | figures/fig4_sensitivity.pdf |
输出文件布局（默认）
project/
├── tables/ table1_balance.{tex,rtf,xlsx,docx} table2_main.{tex,rtf,xlsx,docx}
│ table3_mechanism.{tex,rtf,xlsx,docx} table4_heterogeneity.{tex,rtf,xlsx,docx}
│ table5_robustness.{tex,rtf,xlsx,docx}
└── figures/ fig1_trend.{pdf,png} fig2_event_study.{pdf,png}
fig3_coefplot.{pdf,png} fig4_sensitivity.{pdf,png}
每张表都生成 .tex（LaTeX booktabs）、.rtf（Word）、.xlsx（Excel）和 .docx（Word OOXML）文件。每张图都生成 .pdf（适用于 LaTeX 的矢量格式）和分辨率至少为 300 dpi 的 .png 文件（适用于幻灯片 / 网页）。
何时例外处理
- 单次快速估计 — 只生成相关单元格，但须提醒标准交付内容是上述完整套件，并提出可以执行完整套件。
- 研究设计不支持某张图（横截面研究 → 无事件研究图）— 跳过该图，并用
display打印说明原因的注释；不要默默省略。 - 处理单位为 N=1（
synth） — 用synth_runner生成的 SCM 轨迹图 + 安慰剂分布图替代 F1/F2；T1–T5 仍适用。
必需的软件包
* Run once on a fresh Stata install:
ssc install reghdfe, replace
ssc install ftools, replace // dependency of reghdfe / ivreg2
ssc install ivreg2, replace
ssc install ranktest, replace // dependency of ivreg2
ssc install ivreghdfe, replace // ivreg2 × reghdfe: high-dim FE IV
ssc install ppmlhdfe, replace // Poisson with HD FE
ssc install csdid, replace // Callaway–Sant'Anna (2021)
ssc install drdid, replace // dependency of csdid
ssc install did_imputation, replace // Borusyak–Jaravel–Spiess (2024)
ssc install eventstudyinteract, replace // Sun & Abraham (2021)
ssc install sdid, replace // Synthetic DID (Arkhangelsky et al. 2021)
ssc install did_multiplegt_dyn, replace // de Chaisemartin & D'Haultfœuille
ssc install bacondecomp, replace // Goodman-Bacon (2021)
ssc install honestdid, replace // Rambachan–Roth (2023) PT sensitivity
ssc install rdrobust, replace // Calonico–Cattaneo–Titiunik RD
ssc install rddensity, replace // McCrary / Cattaneo et al. density test
ssc install synth, replace // Abadie–Diamond–Hainmueller SCM
ssc install synth_runner, replace // SCM with placebos + inference
ssc install psmatch2, replace // propensity-score matching
ssc install ebalance, replace // entropy balancing
ssc install coefplot, replace
ssc install estout, replace // provides estout / esttab / eststo
ssc install outreg2, replace
ssc install asdoc, replace // one-click Word/Excel tables
ssc install binscatter, replace
ssc install balancetable, replace
ssc install winsor2, replace
ssc install xtable, replace // better xtreg output tables
ssc install boottest, replace // wild cluster bootstrap (Roodman et al.)
ssc install ritest, replace // randomization inference
ssc install rwolf, replace // Romano–Wolf multiple-testing
ssc install moremata, replace // Mata extensions (dep for several)
ssc install mdesc, replace // missing data description
ssc install missings, replace // missings dropvars, report
ssc install unique, replace // unique IDs in panel
ssc install schemepack, replace // modern publication themes
8 个步骤——标准流程（映射到 AER 论文各章节）
┌──────────────────────────────────────────────────────────────────────┐
│ Step −1 Pre-Analysis Plan (PAP) power/sampsi/clustersampsi/MDE │
│ Step 0 Sample log + data contract sample_log/assert/xtdescribe/JSON │
│ Step 1 Data import & cleaning use/import/destring/misstable/merge│
│ Step 2 Variable construction gen/egen/winsor2/xtile/xtset/L.F.D.│
│ Step 2.5 Empirical strategy equation × ID assumption + pre-reg │
│ Step 3 Descriptive statistics tabstat/balancetable/asdoc/pwcorr │
│ Step 3.5 Identification graphics event-study/1st-stage/McCrary/love │
│ Step 4 Diagnostic tests sktest/hettest/xtserial/vif/dfuller│
│ Step 5 Baseline modeling reghdfe/ivreg2/csdid/rdrobust/synth│
│ Step 6 Robustness battery bacondecomp/honestdid/rwolf/boottest│
│ Step 7 Further analysis triple-diff/subgroup/medsem/margins│
│ Step 8 Tables & figures esttab/outreg2/coefplot/rdplot │
└──────────────────────────────────────────────────────────────────────┘
这 8 个步骤对应应用微观经济学 AER / QJE / AEJ 论文的标准章节。每个步骤对应一个论文章节，并在磁盘上生成一个可直接用于论文的产出：
Paper section Step Stata moves
─────────────────────────── ───── ────────────────────────────────────────────────
Pre-Analysis Plan −1 power/sampsi + freeze protocol.do to disk
§1. Data 0 sample_log + 5-check data contract → JSON
§1. Data 1 use/import/destring/misstable/merge assert/xtset
§1. Data 2 gen/egen/winsor2/xtile/L./F./D./CPI deflation
§1.1 Descriptives (Table 1) 3 tabstat · balancetable · asdoc · pwcorr · twoway
§2. Empirical Strategy 2.5 write equation + ID assumption → strategy.do
§3. Identification graphics 3.5 event-study · 1st-stage F · McCrary · love · SCM
§3.5 Diagnostics 4 swilk · hettest · xtserial · vif · dfuller · hausman
§4. Main Results (Table 2) 5 M1→M6 progressive controls + FE (eststo + esttab)
§5. Heterogeneity (Table 3) 7 margins/marginsplot · subgroup · medsem
§6. Mechanisms / Channels 7 medsem/khb · outcome ladder · DDD interactions
§7. Robustness gauntlet 6 bacondecomp · honestdid · psacalc · boottest · ritest · rwolf
§8. Replication package 8 esttab + outreg2 + coefplot + reproducibility stamp
下面列出每个步骤的标准命令。所有示例都围绕同一个叙述展开——一个劳动经济学面板数据案例，其中 training（处理变量）影响 log_wage（结果变量），协变量为 age、edu、tenure，面板键为 worker_id / firm_id / year。变量名和参数值均为示例；请替换为用户数据集中的实际值。只有命令名称和选项的形式是规范的。
如果某个步骤有多个变体（例如，交错 DID 有 5 种估计量；异方差有 4 种经典检验），SKILL.md 会介绍最常用的一种，并链接到
references/NN-<topic>.md以查看其余变体。如果用户的情形不适用默认方案，请阅读参考文件。
可直接用于论文的图表清单（按章节列出）
现代 AER 论文通常有 5–7 张图和 3–5 张主要表格，另加一张附录稳健性表。以下每个步骤都会在磁盘上生成至少一个带编号的产出。默认文件名假定同时导出 .tex / .rtf 格式（Agent 应同时生成两种格式，以便合作者在 Word 中编辑，并供构建系统使用 LaTeX）：
| § | 产出 | Stata 基本命令 | 文件名 |
|---|---|---|---|
| §1 | 图 1：原始趋势 / 处理推广情况 | collapse + twoway line · 交错推广情况使用 heatplot | figures/fig1_trend.{pdf,png} |
| §1 | 表 1：汇总统计（全部 / 处理组 / 对照组 + Δ + SMD） | balancetable · asdoc sum, by() · tabstat | tables/table1_balance.{tex,rtf,xlsx,docx} |
| §3 | 图 2：识别图示（事件研究 / 第一阶段 / McCrary / RD 散点图 / SCM 轨迹） | 在 eventstudyinteract/csdid 后使用 coefplot · binscatter · rdplot · rddensity, plot · synth | figures/fig2_event_study.{pdf,png} |
| §4 | 表 2：主要结果——逐步增加控制变量 M1→M6 | eststo 运行 6 个规格 → esttab | tables/table2_main.{tex,rtf,xlsx,docx} |
| §4 | 表 2-bis：研究设计比较（OLS / IV / DID / 匹配） | 混合使用 eststo + esttab | tables/table2b_designs.{tex,rtf,xlsx,docx} |
| §4 | 图 3：不同规格的系数图 | coefplot m1 m2 m3 m4 m5 m6, keep(treat) | figures/fig3_coefplot.{pdf,png} |
| §5 | 表 3：按子组划分的异质性 | 每个子样本使用 eststo + esttab + suest Wald | tables/table3_heterogeneity.{tex,rtf,xlsx,docx} |
| §5 | 图 4：剂量-反应 / 按四分位数划分的边际效应 | xtile + margins + marginsplot | figures/fig4_dose.{pdf,png} |
| §6 | 表 4：机制 / 结果变量阶梯 | 循环遍历结果变量，运行 eststo: reghdfe → esttab | tables/table4_mechanism.{tex,rtf,xlsx,docx} |
| §7 | 表 A1：稳健性汇总表（每项检验占一列） | eststo × 各种变体 → esttab | tables/tableA1_robustness.{tex,rtf,xlsx,docx} |
| §7 | 图 5：规格曲线——所有规格的系数 + 95% 置信区间 | 手动编写规格循环 + twoway rcap | figures/fig5_spec_curve.{pdf,png} |
| §7 | 图 6：敏感性分析（HonestDiD / Oster / E-value） | honestdid, coefplot · psacalc plot · evalue 表 | figures/fig6_sensitivity.{pdf,png} |
| §8 | 复现包：将所有表格汇总到一个文档中 | esttab ..., append 输出到一个 .tex / .rtf 文件 · texdoc | replication/paper_tables.{tex,rtf,xlsx,docx} |
上述每种 Stata 估计量都会通过
eststo存储结果，可直接传给esttab/coefplot/outreg2。不要手动编写 LaTeX，也不要用outsheet/putexcel矩阵生成 Word 文档——esttab和outreg2会自动应用书籍式表格边框、AER 风格的星号以及正确的标准误标签。如需了解更深入的导出方法（LaTeX / Word / Markdown 格式、texdoc、frmttable），请参阅references/08-tables-plots.md。
导出手册 — 在一个代码块中导出 LaTeX / Word / RTF
Stata 的导出工具链比 Python 的 StatsPAI 更零散。分为三个层级，按范围选择：
| 层级 | 适用场景 | API | 常用选项 |
|---|---|---|---|
| 1. 单个多列的表格 | 导出带逐步增加列的单个表 2 / 表 3 / 表 A1 | .tex / .rtf：使用 booktabs 的 esttab.xlsx / .docx：outreg2（esttab 原生的 Office 导出只是数据转储工具，而非出版级格式化工具；outreg2 可生成格式规范的 Word/Excel 文件，带边框、对齐的星号和 AER 风格布局） | 用于 tex/rtf 的 esttab：keep()、drop()、mtitles()、stats(N r2 r2_a, labels(...))、star(* 0.10 ** 0.05 *** 0.01)、label、booktabs、addnotes()。用于 xlsx/docx 的 outreg2：label dec(3)，使用 addtext() 添加固定效应指示项；第一列使用 replace，后续列使用 append。始终输出全部四种格式：在 foreach ext in tex rtf 循环中用 esttab 输出 tex/rtf，在单独的 foreach ext in xlsx docx 循环中用 outreg2 输出 xlsx/docx。 |
| 2. 多面板论文格式（表 2 + 表 3 + 表 A1 + 表 A2 放在一个文件中） | 生成论文表格合集——将主结果 + 异质性 + 稳健性 + 安慰剂检验作为单个文档 | 第一个面板使用 esttab ... using "paper.tex", replace；后续每个面板使用 esttab ... using "paper.tex", append；使用 texdoc init "paper.tex" 生成完整的 LaTeX 文档（含正文）。对 .rtf / .xlsx / .docx 文件组重复此操作。 | 第一个面板：replace；后续面板：append；用 texdoc 添加标题 |
3. 完整会话文件组（相当于 Stata 17+ 的 collect） | 生成复现附录，将汇总统计 + 平衡性检验 + 多个回归表 + 标题 + 正文混合在一个文件中 | collect create papercollect get summary, ...collect get est ...collect layout ...collect export "paper.xlsx"（也支持 .docx / .html / .tex / .md） | 仅适用于 Stata 17+；较旧版本的 Stata 使用 texdoc / markdoc |
期刊样式——选择合适的 star 显著性水平和 SE 标签。 AEA 惯例为 * 0.10 ** 0.05 *** 0.01，SE 标签为 “括号内为标准误”；QJE / Econometrica / RES 的不同版本仅在星号 / 注释 / 字体方面有所区别。在 do 文件开头定义一个 esttab 封装程序：
* Top of master.do — journal house-style wrapper
local AER_STAR "* 0.10 ** 0.05 *** 0.01"
local AER_NOTES "Cluster-robust standard errors in parentheses. * p<0.10, ** p<0.05, *** p<0.01."
local AER_STATS stats(N r2_a, labels("N" "Adj. R²"))
* LaTeX + RTF via esttab (booktabs for tex)
* esttab m1 m2 m3 using "tables/table2.tex", replace ///
* se star(`AER_STAR') label booktabs `AER_STATS' addnotes(`AER_NOTES')
* esttab m1 m2 m3 using "tables/table2.rtf", replace ///
* se star(`AER_STAR') label `AER_STATS' addnotes(`AER_NOTES')
* Excel + Word via outreg2 (publication-grade Office formatting)
* First column: replace; subsequent: append
* outreg2 using "tables/table2.xlsx", replace label dec(3) ///
* keep(training age edu tenure) addtext(Worker FE, Yes, Year FE, Yes)
* outreg2 using "tables/table2.xlsx", append label dec(3) ///
* keep(training age edu tenure region) addtext(Worker FE, Yes, Year FE, Yes, Region, Yes)
多格式导出模式——通过 esttab 输出 .tex 和 .rtf，再通过 outreg2 输出 .xlsx 和 .docx。切勿使用 esttab ... using "... .xlsx" 或 esttab ... using "... .docx"——esttab 原生的 Office 导出只是原始数据转储工具，并非出版级格式化工具：
* LaTeX + RTF — use esttab (booktabs for tex, rich text for rtf)
foreach ext in tex rtf {
esttab m1 m2 m3 using "tables/table2_main.`ext'", ///
replace se star(* 0.10 ** 0.05 *** 0.01) ///
label booktabs mtitles("(1)" "(2)" "(3)") ///
stats(N r2_a, labels("N" "Adj. R²")) ///
addnotes("Cluster-robust SE in parentheses. * p<0.10, ** p<0.05, *** p<0.01.")
}
* Excel + Word — use outreg2 (proper Office table formatting, aligned stars, borders)
* First model: replace
eststo m1: qui reghdfe log_wage training age edu tenure, absorb(worker_id year) vce(cluster worker_id)
outreg2 using "tables/table2_main.xlsx", replace label dec(3) ///
keep(training age edu tenure) ///
addtext(Worker FE, Yes, Year FE, Yes)
* Subsequent models: append
eststo m2: qui reghdfe log_wage training age edu tenure region, absorb(worker_id year) vce(cluster worker_id)
outreg2 using "tables/table2_main.xlsx", append label dec(3) ///
keep(training age edu tenure region) ///
addtext(Worker FE, Yes, Year FE, Yes, Region, Yes)
* docx follows the same pattern but with the .doc extension
outreg2 using "tables/table2_main.doc", replace label dec(3) ///
keep(training age edu tenure) ///
addtext(Worker FE, Yes, Year FE, Yes)
eststo m2
outreg2 using "tables/table2_main.doc", append label dec(3) ///
keep(training age edu tenure region) ///
addtext(Worker FE, Yes, Year FE, Yes, Region, Yes)
图表始终同时导出为 .pdf 和 .png，分辨率至少为 300 dpi（LaTeX 使用矢量格式，幻灯片 / 网页 / Word 嵌入使用栅格格式）：
graph export "figures/fig1_trend.pdf", replace
graph export "figures/fig1_trend.png", replace width(2400) height(1800)
关于 collect / texdoc / markdoc 手册（Stata 17+ 多面板论文文件组，以及生成正文加表格的 PDF/Word 文档），请参阅 references/08-tables-plots.md。
步骤 −1 — 预分析计划（数据分析前；AEA RCT Registry 风格）
在接触数据之前，写下 (a) 研究总体，(b) 研究设计，(c) 在计划样本量和 α=0.05、β=0.20 下的最小可检测效应（MDE）。将结果保存为 protocol.do，以便审稿人验证研究设计的统计功效是在查看数据之前而非之后确定的。
* Two-sample MDE (continuous outcome, Cohen's d framing)
power twomeans 0, diff(0.2) sd(1) power(0.80) alpha(0.05)
* → required n per arm
* Cluster-randomized RCT — solve for n_clusters given ICC
power twomeans 0, diff(0.2) sd(1) power(0.80) k1(50) rho(0.05) cluster
* → required clusters per arm under ICC=0.05, cluster size 50
* DID power (Frison–Pocock / Bloom 1995): use -sampsi- + -clustersampsi-
* (no native `power did`; for staggered DID see references/05-modeling.md §5.4)
sampsi 0 0.15, sd(1.0) alpha(0.05) power(0.80) n1(.) n2(.)
* RD power — McCrary-style: solve via Monte Carlo with -simulate-
* (see references/05-modeling.md §5.5 RD power template)
* Persist the protocol — referee will ask whether design was powered ex ante
file open f using "protocol.do", write replace
file write f "* Pre-analysis plan — frozen `c(current_date)' `c(current_time)'" _n
file write f "* Population: manufacturing workers 2010-2020" _n
file write f "* Treatment: training (binary)" _n
file write f "* Outcome: log_wage" _n
file write f "* Design: staggered DID, csdid (Callaway-Sant'Anna 2021)" _n
file write f "* MDE: 0.05 log points at 80% power, α=0.05" _n
file write f "* N planned: ~12,000 worker-years" _n
file close f
在运行步骤 1 之前，将此 .do 文件保存到版本控制中。AEA RCT Registry / OSF 预注册工具接受将其作为分析计划附件。
步骤 0 — 样本构建日志与 5 项数据契约检查
AER 第 1 节 Data 有三项任务：(a) 描述数据来源，(b) 记录所有样本限制（即“脚注 4”样本日志），(c) 锁定面板结构。Stata 用户通常会在单个 do 文件中连续执行 keep / drop 命令，却不记录计数，之后就无法为回复信还原分析样本。数据契约正是解决办法。
0.1 样本构建日志（脚注 4）
use "raw.dta", clear
matrix sample_log = J(0, 2, .)
local row 0
local n0 = _N
local ++row
matrix sample_log = (nullmat(sample_log) \ `row', `n0')
display "Step 0. raw: N = " %12.0fc `n0'
drop if missing(wage)
local n1 = _N
local ++row
matrix sample_log = (nullmat(sample_log) \ `row', `n1')
display "Step 1. drop missing wage: N = " %12.0fc `n1' " (Δ " %10.0fc `n0' - `n1' ")"
drop if !inrange(age, 18, 65)
local n2 = _N
local ++row
matrix sample_log = (nullmat(sample_log) \ `row', `n2')
display "Step 2. drop age outside 18-65: N = " %12.0fc `n2' " (Δ " %10.0fc `n1' - `n2' ")"
keep if inlist(industry, "manuf", "construction", "transport")
local n3 = _N
local ++row
matrix sample_log = (nullmat(sample_log) \ `row', `n3')
display "Step 3. keep target industries: N = " %12.0fc `n3' " (Δ " %10.0fc `n2' - `n3' ")"
* Persist to JSON (for paste into footnote 4)
file open f using "artifacts/sample_construction.json", write replace
file write f "{" _n
file write f `" "step_0_raw": `=sample_log[1,2]',"' _n
file write f `" "step_1_wage": `=sample_log[2,2]',"' _n
file write f `" "step_2_age": `=sample_log[3,2]',"' _n
file write f `" "step_3_indust": `=sample_log[4,2]'"' _n
file write f "}" _n
file close f
将 display 行原样粘贴为论文的脚注 4。
0.2 五项数据契约检查（通过／不通过关卡）
* (1) Shape
display "Check 1. n_obs = " _N
* (2) Dtypes on key vars — assert types are numeric where required
foreach v in wage training worker_id year age edu tenure {
capture confirm numeric variable `v'
if _rc {
display as error "Check 2 FAILED: `v' is not numeric — fix before any panel command"
exit 198
}
}
display "Check 2. dtypes OK on all key vars"
* (3) Missingness pattern on key vars
mdesc wage training worker_id year age edu tenure // ssc install mdesc
foreach v in wage training worker_id year {
qui count if missing(`v')
if r(N) > 0 {
display as error "Check 3 FAILED: `v' has " r(N) " missing — fix before estimation"
exit 198
}
}
* (4) Duplicate (id, time) — fatal for panel methods
duplicates report worker_id year
duplicates tag worker_id year, gen(dup)
qui count if dup > 0
if r(N) > 0 {
display as error "Check 4 FAILED: " r(N) " duplicate (worker_id, year) rows"
exit 198
}
drop dup
* (5) Panel balance
xtset worker_id year
xtdescribe
qui xtdes
local n_balanced = r(N)
display "Check 5. panel n=" _N " (balanced cells if `n_balanced' equals expected unit×period)"
* MCAR sniff test (Rubin) — if missing(y) is associated with covariates,
* listwise deletion biases the estimate. Use `mi` / IPW instead.
gen byte miss_y = missing(wage)
foreach cov in age edu tenure {
qui ttest `cov', by(miss_y)
if r(p) < 0.05 {
display as error "WARNING: y-missingness associates with `cov' (p=" %5.3f r(p) ")."
display as error " -> NOT MCAR, use -mi impute- or IPW, NOT listwise drop."
}
}
drop miss_y
* Persist contract
file open f using "artifacts/data_contract.json", write replace
file write f "{" _n
file write f `" "n_obs": `=_N',"' _n
file write f `" "panel": "worker_id × year","' _n
file write f `" "treatment": "training","' _n
file write f `" "outcome": "log_wage""' _n
file write f "}" _n
file close f
如果任何断言触发，请停止并修复问题。Stata 估计器会静默丢弃 NaN 行，这是回复信中“样本量莫名缩减”问题最常见的来源。
步骤 1 — 数据导入与清理
更深入的模式：references/01-data-cleaning.md — 读取 Excel/CSV/SAS/SPSS，对看似数值的字符串使用 destring，检查 misstable 模式，使用 duplicates 标记重复项，使用 assert(match using) 执行 merge，检查 xtset 平衡性，检查时期／间隔，添加标签。
* 1a. Load + first look
use "raw.dta", clear
describe, short
summarize
misstable summarize
mdesc // missing-data report
* 1b. Dtypes — destring strings-that-should-be-numeric
destring year wage, replace force // force: convert non-numeric to .
gen hire_date = date(hire_date_str, "YMD"); format hire_date %td
* 1c. Missing values — decide PER VARIABLE
local key_vars "wage training worker_id year"
foreach v of local key_vars {
drop if missing(`v')
}
sum tenure, detail
replace tenure = r(p50) if missing(tenure) // median-impute
gen byte tenure_miss = missing(tenure) // keep the flag
* 1d. Outliers — flag first (winsorize in Step 2)
egen wage_z = std(wage)
count if abs(wage_z) > 4
display "Flagged |z|>4 on wage: " r(N)
* 1e. Deduplicate on panel key
duplicates report worker_id year
duplicates tag worker_id year, gen(dup)
assert dup == 0 // hard-fail if panel key not unique
drop dup
* 1f. Merge with assert — never silently lose rows
merge m:1 firm_id using "firm_chars.dta", ///
assert(match using master) keep(master match) nogen
* 1g. Panel structure
xtset worker_id year // declares panel
xtdescribe // balance summary
tab year // observations per year
核心原则：明确记录所有行排除操作，并统计数量、写入日志。步骤 2 及之后的命令都不应静默丢弃行。
步骤 2 — 变量构造与转换
更深入的模式：references/02-data-transformation.md — log/ihs/Box–Cox，组内 winsor2，xtile 和自定义分组边界，egen 用法，时间序列运算符（L.、F.、D.、S.），CPI 平减，分期 DID 时点构造。
* 2a. Log / IHS
gen log_wage = log(max(wage, 1)) // floor at 1
gen ihs_assets = asinh(assets) // handles 0 / negative
* 2b. Winsorize 1/99
winsor2 wage, cuts(1 99) suffix(_w1) by(year) // within-year winsorize
* 2c. Standardize
egen age_std = std(age)
* 2d. Categorical encoding (factor variables — use i. inside regressions)
* Explicit dummies only when needed for export
tab industry, gen(ind_)
drop ind_1 // base category
* 2e. Interactions & polynomials — use c. and i. inline in reg commands
gen age_sq = age^2
gen trt_x_edu = training * edu
* 2f. Panel operators (xtset is required for L./F./D. to work)
xtset worker_id year
gen log_wage_l1 = L.log_wage
gen log_wage_f1 = F.log_wage
gen d_log_wage = D.log_wage
* 2g. Within-unit mean (egen)
egen wage_mean_i = mean(log_wage), by(worker_id)
* 2h. Treatment timing for staggered DID
bysort worker_id (year): egen first_treat = min(cond(training==1, year, .))
gen rel_time = year - first_treat
replace rel_time = . if missing(first_treat) // never-treated = .
gen never_treated = missing(first_treat)
* 2i. Real values (CPI deflation)
merge m:1 year using "cpi.dta", keep(master match) nogen
sum cpi if year == 2010
gen wage_real = wage * r(mean) / cpi
步骤 2.5 — 实证策略（写出方程和识别假设）
这是 AER 论文的核心。编写任何代码之前，先明确写出方程并陈述识别假设。含糊的识别表述是审稿人拒绝应用论文最常见的原因。将策略保存为 strategy.do（或 .md），作为带日期、受版本控制的文档——而不是看过系数后事后编写的合理化解释。
方程 × 识别假设 × Stata 估计器（决策表）
| 设计 | 估计方程 | 识别假设 | Stata 估计器 |
|---|---|---|---|
| 2×2 DID | Y_it = α_i + λ_t + β·D_it + X'γ + ε_it | 给定 X 后满足平行趋势 | reghdfe Y c.D#c.post X, absorb(i t) vce(cluster i) |
| 事件研究（CS / SA） | Y_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t-G_i = e} + ε_it | 无预期效应 + 组别—时间 PT | csdid / eventstudyinteract / did_imputation |
| 2SLS | Y_i = α + β·D_i + X'γ + ε_i; D_i = π·Z_i + X'δ + u_i | 排除限制 + 相关性 + 单调性 | ivreg2 / ivreghdfe |
| Sharp RD | Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i（局部多项式） | E[Y(0)|X] 在 c 处连续，且无操纵 | rdrobust（+ rddensity） |
| SCM | Ŷ_1t(0) = Σ_j ŵ_j Y_jt，τ_t = Y_1t − Ŷ_1t(0)，适用于 t≥T_0 | 处理前拟合 + 插值有效性 | synth / synth_runner / sdid |
| 基于可观测变量的选择（匹配/IPW） | E[Y_i(d) | X_i] = E[Y_i | D=d, X_i] | 无混杂 + 重叠性 | teffects psmatch / ipw / aipw / ipwra，ebalance |
设计选择器（用户不确定时）
┌─ running var + cutoff ───────────────── RDD (rdrobust)
│
├─ exogenous instrument Z ─────────────── IV/2SLS (ivreg2 / ivreghdfe)
data + question ─┤
├─ pre/post × treat/control ─┬ 2 periods ── 2×2 DID (reghdfe / xtreg)
│ └ staggered ── CS / SA / BJS (csdid / eventstudyinteract / did_imputation)
│
├─ 1 treated unit + donor pool + long pre ── SCM (synth / synth_runner / sdid)
│
├─ high-dim X, selection-on-observables ── ML causal (ddml / crforest — see §B)
│
└─ none of the above ──────────────────── matching + sensitivity (teffects + evalue)
预注册 strategy.do 模板
*--- strategy.do — empirical-strategy pre-registration ---
* Frozen `c(current_date)' `c(current_time)' (Git SHA: <paste>)
* Population: manufacturing workers, 2010–2020, balanced panel
* Treatment: training (binary, staggered adoption)
* Outcome: log_wage (CPI-deflated 2010 USD)
* Estimand: ATT on the treated, dynamic horizon -4..+4
*
* Estimating equation (paste from the §2.5 row that matches the design):
* log_wage_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t - G_i = e} + ε_it
*
* Identifying assumption:
* 1. No anticipation: E[Y_it(0) | t < G_i] = E[Y_it(0) | never-treated]
* 2. Group-time PT: Δ E[Y_it(0)] is the same across treatment cohorts
*
* Auto-flagged threats (must defend in §2):
* - Selection of G_i on Y_i(0) -> bacondecomp + honestdid sensitivity
* - Spillover across worker_id within firm -> cluster at firm_id, also try firm_id × year
* - Anticipation in the year before adoption -> include lead in event study
*
* Fallback estimators (Step 6 robustness):
* - eventstudyinteract (Sun-Abraham 2021)
* - did_imputation (Borusyak-Jaravel-Spiess 2024)
* - sdid (Synthetic DID)
在运行步骤 5／步骤 6 之前提交 strategy.do。此文件的 git log 就是分析计划。
步骤 3 — 描述性统计与表 1
更深入的模式：references/03-descriptive-stats.md — 带 SMD 的分层表 1，asdoc / tabstat / balancetable，带显著性星号的相关矩阵（pwcorr, sig star(.05)），按组绘制直方图／kdensity，使用 xtline 展示 DID 动机，使用面板覆盖率 xtdescribe。
* 3a. Full-sample summary
local vars "log_wage age edu tenure training"
tabstat `vars', statistics(n mean sd min p25 p50 p75 max) columns(statistics)
* One-command Word/Excel output:
asdoc sum `vars', stat(N mean sd min median max) ///
save(tables/table1_full.docx) replace
* 3b. Stratified Table 1 (treated vs control + t-tests + SMDs)
balancetable training age edu tenure female ///
using "tables/table1_balance.tex", ///
vce(cluster firm_id) replace ///
varlabels pval
* Manual per-variable t-test + SMD:
foreach v of varlist age edu tenure {
qui sum `v' if training == 1
local m1 = r(mean); local sd1 = r(sd); local n1 = r(N)
qui sum `v' if training == 0
local m0 = r(mean); local sd0 = r(sd); local n0 = r(N)
local smd = (`m1' - `m0') / sqrt((`sd1'^2 + `sd0'^2)/2)
ttest `v', by(training)
display "`v': Δ=" %7.3f (`m1'-`m0') " SMD=" %7.3f `smd' " p=" %6.3f r(p)
}
* 3c. Correlation matrix with significance stars
pwcorr `vars', sig star(.05)
estout using "tables/corr.tex", replace ///
cells("b(star fmt(3))") style(tex) // alternative: estpost correlate
* 3d. Distribution plots
twoway (kdensity log_wage if training==1) ///
(kdensity log_wage if training==0), ///
legend(order(1 "Treated" 2 "Control")) ///
title("Log wage density by treatment") ///
saving(figures/kde_wage, replace)
graph export "figures/kde_wage.pdf", replace
* 3e. Time trends — the DID motivation plot
preserve
collapse (mean) log_wage, by(year training)
twoway (line log_wage year if training==1) ///
(line log_wage year if training==0), ///
xline(`policy_year', lpattern(dash)) ///
legend(order(1 "Treated" 2 "Control"))
graph export "figures/trend_did.pdf", replace
restore
* 3f. Panel coverage
xtdescribe
步骤 3.5 — 识别图形（“Identification, graphical evidence”章节）
AER 惯例：识别图应放在回归表之前。 在要求读者相信点估计之前，应先让他们看到图形证据，证明 PT 成立／第一阶段足够强／RD 跳跃清晰。此部分包含一到两张保存至 figures/ 的图——大部分工作只需一条 Stata 命令和一个 coefplot。
3.5.1 事件研究图 + 数值前趋势检验（DID 识别）
处理前系数 ≈ 0（将 −1 期作为参考期并归一化为零）是平行趋势的图形证据。应将图形与数值前趋势检验结合起来，避免审稿人只能靠目测判断。
* Build relative-time factor variable, base at e = -1
gen rel = year - first_treat
replace rel = . if missing(first_treat) // never-treated dropped from event study
keep if inrange(rel, -4, 4) | missing(first_treat)
gen rel_p = rel + 5 // shift -4..4 → 1..9; ib4 means "base = -1"
* (a) Sun-Abraham via -eventstudyinteract- (preferred for staggered adoption)
forvalues k = 1/9 {
gen rel_d`k' = (rel_p == `k')
}
eventstudyinteract log_wage rel_d1 rel_d2 rel_d3 rel_d5 rel_d6 rel_d7 rel_d8 rel_d9, ///
cohort(first_treat) control_cohort(never_treated) ///
absorb(worker_id year) vce(cluster worker_id)
* (b) Coefficient figure with shaded CI + reference line at e=-1
coefplot, keep(rel_d*) vertical omitted ///
yline(0, lpattern(dash) lcolor(gs10)) xline(4.5, lpattern(dash)) ///
rename(rel_d1="-4" rel_d2="-3" rel_d3="-2" rel_d5="0" ///
rel_d6="1" rel_d7="2" rel_d8="3" rel_d9="4") ///
ciopts(recast(rcap)) levels(95) ///
xtitle("Years relative to treatment") ///
ytitle("Coefficient (ATT, 95% CI)") ///
title("Figure 2a. Event-study coefficients (95% CI; ref. e = -1)") ///
scheme(s2color)
graph export "figures/fig2a_event_study.pdf", replace
graph export "figures/fig2a_event_study.png", replace width(2400)
* (c) Numerical pre-trends F-test (joint zero on the leads e = -4..-2)
test rel_d1 rel_d2 rel_d3
display "Pre-trends F = " %5.2f r(F) " p = " %5.3f r(p)
* (d) Bacon decomposition figure (Goodman-Bacon 2021) — TWFE diagnostic
bacondecomp log_wage training, ddetail
* The -bacondecomp- output names the contaminated 2×2 weights.
* Save the auto-generated figure as figures/fig2a_bacon.pdf.
* (e) Callaway-Sant'Anna dynamic ATT (when -csdid- is the main estimator)
csdid log_wage age edu, ivar(worker_id) time(year) gvar(first_treat) ///
method(dripw) agg(event)
estat event, window(-4 4)
csdid_plot, title("Figure 2a-bis. Dynamic ATT (Callaway-Sant'Anna)")
graph export "figures/fig2a_csdid.pdf", replace
3.5.2 第一阶段 F 统计量 + 散点图（IV 识别）
经验法则：OLS 风格推断要求第一阶段 F ≥ 10；AR 等效推断要求 F ≥ 23（Stock–Yogo / Lee 2022）。ivreg2 默认报告 CD / KP / 弱 IV 统计量，比 ivregress 的精简输出实用得多。
ivreg2 log_wage age edu (training = Z1 Z2), cluster(firm_id) first endog(training)
* Look for the line "Cragg-Donald Wald F" and "Kleibergen-Paap rk Wald F" —
* and "Anderson-Rubin Wald test" for weak-IV-robust CIs.
* First-stage scatter (binscatter — residualizes age + edu)
binscatter training Z1, controls(age edu) nquantiles(20) ///
xtitle("Excluded instrument Z1") ytitle("Pr(training)") ///
title("Figure 2b. First-stage relationship (residualized)") ///
savegraph("figures/fig2b_first_stage.pdf") replace
3.5.3 RD：McCrary 密度检验 + 经典 RD 图
RD 的标志性图表是 rdplot（CCT 风格的分箱散点图，并在两侧分别拟合局部多项式），并配合 McCrary 操纵检验。两者共同回答：（a）是否存在明显的跳跃？（b）密度在截断点处是否连续？
* (a) Canonical RD plot
rdplot outcome running_var, c(0) p(4) kernel(triangular) binselect(esmv) ///
graph_options(title("Figure 2c. RD plot") ///
ytitle("Outcome") xtitle("Running variable") ///
scheme(s2color))
graph export "figures/fig2c_rdplot.pdf", replace
* (b) McCrary density (manipulation test) — Cattaneo–Jansson–Ma 2018
rddensity running_var, c(0) plot ///
plot_options(title("Figure 2c-bis. McCrary density (manipulation test)"))
graph export "figures/fig2c_mccrary.pdf", replace
* (c) Covariate-adjusted continuity test (continuity of *covariates* at c)
foreach v of varlist age edu tenure {
rdrobust `v' running_var, c(0)
}
3.5.4 匹配：love plot（标准化差异的匹配前与匹配后对比）
teffects psmatch (log_wage) (training age edu tenure), atet
tebalance summarize // table of pre/post SMDs
tebalance density ps // density overlap plot — save as figures/fig2d_overlap.pdf
graph export "figures/fig2d_overlap.pdf", replace
* Or use -psmatch2- + -pstest- for the canonical love plot
psmatch2 training age edu tenure, out(log_wage) n(1) common
pstest age edu tenure, both graph // pre/post |std diff| with target |Δ|<10%
graph export "figures/fig2d_loveplot.pdf", replace
3.5.5 SCM：合成控制轨迹图 + 差距图
对于合成控制设计，经典的图 2 是处理组与合成组的时间序列图，并标注处理时间。
synth log_wage age edu tenure, ///
trunit(1) trperiod(2015) fig keep("synth_results.dta", replace)
graph export "figures/fig2e_synth_trajectory.pdf", replace
* Synthetic DID variant
sdid log_wage worker_id year training, vce(bootstrap) graph g1on
graph export "figures/fig2e_sdid.pdf", replace
* Placebo gap distribution
synth_runner log_wage age edu tenure, ///
trunit(1) trperiod(2015) gen_vars
effect_graphs, trlinediff(0)
graph export "figures/fig2e_placebo_gaps.pdf", replace
针对识别的检查（DID 的 PT、弱 IV 的 F、RD 的密度、匹配的共同支撑）也会在第 5 步估计器中自动运行——这里不要重复计算数值，但要生成图表：审稿人会先看图表。
第 4 步 — 诊断性统计检验
更深入的模式请参阅：references/04-statistical-tests.md — Shapiro–Wilk、Jarque–Bera、Breusch–Pagan、White、Cook–Weisberg、Durbin–Watson、Breusch–Godfrey、Wooldridge xtserial、Pesaran CD、xttest3、VIF、estat ovtest（RESET）、增广 Dickey–Fuller、KPSS、Phillips–Perron、Hausman（FE 与 RE）、Sargan–Hansen、estat firststage（弱 IV）。
* Baseline OLS to anchor diagnostics
reg log_wage training age edu tenure
* 4a. Normality of residuals
predict resid, resid
swilk resid // Shapiro–Wilk (N ≤ 5000)
sktest resid // skewness + kurtosis test
* 4b. Heteroskedasticity
estat hettest // Breusch–Pagan / Cook–Weisberg
estat imtest, white // White's general test
* 4c. Autocorrelation (panel)
xtset worker_id year
xtreg log_wage training age edu tenure, fe
xtserial log_wage training age edu tenure // Wooldridge serial-correlation
xttest3 // modified Wald groupwise hetero
* 4d. Cross-sectional dependence (panel)
xtcsd, pesaran abs // Pesaran CD test
* 4e. Multicollinearity
quietly reg log_wage training age edu tenure
estat vif
* 4f. Model specification
estat ovtest // Ramsey RESET
linktest // Stata "linktest"
* 4g. Stationarity (time series)
dfuller log_wage, lags(4) trend // ADF
kpss log_wage, maxlag(4) notrend // KPSS
pperron log_wage, lags(4) // Phillips–Perron
* 4h. Panel unit root (multiple series)
xtunitroot ips log_wage, lags(aic 4) // Im–Pesaran–Shin
xtunitroot llc log_wage, lags(aic 4) // Levin–Lin–Chu
* 4i. Hausman (after running FE and RE)
qui xtreg log_wage training age edu, fe
estimates store fe
qui xtreg log_wage training age edu, re
estimates store re
hausman fe re, sigmamore
决策表：
| 检验 | 原假设 | 若拒绝 |
|---|---|---|
swilk / sktest | 残差服从正态分布 | 大样本：通常忽略；小样本：使用 bootstrap |
estat hettest / imtest, white | 同方差 | 使用 vce(robust) 或 vce(cluster id) |
xtserial / xttest3 | 不存在面板自相关 / 组间异方差 | 按单位聚类 |
xtcsd, pesaran | 不存在截面相关性 | 使用 Driscoll–Kraay 或 xtscc |
estat vif > 10 | — | 删除 / 合并共线回归量 |
estat ovtest | 模型设定正确 | 添加多项式项 / 对数项 |
dfuller 拒绝 + kpss 不拒绝 | 平稳 | 保留水平值 |
dfuller 不拒绝 | 单位根 | 一阶差分或进行协整分析 |
hausman | RE 一致 | 使用 FE |
第 5 步 — 基准实证建模（第 4 节：主要结果）
更深入的模式：references/05-modeling.md — 所有经典估计器及其语法：reg、areg、xtreg、reghdfe、ivreg2 / ivregress / ivreghdfe、logit / probit / ppmlhdfe、csdid / did_imputation / eventstudyinteract / sdid / did_multiplegt_dyn、rdrobust / rddensity / rdmc、synth / synth_runner、psmatch2 / teffects psmatch|ipw|ipwra|aipw、ebalance、heckman、qreg。
这是应用论文中最密集的一节。现代 AER 第 4 节通常包含 2–3 张多回归表和一张系数图：
- 表 2（主要结果）：逐步加入控制变量，4–6 列 — 见下文的 模式 A
- 表 2-bis（设计竞赛）：OLS / IV / DID / 匹配下的同一系数 — 模式 B
- 表 2-ter（多结果）：同一处理变量，多个结果并列呈现 — 模式 C
- 图 3（coefplot）：以图形汇总各规格下的 β̂ 和 95% CI
估计器选择规则（请牢记——选错会悄无声息地得出毫无意义的结果）：
- 无固定效应 / 单个低维固定效应 →
reg y x1 x2, vce(cluster id)- 高维固定效应 →
reghdfe y x1 x2, absorb(fe1 fe2) vce(cluster id)- 双向聚类 →
reghdfe ..., vce(cluster fe1 fe2)- 2SLS / IV →
ivreg2 y x (D = Z), cluster(id) first endog(D)（或在高维固定效应 + IV 时使用ivreghdfe）- DID / 事件研究 →
csdid/eventstudyinteract/did_imputation
根据识别策略选择估计器：
Observational cross-section, selection on obs → reg + controls | teffects psmatch|ipwra
Observational panel, policy shock, parallel trends → csdid / did_imputation / eventstudyinteract / sdid
Exogenous instrument → ivreg2 / ivregress / ivreghdfe
Discontinuity in assignment rule → rdrobust (+ rddensity)
N=1 treated, long panel → synth / synth_runner
Selection on observables + heterogeneity → teffects aipw / ebalance
Binary outcome → logit / probit + margins
Count outcome → poisson / nbreg / ppmlhdfe
规范命令（outreg2 / esttab 的 Stata 对应工具是常用主力——用 eststo 存储 5–6 个规格，再用 esttab 将它们合并为一张表）。关键选项：
keep(...) : list of coefficients to display (e.g. keep(training))
drop(...) : list of coefficients to suppress (controls / intercept)
mtitles("(1)" "(2)" ...) : column labels for the regression table
stats(N r2 r2_a, labels(...)) : footer rows
star(* 0.10 ** 0.05 *** 0.01) : AER stars
addnotes("...") : table footer (cluster level, FE absorbed, sample restrictions)
label booktabs : pretty-print + LaTeX booktabs borders
5.A 模式 A — 逐步加入控制变量（规范的表 2）
各列的 β̂ 保持稳定 ⇒ 不必那么担心观测变量选择驱动了估计结果（Oster 2019 的选择稳定性逻辑；在第 6.j 步中进行量化）。eststo m1...m6 + esttab 是 outreg2 和 R 的 modelsummary 在 Stata 中的对应方案。
eststo clear
eststo m1: qui reg log_wage training, vce(cluster firm_id)
eststo m2: qui reg log_wage training age edu, vce(cluster firm_id)
eststo m3: qui reg log_wage training age edu tenure firm_size, vce(cluster firm_id)
eststo m4: qui reghdfe log_wage training age edu tenure firm_size, absorb(industry year) vce(cluster firm_id)
eststo m5: qui reghdfe log_wage training age edu tenure firm_size, absorb(worker_id year) vce(cluster firm_id)
eststo m6: qui reghdfe log_wage training age edu tenure firm_size, absorb(worker_id year i.industry#i.year) vce(cluster firm_id)
esttab m1 m2 m3 m4 m5 m6 using "tables/table2_main.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) ///
label booktabs ///
mtitles("(1) Baseline" "(2) +Demog" "(3) +Labor-mkt" "(4) Ind×Yr FE" "(5) Worker FE" "(6) Ind×Yr × Worker FE") ///
stats(N r2 r2_a, labels("N" "R²" "Adj. R²")) ///
addnotes("Cluster-robust SE at firm_id in parentheses." ///
"* p<0.10, ** p<0.05, *** p<0.01.")
* Word version: same call with .rtf extension.
AER 惯例：展示所有控制变量——既不传入
keep()，也不传入drop()，这样就能看到每个参数。仅当确实需要只展示核心系数的表格时才使用keep(training)（例如交互形式的异质性分析、IV 第一阶段三列结果）；仅当出于论文排版美观需要隐藏常数项时才使用drop(_cons)。
5.B 模式 B — 设计竞赛（表 2-bis）
展示多个识别策略下同一关注系数的估计值。这是 AER 提升可信度的关键做法：不同设计各自基于不同的识别假设，却得出趋同的证据。
eststo clear
eststo ols: qui reghdfe log_wage training age edu tenure, absorb(industry year) vce(cluster firm_id)
eststo iv: qui ivreg2 log_wage age edu tenure (training = Z1 Z2), cluster(firm_id)
eststo did: qui csdid log_wage age edu tenure, ivar(worker_id) time(year) gvar(first_treat) method(dripw) agg(group)
eststo psm: qui teffects psmatch (log_wage) (training age edu tenure), atet
eststo ebal: qui ebalance training age edu tenure
qui reg log_wage training age edu tenure [pw=_webal], vce(cluster firm_id)
eststo ebal_main
esttab ols iv did psm ebal_main using "tables/table2b_designs.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
keep(training) ///
mtitles("(1) OLS+FE" "(2) 2SLS" "(3) CS-DID" "(4) PSM" "(5) Entropy bal.") ///
stats(N, labels("N")) ///
addnotes("Convergent evidence: same β̂ under five identification strategies.")
5.C 模式 C — 多结果表（相同 X，多个 Y）
一个处理变量，多个结果变量。使用 mtitles，使每列都标注对应的 Y 名称。
eststo clear
foreach y of varlist log_wage weeks_employed left_firm promoted {
eststo `y': qui reghdfe `y' training age edu tenure, ///
absorb(industry year) vce(cluster firm_id)
}
esttab log_wage weeks_employed left_firm promoted using "tables/table2c_multi_outcome.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
keep(training) ///
mtitles("Log wage" "Weeks empl." "Left firm" "Promoted") ///
stats(N r2, labels("N" "R²")) ///
addnotes("Each column is a separate regression on the labelled outcome.")
5.D 模式 D — 堆叠的 Panel A / Panel B 表
相同模型系列，两个时间范围（短期 / 长期）或两个样本（2015 年前 / 2015 年后）。可以通过两次 esttab 调用 + texdoc 拼接后纵向堆叠，也可以使用 esttab ..., refcat() 插入面板标题。
* Panel A — short-run (1 year horizon)
eststo clear
eststo a1: qui reghdfe log_wage_t1 training X, absorb(industry year) vce(cluster firm_id)
eststo a2: qui reghdfe log_wage_t1 training X, absorb(worker_id year) vce(cluster firm_id)
* Panel B — long-run (5 year horizon)
eststo b1: qui reghdfe log_wage_t5 training X, absorb(industry year) vce(cluster firm_id)
eststo b2: qui reghdfe log_wage_t5 training X, absorb(worker_id year) vce(cluster firm_id)
* First panel — write
esttab a1 a2 using "tables/table2d_horizons.tex", replace ///
se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
keep(training) mtitles("(1) Industry FE" "(2) Worker FE") ///
refcat(training "\textbf{Panel A. Short-run (1 year)}", nolabel) ///
stats(N r2, labels("N" "R²"))
* Second panel — append
esttab b1 b2 using "tables/table2d_horizons.tex", append ///
se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
keep(training) mtitles("(1) Industry FE" "(2) Worker FE") ///
refcat(training "\textbf{Panel B. Long-run (5 years)}", nolabel) ///
stats(N r2, labels("N" "R²"))
5.E 模式 E — IV 报告三列（第一阶段 / 简约式 / 2SLS）
教科书式的 AER IV 表会在三列中分别呈现第一阶段、简约式和 2SLS，以便读者验证 Wald 比率 = RF / FS。
eststo clear
eststo fs: qui reghdfe training Z age edu, absorb(industry year) vce(cluster firm_id) // first stage
eststo rf: qui reghdfe log_wage Z age edu, absorb(industry year) vce(cluster firm_id) // reduced form
eststo iv: qui ivreghdfe log_wage age edu (training = Z), absorb(industry year) cluster(firm_id) first
esttab fs rf iv using "tables/table2e_iv_triplet.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
keep(Z training) ///
mtitles("(1) First stage" "(2) Reduced form" "(3) 2SLS") ///
stats(N r2 widstat, labels("N" "R²" "First-stage F (KP)")) ///
addnotes("Wald ratio: $\hat\beta_{2SLS} = \hat\beta_{RF} / \hat\pi_{FS}$.")
**IV 三列结果有意突出核心变量：**仅展示 Z + 内生回归变量，让读者能直观核对 Wald 比率。只有在审稿人要求提供完整系数列表时，才删去
keep(Z training)。
5.F 模式 F — 通过 csdid / synth_runner / teffects 运行因果分析主方法
对于 DID / SCM / 匹配主分析，现代 Stata 估计器会返回自包含的估计结果，并自动提供安慰剂检验 / 前趋势检验 / 重叠诊断。使用估计器自身的报告，然后传入 eststo + esttab。
* CS-DID with pre-trends test
csdid log_wage age edu tenure, ivar(worker_id) time(year) gvar(first_treat) ///
method(dripw) agg(group) saverif(csdid_rif.dta)
estat pretrend // joint zero on pre-period leads
estat simple // ATT(g) summary
estat event, window(-4 4) // dynamic ATT
csdid_estat aggte, type(group) saverif(att_g.dta)
eststo m_csdid
* Synth-runner with 50 placebo treated units → exact p-value on β̂_treated
synth_runner log_wage age edu tenure, trunit(1) trperiod(2015) gen_vars
single_treatment_graphs, do_color(red)
effect_graphs, trlinediff(0)
* Reports exact p-value from the placebo distribution.
* Teffects AIPW with overlap + balance check
teffects aipw (log_wage age edu tenure) (training age edu tenure)
tebalance summarize
eststo m_aipw
5.G 模式 G — 分组 esttab（表 3，参见第 7 步）
每个子组一列。详细代码见 §第 7 步 — 异质性。
5.H 模式 H — 稳健性主表（表 A1，见步骤 6）
将所有稳健性设定与基准设定并列呈现。详细代码见 §步骤 6 — 稳健性主表。
规范估计命令（底层原语）
* 5a. OLS with cluster-robust SE
reg log_wage training age edu tenure, vce(cluster firm_id)
eststo ols_m1
* 5b. Two-way FE — reach for reghdfe (absorb HD FE; clusters any level)
reghdfe log_wage training age edu tenure, ///
absorb(worker_id year) vce(cluster worker_id)
eststo fe_m1
* Multi-way clustering
reghdfe log_wage training, absorb(worker_id year) ///
vce(cluster worker_id firm_id)
* High-dim FE (industry × year, firm × state)
reghdfe log_wage training, absorb(worker_id i.industry#i.year) ///
vce(cluster firm_id)
* 5c. 2×2 DID
reg log_wage i.treated##i.post age edu, vce(cluster worker_id)
* Or with absorbed FE:
reghdfe log_wage c.treated#c.post, absorb(worker_id year) ///
vce(cluster worker_id)
* 5d. Event study (dynamic DID, base period = -1)
* Build relative-time factor var, base at k = -1:
gen rel = rel_time + 10 // shift so -10..10 → 0..20
replace rel = 0 if missing(rel_time) // never-treated → bin 0
* (preferable: use -eventstudyinteract- or -did_imputation- — see reference 5.4)
reghdfe log_wage ib9.rel, absorb(worker_id year) vce(cluster worker_id)
coefplot, keep(*.rel) vertical omitted ///
yline(0) xline(9.5, lpattern(dash)) ///
xtitle("Years relative to treatment") ytitle("Coefficient (ATT)")
graph export "figures/event_study.pdf", replace
* 5e. Staggered DID — modern estimators (see references/05-modeling.md §5.4)
* Callaway–Sant'Anna (2021):
csdid log_wage, ivar(worker_id) time(year) gvar(first_treat) ///
method(dripw) agg(group)
* Sun & Abraham (2021):
eventstudyinteract log_wage rel_dummies_*, ///
cohort(first_treat) control_cohort(never_treated) ///
absorb(worker_id year) vce(cluster worker_id)
* Borusyak–Jaravel–Spiess (2024) imputation:
did_imputation log_wage worker_id year first_treat, ///
allhorizons pretrend(5)
* Synthetic DID (Arkhangelsky et al. 2021):
sdid log_wage worker_id year training, vce(bootstrap) ///
graph g1on g1_opt(ytitle("log wage"))
* 5f. IV / 2SLS — ivreg2 is the de-facto standard (full diagnostics)
ivreg2 log_wage age edu (training = draft_lottery), ///
cluster(firm_id) first endog(training)
* ivreghdfe for HD FE with IV:
ivreghdfe log_wage age (training = draft_lottery), ///
absorb(worker_id year) cluster(worker_id) first
* 5g. Sharp RD
rdrobust outcome running_var, c(0) kernel(triangular) bwselect(mserd)
rdplot outcome running_var, c(0)
* Fuzzy RD
rdrobust outcome running_var, c(0) fuzzy(treatment)
* Density test (McCrary / Cattaneo et al.)
rddensity running_var, c(0)
* 5h. Binary outcome — always follow with margins
logit employed training age edu, vce(cluster firm_id)
margins, dydx(training) atmeans
* 5i. Count outcome with HD FE
ppmlhdfe citations training age, absorb(firm_id year) ///
cluster(firm_id)
* 5j. Quantile regression
qreg log_wage training age edu, quantile(0.5)
sqreg log_wage training age edu, quantile(0.1 0.25 0.5 0.75 0.9) reps(500)
步骤 6 — 稳健性检验组合
深入模式：references/06-robustness.md — 使用 eststo / esttab 逐步构建 M1–M6 设定；使用 bacondecomp 检验 TWFE 偏误；使用 honestdid（Rambachan–Roth）、boottest 野生聚类自助法、ritest 随机化推断、rwolf Romano–Wolf、替代聚类层级、Oster δ*、安慰剂时点、子样本划分，以及通过循环绘制设定曲线。
* 6a. Progressive specifications (M1 → M6)
eststo clear
eststo m1: qui reg log_wage training, vce(cluster firm_id)
eststo m2: qui reg log_wage training age edu, vce(cluster firm_id)
eststo m3: qui reghdfe log_wage training age edu tenure, ///
absorb(worker_id) vce(cluster worker_id)
eststo m4: qui reghdfe log_wage training age edu tenure, ///
absorb(worker_id year) vce(cluster worker_id)
eststo m5: qui reghdfe log_wage training age edu tenure, ///
absorb(worker_id year region) vce(cluster worker_id)
eststo m6: qui reghdfe log_wage training age edu tenure, ///
absorb(worker_id year i.industry#i.year) vce(cluster worker_id)
esttab m1 m2 m3 m4 m5 m6 using "tables/table_main.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) ///
stats(N r2 r2_a, labels("N" "R²" "Adj. R²")) ///
label booktabs
* 6b. Alternative cluster levels
foreach c in worker_id firm_id industry state {
qui reghdfe log_wage training, absorb(worker_id year) vce(cluster `c')
display "cluster=`c' b=" _b[training] " se=" _se[training]
}
* 6c. Wild cluster bootstrap (few clusters)
qui reghdfe log_wage training, absorb(worker_id year) vce(cluster state)
boottest training, cluster(state) reps(9999) seed(42)
* 6d. Subsample splits
foreach mask in "female==0" "female==1" "age<40" "age>=40" {
qui reghdfe log_wage training if `mask', ///
absorb(worker_id year) vce(cluster worker_id)
display "`mask': b=" _b[training] " se=" _se[training] " N=" e(N)
}
* 6e. Placebo timing
gen fake_first = first_treat - 3
gen fake_post = (year >= fake_first)
preserve
keep if year < first_treat // drop real post-period
reghdfe log_wage fake_post, absorb(worker_id year) ///
vce(cluster worker_id)
restore
* 6f. Randomization inference (permutation)
ritest training _b[training], reps(1000) seed(0) ///
strata(worker_id): reghdfe log_wage training, ///
absorb(worker_id year) vce(cluster worker_id)
* 6g. Romano–Wolf multiple testing
rwolf log_wage employed hours_worked, indepvar(training) ///
controls(age edu) reps(500) seed(42) method(reghdfe) ///
fe(worker_id year) cluster(worker_id)
* 6h. TWFE bias diagnosis — Goodman-Bacon decomposition
bacondecomp log_wage training, ddetail
* 6i. Parallel trends sensitivity — HonestDiD
* (after running the event study and saving b/V)
honestdid, pre(1/4) post(5/9) mvec(0(0.1)0.5)
* 6j. Oster (2019) δ*
qui reg log_wage training // short
scalar bs = _b[training]; scalar r2s = e(r2)
qui reghdfe log_wage training age edu tenure, absorb(worker_id year)
scalar bl = _b[training]; scalar r2l = e(r2)
psacalc delta training, mcontrol(age edu tenure) rmax(1.3*r2l)
6.k 模式 H — 稳健性主表（表 A1，每项检验占一列）
标准的 AER 附录表 A1 将所有稳健性设定与基准设定并列呈现，让审稿人一眼就能看出 β̂ 是否稳健。动态构建 eststo 结果列表：
eststo clear
* (1) Baseline
eststo base: qui reghdfe log_wage training age edu tenure, ///
absorb(industry year) vce(cluster firm_id)
* (2) Drop top 1% of wage
qui sum wage, detail
eststo no99: qui reghdfe log_wage training age edu tenure if wage < r(p99), ///
absorb(industry year) vce(cluster firm_id)
* (3) Balanced panel only
preserve
keep if e(sample) // start from baseline sample
bysort worker_id: egen n_obs = count(year)
qui sum n_obs, detail
keep if n_obs == r(max)
eststo balpan: qui reghdfe log_wage training age edu tenure, ///
absorb(industry year) vce(cluster firm_id)
restore
* (4) Drop early adopting cohorts
eststo dropearly: qui reghdfe log_wage training age edu tenure if first_treat > 2008, ///
absorb(industry year) vce(cluster firm_id)
* (5) Worker FE absorbed (selection on FE-removable unobservables)
eststo wfe: qui reghdfe log_wage training age edu tenure, ///
absorb(worker_id year) vce(cluster firm_id)
* (6) Two-way cluster (firm × year)
eststo cl2way: qui reghdfe log_wage training age edu tenure, ///
absorb(industry year) vce(cluster firm_id year)
* (7) Wild cluster bootstrap (after running the baseline)
qui reghdfe log_wage training age edu tenure, absorb(industry year) vce(cluster firm_id)
boottest training, cluster(firm_id) reps(9999) seed(42) nograph // returns CI
* Note bootstrapped p-value in the row footer of the master table.
* (8) Log outcome
eststo logy: qui reghdfe log_log_wage training age edu tenure, ///
absorb(industry year) vce(cluster firm_id)
* (9) IHS outcome (handles zeros)
gen ihs_wage = asinh(wage)
eststo ihsy: qui reghdfe ihs_wage training age edu tenure, ///
absorb(industry year) vce(cluster firm_id)
* (10) PSM-weighted outcome
qui psmatch2 training age edu tenure, out(log_wage) n(1) common
eststo psm: qui reg log_wage training age edu tenure if _support == 1, vce(cluster firm_id)
* (11) Entropy balance
ebalance training age edu tenure
eststo ebal: qui reg log_wage training age edu tenure [pw=_webal], vce(cluster firm_id)
esttab base no99 balpan dropearly wfe cl2way logy ihsy psm ebal ///
using "tables/tableA1_robustness.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
keep(training) ///
mtitles("(1) Baseline" "(2) Drop top 1%" "(3) Balanced" "(4) Drop early" ///
"(5) Worker FE" "(6) 2-way cluster" "(7) log Y" "(8) IHS Y" ///
"(9) PSM" "(10) Entropy bal.") ///
stats(N r2, labels("N" "R²")) ///
addnotes("Each column is one robustness check. β̂ on training is the focal coefficient.")
6.l Stata 中的设定曲线（Simonsohn–Simmons–Nelson 2020）
esttab 没有内置的设定曲线功能；循环遍历 {控制变量 × 样本 × 结果变量变换 × SE 类型}，将每个 β̂ + 95% CI 存入一个 frame，然后绘图。
capture frame drop spec
frame create spec spec_id b se lo hi label
local id 0
foreach controls in "age" "age edu" "age edu tenure" "age edu tenure firm_size" {
foreach trans in "log_wage" "ihs_wage" { // outcome transforms
foreach mask in "1" "industry==\"manuf\"" "wage < r(p99)" { // subsamples
foreach cl in "firm_id" "firm_id year" { // SE types
qui reghdfe `trans' training `controls' if `mask', ///
absorb(industry year) vce(cluster `cl')
local ++id
local b = _b[training]
local se = _se[training]
local lo = `b' - 1.96*`se'
local hi = `b' + 1.96*`se'
frame post spec (`id') (`b') (`se') (`lo') (`hi') ("c=`controls'/y=`trans'/s=`mask'/cl=`cl'")
}
}
}
}
frame change spec
sort b
gen rank = _n
twoway (rcap lo hi rank) (scatter b rank), ///
yline(0, lpattern(dash) lcolor(gs10)) ///
xtitle("Specification (sorted by β̂)") ytitle("Coefficient on training") ///
title("Figure 5. Specification curve") ///
legend(off) scheme(s2color)
graph export "figures/fig5_spec_curve.pdf", replace
graph export "figures/fig5_spec_curve.png", replace width(2400)
frame change default
现代 AER 审稿意见经常会问“设定 X 怎么样？”——设定曲线可以一次性回答所有这类问题。
6.m 敏感性分析仪表板（HonestDiD + Oster + E-value，单个代码块）
对于 DID 主要结果，生成一张图（HonestDiD 部分识别界限）+ 一行敏感性统计量（Oster δ + E-value）。
* (a) HonestDiD — Rambachan-Roth (2023) bound on β̂ under bounded PT violation
* Re-run the event study and store b/V:
qui eventstudyinteract log_wage rel_d1-rel_d9, ///
cohort(first_treat) control_cohort(never_treated) ///
absorb(worker_id year) vce(cluster worker_id)
honestdid, pre(rel_d1 rel_d2 rel_d3) post(rel_d5 rel_d6 rel_d7 rel_d8 rel_d9) ///
mvec(0(0.1)0.5) coefplot
graph export "figures/fig6_honestdid.pdf", replace
* (b) Oster δ (Pattern: how big would selection on unobservables have to be?)
qui psacalc delta training, mcontrol(age edu tenure firm_size) rmax(1.3)
display "Oster δ for β=0: " r(delta)
* (c) E-value (Linden-Mathur 2020) — for binary / risk-ratio outcomes
* Convert OLS coefficient to risk-ratio scale first if applicable.
evalue rr, est(1.45) lcl(1.10) ucl(1.91)
* → reports the minimum strength of unmeasured confounding to nullify the result.
步骤 7 — 进一步分析（机制 / 异质性 / 中介 / 调节）
深入模式：references/07-further-analysis.md — 因子变量交互、margins / marginsplot、三重差分（DDD）、结果指标阶梯、medsem / Baron–Kenny、用于非线性中介分析的 khb（Karlson–Holm–Breen）、通过 xtile 或样条进行剂量–反应分析、亚组事件研究。
* 7a. Heterogeneity via interaction — coefficient IS the heterogeneity test
reghdfe log_wage c.training##i.female age edu, ///
absorb(worker_id year) vce(cluster worker_id)
margins, dydx(training) at(female=(0 1))
marginsplot, title("Marginal effect of training by gender") ///
ytitle("dY/d(training)")
graph export "figures/het_female.pdf", replace
* Continuous moderator
reghdfe log_wage c.training##c.tenure age edu, ///
absorb(worker_id year) vce(cluster worker_id)
margins, dydx(training) at(tenure=(0(2)20))
marginsplot
* 7b. Subgroup estimation + Wald test of equality
eststo clear
eststo m_male: qui reghdfe log_wage training if female==0, ///
absorb(worker_id year) vce(cluster worker_id)
eststo m_female: qui reghdfe log_wage training if female==1, ///
absorb(worker_id year) vce(cluster worker_id)
suest m_male m_female // joint covariance
test [m_male_mean]training = [m_female_mean]training // Wald test
* 7c. Triple difference (DDD)
reghdfe log_wage c.treated##c.post##c.high_exposure, ///
absorb(worker_id year) vce(cluster firm_id)
* coefficient on treated#post#high_exposure IS the differential DID
* 7d. Outcome ladder
eststo clear
foreach y of varlist hours_worked productivity log_wage {
eststo: qui reghdfe `y' training, ///
absorb(worker_id year) vce(cluster worker_id)
}
esttab using "tables/mechanism_ladder.tex", replace ///
se label booktabs
* 7e. Mediation — Baron–Kenny (manual)
qui reghdfe log_wage training age edu, absorb(worker_id year)
scalar c = _b[training]
qui reghdfe hours_worked training age edu, absorb(worker_id year)
scalar a = _b[training]
qui reghdfe log_wage training hours_worked age edu, absorb(worker_id year)
scalar b_m = _b[hours_worked]
scalar cprim = _b[training]
display "Total = " c " Direct = " cprim " Indirect a*b = " a*b_m
* Or via -medsem- (for SEM-based mediation with bootstrap CI)
medsem, indep(training) med(hours_worked) dep(log_wage) ///
mcreps(1000)
* 7f. CATE via Nonparametric (grf) — typically fall back to Python econml
* See references/07-further-analysis.md §7.7 for the rpy2-style bridge.
* 7g. Dose-response (continuous treatment, decile bins)
xtile dose10 = training_hours, nq(10)
reghdfe log_wage i.dose10 age edu, absorb(worker_id year) vce(cluster worker_id)
margins i.dose10
marginsplot, recast(connected) ytitle("Predicted log wage") ///
xtitle("Training-hours decile")
graph export "figures/dose_response.pdf", replace
步骤 8 — 发表用表格与图表
此步骤为必做项——每次分析运行都会生成本技能顶部 默认输出规范 中定义的全部 5 张必需表格（T1–T5）和全部 4 张必需图表（F1–F4）。不要因为“回归已经跑完了”而跳过步骤 8。没有表格和图表的系数，并非应用经济学交流结果的方式。
深入模式：references/08-tables-plots.md — 使用 esttab 和 outreg2 生成回归表；使用 coefplot 绘制系数图、事件研究图和森林图；使用 marginsplot 展示交互效应；使用 binscatter 绘制残差化散点图；使用 rdplot 绘制 RD 图；使用 twoway + 合并图表；Stata 方案 / 图形选项；导出为 LaTeX / Word / Excel。
* ============================================================
* 8a. ★ TABLE 2 — Main results, multi-column regression M1→M6
* (the centerpiece of every economics paper)
* ============================================================
esttab m1 m2 m3 m4 m5 m6 using "tables/table2_main.tex", ///
replace ///
se star(* 0.10 ** 0.05 *** 0.01) ///
stats(N r2 r2_a fe_worker fe_year fe_indyr, ///
labels("N" "R²" "Adj. R²" "Worker FE" "Year FE" "Industry×Year FE")) ///
keep(training age edu tenure) ///
order(training age edu tenure) ///
mtitles("(1) Raw" "(2) +Demog" "(3) +Tenure" "(4) +Unit FE" "(5) +2-way FE" "(6) +Ind×Yr FE") ///
label booktabs nonumbers noobs ///
addnotes("Cluster-robust standard errors at worker_id level in parentheses." ///
"* p<0.10, ** p<0.05, *** p<0.01.")
* Word version
esttab m1 m2 m3 m4 m5 m6 using "tables/table2_main.rtf", replace ///
se star(* 0.10 ** 0.05 *** 0.01) label ///
mtitles("(1)" "(2)" "(3)" "(4)" "(5)" "(6)")
* outreg2 — alternative (one-line, broad compatibility)
* reghdfe log_wage training age edu tenure, absorb(worker_id year)
* outreg2 using "tables/table2_main_outreg2.doc", replace label dec(3) ///
* addtext(Worker FE, YES, Year FE, YES)
* ============================================================
* 8b. TABLE 1 — Summary statistics & balance
* ============================================================
balancetable training age edu tenure female ///
using "tables/table1_balance.tex", ///
vce(cluster firm_id) replace ///
varlabels pval booktabs
asdoc sum log_wage training age edu tenure female, ///
by(training) stat(N mean sd min median max) ///
save(tables/table1_balance.rtf) replace
* ============================================================
* 8c. TABLE 3 — Mechanism / outcome ladder (3+ outcomes)
* ============================================================
eststo clear
foreach y of varlist hours_worked productivity log_wage {
eststo: qui reghdfe `y' training age edu tenure, ///
absorb(worker_id year) vce(cluster worker_id)
}
esttab using "tables/table3_mechanism.tex", replace ///
se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
mtitles("Hours worked" "Productivity" "Log wage") ///
keep(training) ///
addnotes("Each column is a separate regression of the labelled outcome on training," ///
"with worker and year FE and cluster-robust SE at worker_id.")
* ============================================================
* 8d. TABLE 4 — Heterogeneity (subgroup × main coef + Wald)
* ============================================================
eststo clear
eststo all: qui reghdfe log_wage training, ///
absorb(worker_id year) vce(cluster worker_id)
eststo male: qui reghdfe log_wage training if female==0, ///
absorb(worker_id year) vce(cluster worker_id)
eststo female: qui reghdfe log_wage training if female==1, ///
absorb(worker_id year) vce(cluster worker_id)
eststo young: qui reghdfe log_wage training if age<40, ///
absorb(worker_id year) vce(cluster worker_id)
eststo old: qui reghdfe log_wage training if age>=40, ///
absorb(worker_id year) vce(cluster worker_id)
eststo manuf: qui reghdfe log_wage training if industry=="manufacturing", ///
absorb(worker_id year) vce(cluster worker_id)
esttab all male female young old manuf using "tables/table4_heterogeneity.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
mtitles("All" "Female=0" "Female=1" "Age<40" "Age≥40" "Manuf.") ///
keep(training) ///
addnotes("Cluster-robust SE at worker_id. Wald p-values for cross-subgroup equality" ///
"(via -suest-) should accompany this table — see references/07.")
* ============================================================
* 8e. TABLE 5 — Robustness battery (alt SE / cluster / sample / placebo)
* ============================================================
eststo clear
eststo base: qui reghdfe log_wage training, absorb(worker_id year) vce(cluster worker_id)
eststo cl_firm: qui reghdfe log_wage training, absorb(worker_id year) vce(cluster firm_id)
eststo cl_2way: qui reghdfe log_wage training, absorb(worker_id year) vce(cluster worker_id firm_id)
eststo wins: qui reghdfe log_wage_w1 training, absorb(worker_id year) vce(cluster worker_id)
eststo dropmf: qui reghdfe log_wage training if industry!="manufacturing", ///
absorb(worker_id year) vce(cluster worker_id)
eststo placebo: qui reghdfe log_wage fake_post if year < first_treat, ///
absorb(worker_id year) vce(cluster worker_id)
esttab base cl_firm cl_2way wins dropmf placebo using "tables/table5_robustness.tex", ///
replace se star(* 0.10 ** 0.05 *** 0.01) label booktabs ///
mtitles("Baseline" "Cluster=Firm" "2-way Cluster" "Winsor 1/99" "Drop Manuf." "Placebo (-3 yr)") ///
keep(training fake_post)
* ============================================================
* 8f. ★ FIGURE 3 — Coefficient plot across M1→M6
* ============================================================
coefplot m1 m2 m3 m4 m5 m6, keep(training) ///
vertical yline(0, lpattern(dash)) ///
ciopts(recast(rcap)) ///
levels(95) ///
ylabel(, angle(0)) ///
xtitle("Specification") ytitle("Coefficient on training (95% CI)") ///
title("Effect of training across specifications") ///
scheme(s2color)
graph export "figures/fig3_coefplot.pdf", replace
graph export "figures/fig3_coefplot.png", replace width(2400)
* ============================================================
* 8g. FIGURE 2 — Event-study plot (dynamic DID, base period = -1)
* ============================================================
* (after running -eventstudyinteract- or -csdid- or reghdfe with relative-time factor)
coefplot, keep(*.rel) vertical omitted ///
xline(9.5, lpattern(dash)) yline(0, lpattern(dash)) ///
xtitle("Years relative to treatment") ///
ytitle("Coefficient (ATT, 95% CI)") ///
title("Event study: dynamic effect of training") ///
ciopts(recast(rcap)) ///
scheme(s2color)
graph export "figures/fig2_event_study.pdf", replace
graph export "figures/fig2_event_study.png", replace width(2400)
* ============================================================
* 8h. FIGURE 4 — Sensitivity / robustness curve
* (HonestDiD post-DID; or forest plot of robustness battery)
* ============================================================
* HonestDiD sensitivity (after the event study with stored b/V):
honestdid, pre(1/4) post(5/9) mvec(0(0.1)0.5) coefplot
graph export "figures/fig4_sensitivity.pdf", replace
graph export "figures/fig4_sensitivity.png", replace width(2400)
* Alternative: forest plot of subgroup / robustness coefficients
* coefplot (all, label(All)) (male, label(Male)) (female, label(Female)) ///
* (young, label(Age<40)) (old, label(Age≥40)) (manuf, label(Manuf.)), ///
* keep(training) xline(0, lpattern(dash)) ciopts(recast(rcap)) ///
* title("Heterogeneity forest plot")
* graph export "figures/fig4_forest.pdf", replace
* ============================================================
* 8i. FIGURE 1 — Trend / motivation (treated vs control over time)
* ============================================================
preserve
collapse (mean) log_wage, by(year training)
twoway (line log_wage year if training==1, lcolor(navy) lwidth(medthick)) ///
(line log_wage year if training==0, lcolor(cranberry) lwidth(medthick)), ///
xline(`policy_year', lpattern(dash) lcolor(gs8)) ///
legend(order(1 "Treated" 2 "Control") rows(1) position(6)) ///
xtitle("Year") ytitle("Mean log wage") ///
title("Treated vs control trend") ///
scheme(s2color)
graph export "figures/fig1_trend.pdf", replace
graph export "figures/fig1_trend.png", replace width(2400)
restore
* ============================================================
* 8j. Auxiliary plots (optional — produce when relevant to design)
* ============================================================
* Binscatter — residualized scatter
binscatter log_wage tenure, controls(age edu female) nquantiles(20) ///
savegraph("figures/figA_binscatter.pdf") replace
* RD plot (only when running_var exists)
* rdplot outcome running_var, c(0) ///
* graph_options(title("Effect of eligibility") ///
* ytitle("Outcome") xtitle("Running variable"))
* graph export "figures/figA_rdplot.pdf", replace
* Bacon decomposition plot — TWFE diagnostic
* bacondecomp log_wage training, ddetail
* graph export "figures/figA_bacon.pdf", replace
* ============================================================
* 8k. Graph theme (consistent styling across all figures)
* ============================================================
set scheme s2color
* or the modern: set scheme white_tableau (from -schemepack-)
交付物检查清单（在宣布运行完成前核对）：
[ ] tables/table1_balance.tex [ ] figures/fig1_trend.pdf
[ ] tables/table2_main.tex ★ [ ] figures/fig2_event_study.pdf
[ ] tables/table3_mechanism.tex [ ] figures/fig3_coefplot.pdf
[ ] tables/table4_heterogeneity.tex
[ ] tables/table5_robustness.tex [ ] figures/fig4_sensitivity.pdf
[ ] tables/tableA1_robustness.tex [ ] figures/fig5_spec_curve.pdf
[ ] artifacts/sample_construction.json (footnote 4)
[ ] artifacts/data_contract.json
[ ] artifacts/result.json (reproducibility stamp — see 8l)
8l. 可复现性标记
期刊复现办公室（或未来的合著者）复现核心结果所需的单一产物。保存 Stata 版本、随机种子、数据集哈希、基准系数 + CI，以及协议 / 合同的指针：
* After the headline regression has run as `eststo base`:
qui estimates restore base
matrix b = e(b); matrix V = e(V)
local b_hat = b[1, "training"]
local se = sqrt(V["training","training"])
local lo = `b_hat' - 1.96 * `se'
local hi = `b_hat' + 1.96 * `se'
file open f using "artifacts/result.json", write replace
file write f "{" _n
file write f `" "stata_version": "`c(stata_version)'","' _n
file write f `" "current_date": "`c(current_date)' `c(current_time)'","' _n
file write f `" "seed": 42,"' _n
file write f `" "n_obs": `=e(N)',"' _n
file write f `" "estimand": "ATT","' _n
file write f `" "estimator": "`e(cmd)'","' _n
file write f `" "estimate": `b_hat',"' _n
file write f `" "se_cluster": `se',"' _n
file write f `" "ci95_lo": `lo',"' _n
file write f `" "ci95_hi": `hi',"' _n
file write f `" "pre_registration": "strategy.do","' _n
file write f `" "data_contract": "artifacts/data_contract.json","' _n
file write f `" "sample_log": "artifacts/sample_construction.json","' _n
file write f `" "paper_bundle": "tables/table2_main.tex"' _n
file write f "}" _n
file close f
将 artifacts/result.json 与论文 PDF 一并提交。审稿人应能运行 make clean && do main.do，并逐位复现此 JSON。
§A — 流行病学 / 公共卫生模式
当用户措辞表明采用模式 A（目标试验模拟 / IPTW / TMLE / MR / STROBE / 流行病学 / 公共卫生 / RWE / 队列）时，仍适用这 8 个步骤，但步骤 5 会以双重稳健 + 生存分析 + MR 三件套取代 OLS 和固定效应组合，交付物则遵循 STROBE / TRIPOD-AI 惯例。步骤 1–4（清理、构造、表 1、诊断）和步骤 8（表格 / 图表导出）与默认模式相同。
命令配置（在默认 Stata 组合基础上安装）：
ssc install eltmle // TMLE for binary outcome
ssc install mrrobust // Mendelian randomization (IVW / Egger / weighted median)
ssc install mregger // alternative MR-Egger implementation
ssc install mrpresso // outlier-robust MR
ssc install evalue // E-value sensitivity (Linden-Mathur)
ssc install strmst2 // RMST contrast for survival
ssc install gformula // parametric g-formula (time-varying)
* `teffects ipw / ipwra / aipw` are native to Stata 13+; no install needed.
A.0 队列构建 + 目标试验方案
在接触数据之前写好方案。将其保存为 protocol.do，并在论文中引用。
*--- protocol.do — target-trial emulation skeleton ---
* eligibility: age 40-75, no prior MI, ascertained at t0
* treatment: A=1 statin initiation; A=0 no initiation
* assignment: emulated random at t0 via IPTW on baseline covariates
* outcome: incident MI within 5 years
* estimand: ITT ATE on risk difference + hazard ratio
use raw/cohort.dta, clear
keep if inrange(age, 40, 75) & prior_MI == 0
gen t0 = cond(missing(statin_initiation_date), enrollment_date, statin_initiation_date)
gen event_5y = (MI_date - t0 <= 365 * 5) & !missing(MI_date)
gen time_at_risk = min(censor_date - t0, 365 * 5)
stset time_at_risk, failure(event_5y)
A.1 按暴露分组的表 1（与默认模式第 3 步相同）
使用第 3 步中相同的 balancetable / asdoc tabstat，只需添加 by(A)。未测量混杂的 E 值放在页脚。
balancetable A age edu smoke bmi ldl sbp using tables/tableA1_strobe.tex, ///
replace ctitles("Untreated" "Treated" "Diff (SE)") pvalues
A.2 DAG + 倾向评分重叠（正值性检查）
* Estimate PS via logit
logit A age edu smoke bmi ldl sbp
predict ps, pr
* Overlap density (positivity)
twoway (kdensity ps if A == 0) (kdensity ps if A == 1), ///
legend(order(1 "A=0" 2 "A=1")) xtitle("Propensity score")
graph export figures/figA2_ps_overlap.pdf, replace
* Love plot — pre vs post-IPTW SMDs
teffects ipw (event_5y) (A age edu smoke bmi ldl sbp)
tebalance summarize
tebalance density ps
graph export figures/figA2_love.pdf, replace
A.3 IPTW + g 公式 + TMLE 双重稳健三件套（替换第 5 步）
流行病学中的“AER 表 2”：一个 3 行表格，每行对应以下方法之一：{通过 teffects ipw 实现的 IPTW-MSM、通过 teffects aipw 实现的 IPWRA / AIPW、通过 eltmle 实现的 TMLE}，便于读者确认双重稳健方法得出的结果是否一致。
* IPTW-MSM (marginal effect on risk difference)
eststo m_iptw: teffects ipw (event_5y) (A age edu smoke bmi ldl sbp), atet
* AIPW (doubly robust ATE)
eststo m_aipw: teffects aipw (event_5y age edu smoke bmi ldl sbp) ///
(A age edu smoke bmi ldl sbp)
* IPWRA — alternative DR estimator
eststo m_ipwra: teffects ipwra (event_5y age edu smoke bmi ldl sbp) ///
(A age edu smoke bmi ldl sbp)
* TMLE via eltmle (uses SuperLearner under the hood)
preserve
eltmle event_5y A age edu smoke bmi ldl sbp, tmle
restore
* Stack the triplet
esttab m_iptw m_aipw m_ipwra using tables/tableA3_dr_triplet.tex, ///
replace b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("IPTW" "AIPW" "IPWRA") ///
title("Doubly-robust risk-difference triplet")
A.4 生存结局——KM / Cox / AFT / RMST
* KM by treatment (already stset above)
sts graph, by(A) ci risktable
graph export figures/figA4_km.pdf, replace
* Cox HR (covariate-adjusted)
stcox A age edu smoke bmi ldl sbp
estimates store m_cox
* HR + 95% CI for A:
matrix b = e(b); matrix V = e(V)
local HR = exp(b[1, "A"])
local lo = exp(b[1, "A"] - 1.96*sqrt(V["A","A"]))
local hi = exp(b[1, "A"] + 1.96*sqrt(V["A","A"]))
display "HR = " %5.3f `HR' " 95% CI [" %5.3f `lo' ", " %5.3f `hi' "]"
* AFT (Weibull) for time-ratio interpretation
streg A age edu smoke bmi ldl sbp, distribution(weibull) time
* RMST contrast at t = 5 years
strmst2 A, tau(1825) covariates(age edu smoke bmi ldl sbp)
A.5 孟德尔随机化（IVW / Egger / 加权中位数三件套）
* Two-sample MR with pre-extracted instrument betas/SEs
* (BX, BXSE, BY, BYSE in long format, one row per SNP)
mrrobust, beta_outcome(by) se_outcome(byse) beta_exposure(bx) se_exposure(bxse)
mregger by bx [aw=1/byse^2], gxse(bxse)
mrmedian by bx, weighted gxse(bxse) gyse(byse)
* Outlier-robust + heterogeneity sensitivity
mrpresso by bx, gxse(bxse) gyse(byse) outlier distortion seed(1) niter(1000)
* Stack the triplet
esttab using tables/tableA5_mr_triplet.tex, ///
replace mtitles("IVW" "MR-Egger" "Weighted median") ///
title("Mendelian randomization triplet")
A.6 稳健性——E 值 / 界限 / 主分层
* E-value (Linden-Mathur Stata implementation) — required strength of unmeasured confounding
evalue rr, est(1.45) lcl(1.10) ucl(1.91)
A.7 STROBE / TRIPOD-AI 报告清单
保存为 replication/strobe_checklist.do（或 .md），并在提交前逐项勾选：
[ ] Eligibility criteria + dates (target-trial protocol)
[ ] Adjustment set with DAG justification (A.2)
[ ] Positivity / overlap diagnostic (A.2)
[ ] Doubly-robust triplet (IPTW + AIPW + TMLE) (A.3)
[ ] Risk difference + hazard ratio + RMST (A.3, A.4)
[ ] E-value for unmeasured confounding (A.6)
[ ] Loss-to-follow-up rate + censoring assumption (A.0)
[ ] Pre-registered protocol or analysis plan (A.0)
§B — ML 因果推断模式
当用户的措辞表明要使用模式 B（DML / 元学习器 / 因果森林 / CATE / 策略学习 / 因果机器学习）时，流程保留默认模式的步骤 1–4 和步骤 8，将步骤 5 替换为 ML 估计器组合，并在步骤 7 和步骤 8 之间加入 CATE 分布 + 策略价值层。
原生功能与 Python 调用说明：Stata 17+ 附带用于 DML 和因果森林的一流 ML 因果命令。对于神经因果方法（Dragonnet / TARNet / CEVAE）、保形因果方法和公平性审计，该技能会通过 Stata 18 的
python:代码块调用 Python——每一步都会明确标注。
命令组成（在默认 Stata 堆栈之上安装）：
ssc install ddml // double machine learning (Ahrens-Hansen-Schaffer-Wiemann)
ssc install pdslasso // post-double-selection lasso for high-dim controls
ssc install ivlasso // IV variant of pdslasso
ssc install crforest // causal forest (Athey-Tibshirani-Wager) for Stata
ssc install lassopack // covariate selection helpers
* For BART / BCF / Dragonnet / conformal: shell out to Python via `python:` block
* (Stata 18 ships first-class Python integration)
B.0 训练集/留出集划分 + 干扰参数超级学习器
如果干扰参数学习器（结局回归 Q(X,A) = E[Y|X,A] 和倾向评分 g(A|X) = Pr(A=1|X)）收敛速度快于 n^(-1/4)，那么 DML / TMLE / 因果森林就是双重稳健的。使用 pystacked（Stata 18 + Python 调用）在交叉验证下组合 GBM / RF / Lasso / NN——这相当于 StatsPAI 的 SuperLearner。
set seed 42
gen u = runiform()
gen byte holdout = u > 0.7
* pystacked — installed via `python: import pystacked` (requires Stata 18 + Python ≥ 3.8)
ssc install pystacked, replace
* The standard nuisance pair (will be reused inside ddml below):
* ml_g = outcome model pystacked Y X*, type(reg) methods(ols rf gradboost lassocv)
* ml_m = propensity model pystacked D X*, type(class) methods(logit rf gradboost)
* For TMLE-style binary-outcome workflows, use -eltmle- which already wraps SL.
* Standalone holdout is needed for the OFF-policy evaluation in B.4 — DML / GRF
* handle their own cross-fitting internally and don't need a manual split.
B.1 DAG / 估计目标声明（可选用 LLM 辅助）
* Causal discovery is thin in Stata — call out to Python's causal-learn:
python:
import pandas as pd
from causallearn.search.ConstraintBased.PC import pc
from sfi import Data
df = pd.DataFrame({v: Data.get(v) for v in ["A","Y","X1","X2","X3","X4"]})
cg = pc(df.to_numpy())
cg.draw_pydot_graph(labels=list(df.columns)).write_pdf("figures/figB1_dag.pdf")
end
B.2 估计器组合——DML · 元学习器 · 因果森林（替换第 5 步）
ML 因果推断中的“AER 表 2”：一个竞赛表，其中每列代表使用相同 (Y, A, X) 数据的一类估计器——读者希望看到 DML（PLR）、DML（交互式）和因果森林对 ATE 得出一致（或不一致）的结果。
*--- 1. DML — partially linear regression ---
ddml init partial, kfolds(5) reps(5)
ddml E[Y|X]: reg Y X1 X2 X3 X4
ddml E[Y|X]: pystacked Y X1 X2 X3 X4, type(reg) methods(rf gradboost lassocv)
ddml E[D|X]: reg A X1 X2 X3 X4
ddml E[D|X]: pystacked A X1 X2 X3 X4, type(reg) methods(rf gradboost lassocv)
ddml crossfit
ddml estimate, robust
eststo m_dml_plr
*--- 2. DML — interactive (heterogeneous treatment effects) ---
ddml init interactive, kfolds(5) reps(5)
ddml E[Y|X,D]: pystacked Y X1 X2 X3 X4, type(reg) methods(rf gradboost)
ddml E[D|X]: pystacked A X1 X2 X3 X4, type(class) methods(rf gradboost logit)
ddml crossfit
ddml estimate, atet
eststo m_dml_atet
*--- 3. Causal forest — non-parametric CATE ---
crforest Y A X1 X2 X3 X4, ntrees(2000) honesty
predict cate, te
eststo m_cf
*--- 4. (Optional) Neural causal (Dragonnet / TARNet) — Python callout ---
* Stata has no native neural-causal estimator; shell out to econml / causalml:
python:
import numpy as np, pandas as pd
from sfi import Data
from econml.dr import DRLearner
from econml.metalearners import XLearner
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
df = pd.DataFrame({v: Data.get(v) for v in ["Y","A","X1","X2","X3","X4"]})
X = df[["X1","X2","X3","X4"]].values; T = df["A"].values; Y = df["Y"].values
dr = DRLearner(model_propensity=GradientBoostingClassifier(),
model_regression=GradientBoostingRegressor(),
model_final=GradientBoostingRegressor()).fit(Y, T, X=X)
ate, lo, hi = dr.ate(X), *dr.ate_interval(X)
Data.addVarFloat("cate_dr"); Data.store("cate_dr", range(len(df)), dr.effect(X))
end
*--- 5. (Optional) Bayesian Causal Forest (BCF) — R callout ---
* Use the bcf R package via -shell Rscript- to fit the full posterior over CATE.
* The R script writes posterior_summary.csv; we import it back as a frame:
* shell Rscript scripts/bcf.R
* import delimited posterior_summary.csv, clear case(preserve)
*--- 6. Stack the horse-race ---
esttab m_dml_plr m_dml_atet m_cf using tables/tableB2_ml_horserace.tex, ///
replace b(4) se(4) star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("DML (PLR)" "DML (ATET)" "Causal forest") ///
stats(N, labels("N")) ///
title("ML causal horse-race on \$(Y, A, X)\$") ///
addnotes("DR-Learner / Dragonnet / BCF columns added via Python/R callouts when relevant.")
B.3 CATE 分布 + 亚组 CATE 图（扩展第 7 步）
* CATE histogram from causal forest
histogram cate, frequency normal xtitle("CATE") xline(0)
graph export figures/figB3_cate_hist.pdf, replace
* CATE by quartile of a covariate
xtile age_q = X1, n(4)
collapse (mean) cate, by(age_q)
graph bar (asis) cate, over(age_q) ytitle("Mean CATE")
graph export figures/figB3_cate_by_age_q.pdf, replace
B.4 策略学习 + 离线策略评估
* Stata's native policy-tree support is via crforest; for honest depth-3 trees,
* call out to R's policytree from within Stata:
shell Rscript scripts/policy_tree.R // saves policy_pred.csv
import delimited using policy_pred.csv, clear case(preserve)
* DR off-policy evaluation on holdout
keep if holdout
gen DR_match = (policy_pred == A) * Y - (policy_pred != A) * Y
summarize DR_match
display "DR policy value (holdout): " r(mean)
B.5 不确定性（保形因果方法）+ 公平性 + 敏感性分析
python:
# Conformal prediction interval around CATE — distribution-free coverage guarantee
import numpy as np
from sfi import Data
from sklearn.ensemble import GradientBoostingRegressor
from mapie.regression import MapieRegressor
df = ... # same DataFrame as B.1
mapie = MapieRegressor(estimator=GradientBoostingRegressor(), method="plus", cv=10)
# fit on training CATE, predict 90% PI on holdout
end
* Fairness audit — disparate impact (compute via teffects diff across sensitive groups)
foreach g in 0 1 {
quietly mean cate if sensitive_attr == `g'
matrix b = e(b)
display "Mean CATE in group `g': " b[1, 1]
}
B.6 ML 因果推断专用报告清单
保存为 replication/ml_causal_checklist.md：
[ ] Nuisance learners listed (Q model, g model, hyperparameters, CV folds)
[ ] Cross-fitting / sample-splitting documented (ddml kfolds(K) reps(R))
[ ] Overlap / propensity diagnostics (B.0 + A.2-style overlap plot)
[ ] CATE summary (mean, SD, quartiles) + heterogeneity p-value (crforest test)
[ ] Policy value with confidence interval (B.4)
[ ] Conformal coverage rate on holdout (B.5)
[ ] Fairness gaps across sensitive attributes (B.5)
[ ] DAG / adjustment set + sensitivity to unmeasured confounding (E-value or Manski bounds)
库 / 命令速查表
| 步骤 | 任务 | 首选命令 | 备用方案 |
|---|---|---|---|
| 1 | 导入 | use / import excel / import delimited / import sas / import spss | usespss / infile |
| 1 | 缺失值 | misstable summarize / mdesc | missings report |
| 1 | 合并 | merge m:1 ... assert(match using) | — |
| 1 | 面板 | xtset / xtdescribe | 纯时间序列使用 tsset |
| 2 | 缩尾 | winsor2 | 使用 summ, detail + replace 手动处理 |
| 2 | 滞后 / 超前 / 差分 | L. / F. / D.（需要 xtset） | bys id (t): gen lx = x[_n-1] |
| 3 | 表 1 | tabstat + balancetable + asdoc sum | 手动循环 |
| 3 | 相关性 | pwcorr, sig star(.05) | corr |
| 4 | 异方差 | estat hettest / estat imtest, white | — |
| 4 | 面板序列相关 | xtserial / xttest3 | — |
| 4 | 平稳性 | dfuller / kpss / pperron / xtunitroot | — |
| 4 | Hausman 检验 | hausman fe re, sigmamore | — |
| 5 | OLS / 面板固定效应 | reghdfe | areg / xtreg, fe |
| 5 | IV | ivreg2 / ivreghdfe | ivregress 2sls |
| 5 | DID — 2×2 | 使用 c.treated#c.post 的 reg / reghdfe | — |
| 5 | DID — CS | csdid | — |
| 5 | DID — SA | eventstudyinteract | — |
| 5 | DID — 插补法 | did_imputation | did_multiplegt_dyn |
| 5 | DID — SDID | sdid | — |
| 5 | RD | rdrobust / rddensity / rdplot | 多个断点时使用 rdmc |
| 5 | 合成控制 | synth / synth_runner | — |
| 5 | PSM | psmatch2 / teffects psmatch | — |
| 5 | IPW / AIPW | teffects ipw / teffects ipwra / teffects aipw | — |
| 5 | 熵平衡 | ebalance | — |
| 5 | 计数 + 固定效应 | ppmlhdfe | poisson / nbreg |
| 5 | Heckman | heckman | 手动两步法 |
| 5 | 分位数 | qreg / sqreg | bsqreg |
| 6 | 野生聚类自助法 | boottest | bootstrap, cluster() |
| 6 | 随机化推断 | ritest | permute |
| 6 | 多重检验 | rwolf | wyoung |
| 6 | TWFE 诊断 | bacondecomp | — |
| 6 | 平行趋势敏感性分析 | honestdid | — |
| 6 | Oster δ* | psacalc delta | 手动计算公式 |
| 7 | 交互效应边际效应 | margins + marginsplot | — |
| 7 | 中介分析 | medsem / khb | 手动 Baron–Kenny 法 |
| 7 | SEM / 路径分析 | sem / gsem | — |
| 8 | 回归表 | esttab（来自 estout） | outreg2 / xtable |
| 8 | 系数图 / 事件研究 | coefplot | marginsplot |
| 8 | 分箱散点图 | binscatter | 手动使用 twoway |
| 8 | RD 图 | rdplot | — |
常见错误（以及正确做法）
| 错误 | 正确做法 |
|---|---|
在没有 FE 的面板数据上使用 reg y x1 x2 | reghdfe y x1 x2, absorb(unit time) vce(cluster unit) |
| 对聚类数据使用默认的 iid SE | vce(cluster id)；如果聚类数 < 50，则使用 boottest |
在错位实施的情况下使用 TWFE（带处理虚拟变量的 reghdfe） | 使用 csdid / eventstudyinteract / did_imputation |
在大型面板（>10k 个单位）上使用 xtreg, fe | 改用 reghdfe — 速度更快，并支持多维 FE 和多路聚类 |
使用 merge 时不指定 assert() | 始终指定 assert(match using master) 或类似选项 |
未先运行 xtset，就使用 L.x 生成滞后项 | 在进行任何时间序列操作前运行 xtset id time |
使用 ivregress 时不运行 estat firststage | 对于弱 IV 诊断，优先使用 ivreg2 — 它会自动报告 CD、KP 和 AR |
| RD 只使用单个带宽 | rdrobust 默认会报告 MSE 最优带宽和敏感性；还应展示 rdbwselect，并将带宽减半/加倍 |
使用 psmatch2 后不检查平衡性 | 使用 pstest, both 报告匹配前/后的 SMD；目标为 <10% |
直接解读 logit 系数 | 始终运行 margins, dydx(*) — 这是可解读的量 |
只报告 estimates 并手动复制粘贴 | 使用 esttab / outreg2 自动导出为 LaTeX/Word/RTF |
| 只报告核心系数（没有表 2） | 务必提供包含 M1→M6 的多列主表 — 这是经济学论文的核心，而不是摘要中的一句话 |
| 系数表没有任何图表 | 经济学结果至少需要 F1 趋势图 + F2 事件研究图 + F3 系数图 + F4 敏感性分析 — 参见默认输出规范 |
使用 graph save（仅 .gph）保存图表 | 还要使用 graph export ..., replace 导出为 .pdf / .png |
| 不记录会话日志 | 在 do 文件顶部运行 log using main.log, replace |
在回归中使用 gen x2 = x^2，而不是 c.x##c.x | 因子变量表示法与 margins 配合良好 |
典型的 .do 文件框架
* ============================================================
* main.do — reproducible 8-step pipeline
* ============================================================
version 17
clear all
set more off
capture log close
log using "logs/main.log", replace text
cd "/path/to/project"
* 1. Clean -----------------------------------------------------
do "code/01_clean.do" // produces data/analysis.dta
* 2. Transform -------------------------------------------------
do "code/02_transform.do" // adds engineered vars
* 3. Describe --------------------------------------------------
do "code/03_describe.do" // writes tables/table1_*, figures/*_trend.pdf
* 4. Diagnose --------------------------------------------------
do "code/04_diagnose.do" // writes logs/diagnostics.log
* 5. Model -----------------------------------------------------
do "code/05_model.do" // saves e() for m1–m6 via eststo
* 6. Robustify -------------------------------------------------
do "code/06_robust.do"
* 7. Further ---------------------------------------------------
do "code/07_further.do"
* 8. Export ----------------------------------------------------
do "code/08_tables_figures.do" // writes tables/*.tex, figures/*.pdf
log close
每个 .do 文件都会加载 data/analysis.dta，将输出写入 tables/ 或 figures/，并保存日志。运行干净的 make clean && do main.do 即可从原始数据重新生成整篇论文。
回归表（esttab）实用指南（单页配方索引）
eststo + esttab 是 AER 论文中所有多回归表格背后的基础工具。上述八种模式对应如下：
| 模式 | 各列之间的变化项 | 步骤 |
|---|---|---|
| A. 逐步增加控制变量 | 协变量集合 / FE 深度 | 5.A — 表 2 |
| B. 识别策略比较 | 识别策略（OLS / IV / DID / PSM / EB） | 5.B — 表 2-bis |
| C. 多个结果变量 | 因变量 Y | 5.C — 表 2-ter |
| D. 纵向排列面板 A / B | 时间跨度 / 样本（面板行 × 设定列） | 5.D — 表 2-quater |
| E. IV 报告三联表 | 第一阶段 / 简化式 / 2SLS | 5.E — 表 2-quinto |
| F. 因果分析编排 | 1 列，完整诊断（csdid/synth_runner/teffects） | 5.F |
| G. 子组表 | 子样本（全体 / 女性 / 男性 / Q1…Q4） | 7.b — 表 3 |
| H. 稳健性检验总表 | 汇总所有稳健性检验 | 6.k — 表 A1 |
AER 风格的默认 esttab 设置：
esttab m1 ... mN using "tables/tableN.tex", ///
replace ///
se star(* 0.10 ** 0.05 *** 0.01) /// AER stars convention
label booktabs ///
mtitles("(1)" "(2)" ...) /// column labels
stats(N r2 r2_a, labels("N" "R²" "Adj. R²")) ///
addnotes("Cluster-robust SE in parentheses." ///
"* p<0.10, ** p<0.05, *** p<0.01.")
* Word version: same call with .rtf extension.
* For multi-panel paper bundles, chain esttab with `, append` (see 5.D / 8.5.4).
Figure factory（Stata 中的 12 张标准 AER 图）
| # | 图表 | Stata 命令 | 章节 |
|---|---|---|---|
| 1a | 原始趋势（DID Figure 1） | collapse (mean) y, by(year treat) → twoway line | §1（步骤 3.e） |
| 1b | 处理组纳入时间热力图 | heatplot first_treat unit, color(...) | §1 |
| 2a | 事件研究系数 | eventstudyinteract → coefplot, keep(rel_d*) | §3（步骤 3.5.1） |
| 2a' | Bacon 权重 | bacondecomp y treat, ddetail → 内置图表 | §3 |
| 2a'' | CS-DID 动态效应 | csdid → csdid_plot | §3 |
| 2b | 第一阶段散点图 | binscatter D Z, controls(X) | §3（步骤 3.5.2） |
| 2c | RD 经典图 | rdplot Y X, c(0) | §3（步骤 3.5.3） |
| 2c' | McCrary 密度图 | rddensity X, c(0) plot | §3 |
| 2d | 匹配平衡性 Love 图 | pstest X, both graph（在 psmatch2 后运行） | §3（步骤 3.5.4） |
| 2e | SCM 轨迹图 | synth y X, trunit() trperiod() fig · sdid 图表 | §3（步骤 3.5.5） |
| 3 | 主规格系数图 | coefplot m1...m6, keep(D) | §4（步骤 8.f） |
| 4a | 剂量－反应关系 | xtile dose10 = D, nq(10) → margins i.dose10 → marginsplot | §5（步骤 7.g） |
| 4b | 分组边际效应 | margins, dydx(D) at(group=(0 1)) → marginsplot | §5 |
| 5 | 规格曲线 | 手动编写规格循环 → twoway rcap（参见 6.l） | §7 |
| 6 | 敏感性分析仪表板 | honestdid, coefplot · psacalc plot · evalue | §7（步骤 6.m） |
| 7 | 最终主图 | 特定估计量（csdid_plot、synth fig、rdplot） | §8 |
每张 Stata 图均导出为
.png ≥ 300 dpi两种格式（分别用于 LaTeX 和幻灯片／网页）。在master.do开头运行一次set scheme s2color（或使用schemepack中的white_tableau），以保持样式一致。
方法目录
经典 OLS／面板
reg log_wage training age edu, vce(cluster firm_id) // OLS
areg log_wage training age edu, absorb(industry) vce(cluster firm_id) // OLS + 1 FE
xtreg log_wage training age edu, fe // panel FE
xtreg log_wage training age edu, re // panel RE (Hausman)
reghdfe log_wage training age edu, absorb(worker_id year) vce(cluster firm_id) // HD FE workhorse
reghdfe log_wage training, absorb(worker_id year) vce(cluster firm_id year) // 2-way cluster
ppmlhdfe count training, absorb(firm_id year) cluster(firm_id) // Poisson + FE
heckman log_wage training, select(in_lf = age edu marital kids) // Heckman selection
qreg log_wage training age, quantile(0.5) // quantile reg
sqreg log_wage training age, quantile(0.1 0.25 0.5 0.75 0.9) reps(500)
双重差分
* 2×2
reghdfe log_wage c.treat#c.post, absorb(worker_id year) vce(cluster worker_id)
* Staggered — modern stack
csdid log_wage age edu, ivar(worker_id) time(year) gvar(first_treat) method(dripw) agg(group)
eventstudyinteract log_wage rel_d*, cohort(first_treat) control_cohort(never_treated) absorb(worker_id year) vce(cluster worker_id)
did_imputation log_wage worker_id year first_treat, allhorizons pretrend(5)
did_multiplegt_dyn log_wage worker_id year training, effects(5) placebo(3)
sdid log_wage worker_id year training, vce(bootstrap) graph
* Diagnostics + sensitivity
bacondecomp log_wage training, ddetail // TWFE diagnostic
honestdid, pre(rel_d1-rel_d3) post(rel_d5-rel_d9) mvec(0(0.1)0.5) // PT sensitivity
工具变量／2SLS
ivregress 2sls log_wage age edu (training = Z), vce(cluster firm_id) // baseline
ivreg2 log_wage age edu (training = Z), cluster(firm_id) first endog(training) // workhorse: full diagnostics
ivreghdfe log_wage age (training = Z), absorb(industry year) cluster(firm_id) first
* Weak-IV diagnostics
estat firststage // (after ivregress)
* ivreg2 reports CD / KP / AR automatically.
回归不连续
rdrobust outcome running_var, c(0) kernel(triangular) bwselect(mserd) // Sharp RD (CCT)
rdrobust outcome running_var, c(0) fuzzy(treatment) // Fuzzy RD
rddensity running_var, c(0) plot // McCrary density
rdmc outcome running_var, cutoffs(0 5 10) // multi-cutoff
rdplot outcome running_var, c(0)
匹配／重新加权
psmatch2 training age edu tenure, out(log_wage) n(1) common // PSM (canonical)
teffects psmatch (log_wage) (training age edu tenure), atet // PSM (modern)
teffects ipw (log_wage) (training age edu tenure), atet // IPW
teffects ipwra (log_wage age edu tenure) (training age edu tenure) // IPWRA
teffects aipw (log_wage age edu tenure) (training age edu tenure) // AIPW (DR)
ebalance training age edu tenure // entropy balancing
合成控制
synth log_wage age edu, trunit(1) trperiod(2015) fig keep(synth.dta, replace)
synth_runner log_wage age edu, trunit(1) trperiod(2015) gen_vars // + 50 placebos
sdid log_wage worker_id year training, vce(bootstrap) graph
机器学习因果推断（模式 B——完整流程见 §B）
ddml init partial, kfolds(5) reps(5) // DML setup
ddml E[Y|X]: pystacked Y X*, type(reg) methods(rf gradboost lassocv)
ddml E[D|X]: pystacked D X*, type(reg) methods(rf gradboost lassocv)
ddml crossfit
ddml estimate, robust
crforest log_wage training age edu, ntrees(2000) honesty // causal forest
predict cate, te
pdslasso log_wage training (X1-X100), cluster(firm_id) partial // post-double-selection
稳健性、敏感性与推断
boottest training, cluster(state) reps(9999) seed(42) // wild cluster bootstrap
ritest training _b[training], reps(1000) seed(0): reghdfe ... // randomization inf.
rwolf y1 y2 y3, indepvar(training) controls(X) reps(500) ... // Romano-Wolf MTC
psacalc delta training, mcontrol(X) rmax(1.3) // Oster δ
honestdid, pre(...) post(...) mvec(0(0.1)0.5) // RR PT sensitivity
evalue rr, est(1.45) lcl(1.10) ucl(1.91) // E-value
xtcsd, pesaran abs // cross-sec dependence
生存分析／流行病学（模式 A——见 §A）
stset time, failure(event)
sts graph, by(A) ci risktable // Kaplan-Meier
stcox A age edu // Cox
streg A age edu, distribution(weibull) time // AFT
strmst2 A, tau(1825) covariates(age edu) // RMST contrast
gformula ... // parametric g-formula
eltmle event A age edu, tmle // TMLE for binary outcome
mrrobust, beta_outcome(by) se_outcome(byse) beta_exposure(bx) se_exposure(bxse) // MR-IVW
mregger by bx [aw=1/byse^2], gxse(bxse) // MR-Egger
何时移交给其他技能
- Agent 原生的单次导入 Python 工作流（
import statspai as sp）→00-StatsPAI_skill。 - 明确使用 Python 传统技术栈 →
00.1-Full-empirical-analysis-skill。 - Mixtape 风格的跨语言代码模板（Python/R/Stata 并列）→
10-Jill0099-causal-inference-mixtape。 - Stata 语法参考（完整命令目录、Mata 等）→
32-dylantmoore-stata-skill。 - 用于会计研究的 Stata →
18-jusi-aalto-stata-accounting-research。 - 分析完成后的论文撰写 → 本仓库中面向写作的技能（
04-*-scientific-writer、08-*-web-latex等）。
此技能的范围 止于步骤 8——润色完成的 .tex 表格和 .pdf 图表。论文撰写不在范围之内。
步骤 8.5 — 移交给手稿（退出约定）
此技能的职责仍止于步骤 8。不过，“止于”不应意味着“走入死胡同”： 写作阶段需要机器可读的约定，而不是一堆散乱的文件。 停止前，请在表格和图表旁生成两份产物。
1. exhibits_index.md — 每个编号的展示项占一行：
| 展示项 | 文件 | 支持的论点 | 生成者 |
|---|---|---|---|
| 表 2 | tables/table2_main.{tex,docx,xlsx} | 政策对对数输出的主要影响 | analysis.do §baseline |
| 图 2 | figures/fig2_event_study.{png,pdf} | 不存在差异化的事前趋势 | analysis.do §event-study |
2. results_summary.json — 每个主要估计值都作为一条扁平记录：
{"baseline": {"coef": 0.1234, "se": 0.0412, "p": 0.0031, "n": 12345, "spec": "TWFE + region-year FE"}}
为什么两者都要有，以及为什么采用这种格式：
- 手稿阶段通过包含指令引用展示项，绝不重新输入数值 — 在 LaTeX 中使用
\input{tables/table2_main}，在 Markdown 中使用{{ include: table2_main }}。索引可告知作者哪个文件支持哪个论点。 - 论文中打印的每个数字都必须能以其显示精度追溯到结果文件。
results_summary.json让这一点可供核验，而非仅凭信任。 如果某个系数只存在于渲染后的表格中，就无法对照引用它的正文进行核验。
然后移交：
- 完整论文 → 从头到尾生成 Word
.docx/ LaTeX →skills/69-Paper-WorkFlow/。它负责协调第 5–9 阶段（起草 → 润色 → 去 AI 痕迹 → 审阅 → 投稿），其scripts/assemble_manuscript_docx.py会基于这两份产物，组装生成09_submission/main.docx— 将正文、上述每个展示项、图表和参考文献列表汇集到一个文件中。在第 0 阶段选择正文格式：交付物为 Word 时使用 Markdown（Markdown →.docx的保真度高；LaTeX →.docx会有损），投稿平台接收.tex时使用 LaTeX。 - 仅起草正文 → 使用本仓库中的写作技能（
04-*-scientific-writer、35-*-academic-writing-skills）。 - 仅进行 LaTeX / Quarto 排版 →
08-ndpvt-web-latex-document-skill、60-regisely-superpapers。
评论
支持 Markdown