Imported from brycewang-stanford/Auto-Empirical-Research-Skills (plugins/empirical-analysis-stata/skills/pipeline/SKILL.md). Install upstream with npx skills add brycewang-stanford/Auto-Empirical-Research-Skills --skill pipeline. Copyright stays with the author.

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/ 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:

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:

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

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

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

* 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):

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.

Export cookbook — LaTeX / Word / RTF in one block

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

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:

* 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:

* 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):

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.

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.

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

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)

* (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 — 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.

* 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 — 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.

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

*--- 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 — 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.

* 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.

* 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.

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?

* (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)

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.

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 — 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).

* 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:

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

Deeper patterns: 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.

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.

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.

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.

* 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.

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.

* 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_

*Truncated - read the full file at https://github.com/brycewang-stanford/Auto-Empirical-Research-Skills/blob/266adb2ce05d8be2bfa490db70e9d8d2ef53e4a0/plugins/empirical-analysis-stata/skills/pipeline/SKILL.md.*