Full Empirical Analysis — Classical Python Workflow

This skill is the canonical 8-step pipeline an applied economist runs on every empirical paper, written in the traditional Python ecosystem — no opinionated one-stop wrapper. Every step calls libraries directly (pandas, numpy, scipy, statsmodels, linearmodels, pyfixest, rdrobust, econml, causalml, matplotlib, seaborn), so the agent — or the user reading the agent's code — has full visibility and can swap any component.

Companion skill: if the user prefers a single-import agent-native DSL (import statspai as sp), route to 00-StatsPAI_skill instead. This skill is the opposite philosophy: everything explicit, everything inspectable, every diagnostic run by hand, every plot shaped by the user.

Philosophy

1. Traditional stack, no magic. Agents should be able to read every line and know exactly which library / estimator / standard error family is at work.
2. Full pipeline, not just estimation. 80% of the time on a real paper is steps 1–4 and 6–8. This skill treats them as first-class, not an afterthought.
3. Rich outputs. Every step produces at least one table or figure — never a single point estimate in isolation.
4. Progressive disclosure. SKILL.md gives the canonical call at each step; references/ holds variant-specific depth (dozens of tests, estimator-specific diagnostics, plot recipes).
5. Reproducible. Every code block is runnable after pip install -r requirements.txt and df = pd.read_csv(...).

SkillOpt-style execution gate

Use this long playbook as a seed skill, not as a script to exhaustively apply. SkillOpt discipline: treat each local analysis-code 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 analysis code, 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, package stack, and required artifacts.
2. Bounded edit: change one decision at a time (sample rule, estimator, 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 family, or export target as the held-out gate.
4. Held-out gate: define checks before running code: row counts, key uniqueness, 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, code diff, and gate output in analysis_log.md; 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 project template from a single passing run.
7. Promote only after validation: only turn a one-off fix into reusable project 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 / descriptives / diagnostics) and Step 8 (tables & figures) — only Step 5 (estimator) and Step 6/7 (robustness / mechanism) swap libraries:

Mode	Reader convention	Step-5 estimator stack	Reporting stack	Jump to
Default — Applied Econ (AER / QJE / AEJ)	"Show the equation + identifying assumption + design horse-race; controls visible; clustered SE"	DID / IV / RD / SCM / matching / pyfixest.feols HDFE	AER house-style multi-column pf.etable / Stargazer + 8-section paper layout	Steps 1 → 8 (entire playbook below)
Mode A — Epidemiology / Public Health	"STROBE / TRIPOD-AI; target trial protocol; doubly-robust estimand; absolute & relative risk; KM survival"	Target-trial emulation · IPTW (zepid) · g-formula (zepid / hand-rolled) · TMLE (zepid.causal.gformula / econml) · Mendelian randomization (pymr / rpy2+TwoSampleMR) · KM/AFT (lifelines)	Same pf.etable + risk-difference / hazard-ratio / E-value rows	§A. Epidemiology pipeline
Mode B — ML Causal Inference	"DML / meta-learners / causal forest / DR-learner; CATE distribution; policy value"	DML (econml.dml / doubleml) · S/T/X/R/DR-Learner (econml.metalearners / causalml) · GRF causal forest (econml.grf) · Dragonnet/TARNet/CEVAE (causalml) · BCF (pymc-bart/bcf-py) · matrix completion	pf.etable ML horse-race + CATE plot + policy-value table + conformal PI (mapie)	§B. ML causal pipeline

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

User says...	Mode the skill switches to
"Run a DID / IV / RD / event study", "AER table", "applied micro"	Default (AER econ) — Steps 1 → 8
"Target trial emulation", "g-formula", "IPTW", "TMLE", "Mendelian randomization", "STROBE / TRIPOD", "公共健康 / 流行病学", "epi pipeline", "RWE study", "cohort study", "case-control"	Mode A (Epi) — §A
"DML", "double machine learning", "causal forest", "meta-learner", "CATE", "Dragonnet", "BCF", "policy learning", "conformal causal", "fairness audit", "ML causal", "uplift modeling", "因果机器学习"	Mode B (ML causal) — §B
"Mix" (e.g. "estimate DID + then ML CATE on the heterogeneity")	Default + Mode B in sequence — every estimator returns a coefficient + SE pair, drop them all into one pf.etable(...) 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 estimator family you reach for, not the surrounding paper structure. If you only want descriptive stats / Table 1 / a balance check, the AER tableone / gtsummary-style calls in Step 3 work identically across all three modes.

---

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 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	Source / library	Saves to
T1	Summary statistics & balance (treated vs control, with SMD / p-values)	pandas.describe + custom table1() (Step 3)	tables/table1_balance.xlsx + .docx + .tex
T2 ★	Main results — multi-column regression M1→M6 (progressive controls + FE)	pyfixest.feols × 6 specs → pf.etable() / Stargazer (Step 5–6)	tables/table2_main.xlsx + .docx + .tex
T3	Mechanism / outcome ladder — same treatment, 3+ outcomes side-by-side	feols looped over y ∈ {Y1, Y2, Y3, Y_main} → pf.etable	tables/table3_mechanism.xlsx + .docx + .tex
T4	Heterogeneity — subgroup × main coef (gender, age, region, …)	subgroup feols × Wald → pf.etable (Step 7)	tables/table4_heterogeneity.xlsx + .docx + .tex
T5	Robustness battery — alt SE / alt cluster / alt sample / placebo, in one table	feols × variants → pf.etable (Step 6)	tables/table5_robustness.xlsx + .docx + .tex

★ 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 (y ~ treat)
2. M2 + demographics (+ age + edu)
3. M3 + sector controls (+ tenure / firm_size / industry)
4. M4 + unit FE (| worker_id)
5. M5 + 2-way FE (| worker_id + year)
6. M6 + interacted FE (| worker_id + year + industry^year) with cluster-robust SE

Required figures (always produced)

#	Figure	Source / library	Saves to
F1	Trend / motivation — treated vs control over time, with policy line	df.groupby([year, treat])[y].mean().unstack().plot() (Step 3)	figures/fig1_trend.png (300 dpi, 必须导出 PNG) + .pdf
F2	Event-study coefficients with 95% CI, base period at –1	pf.feols(... ~ i(rel_time, ref=-1) ...) → pf.iplot() (Step 5)	figures/fig2_event_study.png (300 dpi, 必须导出 PNG) + .pdf
F3	Coefficient plot across specs M1→M6	matplotlib.errorbar over the 6 fitted models (Step 8)	figures/fig3_coefplot.png (300 dpi, 必须导出 PNG) + .pdf
F4	Robustness / sensitivity curve — spec curve, HonestDiD, or cluster-comparison	spec_curve loop or honest_did plot (Step 6)	figures/fig4_sensitivity.png (300 dpi, 必须导出 PNG) + .pdf

Output file layout (default)

project/
├── tables/    table1_balance.xlsx/.docx/.tex  table2_main.xlsx/.docx/.tex
│              table3_mechanism.xlsx/.docx/.tex table4_heterogeneity.xlsx/.docx/.tex
│              table5_robustness.xlsx/.docx/.tex
└── figures/   fig1_trend.png(300dpi)+.pdf      fig2_event_study.png(300dpi)+.pdf
               fig3_coefplot.png(300dpi)+.pdf   fig4_sensitivity.png(300dpi)+.pdf

关键输出规则（必须遵守）：

• 图片格式：所有图片必须同时导出 PNG 格式（≥300 dpi） 和 PDF 格式（用于 LaTeX 排版）
• 表格格式：所有回归表格必须同时导出 Excel（.xlsx）、Word（.docx） 和 LaTeX（.tex） 三种格式
• PNG 用于幻灯片、Markdown 文档、邮件等场景；PDF 用于学术论文排版

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 note explaining why; do not silently drop.
• N=1 treated unit (synthetic control) — replace F1/F2 with the SCM trajectory + placebo distribution; T1–T5 still apply.

---

Required Libraries

pip install pandas numpy scipy matplotlib seaborn \
            statsmodels linearmodels pyfixest \
            rdrobust rddensity \
            econml causalml \
            stargazer  # publication-ready regression tables
# Optional but commonly needed:
pip install missingno   # missing-data visualization
pip install pyreadstat  # Stata .dta / SPSS .sav import
pip install arch        # GARCH, unit-root tests, HAC
pip install pingouin    # clean stats tests wrapper
pip install pysynth     # synthetic control (N=1 treated)

---

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

┌──────────────────────────────────────────────────────────────────────┐
│ Step −1 Pre-Analysis Plan (PAP)   statsmodels.stats.power / mde     │
│ Step 0  Sample log + data contract sample_log/asserts/JSON dump      │
│ Step 1  Data cleaning              missing / outliers / dtype / join │
│ Step 2  Variable construction      log / winsorize / std / encode    │
│ Step 2.5 Empirical strategy        equation × ID assumption + pre-reg│
│ Step 3  Descriptive statistics     Table 1 / corr / distribution     │
│ Step 3.5 Identification graphics   event-study/1st-stage/McCrary/love│
│ Step 4  Diagnostic tests           normality / hetero / autocorr / VIF│
│ Step 5  Baseline modeling          OLS / panel / IV / DID / RD / SC  │
│ Step 6  Robustness battery         placebo / subsample / spec curve  │
│ Step 7  Further analysis           mechanism / heterogeneity / mediation│
│ Step 8  Tables & figures           stargazer / coefplot / event study│
└──────────────────────────────────────────────────────────────────────┘

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  Python moves
─────────────────────────── ───── ────────────────────────────────────────────────
Pre-Analysis Plan           −1    statsmodels.stats.power + freeze pap.json
§1. Data                     0    sample_log + 5-check data contract → JSON
§1. Data                     1    pandas read_*/dropna/dtype/merge(validate=)
§1. Data                     2    np.log/np.arcsinh/winsorize/groupby.shift/diff
§1.1 Descriptives (Table 1)  3    df.describe() · table1_with_smd · seaborn
§2. Empirical Strategy       2.5  equation × ID assumption table + pre-reg
§3. Identification graphics  3.5  pf.iplot · binscatter · rdplot · mccrary · love
§3.5 Diagnostics             4    statsmodels.diagnostic + scipy.stats
§4. Main Results (Table 2)   5    pf.feols/IV2SLS/CausalForest · pf.etable / Stargazer
§5. Heterogeneity (Table 3)  7    pf.feols(... + i(.):X) · marginaleffects-py
§6. Mechanisms / Channels    7    Baron-Kenny · econml.dml · outcome ladder
§7. Robustness gauntlet      6    placebo · oster · honestdid · spec_curve · 2-way
§8. Replication package      8    Stargazer.render_latex · pf.etable("docx") · result.json

Below is the canonical code 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. Column names and parameter values are illustrative — substitute the real ones from the user's DataFrame. Only library names and call shapes are normative.

When a step has many variants (e.g. staggered DID has five different estimators; heteroskedasticity has four 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 / .docx / .xlsx exports (the agent should produce all three so co-authors can edit in Word, the build system can use LaTeX, and editors can edit raw numbers in Excel). 所有图片必须同时保存 PNG（≥300 dpi）和 PDF 两种格式。

§	Artifact	Python primitive	Filenames
§1	Figure 1: raw trends / treatment rollout	df.groupby([time,treat])[y].mean().unstack().plot() · seaborn.heatmap for staggered rollout	figures/fig1_trend.png(300dpi)+.pdf
§1	Table 1: summary stats (full / treated / control + Δ + SMD)	table1(df, by=, cols=) (Step 3.b) → write LaTeX/Word/Excel	tables/table1_balance.xlsx/.docx/.tex
§3	Figure 2: identification graphic (event-study / first-stage / McCrary / RD scatter / SCM trajectory)	pf.iplot(es) · binscatter · rdplot · rddensity · synthdid	figures/fig2_event_study.png(300dpi)+.pdf
§4	Table 2: main results — progressive controls M1→M6	pf.etable([m1...m6]) · Stargazer([m1.fit ... m6.fit])	tables/table2_main.xlsx/.docx/.tex
§4	Table 2-bis: design horse-race (OLS / IV / DID / DML)	pf.etable([ols, iv, did, dml])	tables/table2b_designs.xlsx/.docx/.tex
§4	Figure 3: coefficient plot across specs	pf.coefplot([m1...m6], coefs=["training"])	figures/fig3_coefplot.png(300dpi)+.pdf
§5	Table 3: heterogeneity by subgroup	pf.etable(g_full, g_male, g_fem, g_q1...q4)	tables/table3_heterogeneity.xlsx/.docx/.tex
§5	Figure 4: dose-response / CATE	econml.CausalForestDML(...).effect() + matplotlib hist	figures/fig4_cate.png(300dpi)+.pdf
§6	Table 4: mechanism / outcome ladder	loop pf.feols over outcomes → pf.etable	tables/table4_mechanism.xlsx/.docx/.tex
§7	Table A1: robustness master (one column per check)	pf.etable([base, no99, balpan, dropearly, wfe, cl2way, logy, ihsy, psm, ebal])	tables/tableA1_robustness.xlsx/.docx/.tex
§7	Figure 5: spec curve	hand-rolled itertools.product + matplotlib errorbar	figures/fig5_spec_curve.png(300dpi)+.pdf
§7	Figure 6: sensitivity (HonestDiD / Oster / E-value)	HonestDiD · oster_bound · evalue	figures/fig6_sensitivity.png(300dpi)+.pdf
§8	Replication bundle: all tables in one document	pf.etable([...] + extra=[...], type="docx") · pylatex / texdoc-style multi-panel	replication/paper_tables.xlsx/.docx/.tex

Every Python estimator above (pf.feols / IV2SLS / att_gt via R-callout / CausalForestDML) returns a result object that can be passed straight into pf.etable(...) / pf.coefplot(...) / Stargazer(...). Don't hand-roll LaTeX from df.to_latex(), and don't render Word via python-docx directly — pf.etable / Stargazer apply book-tab borders, AER stars, and the right SE label automatically. For deeper export recipes (LaTeX / Word / Markdown variants, multi-panel .docx, full gtsummary-style flow), see references/08-tables-plots.md.

---

Export cookbook — LaTeX / Word / Excel in one block

关键规则（必须遵守）：每个表格必须同时导出三种格式——Excel(.xlsx)、Word(.docx)、LaTeX(.tex)。每个图片必须同时保存PNG(≥300dpi)和PDF两种格式。

Three tiers, picked by scope:

Tier	Use when	API	Hot kwargs
1. Single multi-column table	Exporting one Table 2 / Table 3 / Table A1 with progressive columns	pf.etable([m1,...,mN], type="tex"/"docx"/"xlsx", file="...", headers=[...], digits=3, signif_code=[0.1,0.05,0.01]) — or Stargazer([m1.fit,...,mN.fit]).render_latex() for statsmodels-only	keep=, drop=, coef_map=, headers=, digits=, signif_code=, fixef_rm=, notes=
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	Stack via repeat pf.etable(..., extralines=...) calls, or use gtsummary-style chained tables; for true single-file multi-panel, write to a .tex then concat	first panel: write; subsequent: append; surround with LaTeX \section{} headers
3. Full session bundle (the Stata collect / R gt equivalent)	Replication appendix that mixes summary stats + balance + multiple regression tables + headings + prose in one file	Compose programmatically with pylatex / python-docx / quarto — render once. Or use statsmodels.iolib.summary2.summary_col for a quick concat of tables.	journal-style template + per-section heading + footnote macros

Journal styling — pick the right signif_code and SE label. AEA convention is [0.1, 0.05, 0.01] and SE label "Cluster-robust standard errors in parentheses". Define a wrapper once at the top of master.py:

# top of master.py — journal house-style wrapper
AER_SIGNIF = [0.1, 0.05, 0.01]
AER_NOTES  = ("Cluster-robust standard errors in parentheses. "
              "* p<0.10, ** p<0.05, *** p<0.01.")

def aer_table(models, *, file, headers=None, coef_map=None):
    # 同时导出三种格式：.xlsx（用于编辑）、.docx（用于Word）、.tex（用于LaTeX）
    base, ext = os.path.splitext(file)
    for ext, type_ in [(".xlsx", "xlsx"), (".docx", "docx"), (".tex", "tex")]:
        pf.etable(models, type=type_, file=base + ext,
                  headers=headers, coef_map=coef_map,
                  digits=3, signif_code=AER_SIGNIF, notes=AER_NOTES)

For the multi-panel .docx / .xlsx and Markdown / Quarto cookbook (single-file paper-tables bundle), 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 pap.json so a referee can verify the design was powered before, not after, the data were seen.

import json
from statsmodels.stats.power import TTestIndPower, NormalIndPower
from statsmodels.stats.proportion import samplesize_proportions_2indep_onetail

# Two-sample MDE for a continuous outcome (Cohen's d framing)
analysis = TTestIndPower()
n_required = analysis.solve_power(effect_size=0.20, power=0.80, alpha=0.05, ratio=1.0)
print(f"n per arm for d=0.20, 80% power: {n_required:.0f}")

# Solve for MDE given fixed n
mde = analysis.solve_power(nobs1=2000, power=0.80, alpha=0.05, ratio=1.0)
print(f"MDE (Cohen's d) at n=2000 per arm: {mde:.3f}")

# Cluster-randomized RCT — design effect = 1 + (m-1)·ICC
m, icc = 50, 0.05
deff = 1 + (m - 1) * icc
n_eff_required = analysis.solve_power(effect_size=0.20, power=0.80, alpha=0.05) * deff
print(f"n per arm under ICC={icc}, cluster size={m}: {n_eff_required:.0f}")

# DID: use Frison-Pocock / Bloom (1995) — see references/05-modeling.md §5.4
# RD: power via Monte Carlo — see references/05-modeling.md §5.5

# Persist the protocol — the referee will ask whether the design was powered ex ante
pap = {
    "population":        "manufacturing workers, 2010–2020",
    "treatment":         "training (binary, staggered adoption)",
    "outcome":           "log_wage",
    "estimand":          "ATT",
    "design":            "staggered DID, Callaway–Sant'Anna",
    "alpha":             0.05,
    "power_target":      0.80,
    "mde_d":             0.20,
    "n_planned":         12000,
    "frozen_at":         "2026-01-15T09:00:00Z",
    "git_sha":           "<paste>",
}
with open("artifacts/pap.json", "w") as f:
    json.dump(pap, f, indent=2)

Commit artifacts/pap.json in the repo before 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. The data contract is the cure for "mysterious sample-size shrinkage" bugs in the response letter.

0.1 Sample-construction log (footnote 4)

import pandas as pd, json

sample_log = []

df_raw = pd.read_csv("raw.csv")
sample_log.append(("0. raw",                     len(df_raw)))

df1 = df_raw.dropna(subset=["wage"])
sample_log.append(("1. drop missing wage",        len(df1)))

df2 = df1[df1["age"].between(18, 65)]
sample_log.append(("2. drop age outside 18-65",   len(df2)))

df3 = df2[df2["industry"].isin({"manuf","construction","transport"})]
sample_log.append(("3. keep target industries",   len(df3)))

df = df3
for label, n in sample_log:
    print(f"  {label:<30s}  N = {n:>10,d}")

with open("artifacts/sample_construction.json", "w") as f:
    json.dump(sample_log, f, indent=2)

Paste the printed lines verbatim as footnote 4 of the paper.