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","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","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
Stata idioms, not Python-translated.reghdfe, not "statsmodels analogue of reghdfe". esttab, not "Stata's stargazer".
Reproducible .do files. Every code block below is runnable after use data.dta, clear. No Jupyter, no notebooks — just do-files and log files.
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.
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.
Progressive disclosure.SKILL.md gives the canonical command at each step; references/ holds variant-specific depth (dozens of tests, estimator-specific diagnostics, graph recipes).
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>
Route card: record the mode (econ, epi, or ml-causal), estimand, identification design, focal outcome/treatment, Stata command family, and required artifacts.
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.
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.
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.
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.
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.
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.
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:
"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:
M1 raw bivariate (reg y treat)
M2 + demographics (+ age + edu)
M3 + sector controls (+ tenure / firm_size)
M4 + unit FE (reghdfe ..., absorb(unit))
M5 + 2-way FE (absorb(unit year))
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
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:
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
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:
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 .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(). outreg2 for xlsx/docx: label dec(3), addtext() for FE indicators, replace first / append subsequent columns. 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 collect get summary, ... collect get est ... collect layout ... 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:
* 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):
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.
* 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
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)
┌─ 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.dobefore 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
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
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.
* 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)
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_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)
* 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 — 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.
* 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:
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.
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).
* (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.
Deeper patterns: 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.
* 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 — 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.
* ============================================================
* 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" ///
