Stata

Stata Skill

You have access to comprehensive Stata reference files. Do not load all files. Read only the 1-3 files relevant to the user's current task using the routing table below.

---

Critical Gotchas

These are Stata-specific pitfalls that lead to silent bugs. Internalize these before writing any code.

Missing Values Sort to +Infinity

Stata's . (and .a-.z) are greater than all numbers.

* WRONG — includes observations where income is missing!
gen high_income = (income > 50000)

* RIGHT
gen high_income = (income > 50000) if !missing(income)

* WRONG — missing ages appear in this list
list if age > 60

* RIGHT
list if age > 60 & !missing(age)

= vs ==

= is assignment; == is comparison. Mixing them up is a syntax error or silent bug.

* WRONG — syntax error
gen employed = 1 if status = 1

* RIGHT
gen employed = 1 if status == 1

Local Macro Syntax

Locals use ` name' ` (backtick + single-quote). Globals use $name or ${name}`. Forgetting the closing quote is the #1 macro bug.

local controls "age education income"
regress wage `controls'        // correct
regress wage `controls         // WRONG — missing closing quote
regress wage 'controls'        // WRONG — wrong quote characters

by Requires Prior Sort (Use bysort)

* WRONG — error if data not sorted by id
by id: gen first = (_n == 1)

* RIGHT — bysort sorts automatically
bysort id: gen first = (_n == 1)

* Also RIGHT — explicit sort
sort id
by id: gen first = (_n == 1)

Factor Variable Notation (i. and c.)

Use i. for categorical, c. for continuous. Omitting i. treats categories as continuous.

* WRONG — treats race as continuous (e.g., race=3 has 3x effect of race=1)
regress wage race education

* RIGHT — creates dummies automatically
regress wage i.race education

* Interactions
regress wage i.race##c.education    // full interaction
regress wage i.race#c.education     // interaction only (no main effects)

generate vs replace

generate creates new variables; replace modifies existing ones. Using generate on an existing variable name is an error.

gen x = 1
gen x = 2          // ERROR: x already defined
replace x = 2      // correct

String Comparison Is Case-Sensitive

* May miss "Male", "MALE", etc.
keep if gender == "male"

* Safer
keep if lower(gender) == "male"

merge Always Check _merge

Never skip tab _merge — it costs nothing and is the only diagnostic you get when assert fails.

merge 1:1 id using other.dta
tab _merge                      // ALWAYS tab before assert
assert _merge == 3              // fails silently without tab output
drop _merge

preserve / restore + tempfile for Collapse-Merge-Back

The standard pattern for computing group stats and merging them onto the original data:

tempfile stats
preserve
collapse (mean) avg_x=x, by(group)
save `stats'
restore
merge m:1 group using `stats'
tab _merge
assert _merge == 3
drop _merge

For simple group means, bysort group: egen avg_x = mean(x) avoids the round-trip entirely.

Weights Are Not Interchangeable

• fweight — frequency weights (replication)
• aweight — analytic/regression weights (inverse variance)
• pweight — probability/sampling weights (survey data, implies robust SE)
• iweight — importance weights (rarely used)

capture Swallows Errors

capture some_command
if _rc != 0 {
    di as error "Failed with code: " _rc
    exit _rc
}

Line Continuation Uses ///

regress y x1 x2 x3 ///
    x4 x5 x6, ///
    vce(robust)

Stored Results: r() vs e() vs s()

• r() — r-class commands (summarize, tabulate, etc.)
• e() — e-class commands (estimation: regress, logit, etc.)
• s() — s-class commands (parsing)

A new estimation command overwrites previous e() results. Store them first:

regress y x1 x2
estimates store model1

---

Running Stata from the Command Line

Claude can execute Stata code by running .do files in batch mode from the terminal. This is how to run Stata non-interactively.

Finding the Stata Binary

Stata on macOS is a .app bundle. The actual binary is inside it. Common locations:

# Stata 18 / StataNow (most common)
/Applications/Stata/StataMP.app/Contents/MacOS/stata-mp
/Applications/StataNow/StataMP.app/Contents/MacOS/stata-mp

# Other editions (SE, BE)
/Applications/Stata/StataSE.app/Contents/MacOS/stata-se
/Applications/Stata/StataBE.app/Contents/MacOS/stata-be

If Stata isn't on $PATH, find it with: mdfind -name "stata-mp" | grep MacOS

Batch Mode (-b)

# Run a .do file in batch mode — output goes to <filename>.log
/Applications/Stata/StataMP.app/Contents/MacOS/stata-mp -b do analysis.do

# If stata-mp is on PATH (e.g., via symlink or alias):
stata-mp -b do analysis.do

• -b = batch mode (non-interactive, no GUI)
• Output (everything Stata would display) is written to analysis.log in the working directory
• Exit code is 0 on success, non-zero on error
• The log file contains all output, including error messages — check it after execution

Running Inline Stata Code

To run a quick Stata snippet without creating a .do file:

# Write a temp .do file and run it
cat > /tmp/stata_run.do << 'EOF'
sysuse auto, clear
summarize price mpg
EOF
stata-mp -b do /tmp/stata_run.do
cat /tmp/stata_run.log

Checking Results

# Check if it succeeded
stata-mp -b do tests/run_tests.do && echo "SUCCESS" || echo "FAILED"

# Search the log for pass/fail
grep -E "PASS|FAIL|error|r\([0-9]+\)" run_tests.log

Tips

• `clear all` at the top of batch scripts — batch mode starts with a fresh Stata session, but clear all ensures no stale state from prior runs in the same session.
• `set more off` — prevents Stata from pausing for --more-- prompts (fatal in batch mode).
• Log files overwrite silently — analysis.do always writes to analysis.log in the current directory. If you run multiple .do files, check the right log.
• Working directory — Stata's working directory is wherever you run the command from, not where the .do file lives. Use cd in the .do file or absolute paths if needed.

---

Routing Table

Read only the files relevant to the user's task. Paths are relative to this SKILL.md file.

Data Operations

Statistics & Econometrics

Causal Inference

Advanced Methods

Graphics

Programming

Output & Workflow

Community Packages

---

Common Patterns

Regression Table Workflow

* Estimate models
eststo clear
eststo: regress y x1 x2, vce(robust)
eststo: regress y x1 x2 x3, vce(robust)
eststo: regress y x1 x2 x3 x4, vce(cluster id)

* Export table
esttab using "results.tex", replace ///
    se star(* 0.10 ** 0.05 *** 0.01) ///
    label booktabs ///
    title("Main Results") ///
    mtitles("(1)" "(2)" "(3)")

Panel Data Setup

xtset panelid timevar          // declare panel structure
xtdescribe                      // check balance
xtsum outcome                   // within/between variation

* Fixed effects
xtreg y x1 x2, fe vce(cluster panelid)
* Or with reghdfe (preferred for multiple FE)
reghdfe y x1 x2, absorb(panelid timevar) vce(cluster panelid)

Difference-in-Differences

* Classic 2x2 DiD
gen post = (year >= treatment_year)
gen treat_post = treated * post
regress y treated post treat_post, vce(cluster id)

* Event study (uniform timing — must interact with treatment group)
reghdfe y ib(-1).rel_time#1.treated, absorb(id year) vce(cluster id)
testparm *.rel_time#1.treated   // pre-trend test

* Modern staggered DiD (Callaway & Sant'Anna)
csdid y x1 x2, ivar(id) time(year) gvar(first_treat) agg(event)
csdid_plot

Graph Export

* Publication-quality scatter with fit line
twoway (scatter y x, mcolor(navy%50) msize(small)) ///
       (lfit y x, lcolor(cranberry) lwidth(medthick)), ///
    title("Title Here") ///
    xtitle("X Label") ytitle("Y Label") ///
    legend(off) scheme(s2color)
graph export "figure1.pdf", replace as(pdf)
graph export "figure1.png", replace as(png) width(2400)

Data Cleaning Pipeline

* Load and inspect
import delimited "raw_data.csv", clear varnames(1)
describe
codebook, compact

* Clean
rename *, lower                 // lowercase all varnames
destring income, replace force  // convert string to numeric
replace income = . if income < 0

* Label
label variable income "Annual household income (USD)"
label define yesno 0 "No" 1 "Yes"
label values employed yesno

* Save
compress
save "clean_data.dta", replace

Multiple Imputation

mi set mlong
mi register imputed income education
mi impute chained (regress) income (ologit) education = age i.gender, add(20) rseed(12345)
mi estimate: regress wage income education age i.gender

---

Help Improve This Skill

If you produce Stata code with a significant error — wrong syntax, incorrect command usage, or a gotcha you failed to catch — and the issue seems to stem from a gap in these reference files rather than a one-off mistake, consider suggesting to the user that they file an issue on the skill repository. This helps future users.

When to raise this: Only after you've already corrected the error and the user has working code. Frame it as optional: "I made an error with [X] that I think comes from a gap in the Stata skill documentation. If you'd like, I can help you file an issue or a PR so it gets fixed for everyone."

When NOT to raise this: If the user is on Claude Haiku, the error is more likely a model capability issue than a documentation gap. In that case, suggest they try Sonnet or Opus for complex Stata work instead of filing an issue.

If the user agrees, read references/filing-issues.md for instructions on writing a good issue report.

asdoc: Automatic Document Creation in Stata

Overview

asdoc creates formatted tables in Word, LaTeX, HTML, and text by adding a simple prefix to Stata commands. No need to store estimation results -- just prefix your command with asdoc.

Key features:

• Simple prefix syntax for most Stata commands
• Automatic formatting for Word documents (.doc/.docx)
• Supports: summarize, tabstat, correlate, pwcorr, regress, xtreg, logit, probit, tabulate, ttest, anova, and more
• Append multiple tables to one document
• Output to Word (default), LaTeX, HTML, or plain text

---

Installation

ssc install asdoc
help asdoc
which asdoc

Important: Output files are created in your current working directory. Check with pwd.

---

Quick Start

sysuse auto, clear

* Just prefix your command with asdoc
asdoc summarize price mpg weight, replace

* Append additional tables to the same document
asdoc correlate price mpg weight, append
asdoc regress price mpg weight foreign, append

* Default output file is "MyFile.doc" in current directory

---

File Management

* Specify custom filename
asdoc summarize price mpg weight, save(mystats.doc) replace

* Full path
asdoc summarize price mpg, save("C:/Tables/summary.doc") replace

* Pattern: replace for first table, append for the rest
asdoc summarize price mpg, replace
asdoc correlate price mpg, append
asdoc regress price mpg weight, append

Output format is determined by file extension:

asdoc summarize price, save(file.doc) replace   // Word (default)
asdoc summarize price, save(file.tex) replace   // LaTeX
asdoc summarize price, save(file.html) replace  // HTML
asdoc summarize price, save(file.txt) replace   // Plain text

---

Summary Statistics

sysuse auto, clear

* Basic summary
asdoc summarize price mpg weight, replace

* Detailed summary (adds percentiles, skewness, kurtosis)
asdoc summarize price mpg weight, detail replace

* With title and decimal control
asdoc summarize price mpg weight, ///
    title(Table 1: Summary Statistics) dec(2) replace

* tabstat for specific statistics
asdoc tabstat price mpg weight, ///
    statistics(n mean sd min max) replace

* tabstat by groups
asdoc tabstat price mpg weight, ///
    by(foreign) statistics(mean sd) replace

Summary by Groups (Manual Panels)

asdoc, row(Domestic Cars) replace
asdoc summarize price mpg weight if foreign == 0, append

asdoc, row(Foreign Cars) append
asdoc summarize price mpg weight if foreign == 1, append

---

Regression Tables

sysuse auto, clear

* Basic regression
asdoc regress price mpg weight, replace

* With robust standard errors
asdoc regress price mpg weight foreign, robust replace

* Clustered standard errors
asdoc regress price mpg weight, vce(cluster rep78) replace

* Control which statistics to display
asdoc regress price mpg weight foreign, ///
    stat(coef se tstat pvalue) dec(3) replace

Progressive Model Building

asdoc, text(Model 1: Baseline) replace
asdoc regress price mpg, append

asdoc, text(Model 2: Adding Weight) append
asdoc regress price mpg weight, append

asdoc, text(Model 3: Full Model) append
asdoc regress price mpg weight foreign, append

Gotcha: asdoc shows models sequentially (not side-by-side in columns). For side-by-side model comparison, use estout/esttab instead.

---

Panel Data

webuse nlswork, clear
xtset idcode year

* Fixed effects
asdoc xtreg ln_wage age ttl_exp tenure, fe replace

* Random effects
asdoc xtreg ln_wage age ttl_exp tenure, re append

* With robust SE
asdoc xtreg ln_wage age ttl_exp tenure, fe robust append

---

Correlation Matrices

sysuse auto, clear

* Simple correlation
asdoc correlate price mpg weight, replace

* Pairwise with significance stars
asdoc pwcorr price mpg weight foreign, star(all) replace

* With p-values and observation counts
asdoc pwcorr price mpg weight, star(0.05) obs replace

* Custom title and formatting
asdoc pwcorr price mpg weight foreign, ///
    star(0.05) title(Correlation Matrix) dec(3) replace

---

Other Statistical Commands

* Cross-tabulation with chi-square
asdoc tabulate foreign rep78, chi2 replace

* Two-sample t-test
asdoc ttest price, by(foreign) replace

* ANOVA
asdoc anova price foreign rep78, replace

* Logit with odds ratios
asdoc logit foreign mpg weight price, or replace

* Probit
asdoc probit foreign mpg weight price, replace

---

Document Structure: Text, Rows, and Titles

* Section headers with text()
asdoc, text(SECTION 1: DESCRIPTIVE ANALYSIS) replace
asdoc summarize price mpg weight, append

asdoc, text() append   // Blank line for spacing

asdoc, text(SECTION 2: REGRESSION ANALYSIS) append
asdoc regress price mpg weight, append

* Row labels for subsections
asdoc, row(Panel A: Full Sample) replace
asdoc summarize price mpg, append

asdoc, row(Panel B: Domestic Only) append
asdoc summarize price mpg if foreign == 0, append

* Table titles
asdoc summarize price mpg, ///
    title(Table 1: Descriptive Statistics) replace

---

Formatting Options

Helper Commands

asdoc, text(your text here) append   // Add text line
asdoc, row(label text) append        // Add row label
asdoc, text() append                 // Blank line
asdoc, version                       // Check version

---

Loops and Automation

* Loop through variables
asdoc, text(Variable Summaries) replace
foreach var of varlist price mpg weight {
    asdoc, row(`var' Statistics) append
    asdoc summarize `var', append
}

* Loop through groups
levelsof foreign, local(levels)
foreach lev of local levels {
    asdoc, row(Foreign = `lev') append
    asdoc summarize price mpg if foreign == `lev', append
}

* Reusable report program
capture program drop make_report
program make_report
    syntax varlist using/
    asdoc, text(AUTOMATED REPORT) save(`using') replace
    asdoc, text(Generated: `c(current_date)') append
    asdoc summarize `varlist', append
    asdoc correlate `varlist', append
end

make_report price mpg weight using "report.doc"

Dated Filenames

local date: display %tdCYND daily("$S_DATE", "DMY")
asdoc summarize price, save(summary_`date'.doc) replace

---

When to Use asdoc vs. Alternatives

Key asdoc limitations:

• Models displayed sequentially, not in columns
• Limited layout customization compared to estout
• Less LaTeX control than estout

Key asdoc advantages:

• No need to store estimates (eststo not required)
• Works with many command types beyond regression
• Minimal learning curve

---

Common Issues

---

Quick Reference

* Summary statistics
asdoc summarize varlist, replace
asdoc tabstat varlist, statistics(mean sd min max) replace

* Correlation
asdoc correlate varlist, replace
asdoc pwcorr varlist, star(all) replace

* Regression
asdoc regress depvar indepvars, replace
asdoc regress depvar indepvars, robust stat(coef se tstat) dec(3) replace

* Panel data
asdoc xtreg depvar indepvars, fe replace
asdoc xtreg depvar indepvars, re append

* Binary outcome
asdoc logit depvar indepvars, or replace
asdoc probit depvar indepvars, replace

* Cross-tabulation
asdoc tabulate var1 var2, chi2 replace

* T-test
asdoc ttest var, by(group) replace

* Text and structure
asdoc, text(Section Header) replace
asdoc, row(Subsection Label) append
asdoc, text() append   // blank line

BINSREG: Binscatter Estimation and Inference

Overview

The binsreg package (Cattaneo, Crump, Farrell, Feng) implements modern binscatter methods with rigorous statistical foundations: data-driven bin selection, pointwise and uniform inference, hypothesis testing, and support for OLS, logit, probit, and quantile regression.

Key advantage over traditional `binscatter` (Stepner): Proper covariate adjustment (no residualization), formal inference (CI/CB), optimal bin selection, and hypothesis testing.

Package Commands

---

Installation

net install binsreg, from(https://raw.githubusercontent.com/nppackages/binsreg/master/stata) replace
which binsreg
help binsreg

Also available in R (install.packages('binsreg')) and Python (pip install binsreg).

---

Core Syntax and Options

binsreg (Main Command)

binsreg depvar xvar [controls] [if] [in] [weight] [, options]

Binning Options:

• nbins(#) - Number of bins (default: data-driven)
• binspos(es|qs) - Evenly-spaced or quantile-spaced bins
• binsmethod(dpi|rot) - Bin selection method (direct plug-in or rule-of-thumb)

Polynomial Options (all take `(p v)` = polynomial degree, derivative order):

• dots(p v) - Dots (use dots(0,0) for canonical binscatter)
• line(p v) - Smooth line connecting bins
• ci(p v) - Pointwise confidence intervals
• cb(p v) - Uniform confidence bands
• polyreg(#) - Global polynomial overlay of specified degree
• deriv(#) - Plot derivative of specified order

Inference Options:

• level(#) - Confidence level (default: 95)
• vce(robust) - Heteroskedasticity-robust
• vce(cluster clustervar) - Cluster-robust
• vce(bootstrap, reps(#)) - Bootstrap (recommended for quantile regression)

Covariate Options:

• at(mean|median|#) - Evaluation point for control variables
• absorb(varlist) - Fixed effects to absorb

Grouping Options:

• by(groupvar) - Separate binscatter by group
• bycolors(colorlist) - Colors per group
• bysymbols(symbollist) - Symbols per group

Output Options:

• saveplot(filename) - Save graph (.png, .pdf, .gph)
• savedata(filename) - Save bin data for custom plotting
• name(graphname, replace) - Named graph for graph combine

Standard Stata graph options (title(), xtitle(), ytitle(), legend(), scheme(), etc.) all work.

---

Covariate Adjustment (Critical Difference)

WRONG -- traditional residualization approach (biased):

* DON'T DO THIS
reg wage controls
predict wage_resid, resid
reg education controls
predict educ_resid, resid
binscatter wage_resid educ_resid

CORRECT -- binsreg handles adjustment internally:

* DO THIS -- estimates E[wage | education=x, controls=w0]
binsreg wage education controls

---

Binary Outcomes: binslogit / binsprobit

binslogit depvar xvar [controls] [, nolink options]
binsprobit depvar xvar [controls] [, nolink options]

• nolink - Plot fitted probabilities instead of logit/probit index
• All other options same as binsreg

* Employment probability by age
binslogit employed age education, nolink ci(3,3)
binsprobit employed age education, nolink cb(3,3)

---

Quantile Regression: binsqreg

binsqreg depvar xvar [controls] [, quantile(#) options]

• quantile(#) - Quantile to estimate (default: 0.5)
• Bootstrap recommended for inference: vce(bootstrap, reps(100))

binsqreg wage education, quantile(0.5) ci(3,3)
binsqreg wage education experience, quantile(0.75) ci(3,3) vce(bootstrap, reps(100))

---

Hypothesis Testing: binstest

binstest depvar xvar [controls] [, test_options]

Test Options:

• testmodelpoly(#) - Test polynomial specification (1=linear, 2=quadratic, ...)
• testshaper(#) - Shape restriction: 0 for >= (increasing), 1 for <= (decreasing)
• testshapel(#) - Lower bound on derivative
• deriv(#) - Derivative order to test
• lp(#|inf) - Test metric (inf=supremum default, 2=L2, 1=L1)
• estmethod(reg|logit|probit|qreg #) - Estimation method

* Test linearity
binstest wage education, testmodelpoly(1)

* Test quadratic
binstest wage education, testmodelpoly(2)

* Test monotonicity (increasing)
binstest wage education, testshaper(0) deriv(1)

* Test concavity (decreasing returns)
binstest wage education, testshaper(1) deriv(2)

* Test with logit model
binstest employed age, estmethod(logit) testshaper(0) deriv(1)

* Test at specific quantile
binstest wage education, estmethod(qreg 0.5) testmodelpoly(1)

---

Pairwise Comparison: binspwc

binspwc depvar xvar [controls], by(groupvar) [options]

• by(groupvar) - Required grouping variable
• estmethod(method) - Estimation method
• pwc(#) - Comparison type

binspwc wage education, by(gender)
binspwc wage education experience, by(union_status) vce(robust)
binspwc wage education, by(gender) estmethod(qreg 0.4)

---

Bin Selection: binsregselect

binsregselect depvar xvar [controls] [, options]

• binsmethod(dpi|rot) - Selection method
• binspos(es|qs) - Bin positioning
• pselect(# #) - Range of polynomial degrees to search
• sselect(# #) - Range of smoothness degrees

binsregselect wage education
binsregselect wage education, binspos(es) pselect(1/4)

---

Practical Examples

Wage-Education Relationship

sysuse nlsw88, clear
rename grade education
rename ttl_exp experience

* Basic binscatter
binsreg wage education

* With controls, CI, and polynomial overlay
binsreg wage education experience age, ///
    dots(0,0) line(3,3) ci(3,3) polyreg(2) ///
    title("Wage-Education Relationship") ///
    xtitle("Years of Education") ytitle("Hourly Wage ($)")

* Test linearity
binstest wage education experience age, testmodelpoly(1)

* Test monotonicity
binstest wage education experience age, testshaper(0) deriv(1)

Group Comparison

* Visual comparison by union status
binsreg wage education experience, ///
    by(union) dots(0,0) line(3,3) ci(3,3) ///
    legend(order(1 "Non-Union" 2 "Union"))

* Formal pairwise test
binspwc wage education experience, by(union) vce(robust)

Binary Outcome

* Logit binscatter (probability scale)
binslogit employed age education, nolink ci(3,3) ///
    title("Employment Probability by Age") ytitle("Pr(Employed)")

* Test monotonicity
binstest employed age education, estmethod(logit) testshaper(1) deriv(1)

Quantile Comparison

binsqreg wage education experience, quantile(0.25) name(q25, replace)
binsqreg wage education experience, quantile(0.75) name(q75, replace)
graph combine q25 q75, title("Wage Distribution by Education")

* Test linearity at multiple quantiles
foreach q in 0.25 0.5 0.75 {
    display "Testing linearity at quantile `q'"
    binstest wage education experience, estmethod(qreg `q') testmodelpoly(1)
}

Specification Testing Workflow

* Step 1: Visualize with linear overlay
binsreg wage education experience age, ///
    line(3,3) cb(3,3) polyreg(1)

* Step 2: Test linearity
binstest wage education experience age, testmodelpoly(1)
local pval_linear = r(p_val)

* Step 3: If rejected, test quadratic
if `pval_linear' < 0.05 {
    binstest wage education experience age, testmodelpoly(2)
}

* Step 4: Test economic shape restrictions
binstest wage education experience age, testshaper(0) deriv(1)  // Increasing
binstest wage education experience age, testshaper(1) deriv(2)  // Concave

Derivative Plot

* Marginal return to education
binsreg wage education, deriv(1) line(2,2) cb(2,2) ///
    title("Marginal Return to Education")

Saved Data for Custom Plotting

binsreg wage education, savedata(bindata) replace
use bindata, clear
* Variables: dots_x, dots_fit (bin-mean dots), poly_x, poly_fit (polynomial),
*            CI_l, CI_r (pointwise CI), CB_l, CB_r (uniform CB)
* With by(): each group gets suffix: dots_x_1, dots_x_2, etc.
twoway (scatter dots_fit dots_x) ///
       (rcap CI_l CI_r dots_x), ///
       title("Custom Binscatter")

Publication-Ready Plot

binsreg wage education experience age, ///
    nbins(20) dots(0,0) line(3,3) ///
    ci(3,3) cb(3,3) polyreg(2) ///
    title("Returns to Education", size(medium)) ///
    xtitle("Years of Education", size(small)) ///
    ytitle("Hourly Wage ($)", size(small)) ///
    note("95% CI and CB shown. Quadratic overlay.", size(vsmall)) ///
    graphregion(color(white)) scheme(s1color) ///
    saveplot("figure1_wage_education.pdf") replace

---

Migration from Traditional binscatter

* OLD (binscatter):
binscatter wage education, nquantiles(20) ///
    controls(experience age) absorb(state)

* NEW (binsreg) -- with proper covariate adjustment and inference:
binsreg wage education experience age, ///
    nbins(20) binspos(qs) absorb(state) ci(3,3) cb(3,3)

* Or with automatic bin selection:
binsreg wage education experience age, absorb(state) ci(3,3)

---

References

Cattaneo, M. D., Crump, R. K., Farrell, M. H., and Feng, Y. (2024). "On Binscatter." American Economic Review, 114(5): 1488-1514.

Cattaneo, M. D., Crump, R. K., Farrell, M. H., and Feng, Y. (2025). "Binscatter Regressions." The Stata Journal (forthcoming). arXiv:2407.15276.

• Website: https://nppackages.github.io/binsreg/
• GitHub: https://github.com/nppackages/binsreg

Coefplot: Visualizing Regression Coefficients

coefplot creates publication-quality coefficient plots with confidence intervals from stored estimation results.

Installation

ssc install coefplot, replace

// Recommended companions
ssc install estout
ssc install grstyle
ssc install palettes
ssc install colrspace

// Verify
sysuse auto, clear
regress price mpg weight foreign
coefplot

Basic Coefficient Plots

// Basic plot from last estimation
regress price mpg weight foreign length
coefplot

// Specify confidence levels
coefplot, levels(90)
coefplot, levels(90 95 99)

// From stored estimates
regress price mpg weight foreign
estimates store model1
coefplot model1, keep(mpg weight)

// With labels and title
coefplot, ///
    title("Determinants of Automobile Price") ///
    xtitle("Coefficient Estimate") ///
    xlabel(, format(%9.0fc))

Plotting Multiple Models

quietly regress price mpg
estimates store m1
quietly regress price mpg weight
estimates store m2
quietly regress price mpg weight foreign
estimates store m3

// Side-by-side comparison
coefplot m1 m2 m3, drop(_cons)

// With reference line and legend
coefplot (m1) (m2) (m3), drop(_cons) ///
    xline(0, lcolor(red) lpattern(dash)) ///
    legend(order(1 "Model 1" 2 "Model 2" 3 "Model 3"))

// Stacked panels
coefplot (m1, label(Model 1)) ///
         (m2, label(Model 2)) ///
         (m3, label(Model 3)), ///
    drop(_cons) xline(0) byopts(yrescale)

Customization

Markers and Confidence Intervals

// Custom markers
coefplot, ///
    msymbol(square) msize(large) mcolor(navy) mfcolor(navy%50)

// Different markers per model
coefplot m1 m2, drop(_cons) ///
    msymbol(O D) mcolor(navy maroon)

// CI appearance
coefplot, ciopts(lwidth(2) lcolor(navy)) levels(95)

// Capped intervals
coefplot, ciopts(recast(rcap) lwidth(medium))

// Multiple CI levels with different styles
coefplot, levels(90 95) ///
    ciopts1(lwidth(thin) lcolor(gs10)) ///
    ciopts2(lwidth(thick) lcolor(navy))

Colors

// Named colors
coefplot m1 m2 m3, drop(_cons) mcolor(navy maroon forest_green)

// Semi-transparent
coefplot m1 m2 m3, drop(_cons) mcolor(navy%60 maroon%60 forest_green%60)

// Grayscale for print
coefplot m1 m2 m3, drop(_cons) ///
    mcolor(gs2 gs8 gs14) ciopts(lcolor(gs2 gs8 gs14))

Labels

// Use variable labels
label variable mpg "Fuel Efficiency (MPG)"
coefplot, label

// Rename coefficients inline
coefplot, coeflabels(mpg = "Miles per Gallon" ///
                     weight = "Vehicle Weight" ///
                     foreign = "Foreign Indicator")

Horizontal vs Vertical Plots

// Horizontal (default) - better for many variables
coefplot, drop(_cons) xline(0)

// Vertical - better for comparing many models
coefplot m1 m2 m3, drop(_cons) vertical yline(0, lcolor(red) lpattern(dash))

Reordering and Grouping

// Custom order
coefplot, order(foreign mpg weight length) drop(_cons)

// Group headings
coefplot, ///
    headings(mpg = "{bf:Vehicle Characteristics}" ///
             foreign = "{bf:Origin}") ///
    drop(_cons)

// Combined ordering and grouping
coefplot, ///
    headings(mpg = "{bf:Performance}" ///
             foreign = "{bf:Origin & Quality}") ///
    order(mpg weight displacement foreign rep78) ///
    drop(_cons turn)

Omitting and Keeping Variables

// Drop constant (almost always do this)
coefplot, drop(_cons)

// Drop multiple
coefplot, drop(_cons foreign)

// Drop factor variables with wildcards
regress price i.rep78 mpg weight foreign
coefplot, drop(_cons *.rep78)

// Keep only specific variables
coefplot, keep(mpg weight foreign)

// Keep only interactions
coefplot, keep(*#*)

// Drop factor and interaction terms
coefplot, drop(*.foreign *#*)

Interaction Terms

// Basic interaction plot
regress price c.mpg##i.foreign weight
coefplot, drop(_cons weight) xline(0)

// Organized with headings
regress price c.mpg##i.foreign c.weight##i.foreign
coefplot, drop(_cons) ///
    headings(mpg = "{bf:Main Effects}" ///
             1.foreign#c.mpg = "{bf:Interactions}") ///
    xline(0)

// Continuous-by-continuous
regress price c.mpg##c.weight foreign
coefplot, ///
    keep(mpg weight *#*) ///
    coeflabels(mpg#c.weight = "MPG x Weight") ///
    xline(0)

Integration with Margins

// Marginal effects plot
regress price c.mpg##i.foreign weight
margins, dydx(mpg) over(foreign)
coefplot, xline(0) ///
    coeflabels(0.foreign = "Domestic" 1.foreign = "Foreign")

// Average marginal effects from logit
logit foreign price mpg weight
margins, dydx(*)
coefplot, drop(_cons) xline(0) xlabel(, format(%9.3f))

// Marginal effects at representative values
regress price c.mpg##c.weight foreign
margins, dydx(mpg) at(weight=(2000(500)4500))
coefplot, vertical yline(0) ///
    xtitle("Weight (pounds)") ///
    xlabel(1 "2000" 2 "2500" 3 "3000" 4 "3500" 5 "4000" 6 "4500")

// Contrasts
regress price i.rep78 mpg weight
contrast rep78, nowald
estimates store contrasts
coefplot contrasts, xline(0) ///
    coeflabels(2.rep78 = "Fair vs Poor" ///
               3.rep78 = "Average vs Poor" ///
               4.rep78 = "Good vs Poor" ///
               5.rep78 = "Excellent vs Poor")

Publication-Quality Examples

Single Model, Journal Style

set scheme plotplain

regress price mpg weight foreign length, robust

coefplot, ///
    drop(_cons) ///
    xline(0, lcolor(black) lpattern(dash)) ///
    msymbol(D) msize(medium) mcolor(black) ///
    ciopts(lwidth(medthick) lcolor(black)) ///
    levels(95) ///
    title("Figure 1. Determinants of Automobile Prices", ///
          size(medium) position(11) justification(left)) ///
    xtitle("Coefficient Estimate (USD)", size(small)) ///
    xlabel(, format(%9.0fc) labsize(small)) ///
    ylabel(, labsize(small) angle(0)) ///
    graphregion(color(white)) ///
    plotregion(margin(medium) lcolor(black) lwidth(thin)) ///
    note("Note: 95% confidence intervals shown. Robust standard errors.", ///
         size(vsmall) span)

graph export "figure1_coefficients.pdf", replace

Model Comparison with eststo

set scheme plotplain
eststo clear
eststo: quietly regress price mpg, robust
eststo: quietly regress price mpg weight, robust
eststo: quietly regress price mpg weight foreign, robust

coefplot est1 est2 est3, ///
    drop(_cons) ///
    xline(0, lcolor(black) lpattern(dash)) ///
    msymbol(O S D) msize(medium medium medium) ///
    mcolor(gs2 gs6 gs10) ///
    ciopts(lwidth(medium) lcolor(gs2 gs6 gs10)) ///
    legend(order(1 "Model 1" 2 "Model 2" 3 "Model 3") ///
           rows(1) size(small) region(lwidth(none))) ///
    graphregion(color(white))

graph export "figure2_model_comparison.pdf", replace

Vertical with Multiple CI Levels

set scheme plotplain
regress price mpg weight foreign, robust

coefplot, ///
    drop(_cons) vertical ///
    yline(0, lcolor(black) lpattern(dash)) ///
    levels(90 95 99) ///
    msymbol(D) msize(medium) mcolor(black) ///
    ciopts1(lwidth(thin) lcolor(gs12)) ///
    ciopts2(lwidth(medium) lcolor(gs6)) ///
    ciopts3(lwidth(thick) lcolor(black)) ///
    legend(order(2 "90% CI" 4 "95% CI" 6 "99% CI") ///
           rows(1) position(12) size(vsmall))

Advanced Techniques

Using graph combine

// Create named sub-plots with nodraw, then combine
quietly regress price mpg weight foreign
coefplot, drop(_cons) xline(0) title("Full Sample") name(g1, replace)

quietly regress price mpg weight if foreign==0
coefplot, drop(_cons) xline(0) title("Domestic") name(g2, replace)

graph combine g1 g2, rows(1) graphregion(color(white))

Suppress Display During Multi-Plot Construction

coefplot m1, name(g1, replace) nodraw
coefplot m2, name(g2, replace) nodraw
graph combine g1 g2  // display once

Coefficient Plot + esttab Table

eststo clear
eststo: quietly regress price mpg weight
eststo: quietly regress price mpg weight foreign

coefplot est1 est2, drop(_cons) xline(0) bylabel

esttab est1 est2 using "table.tex", ///
    b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) ///
    label booktabs replace

Standardized Coefficients

Approach 1 — Rescale regressors before regression (simplest):

preserve
foreach v of varlist mpg weight length {
    quietly summarize `v'
    replace `v' = `v' / r(sd)
}
regress price mpg weight length foreign, vce(robust)
coefplot, drop(_cons) horizontal xline(0) title("Standardized coefficients")
restore

Approach 2 — `ereturn post` inside eclass program (when you need stored estimates):

program define _post_std, eclass
    args bmat vmat nobs
    ereturn post `bmat' `vmat', obs(`nobs') depname(price)
    ereturn local cmd "regress"
end

// Build b_std (1×k) and V_std (k×k) matrices, then:
_post_std b_std V_std `=_N'
coefplot, horizontal drop(_cons) xline(0)

Note: V_std must be a k×k square matrix. Do NOT use diag(vecdiag(rowvec)) — vecdiag() requires square input. Build with J(k,k,0) and fill the diagonal.

Quick Reference

Common Pitfalls

• Always drop(_cons) -- the constant clutters the plot
• Always add xline(0) (or yline(0) if vertical) for interpretability
• Use label or coeflabels() -- raw variable names are unclear
• Use grayscale (gs2, gs8, etc.) for print publications
• Use keep() to focus on key variables when there are many coefficients
• Wildcard *.varname drops all factor levels; *#* drops all interactions

Resources

• help coefplot
• Ben Jann's coefplot page: https://repec.sowi.unibe.ch/stata/coefplot/
• Related: marginsplot, esttab/estout, grstyle

Data Manipulation Packages in Stata

Overview

Essential user-written packages for data manipulation, especially with large datasets.

Packages covered:

• gtools - Fast replacements for collapse, egen, reshape, etc. (C plugin, 5-20x speedup)
• rangestat - Moving window / rolling statistics over arbitrary intervals
• egenmore - Extended egen functions (strings, outliers, etc.)
• distinct - Count distinct values
• unique - Identify unique value combinations, tag duplicates
• missings - Missing data analysis and management
• carryforward - Forward/backward fill of missing values

---

Installation

ssc install gtools
ssc install ftools       // recommended dependency for gtools
ssc install rangestat
ssc install egenmore
ssc install distinct
ssc install unique
ssc install missings
ssc install carryforward

* Verify
gtools, check
adoupdate, update

---

gtools: Fast Data Manipulation

Drop-in replacements for built-in commands using C plugins. Same syntax, 5-20x faster on large data.

gcollapse: Fast Aggregation

* Identical syntax to collapse
gcollapse (mean) mean_price=price ///
          (median) med_price=price ///
          (sd) sd_price=price ///
          (p25) p25_price=price ///
          (p75) p75_price=price ///
          (count) n=price, by(foreign)

* Weighted
gcollapse (mean) price [aw=weight], by(foreign)

gegen: Fast egen

* Group statistics
gegen mean_price = mean(price), by(foreign)
gegen sd_price = sd(price), by(foreign)
gegen n_obs = count(price), by(foreign rep78)

* Percentiles within groups
gegen p25_price = pctile(price), p(25) by(foreign)
gegen p50_price = pctile(price), p(50) by(foreign)

* Standardize, rank, IQR, MAD
gegen z_price = std(price), by(foreign)
gegen rank_price = rank(price), by(foreign)
gegen iqr_price = iqr(price), by(foreign)
gegen mad_price = mad(price), by(foreign)

* Tag and group
gegen tag = tag(foreign rep78)       // Tag first obs in group
gegen group_id = group(foreign rep78) // Numeric group ID
gegen group_n = count(1), by(foreign rep78)

greshape: Fast Reshape

* Wide to long (same syntax as reshape)
greshape long score, i(id) j(test)

* Long to wide
greshape wide income expenses, i(id) j(year)

Other gtools Commands

* Fast ID check
gisid make              // errors if not unique
gisid id, missok        // allow missing

* Unique levels
glevelsof foreign, local(levels)

* Quantile groups
gquantiles price_quartile = price, nquantiles(4)
gquantiles price_decile = price, nquantiles(10)
gquantiles price_cut = price, cutpoints(5000 10000 15000)

* Top levels by frequency
gtoplevelsof rep78, ntop(5)

Real-World Example: Panel Data Prep

webuse nlswork, clear

gegen mean_wage = mean(ln_wage), by(idcode)
gegen sd_wage = sd(ln_wage), by(idcode)
gegen first_obs = tag(idcode)
gegen n_obs = count(1), by(idcode)

* Within-group standardization
gegen mean_exp = mean(ttl_exp), by(idcode)
gegen sd_exp = sd(ttl_exp), by(idcode)
generate exp_std = (ttl_exp - mean_exp) / sd_exp

keep if n_obs >= 5

gcollapse (mean) avg_wage=ln_wage ///
          (sd) sd_wage=ln_wage ///
          (first) first_year=year ///
          (last) last_year=year ///
          (count) n_years=year, by(idcode)

When to Use gtools

• Use gtools: datasets >100K obs, repeated aggregations, panel data group operations
• Stick with built-in: datasets <10K obs, one-time ops, need egen functions not in gegen

---

rangestat: Moving Window Statistics

Computes statistics on observations within a specified range -- ideal for rolling calculations, event windows, and irregular time intervals.

Syntax

rangestat (statistic) newvar=varname, interval(varname min max) [by(varlist)]

Moving Averages and Rolling Statistics

tsset date

* 7-day and 30-day backward-looking moving averages
rangestat (mean) ma7=value, interval(date -6 0)
rangestat (mean) ma30=value, interval(date -29 0)

* Multiple rolling stats
rangestat (mean) roll_mean=value ///
          (sd) roll_sd=value ///
          (min) roll_min=value ///
          (max) roll_max=value ///
          (count) roll_n=value, ///
          interval(date -9 0)

* Rolling correlation and regression
rangestat (corr) roll_corr=value value2, interval(date -19 0)
rangestat (reg) roll_beta=value value2, interval(date -29 0)

* Forward-looking, centered, or exact lag/lead
rangestat (mean) forward_ma5=value, interval(date 0 5)
rangestat (mean) centered_ma7=value, interval(date -3 3)
rangestat (mean) lag10=value, interval(date -10 -10)

Event Study Windows

* Pre- and post-event averages by group
rangestat (mean) pre_avg=outcome, interval(date -5 -1) by(firm_id)
rangestat (mean) post_avg=outcome, interval(date 0 5) by(firm_id)

Panel Data Applications

webuse nlswork, clear
xtset idcode year

rangestat (mean) wage_trend=ln_wage, interval(year -2 2) by(idcode)
rangestat (mean) lag3_wage=ln_wage, interval(year -3 -1) by(idcode)

Performance Tips

* rangestat can be slow with wide windows on large data
* For simple MA on regular time series, tssmooth or lag operators are faster:
tssmooth ma value_ma5 = value, window(5)
generate ma5_manual = (value + L1.value + L2.value + L3.value + L4.value)/5

---

egenmore: Extended egen Functions

Adds specialized egen functions. Note: many have been incorporated into modern Stata -- check built-in egen first.

Key Functions

ssc install egenmore

* String functions
egen first_word = nss(make), find(1)
egen n_words = nwords(make)
egen combined = concat(var1 var2), punct(" ")
egen clean_text = sieve(text), omit(":$,.")

* Outlier detection (1.5*IQR rule)
egen price_outlier = outside(price)
egen price_out_grp = outside(price), by(foreign)

* Rounding
egen price_rounded = roundi(price), nearest(100)

---

distinct: Count Distinct Values

More flexible than codebook for counting unique values.

ssc install distinct

distinct rep78                    // One variable
distinct foreign rep78, joint     // Joint combinations
distinct mpg, by(foreign)         // By groups

* Stored results
distinct foreign
local n = r(ndistinct)

* Find low-variation variables
foreach var of varlist _all {
    quietly distinct `var'
    if r(ndistinct) < 5 {
        display "`var' has only " r(ndistinct) " distinct values"
    }
}

---

unique: Unique Value Combinations

Like distinct but can tag observations and identify duplicates.

ssc install unique

* Check uniqueness
unique make
* If r(unique) == r(N), variable is a unique ID

* Tag first occurrence
unique foreign rep78, gen(tag)
list foreign rep78 if tag == 1

* Find duplicates in panel
unique idcode year, gen(unique_obs)
list idcode year if unique_obs == 0   // These are duplicates

* Programmatic duplicate check
unique idcode year
if r(unique) != r(N) {
    display "WARNING: Duplicate observations found!"
}

distinct vs unique: Use distinct for quick counts. Use unique when you need to tag/identify specific observations.

---

missings: Missing Data Utilities

ssc install missings

Commands

* Report missing counts
missings report
missings report price mpg rep78
missings report price mpg, by(foreign)

* Cross-tabulation of missing patterns
missings table price mpg rep78

* Tag observations with missing
missings tag price mpg rep78, generate(any_miss)
missings tag price mpg rep78, generate(n_miss) count

* List observations with missing
missings list price mpg rep78

* Drop observations with any missing in specified vars
missings dropobs price mpg rep78
* Drop only if ALL specified vars are missing
missings dropobs price mpg rep78, all

* Drop variables with >50% missing
missings dropvars, min(50)

Real-World Example

webuse nlswork, clear
missings report ln_wage ttl_exp tenure union, by(year)
missings tag ln_wage ttl_exp tenure, generate(n_miss) count

* Compare complete vs incomplete cases
summarize ln_wage if n_miss == 0
summarize ln_wage if n_miss > 0

---

carryforward: Fill Missing Values

Fills missing values by carrying forward (or backward) the last non-missing value within groups.

ssc install carryforward

Usage

* Forward fill within groups
bysort id (year): carryforward value, gen(value_filled)

* Replace in place (destructive -- prefer gen())
carryforward value, replace

* Backward fill: reverse sort, carry forward, re-sort
gsort id -year
by id: carryforward value, gen(value_back)
sort id year

* Forward then backward for max coverage
bysort id (year): carryforward value, gen(value_fwd)
gsort id -year
by id: carryforward value_fwd, replace
sort id year

Panel Data Example

* Fill time-invariant variables with occasional gaps
bysort id (year): carryforward city, gen(city_filled)
bysort id (year): carryforward race, gen(race_filled)

Caution: Only appropriate for time-invariant or slowly-changing variables. Don't use for truly missing data (introduces bias). Always use gen() over replace to preserve originals.

---

Choosing the Right Tool

---

Related Packages

ssc install ftools        // Fast fixed effects (works with gtools)
ssc install reghdfe       // Fast fixed effects regression
ssc install fastxtile     // Fast quantile groups

Diagnostic Packages for Stata

Overview

Packages covered:

• bacondecomp: Bacon decomposition for TWFE DiD -- reveals how staggered treatment timing creates weighted averages of 2x2 comparisons, some potentially biased
• xttest3: Modified Wald test for groupwise heteroskedasticity in panel FE models
• xtcsd: Cross-sectional dependence tests (Pesaran CD, Friedman, Frees) for panel data
• xtserial: Wooldridge test for first-order serial correlation in panels
• vif/collin: Multicollinearity diagnostics (VIF, condition number, variance decomposition)
• hettest/imtest: Breusch-Pagan and White heteroskedasticity tests (built-in)
• ovtest/linktest: Ramsey RESET and link specification tests (built-in)

Installation

// Panel data diagnostics
ssc install xttest3
ssc install xtcsd
ssc install xtserial
ssc install bacondecomp
ssc install collin

// Built-in (no install needed): vif, hettest, imtest, ovtest, linktest

// Verify
which xttest3
which xtcsd
which xtserial
which bacondecomp
which collin

---

Bacon Decomposition for TWFE DiD

The Bacon decomposition (Goodman-Bacon 2021) shows how a TWFE DiD estimate is a weighted average of all 2x2 comparisons. Critical for detecting negative weights and bias from staggered adoption.

Usage

ssc install bacondecomp

bacondecomp outcome treatment_var, ddetail

// With plot options
bacondecomp outcome treatment_var, ///
    ddetail stub(bacon_) ///
    gropt(title("Bacon Decomposition") ///
          xtitle("Weight") ytitle("2x2 DD Estimate"))

Interpreting Results

Comparison types:
  "Earlier T vs. Later C": Early treated vs. not-yet-treated (GOOD)
  "Later T vs. Earlier C": Late treated vs. already-treated (PROBLEMATIC)
  "T vs. Never treated": Treated vs. never-treated (GOOD)

Warning signs:
  - Large weight on "Later T vs. Earlier C" (uses already-treated as control)
  - 2x2 estimates with opposite signs across comparison types
  - Any negative weights

Workflow

// Step 1: TWFE estimate
reghdfe outcome treatment, absorb(unit time) cluster(unit)

// Step 2: Decompose
bacondecomp outcome treatment, ddetail stub(bacon_)

// Step 3: Check problematic weight
preserve
collapse (sum) bacon_weight (mean) bacon_dd_estimate, by(bacon_type)
list bacon_type bacon_weight bacon_dd_estimate
restore

// Step 4: If >20% weight on "Later T vs. Earlier C", use robust estimators:
//   - Sun & Abraham (eventstudyinteract)
//   - Callaway & Sant'Anna (csdid)
//   - De Chaisemartin & D'Haultfoeuille (did_multiplegt)

---

Panel Data Heteroskedasticity: xttest3

Tests whether error variances differ across panel groups (H0: sigma_i^2 = sigma^2 for all i).

Usage

// MUST run after xtreg, fe
xtreg depvar indepvars, fe
xttest3

// No options needed -- automatic after xtreg, fe

Remedies

// If xttest3 rejects (heteroskedasticity detected):

// SOLUTION 1: Robust SE (most common)
xtreg y x1 x2, fe vce(robust)

// SOLUTION 2: Cluster-robust SE (also handles serial correlation)
xtreg y x1 x2, fe vce(cluster id)

// SOLUTION 3: reghdfe (automatically robust)
reghdfe y x1 x2, absorb(id) vce(robust)
// or
reghdfe y x1 x2, absorb(id) vce(cluster id)

// SOLUTION 4: Bootstrap
xtreg y x1 x2, fe vce(bootstrap, reps(1000))

/*
Decision tree:
  Just heteroskedasticity     -> vce(robust)
  + clustering structure      -> vce(cluster id)
  Few clusters (<30)          -> Bootstrap or wild bootstrap
  Need efficiency             -> FGLS (if specification correct)
*/

---

Cross-Sectional Dependence: xtcsd

Tests whether errors are correlated across panel units (common shocks, spatial correlation, network effects).

Usage

// Run after xtreg, fe
xtreg depvar indepvars, fe

xtcsd, pesaran          // Pesaran CD test (recommended, large N large T)
xtcsd, pesaran abs      // With absolute correlation
xtcsd, frees            // Frees' Q test (unbalanced panels)
xtcsd, friedman         // Friedman test (small N)

Remedies

// SOLUTION 1: Driscoll-Kraay SE (recommended for general CSD)
ssc install xtscc
xtscc y x1 x2, fe lag(2)

// SOLUTION 2: Add time fixed effects (absorbs common shocks)
reghdfe y x1 x2, absorb(unit year) vce(cluster unit)

// SOLUTION 3: Common correlated effects (Pesaran 2006)
foreach var of varlist y x1 x2 {
    bysort year: egen `var'_avg = mean(`var')
}
xtreg y x1 x2 y_avg x1_avg x2_avg, fe

// SOLUTION 4: Two-way clustering
reghdfe y x1 x2, absorb(unit) vce(cluster unit year)

/*
Decision tree:
  Large N, large T        -> Driscoll-Kraay SE (xtscc)
  Common shocks           -> Add time fixed effects
  Factor structure        -> CCE or factor models
  Spatial correlation     -> Spatial econometrics
  General case            -> Two-way clustering
*/

---

Serial Correlation: xtserial

Implements Wooldridge's (2002) test for first-order autocorrelation in panel data (H0: no AR(1) serial correlation).

Usage

xtserial depvar indepvars [if] [in], [output]

Example

webuse nlswork, clear
xtset idcode year

xtserial ln_wage tenure ttl_exp age grade
// Reject H0 (p < 0.05): Serial correlation present

Remedies

// SOLUTION 1: Cluster SE by panel unit (most common)
xtreg y x1 x2, fe vce(cluster id)

// SOLUTION 2: AR(1) model (xtregar)
xtregar y x1 x2, fe

// SOLUTION 3: Driscoll-Kraay SE (handles serial + CSD)
xtscc y x1 x2, fe lag(4)

// SOLUTION 4: Dynamic panel GMM (Arellano-Bond)
ssc install xtabond2
xtabond2 y L.y x1 x2, gmm(L.y) iv(x1 x2) robust small

/*
Decision tree:
  Short T          -> Cluster by id
  Long T           -> HAC (Newey-West) or Driscoll-Kraay
  Dynamic model    -> GMM (xtabond2)
  AR(1) structure  -> xtregar
  CSD + serial     -> Driscoll-Kraay (xtscc)
*/

---

Multicollinearity: VIF and collin

Built-in VIF

regress y x1 x2 x3 x4
vif

// Interpretation:
// VIF = 1: No collinearity
// VIF < 5: Acceptable
// 5 <= VIF < 10: Moderate (caution)
// VIF >= 10: Problematic

Enhanced collin

ssc install collin

collin x1 x2 x3 x4

// Additional output:
// - Condition number (<30 OK, 30-100 moderate, >100 severe)
// - Condition indices (>30 indicates problem)
// - Variance decomposition proportions (>0.5 on high index = problem)

Remedies

// SOLUTION 1: Drop redundant variables
correlate weight displacement length
// If r > 0.9, drop one
regress y weight length turn foreign
vif

// SOLUTION 2: Combine correlated variables
egen size = rowmean(weight length)
// Or factor analysis
factor weight length displacement
predict size_factor

// SOLUTION 3: Center variables (especially for interactions)
summarize weight
gen weight_c = weight - r(mean)
summarize length
gen length_c = length - r(mean)
gen weight_X_length = weight_c * length_c
regress y weight_c length_c weight_X_length
vif

// SOLUTION 4: Ridge regression
ssc install ridgereg
ridgereg y x1 x2 x3, model(orr) ridge(0.001 0.01 0.1 1 10)

// SOLUTION 5: Principal components regression
pca weight length turn displacement
predict pc1 pc2
regress y pc1 pc2 foreign
vif  // VIF will be low

// SOLUTION 6: Use joint tests instead of individual tests
regress y weight length turn displacement foreign
testparm weight displacement

/*
Decision tree:
  VIF < 5                           -> No action needed
  Perfect collinearity              -> Drop redundant variable
  Variables measure same concept    -> Combine or choose one
  Polynomial terms causing high VIF -> Center variables
  Prediction focus                  -> Ridge or PCA
  Inference focus                   -> Drop or collect more data
*/

---

Heteroskedasticity Tests: hettest and imtest

Usage

regress y x1 x2 x3

// Breusch-Pagan / Cook-Weisberg (H0: constant variance)
hettest              // Test against fitted values
hettest x1 x2       // Test against specific variables
hettest, fstat       // F-test version

// White's General Test (more general, no functional form assumption)
imtest, white

Remedies

// SOLUTION 1: Robust (Huber-White) SE -- most common
regress y x1 x2, robust

// SOLUTION 2: Log transformation (for positive Y)
gen ln_y = ln(y)
regress ln_y x1 x2
hettest

// SOLUTION 3: WLS (if variance structure known)
quietly regress y x1 x2
predict resid, residuals
gen resid2 = resid^2
regress resid2 x1 x2
predict var_hat
gen wt = 1/var_hat
regress y x1 x2 [aweight=wt]

// SOLUTION 4: Bootstrap SE
regress y x1 x2, vce(bootstrap, reps(1000))

// SOLUTION 5: Cluster-robust SE
regress y x1 x2, cluster(group_var)

/*
Decision tree:
  Just need valid inference   -> robust SE
  Positive Y, wide range     -> log transformation
  Variance structure known    -> WLS
  Clustered data              -> cluster-robust SE
  Unknown structure           -> robust SE or bootstrap
*/

---

Specification Tests: ovtest and linktest

Ramsey RESET Test

Tests for omitted variables/nonlinearity by adding powers of fitted values.

regress y x1 x2 x3
ovtest           // H0: no omitted variables
ovtest, rhs      // Use powers of RHS variables instead

Link Test

Tests specification by regressing y on y-hat and y-hat-squared.

regress y x1 x2 x3
linktest
// _hat should be significant (model has explanatory power)
// _hatsq should NOT be significant (no misspecification)

Remedies for Specification Errors

// SOLUTION 1: Add polynomial terms
gen x1_sq = x1^2
regress y x1 x1_sq x2
ovtest

// SOLUTION 2: Add interactions
gen x1_X_x2 = x1 * x2
regress y x1 x2 x1_X_x2
ovtest

// SOLUTION 3: Fractional polynomials
fracpoly: regress y x1 x2

// SOLUTION 4: Log transformations
gen ln_y = ln(y)
gen ln_x1 = ln(x1)
regress ln_y ln_x1 x2
ovtest

// SOLUTION 5: Splines
mkspline x1_sp = x1, nknots(4) cubic
regress y x1_sp* x2
ovtest

// Compare alternatives with AIC/BIC
estimates stats baseline alt1 alt2 alt3

/*
Decision tree:
  Single variable nonlinear   -> log(Y)/log(X) or fractional polynomial
  Interaction suspected       -> Add interaction terms
  Omitted variables           -> Add controls based on theory
  Very nonlinear, prediction  -> Splines or GAM
  Very nonlinear, inference   -> Fractional polynomial
*/

---

Complete Diagnostic Workflow: Cross-Sectional

capture program drop diagnostic_suite
program define diagnostic_suite
    syntax varlist(min=2) [if] [in]
    gettoken depvar indepvars : varlist

    display "{hline 70}"
    display "COMPREHENSIVE REGRESSION DIAGNOSTICS"
    display "{hline 70}"

    quietly regress `depvar' `indepvars' `if' `in'

    // 1. Heteroskedasticity
    display ""
    display "1. HETEROSKEDASTICITY"
    hettest
    scalar bp_p = r(p)
    quietly imtest, white
    display "White test: chi2=" %7.2f r(chi2) ", p=" %6.4f r(p)
    scalar white_p = r(p)

    if bp_p < 0.05 | white_p < 0.05 {
        display "-> Heteroskedasticity detected: Use robust SE"
    }

    // 2. Specification
    display ""
    display "2. SPECIFICATION"
    quietly regress `depvar' `indepvars' `if' `in'
    ovtest
    scalar reset_p = r(p)
    quietly linktest
    test _hatsq
    scalar link_p = r(p)

    if reset_p < 0.05 | link_p < 0.05 {
        display "-> Misspecification detected: Check functional form"
    }

    // 3. Normality
    display ""
    display "3. NORMALITY"
    quietly regress `depvar' `indepvars' `if' `in'
    predict resid_temp, residuals
    swilk resid_temp
    scalar norm_p = r(p)
    if norm_p < 0.05 {
        display "-> Non-normality: Consider robust inference"
    }

    // 4. Multicollinearity
    display ""
    display "4. MULTICOLLINEARITY"
    quietly regress `depvar' `indepvars' `if' `in'
    vif

    // Recommendation
    display ""
    display "{hline 70}"
    display "RECOMMENDED SPECIFICATION"
    if bp_p < 0.05 | white_p < 0.05 {
        display "-> Use robust SE"
        regress `depvar' `indepvars' `if' `in', robust
    }
    if reset_p < 0.05 | link_p < 0.05 {
        display "-> Check functional form (quadratic, log, interactions)"
    }

    capture drop resid_temp
end

// Usage:
sysuse auto, clear
diagnostic_suite mpg weight length foreign

---

Complete Diagnostic Workflow: Panel Data

capture program drop panel_diagnostics
program define panel_diagnostics
    syntax varlist [if] [in], Panel(varname) Time(varname)
    gettoken depvar indepvars : varlist

    display "{hline 70}"
    display "PANEL DATA DIAGNOSTICS"
    display "{hline 70}"

    xtset `panel' `time'
    quietly xtreg `depvar' `indepvars' `if' `in', fe

    // 1. Groupwise heteroskedasticity
    display ""
    display "1. HETEROSKEDASTICITY (xttest3)"
    xttest3
    scalar het_p = r(p)

    // 2. Serial correlation
    display ""
    display "2. SERIAL CORRELATION (xtserial)"
    xtserial `depvar' `indepvars' `if' `in'
    scalar serial_p = r(p)

    // 3. Cross-sectional dependence
    display ""
    display "3. CROSS-SECTIONAL DEPENDENCE (xtcsd)"
    quietly xtreg `depvar' `indepvars' `if' `in', fe
    xtcsd, pesaran abs
    scalar csd_p = r(p)

    // Recommendation
    display ""
    display "{hline 70}"
    display "RECOMMENDED SPECIFICATION"
    display "{hline 70}"

    if het_p < 0.05 & serial_p < 0.05 & csd_p < 0.05 {
        display "All three violations -> Driscoll-Kraay SE"
        xtscc `depvar' `indepvars' `if' `in', fe lag(2)
    }
    else if het_p < 0.05 & serial_p < 0.05 {
        display "Heteroskedasticity + serial correlation -> Clustered SE"
        xtreg `depvar' `indepvars' `if' `in', fe vce(cluster `panel')
    }
    else if het_p < 0.05 {
        display "Heteroskedasticity only -> Robust SE"
        xtreg `depvar' `indepvars' `if' `in', fe vce(robust)
    }
    else if serial_p < 0.05 {
        display "Serial correlation only -> Clustered SE"
        xtreg `depvar' `indepvars' `if' `in', fe vce(cluster `panel')
    }
    else if csd_p < 0.05 {
        display "CSD only -> Time FE or Driscoll-Kraay"
    }
    else {
        display "No violations -> Conventional SE valid"
    }
end

// Usage:
webuse nlswork, clear
panel_diagnostics ln_wage tenure ttl_exp age grade, panel(idcode) time(year)

---

Quick Reference

// Install all
ssc install bacondecomp xttest3 xtcsd xtserial collin xtscc

// Panel diagnostics
xtreg y x1 x2, fe
xttest3              // Groupwise heteroskedasticity
xtserial y x1 x2    // Serial correlation
xtcsd, pesaran abs   // Cross-sectional dependence

// Cross-sectional diagnostics
regress y x1 x2 x3
hettest              // Breusch-Pagan
imtest, white        // White's test
ovtest               // RESET specification
linktest             // Link specification
vif                  // VIF
collin x1 x2 x3     // Enhanced multicollinearity

// DiD diagnostics
bacondecomp y treatment, ddetail

// Remedies
regress y x1 x2, robust                    // Robust SE
regress y x1 x2, cluster(id)               // Cluster SE
xtscc y x1 x2, fe lag(2)                   // Driscoll-Kraay
xtreg y x1 x2, fe vce(cluster id)          // Panel cluster

References

• Goodman-Bacon (2021). "Difference-in-differences with variation in treatment timing." Journal of Econometrics.
• Hoechle (2007). "Robust standard errors for panel regressions with cross-sectional dependence." Stata Journal.
• Pesaran (2004). "General diagnostic tests for cross section dependence in panels."
• White (1980). "A heteroskedasticity-consistent covariance matrix estimator." Econometrica.
• Wooldridge (2002). Econometric Analysis of Cross Section and Panel Data. MIT Press.

---

Stata version: 17.0+ Required packages: bacondecomp, xttest3, xtcsd, xtserial, collin, xtscc

Modern Difference-in-Differences Packages in Stata

Introduction

Traditional two-way fixed effects (TWFE) regression can produce biased estimates under staggered treatment adoption with heterogeneous effects. TWFE uses already-treated units as controls for newly-treated units ("forbidden comparisons"), creating negative weights.

Modern packages solve this by:

• Using only valid comparisons (never-treated or not-yet-treated as controls)
• Explicitly modeling treatment effect heterogeneity
• Providing transparent aggregation schemes

Packages covered:

• csdid: Callaway-Sant'Anna group-time ATTs with flexible aggregation
• did_multiplegt: de Chaisemartin-D'Haultfoeuille with built-in placebos and dynamics
• did_imputation: Borusyak-Jaravel-Spiess imputation-based estimator (most efficient)
• didregress: Stata 18 built-in (AIPW, IPW, RA, TWFE)
• eventstudyinteract: Sun-Abraham interaction-weighted estimator
• stackedev: Stacked event study approach

Installation

* Core packages
ssc install csdid
ssc install drdid
ssc install did_multiplegt
ssc install did_imputation
ssc install eventstudyinteract
ssc install bacondecomp

* Supporting packages
ssc install reghdfe
ssc install ftools
ssc install avar
ssc install boottest
ssc install coefplot
ssc install event_plot

Diagnosing the Problem: Goodman-Bacon Decomposition

ssc install bacondecomp

* Decompose TWFE into component 2x2 DiD comparisons
bacondecomp outcome, ddetail stub(bacon_)

* Visualize: look for negative weights and heterogeneous estimates
scatter bacon_estimate bacon_weight, ///
    yline(`twfe_coef', lcolor(red) lpattern(dash)) yline(0)

When is TWFE valid? Only when treatment effects are homogeneous across units and time, and there are no dynamic effects.

Data Structure Requirements

All modern DiD packages need:

* Key variables:
* - unit_id: Panel identifier
* - year: Time period
* - treatment_year: Year unit FIRST receives treatment
*   (0 or missing for never-treated units)
* - outcome: Outcome variable

* Create binary treatment indicator
generate treated = (year >= treatment_year & treatment_year > 0)

---

csdid: Callaway and Sant'Anna (2021)

Estimates group-time ATTs (ATT(g,t)) and aggregates them flexibly.

Syntax

csdid outcome [covariates], ///
    ivar(unit_id) ///           // Panel identifier
    time(time_var) ///          // Time variable
    gvar(first_treat_var) ///   // First treatment period (0 or . for never-treated)
    [method(dripw|reg|ipw)] /// // Estimation method
    [notyet | long2]            // Control group choice

Key options:

• method(dripw): Doubly-robust (default, most robust)
• method(reg): Regression-based (faster)
• method(ipw): Inverse probability weighting
• notyet: Use not-yet-treated as controls (default)
• long2: Use only never-treated as controls
• wboot rseed(#): Wild bootstrap SEs
• wboot_reps(#): Number of bootstrap reps

Basic Usage

csdid wage_outcome age grade, ///
    ivar(idcode) time(year) gvar(treatment_year) ///
    method(dripw) notyet

Aggregation (estat commands)

* Overall ATT
estat simple
matrix simple_att = r(table)
local att = simple_att[1,1]
local se = simple_att[2,1]

* By treatment cohort
estat group
csdid_plot, group

* By calendar time
estat calendar
csdid_plot, calendar

* Event study
estat event, window(-5 10) estore(cs_event)
csdid_plot, name(event_study, replace)

* Pre-trend test
estat pretrend, pre(5)

Sensitivity to Control Group

* Compare not-yet-treated vs never-treated
csdid outcome, ivar(id) time(year) gvar(treat_year) method(dripw) notyet
estat simple
local att_notyet = r(table)[1,1]

csdid outcome, ivar(id) time(year) gvar(treat_year) method(dripw) long2
estat simple
local att_never = r(table)[1,1]

display "Not-yet: " `att_notyet' "  Never: " `att_never'

---

did_multiplegt: de Chaisemartin and D'Haultfoeuille (2020)

Heterogeneity-robust estimator with built-in placebo tests and dynamic effects. Can handle treatment switchers.

Syntax

did_multiplegt outcome unit_id time_var treatment_var, ///
    [robust_dynamic] ///        // Robust SEs with dynamic effects
    [dynamic(#)] ///            // Number of post-treatment periods
    [placebo(#)] ///            // Number of pre-treatment placebo periods
    [breps(#)] ///              // Bootstrap replications for SEs
    [cluster(varname)] ///      // Cluster variable
    [controls(varlist)] ///     // Time-varying controls
    [trends_nonparam(varname)] ///  // Unit-specific nonparametric trends
    [jointtestplacebo] ///      // Joint test of all placebos
    [only_never_switchers]      // Restrict to absorbing treatment

Basic Usage

did_multiplegt price hospital_id year merged, ///
    robust_dynamic dynamic(5) placebo(3) breps(100) cluster(hospital_id)

* Results stored in e()
display "Instantaneous: " e(effect) " (" e(se_effect) ")"
forvalues k = 1/5 {
    display "Dynamic +`k': " e(dynamic_`k') " (" e(se_dynamic_`k') ")"
}
forvalues k = 1/3 {
    display "Placebo -`k': " e(placebo_`k') " (" e(se_placebo_`k') ")"
}

Joint Placebo Test

did_multiplegt price hospital_id year merged, ///
    robust_dynamic dynamic(5) placebo(3) breps(200) ///
    cluster(hospital_id) jointtestplacebo

display "Joint placebo p-value: " e(p_jointplacebo)

Plotting Dynamic Effects

* Extract results and build dataset manually
preserve
clear
local n_placebo = 3
local n_dynamic = 5
set obs `= `n_placebo' + 1 + `n_dynamic''
generate period = _n - `n_placebo' - 1
generate estimate = .
generate se = .

forvalues k = 1/`n_placebo' {
    local row = `n_placebo' - `k' + 1
    replace estimate = e(placebo_`k') in `row'
    replace se = e(se_placebo_`k') in `row'
}
local row = `n_placebo' + 1
replace estimate = e(effect) in `row'
replace se = e(se_effect) in `row'
forvalues k = 1/`n_dynamic' {
    local row = `n_placebo' + 1 + `k'
    replace estimate = e(dynamic_`k') in `row'
    replace se = e(se_dynamic_`k') in `row'
}

generate ci_lower = estimate - 1.96 * se
generate ci_upper = estimate + 1.96 * se

twoway (scatter estimate period) (rcap ci_lower ci_upper period), ///
    xline(-0.5, lcolor(red) lpattern(dash)) yline(0) ///
    xtitle("Periods Relative to Treatment") ytitle("Effect") legend(off)
restore

Saving Results

did_multiplegt price hospital_id year merged, ///
    robust_dynamic dynamic(5) placebo(3) breps(100) ///
    cluster(hospital_id) save_results("did_multiplegt_output")
* Creates did_multiplegt_output.dta

---

did_imputation: Borusyak, Jaravel, and Spiess (2021)

Imputation-based approach: imputes counterfactual outcomes for treated units using untreated observations. Most efficient estimator under homoskedasticity.

Syntax

did_imputation outcome unit_id time_var first_treat_var, ///
    [horizons(numlist)] ///     // Event times (e.g., 0/10)
    [pretrend(#)] ///           // Pre-treatment periods to test
    [autosample] ///            // Automatically balance sample
    [minn(#)] ///               // Min obs per group-time
    [controls(varlist)] ///     // Control variables (partialled out)
    [cluster(varname)] ///      // Cluster SEs
    [agg(simple|event|cohort|time)] // Aggregation type

Basic Usage

did_imputation employed idcode year first_treatment, autosample minn(0)

Event Study

did_imputation employed idcode year first_treatment, ///
    horizons(0/10) pretrend(5) autosample cluster(idcode)

* Plot
event_plot, ///
    default_look stub_lag(tau#) stub_lead(pre#) together ///
    graph_opt(xtitle("Periods Since Treatment") ytitle("Effect") xlabel(-5(1)10))

Pre-trend Testing

did_imputation employed idcode year first_treatment, ///
    horizons(0/10) pretrend(5) autosample cluster(idcode)

* Joint test of pre-treatment coefficients
test pre1 pre2 pre3 pre4 pre5

Aggregation Types

* Cohort-specific effects
did_imputation employed idcode year first_treatment, horizons(0/5) agg(cohort) autosample

* Calendar-time effects
did_imputation employed idcode year first_treatment, horizons(0/5) agg(time) autosample

---

didregress: Stata 18 Built-in Command

Requires Stata 18+. Multiple estimators in one command.

Syntax

didregress (outcome [covariates]) (treatment_var [covariates]), ///
    group(unit_id) time(time_var) ///
    [aipw | ipw | ra | twfe] ///
    [vce(cluster varname)]

• Covariates in outcome equation: for RA and AIPW
• Covariates in treatment equation: for IPW and AIPW

Usage

* AIPW (default, doubly robust)
didregress (outcome gdp_growth) (treated gdp_growth), ///
    group(state) time(year) aipw vce(cluster state)

* Post-estimation aggregations
estat aggregation, event window(-5 10)
estat aggregation, cohort
estat aggregation, time
estat grlevelsof

Advantages: Official support, multiple estimators, clean post-estimation. Limitations: Stata 18 only, fewer diagnostics than specialized packages, no built-in event study plots.

---

eventstudyinteract: Sun and Abraham (2021)

Interaction-weighted estimator avoiding negative weights.

ssc install eventstudyinteract

* Create interaction terms for each cohort x relative-time
levelsof treatment_year if treatment_year > 0, local(cohorts)
foreach g of local cohorts {
    forvalues k = -5/10 {
        if `k' != -1 {
            generate treat_`g'_`k' = (treatment_year == `g' & rel_time == `k')
        }
    }
}

eventstudyinteract outcome treat_*, ///
    absorb(unit_id year) ///
    cohort(treatment_year) ///
    control_cohort(never_treated) ///
    vce(cluster unit_id)

event_plot, default_look ///
    graph_opt(xtitle("Event Time") ytitle("Effect"))

---

Stacked Event Study (Manual)

Stack separate datasets per cohort, each with its own never-treated control group.

levelsof treatment_year if !missing(treatment_year), local(cohorts)
tempfile stacked
local first = 1

foreach cohort of local cohorts {
    preserve
    keep if treatment_year == `cohort' | missing(treatment_year)
    generate stack_cohort = `cohort'
    generate stack_id = unit_id * 10000 + `cohort'
    generate stack_reltime = year - `cohort'

    if `first' {
        save `stacked', replace
        local first = 0
    }
    else {
        append using `stacked'
        save `stacked', replace
    }
    restore
}

use `stacked', clear

* Create event time dummies (omit -1)
forvalues k = 5(-1)1 {
    generate lead`k' = (stack_reltime == -`k')
}
forvalues k = 0/10 {
    generate lag`k' = (stack_reltime == `k')
}

reghdfe outcome lead* lag*, absorb(stack_id stack_cohort#year) cluster(unit_id)

---

When to Use Each Package

Decision rule: If staggered adoption, always use at least one robust method. Report TWFE alongside for comparison. If estimates diverge, heterogeneity is present and TWFE is biased.

---

Event Study Setup

* Relative time variable
generate event_time = year - treatment_year
replace event_time = -1000 if missing(treatment_year) | treatment_year == 0

* Event time dummies (omit -1 as reference)
forvalues k = 10(-1)2 {
    generate lead`k' = (event_time == -`k')
}
forvalues k = 0/10 {
    generate lag`k' = (event_time == `k')
}

GOTCHA — event-study interaction syntax:

• Correct: ib(-1).rel_time#1.treated — base prefix goes on the factored variable
• Wrong: treated#ib(-1).rel_time — produces wrong model or error
• NEVER use hyphenated variable names for negative event times (e.g., lead-3). Stata parses event-3 as event MINUS 3, not a variable name. Use: local pre = abs(k'); gen leadpre' = (rel_time == -k')`
• Use testparm for testing negative factor levels, not test (which fails with negative levels)
• Prefer factor notation (i.rel_time#i.treated) over manual dummies to avoid these issues entirely

Binning Distant Endpoints

generate event_time_binned = event_time
replace event_time_binned = -10 if event_time < -10 & event_time != -1000
replace event_time_binned = 10 if event_time > 10

Testing Pre-Trends

* Joint test of all pre-treatment coefficients
test lead10 lead9 lead8 lead7 lead6 lead5 lead4 lead3 lead2

* Or after csdid:
estat pretrend, pre(5)

---

Standard Errors and Clustering

General rule: Cluster at the level of treatment assignment or higher.

By Package

* csdid: auto-clusters at ivar() level; wild bootstrap available
csdid outcome, ivar(id) time(year) gvar(g) wboot rseed(12345)

* did_multiplegt: explicit cluster + bootstrap reps
did_multiplegt outcome id year treat, breps(200) cluster(state_id)

* did_imputation: cluster option
did_imputation outcome id year g, autosample cluster(id)

* didregress: standard vce()
didregress (outcome) (treat), group(id) time(year) vce(cluster state_id)

Few Clusters (< 30)

* Wild cluster bootstrap (install boottest)
ssc install boottest
reghdfe outcome treatment, absorb(unit_id year) cluster(state_id)
boottest treatment, cluster(state_id) boottype(wild) reps(9999) seed(12345)

Two-Way Clustering

reghdfe outcome treatment, absorb(unit_id year) vce(cluster state_id year)

---

Side-by-Side Comparison of All Estimators

* 1. TWFE
reghdfe outcome treated x1 x2, absorb(unit_id year) cluster(unit_id)
local est1 = _b[treated]

* 2. Callaway-Sant'Anna
csdid outcome x1 x2, ivar(unit_id) time(year) gvar(treatment_year) method(dripw) notyet
estat simple
local est2 = r(table)[1,1]

* 3. did_multiplegt
did_multiplegt outcome unit_id year treated, robust_dynamic breps(100) cluster(unit_id)
local est3 = e(effect)

* 4. did_imputation
did_imputation outcome unit_id year treatment_year, autosample cluster(unit_id)
local est4 = _b[tau]

* 5. didregress (Stata 18)
capture didregress (outcome x1 x2) (treated x1 x2), group(unit_id) time(year) aipw
capture local est5 = _b[ATET:r1vs0.treated]

* Compare
display "TWFE: " `est1'
display "CS:   " `est2'
display "DM:   " `est3'
display "DI:   " `est4'

---

Complete Workflow Template (csdid)

clear all
set seed 98765

use "your_data.dta", clear

* 1. TWFE baseline
reghdfe outcome treated covars, absorb(unit_id year) vce(cluster unit_id)
local twfe = _b[treated]

* 2. Bacon decomposition
bacondecomp outcome, ddetail

* 3. Main estimation
csdid outcome covars, ivar(unit_id) time(year) gvar(treatment_year) ///
    method(dripw) notyet

* 4. Aggregations
estat simple                              // Overall ATT
estat group                               // By cohort
estat event, window(-5 10) estore(es)     // Event study
estat pretrend, pre(5)                    // Pre-trend test

* 5. Plots
csdid_plot, name(event_study, replace)

* 6. Robustness: never-treated controls
csdid outcome covars, ivar(unit_id) time(year) gvar(treatment_year) ///
    method(dripw) long2
estat simple

* 7. Robustness: alternative methods
csdid outcome covars, ivar(unit_id) time(year) gvar(treatment_year) ///
    method(reg) notyet
estat simple

csdid outcome covars, ivar(unit_id) time(year) gvar(treatment_year) ///
    method(ipw) notyet
estat simple

Pre-Publication Checklist

• [ ] Test parallel trends (visual + formal estat pretrend)
• [ ] Show event study plot
• [ ] Report Goodman-Bacon decomposition (if staggered)
• [ ] Use at least one heterogeneity-robust estimator
• [ ] Report treatment effect dynamics
• [ ] Cluster standard errors at treatment level
• [ ] Check robustness to control group choice (not-yet vs never)
• [ ] Report number of treated/control units and cohort sizes

estout/esttab: Publication-Quality Tables in Stata

Overview

The estout package is the most widely used tool for creating publication-quality regression tables in Stata (by Ben Jann). It stores estimation results and exports formatted side-by-side model comparisons to LaTeX, HTML, RTF (Word), and CSV.

Main commands:

• eststo - Store estimation results
• esttab - User-friendly table creation and export
• estout - Low-level flexible table creation
• estadd - Add custom statistics to stored results

---

Contents

• Installation
• Quick Start
• eststo: Storing Estimation Results
• esttab: Creating Tables
• estadd: Adding Custom Statistics
• estout: Low-Level Flexible Tables
• Output Formats
• Multi-Column Model Comparisons
• Summary Statistics Tables
• Advanced Features
• Publication Template: AER Style
• Publication Template: Multi-Panel
• Tips and Best Practices
• Common Issues
• Quick Reference

---

Installation

ssc install estout
help esttab
help estout
help eststo
help estadd

---

Quick Start

sysuse auto, clear

eststo clear
eststo: regress price mpg
eststo: regress price mpg weight
eststo: regress price mpg weight foreign

* Display in Stata
esttab

* Export to LaTeX / CSV / Word
esttab using results.tex, replace
esttab using results.csv, replace csv
esttab using results.rtf, replace

---

eststo: Storing Estimation Results

eststo clear                                    // Clear all stored estimates

eststo: regress price mpg weight                // Auto-named (est1, est2, ...)
eststo model1: regress price mpg weight foreign // Custom name
quietly eststo: regress price mpg weight        // Store without displaying

* Manage stored estimates
estimates dir                   // List stored estimates
estimates replay model1         // Display a stored estimate
eststo drop model1              // Drop specific estimate
estimates save mymodel, replace // Save to disk
estimates use mymodel           // Load from disk

* Store with title
eststo model1, title("Basic Model"): regress price mpg weight

* Works with prefixed commands
eststo: xi: regress price i.rep78 mpg weight
eststo: ivregress 2sls price (mpg = gear_ratio) weight

---

esttab: Creating Tables

Basic Options

eststo clear
eststo: regress price mpg
eststo: regress price mpg weight
eststo: regress price mpg weight foreign

esttab               // Default output
esttab, se           // Standard errors in parentheses
esttab, t            // t-statistics
esttab, p            // p-values
esttab, ci           // Confidence intervals
esttab, not          // Suppress secondary stats
esttab, b(3) se(3)   // 3 decimal places
esttab, compress     // Compact output
esttab, wide         // Wide format
esttab, plain        // No formatting

Titles, Labels, and Stars

esttab, label                                        // Use variable labels
esttab, mtitles("Model 1" "Model 2" "Model 3")     // Column titles
esttab, nomtitles                                    // No column titles
esttab, title("Table 1: Price Regressions")         // Table title
esttab, star(* 0.10 ** 0.05 *** 0.01)              // Significance stars
esttab, nostar                                       // No stars

* Model groups (for panel headers over columns)
esttab, mgroups("OLS Models" "IV Models", pattern(1 0 0 1 0 0))

Statistics and Scalars

* Add model statistics below coefficients
esttab, stats(N r2 r2_a F, ///
    labels("Observations" "R-squared" "Adjusted R-squared" "F-statistic") ///
    fmt(%9.0fc %9.3f %9.3f %9.2f))

* Notes
esttab, addnote("Standard errors in parentheses." ///
                "* p<0.10, ** p<0.05, *** p<0.01")

* Replace default note
esttab, nonotes addnote("Custom note here")

Variable Selection and Ordering

esttab, keep(mpg weight foreign)          // Keep specific variables
esttab, drop(_cons)                        // Drop variables
esttab, order(foreign mpg weight)          // Reorder variables

* Custom variable labels in table (overrides label variable)
esttab, varlabels(mpg "Miles per Gallon" ///
                  weight "Weight (lbs)" ///
                  foreign "Foreign Car" ///
                  _cons "Constant")

* Equivalent alternative
esttab, coeflabels(mpg "Miles per Gallon" ///
                   weight "Vehicle Weight")

---

estadd: Adding Custom Statistics

eststo clear
eststo: regress price mpg weight foreign

* Add scalar
estadd scalar mystat = 1.234

* Add string (e.g., fixed effects indicators)
estadd local controls "Yes"

* Compute and add
quietly summarize price
estadd scalar mean_y = r(mean)
estadd scalar rmse = e(rmse)

* AIC/BIC
estat ic
matrix ic = r(S)
estadd scalar aic = ic[1,5]
estadd scalar bic = ic[1,6]

* Display with custom stats
esttab, scalars("mean_y Mean of Dep. Var." "rmse RMSE" "aic AIC" "bic BIC")

Fixed Effects Indicators Pattern

webuse nlswork, clear

eststo clear
eststo: regress ln_wage age ttl_exp tenure
estadd local fe_ind "No"
estadd local fe_year "No"

eststo: areg ln_wage age ttl_exp tenure, absorb(idcode)
estadd local fe_ind "Yes"
estadd local fe_year "No"

eststo: xtreg ln_wage age ttl_exp tenure i.year, fe
estadd local fe_ind "Yes"
estadd local fe_year "Yes"

esttab, scalars("fe_ind Individual FE" "fe_year Year FE")

Adding Test Statistics

eststo clear
eststo: regress price mpg weight foreign length headroom trunk
test mpg weight foreign
estadd scalar F_joint = r(F)
estadd scalar p_joint = r(p)

esttab, scalars("F_joint F-test (joint)" "p_joint p-value (joint)")

---

estout: Low-Level Flexible Tables

estout gives more control over cell formatting than esttab.

eststo clear
eststo: regress price mpg weight
eststo: regress price mpg weight foreign

* Custom cell layout
estout, cells(b(star fmt(3)) se(par fmt(3)))

* Multiple rows per coefficient
estout, cells("b(star fmt(3) label(Coef.))" ///
              "se(par fmt(3) label(Std. Err.))" ///
              "t(fmt(2) label(t-stat))")

* Confidence intervals
estout, cells(b(star) ci(par))

* Standardized (beta) coefficients
estout, cells(b(star fmt(3)) beta(fmt(3)))

---

Output Formats

LaTeX

* Basic LaTeX with booktabs
esttab using table.tex, replace ///
    booktabs label ///
    b(%9.3f) se(%9.3f) ///
    star(* 0.10 ** 0.05 *** 0.01) ///
    mtitles("Basic" "Full Model") ///
    stats(N r2 r2_a, labels("Observations" "R-squared" "Adj. R-squared")) ///
    title("Determinants of Price\label{tab:price}") ///
    addnote("Standard errors in parentheses." ///
            "* \(p<0.10\), ** \(p<0.05\), *** \(p<0.01\)")

* Longtable (multi-page)
esttab using table.tex, replace booktabs longtable label

* Custom prehead/postfoot for full table environment control
esttab using table.tex, replace ///
    booktabs b(3) se(3) label ///
    prehead("\begin{table}[htbp]\centering" ///
            "\caption{Results}\label{tab:main}" ///
            "\begin{tabular}{l*{@M}{c}}" ///
            "\toprule") ///
    postfoot("\bottomrule" ///
             "\end{tabular}" ///
             "\end{table}")

Gotcha: In LaTeX, escape special characters: \%, \$, \_, \&. In label text use \(p<0.05\) to get math mode.

HTML

esttab using table.html, replace html label ///
    b(%9.3f) se(%9.3f) ///
    star(* 0.10 ** 0.05 *** 0.01) ///
    mtitles("Model 1" "Model 2") ///
    title("Table 1: Regression Results")

RTF (Word)

esttab using table.rtf, replace label ///
    b(%9.3f) se(%9.3f) ///
    star(* 0.10 ** 0.05 *** 0.01) ///
    mtitles("Model 1" "Model 2") ///
    note("Standard errors in parentheses.")

CSV

esttab using table.csv, replace csv label ///
    b(3) se(3) star(* 0.10 ** 0.05 *** 0.01)

* Plain CSV for Excel
esttab using table.csv, replace csv plain label b(3) se(3)

---

Multi-Column Model Comparisons

Progressive Specifications

sysuse auto, clear
eststo clear
eststo m1: regress price mpg
eststo m2: regress price mpg weight
eststo m3: regress price mpg weight foreign
eststo m4: regress price mpg weight foreign length headroom

esttab m1 m2 m3 m4, ///
    b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) label ///
    mtitles("(1)" "(2)" "(3)" "(4)") ///
    stats(N r2 r2_a, labels("Observations" "R-squared" "Adj. R-squared"))

Different Estimators

eststo clear
eststo ols: regress price mpg weight foreign
eststo robust: regress price mpg weight foreign, robust
eststo cluster: regress price mpg weight foreign, vce(cluster rep78)

esttab ols robust cluster, ///
    b(3) se(3) label ///
    mtitles("OLS" "Robust SE" "Clustered SE")

Different Subsamples

eststo clear
eststo full: regress price mpg weight foreign
eststo domestic: regress price mpg weight if foreign == 0
eststo foreign_only: regress price mpg weight if foreign == 1

esttab full domestic foreign_only, ///
    b(3) se(3) label ///
    mtitles("Full Sample" "Domestic" "Foreign")

Panel Groupings with mgroups

eststo clear
eststo pa1: regress price mpg weight
eststo pa2: regress price mpg weight foreign
eststo pb1: ivregress 2sls price weight (mpg = gear_ratio)
eststo pb2: ivregress 2sls price weight foreign (mpg = gear_ratio)

esttab pa1 pa2 pb1 pb2 using multipanel.tex, replace booktabs ///
    b(3) se(3) label ///
    mgroups("Panel A: OLS" "Panel B: 2SLS", pattern(1 0 1 0)) ///
    mtitles("(1)" "(2)" "(3)" "(4)")

---

Summary Statistics Tables

sysuse auto, clear

* Basic summary statistics
estpost summarize price mpg weight length
esttab, cells("count mean(fmt(2)) sd(fmt(2)) min max") ///
    label title("Summary Statistics") nomtitles

* Formatted for LaTeX
estpost summarize price mpg weight length
esttab using summary.tex, replace booktabs ///
    cells("count(fmt(%9.0fc)) mean(fmt(%9.2fc)) sd(fmt(%9.2fc)) min(fmt(%9.2fc)) max(fmt(%9.2fc))") ///
    collabels("N" "Mean" "SD" "Min" "Max") ///
    label nomtitles ///
    title("Table 1: Descriptive Statistics\label{tab:summary}")

* Detailed (with percentiles)
estpost summarize price mpg weight, detail
esttab, cells("count mean sd min p25 p50 p75 max") label nomtitles

Summary by Groups with T-Test

eststo clear
eststo domestic: estpost summarize price mpg weight if foreign == 0
eststo foreign: estpost summarize price mpg weight if foreign == 1
eststo diff: estpost ttest price mpg weight, by(foreign) unequal

esttab domestic foreign diff, ///
    cells("mean(fmt(2)) sd(fmt(2) par)") ///
    label mtitles("Domestic" "Foreign" "Difference")

---

Advanced Features

Factor Variables and Interactions

eststo clear
eststo: regress price i.rep78 c.mpg##c.weight i.foreign#c.mpg

* Hide FE coefficients but indicate their presence
esttab, label ///
    drop(*.rep78) ///
    indicate("Repair Record FE = *.rep78")

* With Yes/No labels
esttab, label ///
    indicate("Repair FE = *.rep78" "Headroom FE = *.headroom", labels("Yes" "No"))

Margins and Marginal Effects

sysuse auto, clear
eststo clear

* Logit coefficients
quietly logit foreign mpg weight price
estimates store logit_model

* Average marginal effects
margins, dydx(*) post
eststo ame

* Marginal effects at means
quietly logit foreign mpg weight price
margins, dydx(*) atmeans post
eststo mem

esttab logit_model ame mem, ///
    b(3) se(3) label ///
    mtitles("Logit Coef." "AME" "MEM")

Multiple Equation Models

* Heckman selection model
webuse womenwk, clear
eststo clear
eststo: heckman wage education age, ///
    select(married children education age)

esttab, label b(3) se(3) ///
    eqlabels("Wage Equation" "Selection Equation")

Instrumental Variables

sysuse auto, clear
eststo clear
eststo ols: regress price mpg weight
eststo iv: ivregress 2sls price weight (mpg = gear_ratio turn)
eststo first: regress mpg gear_ratio turn weight

esttab ols iv first, ///
    b(3) se(3) label ///
    mtitles("OLS" "2SLS" "First Stage") ///
    stats(N r2, labels("Observations" "R-squared"))

Fixed Effects with Indicators

webuse nlswork, clear
eststo clear

eststo pooled: regress ln_wage age ttl_exp tenure
estadd local fe "No"

eststo fe: xtreg ln_wage age ttl_exp tenure, fe
estadd local fe "Yes"

eststo fe_time: xtreg ln_wage age ttl_exp tenure i.year, fe
estadd local fe "Yes"
estadd local time_fe "Yes"

esttab pooled fe fe_time, ///
    b(3) se(3) label ///
    drop(*.year) ///
    scalars("fe Individual FE" "time_fe Year FE") ///
    mtitles("Pooled" "FE" "FE+Year")

---

Publication Template: AER Style

sysuse auto, clear
eststo clear
eststo: regress price mpg weight
eststo: regress price mpg weight foreign
eststo: regress price mpg weight foreign length

esttab using "table_aer.tex", replace ///
    booktabs b(%9.3f) se(%9.3f) ///
    star(* 0.10 ** 0.05 *** 0.01) label nomtitles ///
    mgroups("Dependent Variable: Price", pattern(1 0 0) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) ///
            span erepeat(\cmidrule(lr){@span})) ///
    stats(N r2, fmt(%9.0fc %9.3f) labels("Observations" "R-squared")) ///
    prehead("\begin{table}[htbp]\centering" ///
            "\caption{Determinants of Automobile Prices}" ///
            "\label{tab:main}" ///
            "\begin{tabular}{l*{3}{c}}" ///
            "\toprule") ///
    posthead("\midrule") ///
    prefoot("\midrule") ///
    postfoot("\bottomrule" ///
             "\multicolumn{4}{p{0.9\textwidth}}{\footnotesize \textit{Notes:} Standard errors in parentheses. * \(p<0.10\), ** \(p<0.05\), *** \(p<0.01\).}" ///
             "\\" ///
             "\end{tabular}" ///
             "\end{table}")

Publication Template: Multi-Panel

eststo clear
eststo pa1: regress price mpg weight
eststo pa2: regress price mpg weight foreign
eststo pb1: regress price mpg weight, robust
eststo pb2: regress price mpg weight foreign, robust

esttab pa1 pa2 pb1 pb2 using "table_panels.tex", replace ///
    booktabs b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) label ///
    mgroups("Panel A: Standard SE" "Panel B: Robust SE", ///
            pattern(1 0 1 0) ///
            prefix(\multicolumn{@span}{c}{) suffix(}) ///
            span erepeat(\cmidrule(lr){@span})) ///
    mtitles("(1)" "(2)" "(3)" "(4)") ///
    stats(N r2, fmt(%9.0fc %9.3f) labels("Observations" "R-squared"))

---

Tips and Best Practices

Workflow Pattern

* 1. Clear at start
eststo clear

* 2. Store with descriptive names
eststo baseline: regress y x1 x2

* 3. Add custom stats immediately after estimation
estadd local fe "Yes"

* 4. Preview in Stata window first
esttab, label

* 5. Then export
esttab using "table.tex", replace booktabs label ...

Debugging

• Preview first: Always esttab, label before exporting to file
• LaTeX errors: Check for unescaped special characters (%, $, &, _, #)
• Missing packages: LaTeX needs \usepackage{booktabs} for booktabs option
• Wide tables: Use compress, smaller font (\small), or landscape mode

Journal-Specific Choices

* Standard errors vs t-stats vs CI
esttab, se    // Most journals
esttab, t     // Some journals
esttab, ci    // Medical journals

* Star conventions vary
esttab, star(* 0.10 ** 0.05 *** 0.01)   // Economics
esttab, star(+ 0.10 * 0.05 ** 0.01)     // Political science
esttab, star(* 0.05 ** 0.01 *** 0.001)  // Some social sciences

---

Common Issues

---

Quick Reference

Key Options

Common Patterns

* Store + display
eststo clear
eststo: regress y x1 x2
esttab, label se r2

* Export to LaTeX
esttab using "table.tex", replace booktabs label b(3) se(3) ///
    star(* 0.10 ** 0.05 *** 0.01) stats(N r2)

* Export to Word
esttab using "table.rtf", replace label b(3) se(3)

* Summary statistics
estpost summarize var1 var2 var3
esttab, cells("mean sd min max") label

* Clear stored estimates
eststo clear

Help Files

help esttab       // Main table command
help estout       // Low-level table command
help eststo       // Store estimates
help estadd       // Add custom statistics

Official documentation: http://repec.sowi.unibe.ch/stata/estout/

Graph Schemes: Professional Visualization Design

Table of Contents

• Installation and Setup
• grstyle: Dynamic Graph Styling
• schemepack: Professional Scheme Collection
• blindschemes: Colorblind-Friendly Palettes
• Activating Schemes
• Color Palette Selection
• Font and Size Customization
• Creating Custom Schemes
• Publication Requirements
• Quick Reference

---

Installation and Setup

Install All Major Packages

ssc install grstyle, replace
ssc install palettes, replace
ssc install colrspace, replace
ssc install schemepack, replace
ssc install blindschemes, replace

Built-in Stata Schemes

graph query, schemes       // View all available schemes

set scheme s2color         // Default colored scheme
set scheme s1mono          // Monochrome scheme
set scheme s1color         // Alternative color scheme

Setting Default Scheme

set scheme plotplain              // Current session only
set scheme plotplain, permanently // Persists across sessions

query graphics                    // Check current scheme

Profile Setup

Create a profile.do to set preferences automatically:

• Windows: C:\Users\[username]\ado\profile.do
• Mac: ~/Library/Application Support/Stata/ado/profile.do
• Linux: ~/.stata/ado/profile.do

// Contents of profile.do:
set scheme plotplain, permanently

Workspace Setup (Top of Do-File)

clear all
set more off
set scheme plotplain

grstyle clear
grstyle init
grstyle set plain, box
grstyle set color tableau

---

grstyle: Dynamic Graph Styling

On-the-fly customization without creating permanent scheme files. Requires palettes and colrspace as dependencies.

Basic Usage

grstyle init          // Must run first
grstyle set plain     // Clean graph style

scatter price mpg     // Uses new style

Color Palettes

grstyle set color tableau             // Tableau palette
grstyle set color viridis             // Colorblind-friendly
grstyle set color Set1                // ColorBrewer
grstyle set color tol bright          // Paul Tol palettes
grstyle set color tol muted
grstyle set color okabe               // Okabe & Ito

// RGB values
grstyle set color "0 114 178" "213 94 0" "0 158 115"

// Named colors
grstyle set color navy maroon forest_green

Symbols, Lines, Grids

grstyle set symbol O D T S            // Marker symbols
grstyle set lpattern solid dash dot   // Line patterns
grstyle set linewidth thin medium thick

grstyle set plain, box                // White bg with box
grstyle set grid, horizontal          // Horizontal gridlines
grstyle set nogrid                    // Remove grids

Font Sizes

grstyle set size 11pt                 // Overall
grstyle set size axis_title: 12pt     // Specific elements
grstyle set size tick_label: 10pt
grstyle set size title: 14pt

Comprehensive Setup

grstyle clear
grstyle init
grstyle set plain, box horizontal
grstyle set color tableau
grstyle set symbol
grstyle set lpattern
grstyle set linewidth thin medium thick
grstyle set size 11pt

Resetting

grstyle clear           // Clear all grstyle settings
set scheme s2color      // Return to default

---

schemepack: Professional Scheme Collection

Available Schemes

Usage

set scheme 538
scatter price mpg, title("FiveThirtyEight Style")

set scheme economist
scatter price mpg, title("Economist Style")

set scheme tableau
scatter price mpg, title("Tableau Style")

set scheme tufte
scatter price mpg, title("Tufte Minimalist")

---

blindschemes: Colorblind-Friendly Palettes

Installed via ssc install blindschemes. Includes:

plotplain for Publications

set scheme plotplain

scatter price mpg, ///
    title("Automobile Prices and Fuel Efficiency") ///
    ytitle("Price (USD)") xtitle("Miles per Gallon")

Colorblind-Safe Multi-Series

set scheme plotplainblind

twoway (scatter price mpg if rep78<=2, msymbol(O)) ///
       (scatter price mpg if rep78==3, msymbol(D)) ///
       (scatter price mpg if rep78>=4, msymbol(T)), ///
    legend(order(1 "Poor" 2 "Average" 3 "Good"))

Extra Accessibility: Color + Patterns

set scheme plotplainblind

graph bar (mean) price, over(rep78) ///
    asyvars showyvars legend(off) blabel(group)

---

Activating Schemes

Three Levels of Scope

// 1. Permanent (persists across sessions)
set scheme plotplain, permanently

// 2. Session-level (all graphs in session)
set scheme plotplain

// 3. Graph-specific (single graph only)
scatter price mpg, scheme(economist)
histogram mpg  // Uses session default, not economist

grstyle Scope

// grstyle is global within session until cleared
grstyle init
grstyle set color tableau    // Affects all graphs

grstyle clear                // Returns to base scheme

Comparing Schemes Side-by-Side

local schemes "plotplain 538 economist tableau"
local i = 1
foreach scheme of local schemes {
    set scheme `scheme'
    scatter price mpg, title("`scheme'") name(g`i', replace)
    local i = `i' + 1
}
graph combine g1 g2 g3 g4, cols(2)

---

Color Palette Selection

colorpalette Command

colorpalette tableau           // Display palette visually
colorpalette tableau, nograph  // Get colors programmatically
local colors `r(p)'

Palette Types

// Sequential (continuous data): Blues, Reds, viridis
// Diverging (data with midpoint): RdBu, PiYG
// Qualitative (categories): Set1, Paired, tableau, okabe

Colorblind-Safe Palettes

colorpalette tol bright     // Paul Tol
colorpalette okabe          // Okabe & Ito
colorpalette Set2           // ColorBrewer (many are safe)
colorpalette Dark2

Manual Color Specification

// RGB values in Stata
twoway (scatter price mpg if foreign==0, mcolor("0 114 178")) ///
       (scatter price mpg if foreign==1, mcolor("213 94 0"))

// Named locals for readability
local blue "0 114 178"
local orange "213 94 0"
scatter price mpg, mcolor("`blue'")

Grayscale

grstyle set color gs2 gs6 gs10 gs14

// Or manually
twoway (scatter price mpg if rep78<=2, mcolor(gs2) msymbol(O)) ///
       (scatter price mpg if rep78==3, mcolor(gs8) msymbol(D)) ///
       (scatter price mpg if rep78>=4, mcolor(gs14) msymbol(T))

Dynamic Palette from Data

levelsof rep78, local(levels)
local n : word count `levels'
colorpalette viridis, n(`n') nograph
local colors `r(p)'

local i = 1
local plots ""
foreach lev of local levels {
    local color : word `i' of `colors'
    local plots `plots' (scatter price mpg if rep78==`lev', mcolor("`color'"))
    local i = `i' + 1
}
twoway `plots', legend(order(1 "Poor" 2 "Fair" 3 "Average" 4 "Good" 5 "Excellent"))

---

Font and Size Customization

Presentation vs Publication Sizes

// Publication (journal article)
grstyle clear
grstyle init
grstyle set size 9pt
grstyle set size title: 11pt
grstyle set size axis_title: 10pt

// Presentation (slides)
grstyle clear
grstyle init
grstyle set size 14pt
grstyle set size title: 18pt
grstyle set size axis_title: 16pt

Size Hierarchy in a Single Graph

scatter price mpg, ///
    title("Main Title", size(large) color(black)) ///
    subtitle("Context", size(medium) color(gs4)) ///
    ytitle("Y-Axis", size(medsmall)) ///
    xtitle("X-Axis", size(medsmall)) ///
    ylabel(, labsize(small)) ///
    xlabel(, labsize(small)) ///
    note("Source info", size(vsmall) color(gs8))

Legend Customization

twoway (scatter price mpg if foreign==0) ///
       (scatter price mpg if foreign==1), ///
    legend(order(1 "Domestic" 2 "Foreign") ///
           size(small) region(lwidth(thin)) ///
           position(6) rows(1))

---

Creating Custom Schemes

Method 1: Scheme File

Scheme files are .scheme text files in ~/ado/plus/s/. They inherit from existing schemes:

#include s2color

// Background
color background white
color plotregion white

// Lines
color p1line "0 82 155"
color p2line "232 119 34"
color p3line "0 150 136"
linewidth p medthick

// Markers
symbol p1 O
symbol p2 D
symbol p3 T
color p1markfill "0 82 155"
color p2markfill "232 119 34"

// Text sizes
gsize heading medlarge
gsize axis_title medsmall
gsize tick_label small

// Grid
color grid gs12
linestyle grid dot
linewidth grid thin

// Install: place myscheme.scheme in ado/plus/s/
graph query, schemes   // Verify it appears
set scheme myscheme

Method 2: grstyle-Generated Scheme

grstyle clear
grstyle init myscheme, replace
grstyle set plain, box
grstyle set color navy maroon forest_green
grstyle set symbol O D T
grstyle set lpattern solid dash dot

// Creates myscheme.scheme file
set scheme myscheme

Project Template (Reusable Do-File)

// graph_setup.do
program define setup_project_graphs
    grstyle clear
    set scheme plotplain
    grstyle init
    grstyle set color "0 75 135" "242 101 34" "0 128 128"
    grstyle set plain, box horizontal
    grstyle set symbol O D T S
    grstyle set size 10pt
    grstyle set size axis_title: 11pt
    grstyle set size title: 12pt
end

// In analysis do-files:
do graph_setup.do
setup_project_graphs

---

Publication Requirements

Journal-Style Figure

set scheme plotplain

scatter price mpg, ///
    title("Figure 1. Price and Fuel Efficiency", ///
          position(11) justification(left) size(medium)) ///
    ytitle("Price (USD)", size(small)) ///
    xtitle("Miles per Gallon", size(small)) ///
    ylabel(, labsize(vsmall) angle(0)) ///
    xlabel(, labsize(vsmall)) ///
    graphregion(color(white) margin(medium)) ///
    plotregion(lcolor(black) margin(zero)) ///
    note("Notes: Sample of 74 automobiles from 1978.", ///
         size(vsmall) span)

graph export "figure1.pdf", replace

Multi-Panel Figure

set scheme plotplain

scatter price mpg, title("Panel A.", position(11) size(medsmall)) ///
    name(panelA, replace) nodraw

histogram price, title("Panel B.", position(11) size(medsmall)) ///
    name(panelB, replace) nodraw

graph combine panelA panelB, cols(2) ///
    title("Figure 2.", position(11) size(medium)) ///
    graphregion(color(white)) ///
    note("Notes: 1978 automobile data (N=74).", size(vsmall) span)

graph export "figure2.pdf", replace

Export Formats

graph export "figure.pdf", replace                // PDF vector (best for LaTeX)
graph export "figure.eps", replace fontface(Times) // EPS vector
graph export "figure.png", width(3000) replace     // PNG at ~300 DPI

Size Specifications

// Full page width (6.5 inches typical)
scatter price mpg, xsize(6.5) ysize(4) graphregion(color(white))

// Two-column format (3.25 inches typical)
scatter price mpg, xsize(3.25) ysize(3) graphregion(color(white))

Common Journal Styles

// NBER Working Paper
set scheme plotplain
grstyle init
grstyle set color navy maroon forest_green
grstyle set size 10pt

// Econometrica (very minimal)
set scheme plotplain
grstyle init
grstyle set plain
grstyle set size 9pt

// PLOS (colorblind-safe required)
set scheme plotplain
grstyle init
grstyle set color "0 114 178" "213 94 0" "0 158 115"

Grayscale for Print

set scheme plotplain
grstyle init
grstyle set color gs2 gs6 gs10 gs14
grstyle set lpattern solid dash dot dash_dot

---

Quick Reference

Common Schemes

Quick Setup Templates

// Academic Paper
set scheme plotplain
grstyle init
grstyle set plain, box
grstyle set color navy maroon forest_green
grstyle set size 10pt

// Presentation
set scheme 538
grstyle init
grstyle set size 14pt

// Colorblind-Safe
set scheme plotplainblind
grstyle init
grstyle set plain, box horizontal

// Grayscale
set scheme s1mono
grstyle init
grstyle set color gs2 gs6 gs10 gs14

Troubleshooting

// Scheme not found
graph query, schemes       // List available
ssc install schemepack     // Install package
findfile xyz.scheme        // Check file exists

// grstyle not taking effect
grstyle clear              // Clear first
grstyle init               // Then reinitialize
grstyle set color tableau

// Export looks different
graph export "fig.pdf", replace           // Use vector format
graph export "fig.png", width(3000) replace // Or high-res raster

// Scheme reverts on restart
set scheme plotplain, permanently

Preserving and Restoring

local savedscheme "`c(scheme)'"
set scheme economist
scatter price mpg
set scheme `savedscheme'

Cleanup at End of Do-File

grstyle clear
graph drop _all
set scheme s2color  // Return to default

Winsorizing in Stata: winsor and winsor2

Overview

Winsorizing caps extreme values at specified percentiles rather than removing observations. Two packages are available:

• winsor - Basic command by Nicholas Cox (single variable at a time)
• winsor2 - Enhanced version by Yujun Lian (multiple variables, by-group, replace option)

Winsorizing vs. trimming: Winsorizing replaces extremes with percentile values (preserves N). Trimming deletes extremes (reduces N).

---

Installation

ssc install winsor
ssc install winsor2
help winsor
help winsor2

---

winsor (Basic Command)

* Syntax
winsor varname [if] [in], gen(newvar) p(#)

* Options:
*   p(#)   - Fraction for each tail (0.01 = 1%/99%, 0.05 = 5%/95%)
*   gen()  - Name for new variable (required)
*   h(#)   - High-tail fraction (asymmetric)
*   l(#)   - Low-tail fraction (asymmetric)

Examples

sysuse auto, clear

* Symmetric: 1% each tail
winsor price, gen(price_w1) p(0.01)

* Symmetric: 5% each tail
winsor price, gen(price_w5) p(0.05)

* Asymmetric cutoffs
winsor price, gen(price_asym) l(0.01) h(0.05)

* Only upper tail
winsor price, gen(price_upper) l(0) h(0.05)

* Only lower tail
winsor price, gen(price_lower) l(0.05) h(1)

Limitation: Only handles one variable at a time. Use a loop for multiple:

local varlist "price mpg weight length"
foreach var of local varlist {
    winsor `var', gen(`var'_w) p(0.05)
}

---

winsor2 (Enhanced Command)

* Syntax
winsor2 varlist [if] [in] [, options]

* Key options:
*   cuts(# #)    - Lower and upper percentile cutoffs (e.g., cuts(1 99))
*   suffix(str)  - Suffix for new variables (default: _w)
*   replace      - Overwrite original variables (destructive!)
*   trim         - Trim (delete) instead of winsorize
*   by(varlist)  - Winsorize within groups
*   label        - Add labels to new variables

Basic Usage

sysuse auto, clear

* Single variable
winsor2 price, cuts(1 99) suffix(_w)

* Multiple variables at once
winsor2 price mpg weight, cuts(1 99) suffix(_wins)

* Variable range
winsor2 price-trunk, cuts(5 95) suffix(_w)

* Replace original (use with caution)
preserve
winsor2 price mpg weight, cuts(1 99) replace
restore

By-Group Winsorizing

When groups have different distributions, winsorize within groups:

sysuse auto, clear

* By a single grouping variable
winsor2 price, cuts(5 95) suffix(_bygroup) by(foreign)

* Panel data: by firm
webuse grunfeld, clear
winsor2 invest, cuts(1 99) suffix(_w) by(company)

* By industry-year (common in finance research)
* egen industry_year = group(industry year)
* winsor2 roa roe leverage, cuts(1 99) suffix(_w) by(industry_year)

Trimming with winsor2

* trim option deletes observations instead of capping
winsor2 price, cuts(1 99) trim
* WARNING: This drops observations, reducing N

---

Replace vs. Suffix

* suffix (safer -- keeps original)
winsor2 price mpg, cuts(1 99) suffix(_w)
* Creates price_w, mpg_w; originals unchanged

* replace (destructive -- overwrites original)
preserve
winsor2 price mpg, cuts(1 99) replace
* Original values are gone
restore

Best practice: Use suffix() and keep originals for robustness checks. If you must use replace, wrap in preserve/restore or back up data first.

---

Common Percentile Choices

* Common specifications
winsor2 price, cuts(1 99) suffix(_w1)       // Conservative
winsor2 price, cuts(2.5 97.5) suffix(_w2p5) // Moderate
winsor2 price, cuts(5 95) suffix(_w5)       // Aggressive

Rule of thumb: At 1% with N=100, you are capping ~1 observation per tail. With N < 50, winsorizing removes too much information -- use rreg (robust regression) or qreg (quantile regression) instead.

---

Quality Control After Winsorizing

winsor2 roa, cuts(1 99) suffix(_w)

* Check how many observations changed
gen changed = (roa != roa_w)
tab changed
* Should be ~2% (1% each tail)

* Verify percentile boundaries
_pctile roa, p(1 99)
summarize roa_w
* Min should equal 1st percentile, max should equal 99th

* Visual check
scatter roa_w roa
* Should see diagonal line with horizontal flats at tails

---

Workflow: Order of Operations

* CORRECT order:
use firm_data, clear

* 1. Construct variables
gen roa = net_income / total_assets

* 2. Apply sample restrictions
drop if total_assets < 10
keep if year >= 2000

* 3. Winsorize within final sample
winsor2 roa leverage, cuts(1 99) suffix(_w)

* 4. Run analysis
regress y roa_w leverage_w

Important: Winsorize AFTER sample restrictions so percentiles are based on the observations actually used in analysis.

---

Robustness: Comparing Outlier Treatments

sysuse auto, clear
estimates clear

* No treatment
regress mpg price weight
estimates store m1

* Winsorize 1%
winsor2 price weight, cuts(1 99) suffix(_w1)
regress mpg price_w1 weight_w1
estimates store m2

* Winsorize 5%
winsor2 price weight, cuts(5 95) suffix(_w5)
regress mpg price_w5 weight_w5
estimates store m3

* Log transformation
gen log_price = log(price)
gen log_weight = log(weight)
regress mpg log_price log_weight
estimates store m4

* Robust regression
rreg mpg price weight
estimates store m5

esttab m1 m2 m3 m4 m5, ///
    mtitles("Original" "Wins1%" "Wins5%" "Log" "Robust") ///
    b(3) se(3) r2 star(* 0.10 ** 0.05 *** 0.01)

---

Quick Reference

* Install
ssc install winsor
ssc install winsor2

* winsor: single variable, fraction-based
winsor price, gen(price_w) p(0.01)           // 1%/99%
winsor price, gen(price_w) p(0.05)           // 5%/95%
winsor price, gen(price_w) l(0.01) h(0.05)   // Asymmetric

* winsor2: multiple variables, percentile-based
winsor2 price mpg weight, cuts(1 99) suffix(_w)
winsor2 price mpg weight, cuts(5 95) replace
winsor2 roa roe, cuts(1 99) suffix(_w) by(industry)
winsor2 price, cuts(1 99) trim               // Trim instead

* Related commands
help _pctile   // Calculate percentiles manually
help rreg      // Robust regression (alternative to winsorizing)
help qreg      // Quantile regression (alternative)

When NOT to Winsorize

• Extreme values are your research focus (e.g., fraud detection, tail risk)
• Very small samples (N < 50)
• Outliers are clearly data errors (fix or delete instead)
• Theory predicts heavy tails (e.g., income, firm size -- use log transform)
• Distribution is already well-behaved (check with summarize, detail)