00.2 Full Empirical Analysis Skill Stata

---
name: Full-empirical-analysis-skill-Stata
description: Classical end-to-end empirical analysis workflow in the traditional Stata ecosystem — native Stata + reghdfe + ivreg2 + csdid + did_imputation + eventstudyinteract + sdid + rdrobust + rddensity + synth + synth_runner + psmatch2 + teffects + ebalance + coefplot + esttab + asdoc + binscatter. **Defaults to economics empirical-paper style** (AER / QJE / AEJ) — every run produces a publication-ready output set with a multi-column regression table (M1→M6 progressive controls/FE) as the centerpiece, plus Table 1 (descriptives), mechanism / heterogeneity / robustness tables, and event-study + coefficient + trend figures. Covers the full 8-step Stata pipeline an applied economist runs on every paper — (1) data import & cleaning (use/import, destring, misstable, duplicates, merge assert), (2) variable construction (gen/egen/winsor2/xtile/xtset with L./F./D.), (3) descriptive statistics & Table 1 (tabstat/balancetable/asdoc), (4) classical diagnostic tests (sktest/swilk/hettest/imtest/xtserial/xttest3/vif/dfuller/kpss/hausman/estat overid), (5) baseline modeling (reg/xtreg/reghdfe/ivreg2/ivregress/csdid/did_imputation/eventstudyinteract/sdid/rdrobust/synth/psmatch2/teffects/heckman/qreg/ppmlhdfe), (6) robustness battery (bacondecomp/honestdid/rwolf/ritest/wildbootstrap/oster), (7) further analysis (subgroup/triple-diff/interactions/medsem/marginsplot/binscatter by group), (8) publication-ready tables & figures (esttab/outreg2/estout/coefplot/marginsplot/rdplot/twoway combined). **Also covers two parallel domain modes that share the same 8-step scaffolding** — **Mode A — Epidemiology / public health** (target-trial emulation, IPTW + g-formula + TMLE doubly-robust triplet via `teffects ipw` / `teffects ipwra` / `teffects aipw` / `eltmle`, Mendelian randomization via `mrrobust` (IVW / Egger / weighted median) and `mregger` / `mrpresso`, KM / Cox / AFT / RMST survival via `sts` / `stcox` / `streg` / `strmst2`, E-value sensitivity via `evalue` (Linden-Mathur), principal stratification — STROBE / TRIPOD reporting), and **Mode B — ML causal inference** (DML via `ddml` / `pdslasso`, S/T/X/R/DR meta-learners via `crforest` and `ddml interactive`, causal forest via `crforest` / `cforest`, BART/BCF via `bart` / `bartCause`-style externals, CATE distribution + policy tree via `crforest`, off-policy evaluation, conformal causal externals, fairness audit, DAG learning via `pcalg` / external Python callouts). Use when the user asks for a complete Stata empirical analysis, wants a reproducible .do-file pipeline, needs a Stata counterpart to the Python StatsPAI / Full-empirical-analysis-skill, or names a specific Stata step in isolation ("run reghdfe with two-way clustering", "csdid event study", "winsor2 at 1%", "esttab to LaTeX", "coefplot with CI", "ivreg2 weak-IV test", "synth_runner placebos", "teffects psmatch balance check"). Mode A triggers on "target trial emulation Stata", "teffects ipw aipw", "eltmle", "mrrobust", "mregger weighted median", "stcox AFT survival", "strmst2", "evalue Stata", "STROBE Stata", "公共健康 Stata", "流行病学 Stata". Mode B triggers on "ddml Stata", "pdslasso", "crforest causal forest Stata", "policy tree Stata", "因果机器学习 Stata".
triggers:
  - Stata empirical analysis
  - full Stata pipeline
  - reproducible do-file
  - Stata do-file workflow
  - reghdfe two-way FE
  - high-dimensional fixed effects Stata
  - ivreg2 weak instruments
  - ivregress 2sls liml gmm
  - csdid Callaway SantAnna
  - did_imputation Borusyak
  - eventstudyinteract Sun Abraham
  - sdid synthetic DID Stata
  - rdrobust Stata
  - rddensity manipulation test
  - synth synthetic control
  - synth_runner placebo
  - psmatch2 propensity score
  - teffects psmatch
  - teffects ipwra AIPW Stata
  - ebalance entropy balancing
  - xtreg fe re hausman
  - ppmlhdfe Poisson
  - quantile regression qreg
  - heckman selection model
  - esttab publication table
  - outreg2 LaTeX
  - estout coefplot
  - marginsplot interaction
  - rdplot binned scatter
  - binscatter Stata
  - bacondecomp Goodman Bacon
  - honestdid Rambachan Roth
  - wild cluster bootstrap boottest
  - ritest randomization inference
  - rwolf Romano-Wolf
  - oster delta
  - winsor2 winsorize Stata
  - tabstat table 1
  - balancetable Stata
  - misstable patterns
  - destring dates
  - xtset panel
  # Mode A — Epidemiology / public health
  - epidemiology pipeline Stata
  - public health causal inference Stata
  - target trial emulation Stata
  - teffects ipw aipw ipwra
  - g-formula Stata
  - eltmle TMLE Stata
  - HAL-TMLE Stata
  - Mendelian randomization Stata
  - mrrobust mregger
  - MR-PRESSO Stata
  - MR-Egger weighted median Stata
  - STROBE TRIPOD reporting Stata
  - evalue sensitivity Stata
  - Kaplan-Meier AFT survival Stata
  - sts stcox streg strmst2
  - 流行病学 Stata
  - 公共健康 Stata
  # Mode B — ML causal inference
  - ML causal inference Stata
  - ddml double machine learning Stata
  - pdslasso ivlasso
  - crforest causal forest Stata
  - cforest Stata
  - meta-learner S T X R DR Stata
  - CATE distribution Stata
  - policy tree Stata
  - off-policy evaluation Stata
  - conformal causal Stata
  - causal discovery PC Stata
  - 因果机器学习 Stata
---

# Full Empirical Analysis — Classical Stata Workflow

This skill is the *canonical* 8-step pipeline an applied economist runs on every empirical paper, written in the **traditional Stata ecosystem** — native Stata + the 20+ community commands that have become de-facto standards (`reghdfe`, `ivreg2`, `csdid`, `did_imputation`, `eventstudyinteract`, `sdid`, `rdrobust`, `rddensity`, `synth`, `synth_runner`, `psmatch2`, `teffects`, `ebalance`, `coefplot`, `esttab`, `outreg2`, `boottest`, `ritest`, `rwolf`, `bacondecomp`, `honestdid`, `binscatter`).

**Companion skills**: if the user wants the same pipeline in Python, route to `00-StatsPAI_skill` (agent-native DSL) or `00.1-Full-empirical-analysis-skill` (explicit Python stack). **This skill is the Stata counterpart** — every step produces a `.do` file you can hand to a journal's replication office or a co-author who refuses to leave Stata.

## Philosophy

1. **Stata idioms, not Python-translated.** `reghdfe`, not "statsmodels analogue of reghdfe". `esttab`, not "Stata's stargazer".
2. **Reproducible .do files.** Every code block below is runnable after `use data.dta, clear`. No Jupyter, no notebooks — just do-files and log files.
3. **Full pipeline, not just regressions.** Stata users historically over-invest in Step 5 (modeling) and under-invest in Steps 1–4 and 6–8. This skill treats them as first-class.
4. **Rich outputs.** Every step yields at least one table (`.tex`/`.rtf`) or figure (`.pdf`/`.png`) — never a coefficient printed to the Results window and forgotten.
5. **Progressive disclosure.** `SKILL.md` gives the canonical command at each step; [`references/`](references/) holds variant-specific depth (dozens of tests, estimator-specific diagnostics, graph recipes).

## SkillOpt-style execution gate

Use this long playbook as a seed skill, not as a script to exhaustively apply. SkillOpt discipline: treat each local `.do`-file change as a candidate patch that must beat a selection check and survive a held-out check before it becomes reusable boilerplate. Before writing or revising a `.do` file, compress the user's request into a task-local `best_skill` card:

```text
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>
```

1. **Route card**: record the mode (`econ`, `epi`, or `ml-causal`), estimand, identification design, focal outcome/treatment, Stata command family, and required artifacts.
2. **Bounded edit**: change one decision at a time (sample rule, estimator command, fixed effects, clustering, export format, or robustness check). Prefer the smallest patch that can pass validation.
3. **Selection split discipline**: treat the user's immediate failure or requested artifact as the selection split. Reserve at least one alternate outcome, sample window, estimator command family, or export target as the held-out gate.
4. **Held-out gate**: define checks before running code: `isid`/duplicate status, merge assertions, treatment support, missingness thresholds, expected table/figure files, and one non-focal robustness/specification that the change must not break.
5. **Reject buffer**: if a candidate spec fails the gate, log the failure, `.do` diff, and gate output in `analysis_log.md` or the Stata log; revert to the last passing spec and do not retry the same unchecked pattern.
6. **Slow/meta update**: at the end of the task, write down `accepted_patterns` and `rejected_patterns` from the trajectory. Do not widen the canonical `.do` template from a single passing run.
7. **Promote only after validation**: only turn a one-off fix into reusable `.do` boilerplate after it passes the current data and at least one alternate outcome/sample/specification.

---

## Three domain modes (default = AER econ; alternates = epi & ML-causal)

The default playbook above is **AER-style applied econometrics** — the AEA convention: written-out estimating equation, identifying assumption, design horse-race, full robustness gauntlet. The skill **also** ships two parallel sub-pipelines for the other two big causal-inference traditions, each reusing the same Steps 1–4 (cleaning / construction / Table 1 / diagnostics) and Step 8 (tables/figures) — only Step 5 (estimator) and Step 6/7 swap commands:

| Mode | Reader convention | Step-5 estimator stack | Reporting stack | Jump to |
|---|---|---|---|---|
| **Default — Applied Econ (AER / QJE / AEJ)** | "Show the equation + identifying assumption + design horse-race; controls visible; clustered SE" | DID / IV / RD / SCM / matching / `reghdfe` HDFE | AER house-style multi-column `esttab` / `outreg2` / `coefplot` + 8-section paper layout | Steps 1 → 8 (entire playbook below) |
| **Mode A — Epidemiology / Public Health** | "STROBE / TRIPOD-AI; target trial protocol; doubly-robust estimand; absolute & relative risk; KM survival" | Target-trial emulation · IPTW (`teffects ipw`) · IPWRA / AIPW (`teffects ipwra` / `teffects aipw`) · g-formula (`gformula`) · TMLE (`eltmle`) · Mendelian randomization (`mrrobust` IVW / `mregger` / `mrpresso`) · KM/Cox/AFT (`sts`/`stcox`/`streg`/`strmst2`) | Same `esttab` + risk-difference / hazard-ratio / E-value rows | §A. Epidemiology pipeline |
| **Mode B — ML Causal Inference** | "DML / meta-learners / causal forest / DR-learner; CATE distribution; policy value" | DML (`ddml` / `pdslasso`) · S/T/X/R/DR-Learner (`ddml interactive`) · GRF causal forest (`crforest` / `cforest`) · BART / BCF (external Python via `python_user`) · matrix completion (external) | `esttab` ML horse-race + `crforest` CATE plot + policy-value table | §B. ML causal pipeline |

**How to invoke a non-default mode** (Claude / agent picks this up from the user's wording):

| User says... | Mode the skill switches to |
|---|---|
| "Run a DID / IV / RD / event study", "AER table", "applied micro" | Default (AER econ) — Steps 1 → 8 |
| "Target trial emulation", "g-formula", "IPTW", "TMLE", "Mendelian randomization", "STROBE / TRIPOD", "公共健康 / 流行病学", "epi pipeline", "RWE study", "cohort study", "case-control" | Mode A (Epi) — §A |
| "DML", "double machine learning", "ddml", "causal forest", "crforest", "meta-learner", "CATE", "policy learning", "ML causal", "因果机器学习" | Mode B (ML causal) — §B |
| "Mix" (e.g. "estimate DID + then ML CATE on the heterogeneity") | Default + Mode B in sequence — every estimator stores results via `eststo`, drop them all into one `esttab` for the horse-race column |

The three modes share **the same Step 1–4 cleaning / Table 1 / diagnostics scaffolding, the same Step 8 export stack, and the same DAG-first identification logic** — switching modes only changes which Step-5 command family you reach for. If you only want descriptive stats / Table 1 / a balance check, the AER `tabstat` / `balancetable` / `asdoc` calls in Step 3 work identically across all three modes.

> **Stata-specific caveat for Mode B**: Stata's first-party ML-causal coverage is thinner than Python/R. For Dragonnet / TARNet / CEVAE / cfcausal / fairness audit, call out to Python via Stata 18's `python:` / `python script` block (or shell out to a sister `.py`) and read the result back via `frame` or `import delimited`. The skill prefers native Stata commands (`ddml`, `pdslasso`, `crforest`) where they exist, and explicitly marks the Python callouts in §B.

---

## Default Output Spec — Economics Empirical Paper

This skill defaults to the **applied-economics paper convention**. Unless the user explicitly asks for a single point estimate, every `.do`-file run produces the full publication-ready output set below. Treat it as the contract of Step 8 — **mandatory**, not opt-in.

### Required tables (always produced)

| # | Table | Stata source | Saves to |
|---|---|---|---|
| **T1** | Summary statistics & balance (treated vs control, with SMD / p-values) | `balancetable` + `asdoc sum` (Step 3) | `tables/table1_balance.{tex,rtf,xlsx,docx}` |
| **T2** ★ | **Main results — multi-column regression M1→M6** (progressive controls + FE) | `eststo` 6 specs → `esttab` (Step 5–6) | `tables/table2_main.{tex,rtf,xlsx,docx}` |
| **T3** | Mechanism / outcome ladder — same treatment, 3+ outcomes side-by-side | loop `eststo: reghdfe` over outcomes → `esttab` (Step 7) | `tables/table3_mechanism.{tex,rtf,xlsx,docx}` |
| **T4** | Heterogeneity — subgroup × main coef (gender, age, region, …) | subgroup `eststo` + `suest` Wald → `esttab` (Step 7) | `tables/table4_heterogeneity.{tex,rtf,xlsx,docx}` |
| **T5** | Robustness battery — alt SE / cluster / sample / placebo, in **one** table | `eststo` × variants → `esttab` (Step 6) | `tables/table5_robustness.{tex,rtf,xlsx,docx}` |

> **★ Table 2 is the centerpiece of every economics paper.** It is the multi-column regression table that walks the reader from raw correlation (M1) to the fully-specified design (M6: 2-way FE + interacted FE + cluster-robust SE). Do **not** collapse it into a single column. Do **not** report only the headline coefficient. The progression *is* the credibility argument: if M1→M6 is monotone and stable, the design is plausibly identifying; if it collapses on adding FE, that *is* the result.
>
> **Canonical 6 columns, in order:**
> 1. **M1** raw bivariate (`reg y treat`)
> 2. **M2** + demographics (`+ age + edu`)
> 3. **M3** + sector controls (`+ tenure / firm_size`)
> 4. **M4** + unit FE (`reghdfe ..., absorb(unit)`)
> 5. **M5** + 2-way FE (`absorb(unit year)`)
> 6. **M6** + interacted FE (`absorb(unit year i.industry#i.year)`) with `vce(cluster unit)`

### Required figures (always produced)

| # | Figure | Stata source | Saves to |
|---|---|---|---|
| **F1** | Trend / motivation — treated vs control over time, with policy line | `collapse (mean) y, by(year treat)` → `twoway line` (Step 3) | `figures/fig1_trend.pdf` (+ `.png`) |
| **F2** | Event-study coefficients with 95% CI, base period at –1 | `eventstudyinteract` / `csdid` / `coefplot, keep(*.rel)` (Step 5) | `figures/fig2_event_study.pdf` |
| **F3** | Coefficient plot across specs M1→M6 | `coefplot m1 m2 m3 m4 m5 m6, keep(treat) vertical` (Step 8) | `figures/fig3_coefplot.pdf` |
| **F4** | Robustness / sensitivity curve — `bacondecomp` plot, `honestdid` plot, or cluster-comparison forest | scenario-specific (Step 6) | `figures/fig4_sensitivity.pdf` |

### Output file layout (default)

```
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}
```

Every table → `.tex` (LaTeX `booktabs`) **and** `.rtf` (Word) **and** `.xlsx` (Excel) **and** `.docx` (Word OOXML). Every figure → `.pdf` (vector for LaTeX) **and** `.png` at ≥300 dpi (slides / web).

### When to deviate

- **Single quick estimate** — produce only the relevant cell, but warn that the standard deliverable is the full set above and offer to run it.
- **Design does not support a figure** (cross-section → no event study) — skip with a printed `display` note explaining why; do **not** silently drop.
- **N=1 treated unit (`synth`)** — replace F1/F2 with the SCM trajectory + placebo distribution from `synth_runner`; T1–T5 still apply.

---

## Required packages

```stata
* 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
```

---

## The 8 Steps — Canonical Pipeline (mapped to AER paper sections)

```
┌──────────────────────────────────────────────────────────────────────┐
│ 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     │
└──────────────────────────────────────────────────────────────────────┘
```

The 8 steps mirror the canonical sections of an applied AER / QJE / AEJ paper. Each step is one paper section and emits a paper-ready artifact on disk:

```
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
```

Below is the canonical command at each step. **All examples share one running narrative** — a labor-economics panel where `training` (treatment) affects `log_wage` (outcome), with covariates `age`, `edu`, `tenure`, panel keys `worker_id` / `firm_id` / `year`. Variable names and parameter values are **illustrative**; substitute the real ones from the user's dataset. Only command names and option *shapes* are normative.

> **When a step has many variants** (e.g. staggered DID has 5 estimators; heteroskedasticity has 4 classic tests), SKILL.md shows the one you reach for first and links to `references/NN-<topic>.md` for the rest. **Read the reference file when the user's case doesn't fit the default.**

---

## Paper-ready figure & table inventory (what to produce by section)

A modern AER paper has **5–7 figures** and **3–5 main tables** + an appendix robustness table. Every step below leaves at least one numbered artifact on disk. Default file names assume parallel `.tex` / `.rtf` exports (the agent should produce both so co-authors can edit in Word, and the build system can use LaTeX):

| § | Artifact | Stata primitive | Filenames |
|---|---|---|---|
| §1 | **Figure 1**: raw trends / treatment rollout | `collapse` + `twoway line` · `heatplot` for staggered rollout | `figures/fig1_trend.{pdf,png}` |
| §1 | **Table 1**: summary stats (full / treated / control + Δ + SMD) | `balancetable` · `asdoc sum, by()` · `tabstat` | `tables/table1_balance.{tex,rtf,xlsx,docx}` |
| §3 | **Figure 2**: identification graphic (event-study / first-stage / McCrary / RD scatter / SCM trajectory) | `coefplot` after `eventstudyinteract`/`csdid` · `binscatter` · `rdplot` · `rddensity, plot` · `synth` | `figures/fig2_event_study.{pdf,png}` |
| §4 | **Table 2**: main results — progressive controls M1→M6 | `eststo` 6 specs → `esttab` | `tables/table2_main.{tex,rtf,xlsx,docx}` |
| §4 | **Table 2-bis**: design horse-race (OLS / IV / DID / matching) | `eststo` mix + `esttab` | `tables/table2b_designs.{tex,rtf,xlsx,docx}` |
| §4 | **Figure 3**: coefficient plot across specs | `coefplot m1 m2 m3 m4 m5 m6, keep(treat)` | `figures/fig3_coefplot.{pdf,png}` |
| §5 | **Table 3**: heterogeneity by subgroup | `eststo` per slice + `esttab` + `suest` Wald | `tables/table3_heterogeneity.{tex,rtf,xlsx,docx}` |
| §5 | **Figure 4**: dose-response / margins-by-quartile | `xtile` + `margins` + `marginsplot` | `figures/fig4_dose.{pdf,png}` |
| §6 | **Table 4**: mechanism / outcome ladder | loop `eststo: reghdfe` over outcomes → `esttab` | `tables/table4_mechanism.{tex,rtf,xlsx,docx}` |
| §7 | **Table A1**: robustness master (one column per check) | `eststo` × variants → `esttab` | `tables/tableA1_robustness.{tex,rtf,xlsx,docx}` |
| §7 | **Figure 5**: spec curve — coefficient + 95% CI across all specs | hand-rolled spec loop + `twoway rcap` | `figures/fig5_spec_curve.{pdf,png}` |
| §7 | **Figure 6**: sensitivity (HonestDiD / Oster / E-value) | `honestdid, coefplot` · `psacalc plot` · `evalue` table | `figures/fig6_sensitivity.{pdf,png}` |
| §8 | **Replication bundle**: all tables in one document | `esttab ..., append` to one `.tex` / `.rtf` · `texdoc` | `replication/paper_tables.{tex,rtf,xlsx,docx}` |

> Every Stata estimator above stores results via `eststo` and can be passed straight into `esttab` / `coefplot` / `outreg2`. Don't hand-roll LaTeX, and don't render Word from `outsheet`/`putexcel` matrices — `esttab` and `outreg2` apply book-tab borders, AER-style stars, and the right SE label automatically. For deeper export recipes (LaTeX / Word / Markdown variants, `texdoc`, `frmttable`), see [`references/08-tables-plots.md`](references/08-tables-plots.md).

---

## Export cookbook — LaTeX / Word / RTF in one block

Stata's export stack is more fragmented than Python's StatsPAI. Three tiers, picked by **scope**:

| Tier | Use when | API | Hot options |
|---|---|---|---|
| **1. Single multi-column table** | Exporting *one* Table 2 / Table 3 / Table A1 with progressive columns | **`.tex` / `.rtf`**: `esttab` with `booktabs`<br>**`.xlsx` / `.docx`**: `outreg2` (esttab's native Office export is a data dumper, not a publication formatter — outreg2 produces properly formatted Word/Excel with borders, aligned stars, and AER-style layout) | **`esttab`** for tex/rtf: `keep()`, `drop()`, `mtitles()`, `stats(N r2 r2_a, labels(...))`, `star(* 0.10 ** 0.05 *** 0.01)`, `label`, `booktabs`, `addnotes()`.<br>**`outreg2`** for xlsx/docx: `label dec(3)`, `addtext()` for FE indicators, `replace` first / `append` subsequent columns.<br>Always emit all four formats — esttab for tex/rtf in a `foreach ext in tex rtf` loop, outreg2 for xlsx/docx in a separate `foreach ext in xlsx docx` loop. |
| **2. Multi-panel paper format** (Tables 2 + 3 + A1 + A2 in one file) | Producing the paper-tables block — main + heterogeneity + robustness + placebo as a single document | `esttab ... using "paper.tex", replace` for first panel; subsequent `esttab ... using "paper.tex", append` for each next panel; `texdoc init "paper.tex"` for full LaTeX with prose. Repeat for `.rtf`/`.xlsx`/`.docx` bundles. | first panel: `replace`; subsequent: `append`; surround with `texdoc` for headings |
| **3. Full session bundle** (the Stata 17+ `collect` equivalent) | Replication appendix that mixes summary stats + balance + multiple regression tables + headings + prose in **one** file | `collect create paper`<br>`collect get summary, ...`<br>`collect get est ...`<br>`collect layout ...`<br>`collect export "paper.xlsx"` (also `.docx`/`.html`/`.tex`/`.md`) | Stata 17+ only; for older Stata use `texdoc` / `markdoc` |

**Journal styling — pick the right `star` levels and SE label.** The AEA convention is `* 0.10 ** 0.05 *** 0.01` and SE label "Standard errors in parentheses"; QJE / Econometrica / RES variants only differ in stars / notes / fonts. Define an `esttab` wrapper once at the top of the do-file:

```stata
* 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)
```

**Multi-format export pattern** — emit `.tex` and `.rtf` via `esttab`, then `.xlsx` and `.docx` via `outreg2`. Never use `esttab ... using "... .xlsx"` or `esttab ... using "... .docx"` — esttab's native Office export is a raw data dumper, not a publication formatter:

```stata
* 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)
```

**Figures always export to both `.pdf` and `.png` at ≥300 dpi** (vector for LaTeX, raster for slides/web/Word embedding):

```stata
graph export "figures/fig1_trend.pdf", replace
graph export "figures/fig1_trend.png", replace width(2400) height(1800)
```

For the `collect` / `texdoc` / `markdoc` cookbook (Stata 17+ multi-panel paper bundle and prose+tables PDF/Word generation), see [`references/08-tables-plots.md`](references/08-tables-plots.md).

---

## Step −1 — Pre-Analysis Plan (pre-data; AEA RCT Registry style)

Before touching the data, write down (a) the population, (b) the design, (c) the **minimum detectable effect (MDE)** under the planned sample size and α=0.05, β=0.20. Persist the result as `protocol.do` so a referee can verify the design was powered before, not after, the data were seen.

```stata
* 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
```

Save this `.do` file in version control **before** running Step 1. AEA RCT Registry / OSF preregistration tools accept it as the analysis-plan exhibit.

---

## Step 0 — Sample-construction log & 5-check data contract

An AER §1 *Data* section has three jobs: (a) describe sources, (b) document **every** sample restriction (the "footnote 4" sample log), (c) lock the panel structure. Stata users typically chain `keep` / `drop` commands in a single do-file with no count log, then can't reconstruct the analysis sample for the response letter. The data contract is the cure.

### 0.1 Sample-construction log (footnote 4)

```stata
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
```

Paste the `display` lines verbatim as footnote 4 of the paper.

### 0.2 Five-check data contract (go / no-go gate)

```stata
* (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
```

If any assertion fires, **stop** and fix it. Stata estimators silently drop NaN rows, the most common source of "mysterious sample-size shrinkage" bugs in the response letter.

---

### Step 1 — Data import & cleaning

Deeper patterns: [references/01-data-cleaning.md](references/01-data-cleaning.md) — reading Excel/CSV/SAS/SPSS, `destring` on numeric-looking strings, `misstable` patterns, `duplicates` tagging, `merge` with `assert(match using)`, `xtset` balance checks, spells / gaps, labels.

```stata
* 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
```

**Key principle**: all row exclusions are explicit, counted, logged. No command in Steps 2+ should silently drop rows.

---

### Step 2 — Variable construction & transformation

Deeper patterns: [references/02-data-transformation.md](references/02-data-transformation.md) — log / ihs / Box–Cox, within-group `winsor2`, `xtile` and custom cuts, `egen` recipes, time-series operators (`L.`, `F.`, `D.`, `S.`), CPI deflation, staggered-DID timing construction.

```stata
* 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
```

---

### Step 2.5 — Empirical strategy (write the equation + identifying assumption)

This is the heart of an AER paper. **Before any code**, write down the equation explicitly and state the identifying assumption. Vague identification language is the single most common reason a referee rejects an applied paper. Persist the strategy as `strategy.do` (or `.md`) so it is a dated, version-controlled artifact — *not* a post-hoc rationalization written after seeing the coefficient.

#### Equation × identifying assumption × Stata estimator (decision table)

| Design | Estimating equation | Identifying assumption | Stata estimator |
|---|---|---|---|
| 2×2 DID | `Y_it = α_i + λ_t + β·D_it + X'γ + ε_it` | parallel trends conditional on X | `reghdfe Y c.D#c.post X, absorb(i t) vce(cluster i)` |
| Event-study (CS / SA) | `Y_it = α_i + λ_t + Σ_{e≠-1} β_e · 1{t-G_i = e} + ε_it` | no anticipation + group-time PT | `csdid` / `eventstudyinteract` / `did_imputation` |
| 2SLS | `Y_i = α + β·D_i + X'γ + ε_i;  D_i = π·Z_i + X'δ + u_i` | exclusion + relevance + monotonicity | `ivreg2` / `ivreghdfe` |
| Sharp RD | `Y_i = α + β·1{X_i ≥ c} + f(X_i) + ε_i` (local poly) | continuity of E[Y(0)\|X] at c, no manipulation | `rdrobust` (+ `rddensity`) |
| SCM | `Ŷ_1t(0) = Σ_j ŵ_j Y_jt`, τ_t = `Y_1t − Ŷ_1t(0)` for t≥T_0 | pre-period fit + interpolation validity | `synth` / `synth_runner` / `sdid` |
| Selection-on-observables (matching/IPW) | `E[Y_i(d) \| X_i] = E[Y_i \| D=d, X_i]` | unconfoundedness + overlap | `teffects psmatch / ipw / aipw / ipwra`, `ebalance` |

#### Design picker (when the user is unsure)

```
                 ┌─ 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)
```

#### Pre-registration `strategy.do` template

```stata
*--- 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)
```

Commit `strategy.do` **before** running Step 5 / Step 6. The `git log` of this file *is* the analysis plan.

---

### Step 3 — Descriptive statistics & Table 1

Deeper patterns: [references/03-descriptive-stats.md](references/03-descriptive-stats.md) — stratified Table 1 with SMDs, `asdoc` / `tabstat` / `balancetable`, correlation matrix with significance stars (`pwcorr, sig star(.05)`), histograms / kdensity by group, `xtline` for DID motivation, panel-coverage `xtdescribe`.

```stata
* 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
```

---

### Step 3.5 — Identification graphics (Section "Identification, graphical evidence")

**AER convention: the identification figure precedes the regression table.** The reader should see graphical evidence that PT holds / first stage is strong / RD jumps cleanly *before* you ask them to trust your point estimate. This block is one or two figures saved to `figures/` — most of the heavy lifting is one Stata command + one `coefplot`.

#### 3.5.1 Event-study figure + numerical pre-trends test (DID identification)

Pre-period coefficients ≈ 0 (with the −1 reference period normalized to zero) is the visual evidence for parallel trends. Pair the **figure** with a **numerical** pre-trends test so reviewers don't have to eyeball it.

```stata
* 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 First-stage F-statistic + scatter (IV identification)

Rule of thumb: first-stage F ≥ 10 for OLS-style inference; F ≥ 23 for AR-equivalent inference (Stock–Yogo / Lee 2022). `ivreg2` reports CD / KP / weak-IV statistics by default — far more useful than `ivregress`'s minimal output.

```stata
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 density + canonical RD plot

The signature RD figure is `rdplot` (CCT-style binned scatter with local-polynomial fit on each side), paired with the McCrary manipulation test. Together they answer: (a) is there a visual jump? (b) is the density continuous at the cutoff?

```stata
* (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 Matching: love plot (standardized differences pre vs post)

```stata
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: synthetic-control trajectory + gap plot

For synthetic-control designs the canonical Figure 2 is the treated-vs-synthetic time series with treatment time annotated.

```stata
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
```

> Identification-specific checks (PT for DID, weak-IV F, density for RD, common support for matching) **are also auto-run inside the Step-5 estimators** — don't duplicate the numerics here, but DO produce the figures: a referee scans the figures first.

---

### Step 4 — Diagnostic statistical tests

Deeper patterns: [references/04-statistical-tests.md](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), Augmented Dickey–Fuller, KPSS, Phillips–Perron, Hausman (FE vs RE), Sargan–Hansen, `estat firststage` (weak IV).

```stata
* 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
```

**Decision table**:

| Test | Null | If rejected |
|------|------|-------------|
| `swilk` / `sktest` | residuals Normal | large N: usually ignore; small N: bootstrap |
| `estat hettest` / `imtest, white` | homoskedastic | use `vce(robust)` or `vce(cluster id)` |
| `xtserial` / `xttest3` | no panel autocorr / no groupwise hetero | cluster by unit |
| `xtcsd, pesaran` | no cross-sectional dependence | Driscoll–Kraay or `xtscc` |
| `estat vif` > 10 | — | drop/combine collinear regressors |
| `estat ovtest` | specification OK | add polynomials / logs |
| `dfuller` reject + `kpss` fail to reject | stationary | keep levels |
| `dfuller` fail to reject | unit root | first-difference or cointegrate |
| `hausman` | RE consistent | use FE |

---

### Step 5 — Baseline empirical modeling (Section 4: Main Results)

Deeper patterns: [references/05-modeling.md](references/05-modeling.md) — every classical estimator with syntax: `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`.

This is the densest section of an applied paper. A modern AER §4 typically contains **2–3 multi-regression tables and one coefficient plot**:

- **Table 2** (main): progressive controls, 4–6 columns — **Pattern A** below
- **Table 2-bis** (design horse race): same coefficient under OLS / IV / DID / matching — **Pattern B**
- **Table 2-ter** (multi-outcome): same treatment, several outcomes side-by-side — **Pattern C**
- **Figure 3** (coefplot): visual summary of β̂ and 95% CI across specs

> **Estimator routing** (memorize this — getting it wrong silently produces nonsense):
> - **No FE / single low-card FE** → `reg y x1 x2, vce(cluster id)`
> - **High-dim FE** → `reghdfe y x1 x2, absorb(fe1 fe2) vce(cluster id)`
> - **Two-way cluster** → `reghdfe ..., vce(cluster fe1 fe2)`
> - **2SLS / IV** → `ivreg2 y x (D = Z), cluster(id) first endog(D)` (or `ivreghdfe` for HD FE + IV)
> - **DID / event-study** → `csdid` / `eventstudyinteract` / `did_imputation`

**Pick the estimator by identification strategy**:

```
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
```

Canonical commands (a Stata equivalent of `outreg2` / `esttab` is the workhorse — `eststo` 5–6 specs, then `esttab` consolidates them into one table). Key options:

```
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 Pattern A — Progressive controls (the canonical Table 2)
Stable β̂ across columns ⇒ less concern that selection on observables is driving the estimate (Oster 2019 selection-stability logic; quantified in Step 6.j). **`eststo m1...m6` + `esttab` is the Stata equivalent of `outreg2` and R's `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 convention: show all controls — pass NEITHER `keep()` NOR `drop()`** so every parameter is visible. Use `keep(training)` only when a focal-coefficient-only table is intentional (e.g. interaction-form heterogeneity, IV first-stage triplet); use `drop(_cons)` only when you want to suppress the constant for paper aesthetics.

#### 5.B Pattern B — Design horse race (Table 2-bis)
Show the same coefficient of interest under multiple identification strategies. This is *the* AER credibility move: convergent evidence across designs each making different identifying assumptions.

```stata
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 Pattern C — Multi-outcome table (same X, several Y's)
A single treatment, several outcomes. Use `mtitles` so each column carries the Y name.

```stata
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 Pattern D — Stacked Panel A / Panel B table
Same model family, two horizons (short-run / long-run) or two samples (pre-2015 / post-2015). Stack vertically with two `esttab` calls + `texdoc` glue, OR use `esttab ..., refcat()` to inject panel headers.

```stata
* 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 Pattern E — IV reporting triplet (first-stage / reduced-form / 2SLS)
The textbook AER IV table presents the **first stage**, the **reduced form**, and the **2SLS** in three columns so the reader can verify Wald-ratio = RF / FS.

```stata
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 triplet is intentionally focal:** show only Z + endogenous regressor so the reader can eyeball the Wald ratio. Drop `keep(Z training)` only if a referee asks for the full coefficient list.

#### 5.F Pattern F — Causal-orchestrator main via `csdid` / `synth_runner` / `teffects`
For DID / SCM / matching mains, the modern Stata estimator returns a self-contained estimate + automatic placebos / pre-trends / overlap diagnostics. Use the estimator's own report, then pipe into `eststo` + `esttab`.

```stata
* 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 Pattern G — Subgroup `esttab` (Table 3, see Step 7)
One column per subgroup. Detailed code in §Step 7 — Heterogeneity.

#### 5.H Pattern H — Robustness master (Table A1, see Step 6)
Stack every robustness specification next to the baseline. Detailed code in §Step 6 — Robustness master.

---

#### Canonical estimator commands (the underlying primitives)

```stata
* 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)
```

---

### Step 6 — Robustness battery

Deeper patterns: [references/06-robustness.md](references/06-robustness.md) — progressive specs M1–M6 with `eststo` / `esttab`, `bacondecomp` for TWFE bias, `honestdid` (Rambachan–Roth), `boottest` wild-cluster bootstrap, `ritest` randomization inference, `rwolf` Romano–Wolf, alternative cluster levels, Oster δ*, placebo timing, subsample splits, specification curve via loops.

```stata
* 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 Pattern H — Robustness master table (Table A1, one column per check)

The canonical AER appendix Table A1 stacks every robustness specification next to the baseline so reviewers see at a glance that β̂ survives. Build the list of `eststo` results dynamically:

```stata
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 Specification curve (Simonsohn–Simmons–Nelson 2020) in Stata

`esttab` doesn't have a native spec-curve; loop over {controls × samples × outcome transforms × SE types}, store each β̂ + 95% CI in a frame, then plot.

```stata
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
```

A modern AER referee letter often asks "what about specification X?" — the spec curve answers all such asks at once.

#### 6.m Sensitivity dashboard (HonestDiD + Oster + E-value, one block)

For DID main results, produce one figure (`HonestDiD` partial-identification bound) + one row of sensitivity stats (Oster δ + E-value).

```stata
* (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.
```

---

### Step 7 — Further analysis (mechanism / heterogeneity / mediation / moderation)

Deeper patterns: [references/07-further-analysis.md](references/07-further-analysis.md) — factor-variable interactions, `margins` / `marginsplot`, triple-difference (DDD), outcome ladder, `medsem` / Baron–Kenny, `khb` (Karlson–Holm–Breen) for non-linear mediation, dose-response via `xtile` or splines, subgroup event studies.

```stata
* 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
```

---

### Step 8 — Publication tables & figures

> **This step is mandatory** — every analysis run produces all 5 required tables (T1–T5) and all 4 required figures (F1–F4) defined in the *Default Output Spec* at the top of this skill. Do not skip Step 8 because "the regression already ran". A coefficient without a table and a figure is not how applied economics communicates a result.

Deeper patterns: [references/08-tables-plots.md](references/08-tables-plots.md) — `esttab` and `outreg2` for regression tables; `coefplot` for coefficient plots, event studies, forest plots; `marginsplot` for interactions; `binscatter` for residualized scatter; `rdplot` for RD; `twoway` + combined graphs; Stata scheme / graph options; LaTeX / Word / Excel export.

```stata
* ============================================================
* 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-)
```

**Deliverables checklist** (verify before declaring the run complete):

```
[ ] 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. Reproducibility stamp

The single artifact a journal's replication office (or a future co-author) needs to reproduce the headline number. Persist Stata version, seed, dataset hash, baseline coefficient + CI, and pointers to the protocol/contract:

```stata
* 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
```

Commit `artifacts/result.json` alongside the paper PDF. A referee should be able to run `make clean && do main.do` and bit-identically reproduce this JSON.

---

## §A — Epidemiology / Public Health Mode

When the user's wording flags Mode A (target-trial emulation / IPTW / TMLE / MR / STROBE / 流行病学 / 公共健康 / RWE / cohort), the 8 steps still apply — but Step 5 swaps the OLS-and-FE stack for the doubly-robust + survival + MR triplet, and the deliverables follow STROBE / TRIPOD-AI conventions. **Steps 1–4 (cleaning, construction, Table 1, diagnostics) and Step 8 (tables/figures export) are identical to the Default mode.**

**Command footprint** (install on top of the Default Stata stack):

```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 Cohort construction + target-trial protocol

Write the protocol **before** touching the data. Save it as `protocol.do` and quote it in the paper.

```stata
*--- 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 Table 1 by exposure (identical to Default Step 3)

Use the same `balancetable` / `asdoc tabstat` from Step 3, just `by(A)`. E-values for unmeasured confounding go in the footer.

```stata
balancetable A age edu smoke bmi ldl sbp using tables/tableA1_strobe.tex, ///
    replace ctitles("Untreated" "Treated" "Diff (SE)") pvalues
```

### A.2 DAG + propensity-score overlap (positivity check)

```stata
* 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-formula + TMLE doubly-robust triplet (Step 5 swap)

The "AER Table 2" of epi: a 3-row table where each row is one of {IPTW-MSM via `teffects ipw`, IPWRA / AIPW via `teffects aipw`, TMLE via `eltmle`}, so the reader can confirm doubly-robust agreement.

```stata
* 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 Survival outcomes — KM / Cox / AFT / RMST

```stata
* 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 Mendelian randomization (IVW / Egger / weighted-median triplet)

```stata
* 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 Robustness — E-value / bounds / principal stratification

```stata
* 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 reporting checklist

Save as `replication/strobe_checklist.do` (or `.md`) and tick before submission:

```
[ ] 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 Causal Inference Mode

When the user's wording flags Mode B (DML / meta-learner / causal forest / CATE / policy learning / 因果机器学习), the pipeline keeps Steps 1–4 and Step 8 from the Default mode, swaps Step 5 for the ML estimator stack, and adds a CATE-distribution + policy-value layer between Step 7 and Step 8.

> **Native vs Python callout**: Stata 17+ ships first-class ML-causal commands for DML and causal forests. For neural causal (Dragonnet / TARNet / CEVAE), conformal causal, and fairness audit, the skill calls out to Python via Stata 18's `python:` block — clearly marked at each step.

**Command footprint** (install on top of the Default Stata stack):

```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 Train/holdout split + nuisance super-learner

DML / TMLE / Causal Forest are doubly-robust *if* the nuisance learners (outcome regression `Q(X,A) = E[Y|X,A]` and propensity `g(A|X) = Pr(A=1|X)`) converge faster than `n^(-1/4)`. Use `pystacked` (Stata 18 + Python callout) to stack GBM / RF / Lasso / NN under cross-validation — that is the Stata equivalent of StatsPAI's SuperLearner.

```stata
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 / estimand declaration (optionally LLM-assisted)

```stata
* 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 Estimator stack — DML · meta-learners · causal forest (Step 5 swap)

The "AER Table 2" of ML causal: a horse-race table where each column is one estimator family on the same `(Y, A, X)` data — readers want to see DML (PLR), DML (interactive), and causal forest all agree (or disagree) on the ATE.

```stata
*--- 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 distribution + subgroup CATE plot (Step 7 extension)

```stata
* 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 Policy learning + off-policy evaluation

```stata
* 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 Uncertainty (conformal causal) + fairness + sensitivity

```stata
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-causal-specific reporting checklist

Save as `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)
```

---

## Library / command cheat-sheet

| Step | Task | Go-to command | Fallback |
|------|------|---------------|----------|
| 1 | Import | `use` / `import excel` / `import delimited` / `import sas` / `import spss` | `usespss` / `infile` |
| 1 | Missing | `misstable summarize` / `mdesc` | `missings report` |
| 1 | Merge | `merge m:1 ... assert(match using)` | — |
| 1 | Panel | `xtset` / `xtdescribe` | `tsset` for pure TS |
| 2 | Winsorize | `winsor2` | hand-roll via `summ, detail` + `replace` |
| 2 | Lag / lead / diff | `L.` / `F.` / `D.` (require `xtset`) | `bys id (t): gen lx = x[_n-1]` |
| 3 | Table 1 | `tabstat` + `balancetable` + `asdoc sum` | manual loop |
| 3 | Correlations | `pwcorr, sig star(.05)` | `corr` |
| 4 | Hetero | `estat hettest` / `estat imtest, white` | — |
| 4 | Panel serial | `xtserial` / `xttest3` | — |
| 4 | Stationarity | `dfuller` / `kpss` / `pperron` / `xtunitroot` | — |
| 4 | Hausman | `hausman fe re, sigmamore` | — |
| 5 | OLS / panel FE | `reghdfe` | `areg` / `xtreg, fe` |
| 5 | IV | `ivreg2` / `ivreghdfe` | `ivregress 2sls` |
| 5 | DID — 2×2 | `reg` / `reghdfe` with `c.treated#c.post` | — |
| 5 | DID — CS | `csdid` | — |
| 5 | DID — SA | `eventstudyinteract` | — |
| 5 | DID — imputation | `did_imputation` | `did_multiplegt_dyn` |
| 5 | DID — SDID | `sdid` | — |
| 5 | RD | `rdrobust` / `rddensity` / `rdplot` | `rdmc` for multi-cutoff |
| 5 | Synthetic control | `synth` / `synth_runner` | — |
| 5 | PSM | `psmatch2` / `teffects psmatch` | — |
| 5 | IPW / AIPW | `teffects ipw` / `teffects ipwra` / `teffects aipw` | — |
| 5 | Entropy balancing | `ebalance` | — |
| 5 | Count + FE | `ppmlhdfe` | `poisson` / `nbreg` |
| 5 | Heckman | `heckman` | manual 2-step |
| 5 | Quantile | `qreg` / `sqreg` | `bsqreg` |
| 6 | Wild cluster boot | `boottest` | `bootstrap, cluster()` |
| 6 | Randomization inf | `ritest` | `permute` |
| 6 | Multiple testing | `rwolf` | `wyoung` |
| 6 | TWFE diagnosis | `bacondecomp` | — |
| 6 | PT sensitivity | `honestdid` | — |
| 6 | Oster δ* | `psacalc delta` | manual formula |
| 7 | Interaction margins | `margins` + `marginsplot` | — |
| 7 | Mediation | `medsem` / `khb` | manual Baron–Kenny |
| 7 | SEM / path | `sem` / `gsem` | — |
| 8 | Regression table | `esttab` (from `estout`) | `outreg2` / `xtable` |
| 8 | Coefplot / eventstudy | `coefplot` | `marginsplot` |
| 8 | Binscatter | `binscatter` | `twoway` manual |
| 8 | RD plot | `rdplot` | — |

---

## Common mistakes (and what to do instead)

| Mistake | Correct approach |
|---------|------------------|
| `reg y x1 x2` on panel data with no FE | `reghdfe y x1 x2, absorb(unit time) vce(cluster unit)` |
| Default iid SEs on clustered data | `vce(cluster id)`; `boottest` if clusters < 50 |
| TWFE (`reghdfe` with treatment dummy) on staggered adoption | Use `csdid` / `eventstudyinteract` / `did_imputation` |
| `xtreg, fe` with large panel (>10k units) | Switch to `reghdfe` — faster and allows multi-dim FE + multi-way cluster |
| `merge` without `assert()` | Always specify `assert(match using master)` or similar |
| Generating lag as `L.x` without `xtset` first | Run `xtset id time` before any time-series op |
| `ivregress` without `estat firststage` | For weak-IV diagnostics, prefer `ivreg2` — reports CD, KP, AR automatically |
| RD with a single bandwidth | `rdrobust` default reports MSE-optimal + sensitivity; also show `rdbwselect` and halve/double |
| `psmatch2` without balance check | Use `pstest, both` to report pre/post SMDs; target <10% |
| Interpreting `logit` coefficients directly | Always run `margins, dydx(*)` — that's the interpretable quantity |
| Reporting only `estimates` and manual copy-paste | `esttab` / `outreg2` → LaTeX/Word/RTF automatically |
| Reporting only the headline coefficient (no Table 2) | **Always** ship the multi-column M1→M6 main table — that is the centerpiece of an economics paper, not the abstract sentence |
| Coefficient table without any figures | An economics result needs **at least** F1 trend + F2 event study + F3 coefplot + F4 sensitivity — see the Default Output Spec |
| Saving graphs with `graph save` (`.gph` only) | Also `graph export ..., replace` to `.pdf` / `.png` |
| Not logging the session | `log using main.log, replace` at top of do-file |
| `gen x2 = x^2` instead of `c.x##c.x` inside regressions | Factor-variable notation plays nicely with `margins` |

---

## Typical .do-file skeleton

```stata
* ============================================================
* 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
```

Every `.do` file loads `data/analysis.dta`, writes its outputs to `tables/` or `figures/`, and saves a log. A clean `make clean && do main.do` regenerates the entire paper from raw data.

---

## Regtable (esttab) cookbook (one-page recipe index)

`eststo` + `esttab` is the single primitive behind every multi-regression table in an AER paper. The eight patterns above map to:

| Pattern | What varies across columns | Step |
|---|---|---|
| **A. Progressive controls** | covariate set / FE depth | 5.A — Table 2 |
| **B. Design horse race** | identification strategy (OLS / IV / DID / PSM / EB) | 5.B — Table 2-bis |
| **C. Multi-outcome** | dependent variable Y | 5.C — Table 2-ter |
| **D. Stacked Panel A / B** | horizon / sample (panel rows × spec columns) | 5.D — Table 2-quater |
| **E. IV reporting triplet** | first stage / reduced form / 2SLS | 5.E — Table 2-quinto |
| **F. Causal-orchestrator** | 1 column, full diagnostics (`csdid`/`synth_runner`/`teffects`) | 5.F |
| **G. Subgroup table** | subsample (full / female / male / Q1…Q4) | 7.b — Table 3 |
| **H. Robustness master** | every robustness check stacked | 6.k — Table A1 |

Default `esttab` settings for AER house style:

```stata
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 (the 12 standard AER figures in Stata)

| # | Figure | Stata commands | Section |
|---|---|---|---|
| 1a | Raw trends (DID Figure 1) | `collapse (mean) y, by(year treat)` → `twoway line` | §1 (Step 3.e) |
| 1b | Treatment rollout heatmap | `heatplot first_treat unit, color(...)` | §1 |
| 2a | Event-study coefficients | `eventstudyinteract` → `coefplot, keep(rel_d*)` | §3 (Step 3.5.1) |
| 2a' | Bacon weights | `bacondecomp y treat, ddetail` → built-in graph | §3 |
| 2a'' | CS-DID dynamic effects | `csdid` → `csdid_plot` | §3 |
| 2b | First-stage scatter | `binscatter D Z, controls(X)` | §3 (Step 3.5.2) |
| 2c | RD canonical plot | `rdplot Y X, c(0)` | §3 (Step 3.5.3) |
| 2c' | McCrary density | `rddensity X, c(0) plot` | §3 |
| 2d | Matching love plot | `pstest X, both graph` (after `psmatch2`) | §3 (Step 3.5.4) |
| 2e | SCM trajectory | `synth y X, trunit() trperiod() fig` · `sdid` graph | §3 (Step 3.5.5) |
| 3 | Coefficient plot of main specs | `coefplot m1...m6, keep(D)` | §4 (Step 8.f) |
| 4a | Dose-response | `xtile dose10 = D, nq(10)` → `margins i.dose10` → `marginsplot` | §5 (Step 7.g) |
| 4b | Margins-by-subgroup | `margins, dydx(D) at(group=(0 1))` → `marginsplot` | §5 |
| 5 | Specification curve | hand-rolled spec loop → `twoway rcap` (see 6.l) | §7 |
| 6 | Sensitivity dashboard | `honestdid, coefplot` · `psacalc plot` · `evalue` | §7 (Step 6.m) |
| 7 | Final main figure | estimator-specific (`csdid_plot`, `synth` fig, `rdplot`) | §8 |

> Every Stata figure is exported as **both** `.pdf` (for LaTeX) and `.png ≥ 300 dpi` (for slides / web). Use `set scheme s2color` (or `white_tableau` from `schemepack`) once at the top of `master.do` for consistent styling.

---

## Method Catalog

### Classical OLS / Panel
```stata
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)
```

### Difference-in-Differences
```stata
* 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
```

### Instrumental Variables / 2SLS
```stata
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.
```

### Regression Discontinuity
```stata
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)
```

### Matching / Reweighting
```stata
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
```

### Synthetic Control
```stata
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
```

### ML Causal (Mode B — see §B for the full pipeline)
```stata
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
```

### Robustness, Sensitivity & Inference
```stata
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
```

### Survival / Epi (Mode A — see §A)
```stata
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
```

---

## When to hand off to other skills

- **Agent-native single-import Python workflow** (`import statspai as sp`) → `00-StatsPAI_skill`.
- **Explicit Python traditional stack** → `00.1-Full-empirical-analysis-skill`.
- **Mixtape-style cross-language code templates** (Python/R/Stata side-by-side) → `10-Jill0099-causal-inference-mixtape`.
- **Stata syntax reference** (complete command catalog, Mata, etc.) → `32-dylantmoore-stata-skill`.
- **Stata for accounting research** → `18-jusi-aalto-stata-accounting-research`.
- **Paper drafting** after the analysis is done → the writing-oriented skills in this repo (`04-*-scientific-writer`, `08-*-web-latex`, etc.).

This skill's remit **ends at Step 8** — polished `.tex` tables and `.pdf` figures. Paper drafting is out of scope.

---

## Step 8.5 — Handoff to the manuscript (exit contract)

This skill's remit still ends at Step 8. "Ends" should not mean "dead-ends", though:
the writing stage needs a machine-readable contract, not a folder of loose files.
Before you stop, emit two artifacts next to the tables and figures.

**1. `exhibits_index.md`** — one row per numbered exhibit:

| Exhibit | Files | Claim it supports | Produced by |
|---|---|---|---|
| Table 2 | `tables/table2_main.{tex,docx,xlsx}` | Main effect of the policy on log output | `analysis.do` §baseline |
| Figure 2 | `figures/fig2_event_study.{png,pdf}` | No differential pre-trend | `analysis.do` §event-study |

**2. `results_summary.json`** — every headline estimate as a flat record:

```json
{"baseline": {"coef": 0.1234, "se": 0.0412, "p": 0.0031, "n": 12345, "spec": "TWFE + region-year FE"}}
```

Why both, and why in these shapes:

- The manuscript stage cites exhibits by **include directive**, never by retyped
  numbers — `\input{tables/table2_main}` in LaTeX, `{{ include: table2_main }}` in
  Markdown. The index is what tells the writer which file backs which claim.
- Every number a paper prints must trace back to a result file at the precision it is
  displayed. `results_summary.json` is what makes that checkable instead of trusted.
  A coefficient that exists only inside a rendered table cannot be verified against
  the prose that cites it.

**Then hand off:**

- **Full paper → Word `.docx` / LaTeX, end to end** → `skills/69-Paper-WorkFlow/`. It
  orchestrates Stage 5–9 (draft → polish → de-AIGC → review → submission) and its
  `scripts/assemble_manuscript_docx.py` assembles `09_submission/main.docx` — body,
  every exhibit above, figures and reference list in one file — from exactly these two
  artifacts. Pick the body format at its Stage 0: **Markdown when the deliverable is
  Word** (Markdown → `.docx` is high fidelity; LaTeX → `.docx` is lossy), LaTeX when
  the venue takes `.tex`.
- **Prose drafting only** → the writing skills in this repo (`04-*-scientific-writer`,
  `35-*-academic-writing-skills`).
- **LaTeX / Quarto typesetting only** → `08-ndpvt-web-latex-document-skill`,
  `60-regisely-superpapers`.