westata

Stata-conformant econometrics with a statsmodels-style API.

Scope. westata reimplements a subset of Stata's estimation commands in Python. It aims to match Stata's numbers, and where it cannot it says so rather than rounding the gap away. Every command carries an evidence grade, and the commands that deviate are listed in Known deviations with the measured gap.

westata is a personal research library that answers one question: if I run this in Python, do I get the same numbers Stata would print?

• API follows statsmodels — Model.from_formula(...).fit(), params, bse, pvalues, conf_int(), cov_params(). Results are duck-typed, so existing statsmodels-style downstream code keeps working.
• Numerics follow Stata — small-sample corrections, weight semantics, factor-variable coding, degrees of freedom, and degrees-of-freedom splits are implemented per the corresponding Stata manual entry.
• Claims are auditable — every command carries an evidence level in tests/conformance/registry.yml. Nothing claims more than it can show.

Status

Tests	4167 passing
Commands covered	93 registry entries
Evidence	52 × E3, 26 × E2, 15 × E1
Independent sweep	386 comparison points against live Stata (retest/)

Verified against StataNow 19.5 (MP) through pystata; the captures are committed so the comparison replays without a licence.

Install

westata is on PyPI:

python -m venv .venv
.venv/Scripts/python.exe -m pip install westata

On Linux or macOS the second line is ./.venv/bin/python -m pip install westata.

Requires Python ≥ 3.10. numpy, pandas, scipy, statsmodels, patsy, pyfixest (HDFE), linearmodels (IV, panel) and pyyaml come in automatically — the package does not import without the last three.

arch is optional and only needed for GARCH, EGARCH and the Phillips–Perron test:

.venv/Scripts/python.exe -m pip install "westata[arch]"

Without it those commands raise an UnsupportedError naming the missing extra, and their tests skip rather than fail.

Working on the library itself

The released version is what pip install westata gives you. To work on the source, clone and install it in editable mode with the test dependencies:

git clone https://github.com/bigchill-z/westata.git
cd westata
python -m venv .venv
.venv/Scripts/python.exe -m pip install -e ".[test]"

python -m pytest -q then runs the 1455-test suite. It needs no Stata licence: the conformance captures are committed and replayed. Re-running the comparisons against live Stata (retest/, docs/collect_comparison.py) does need one.

Quick start

Every line below is executed as written against tests/data/auto.dta.

import pandas as pd
import westata as ws

df = pd.read_stata("tests/data/auto.dta")

# OLS with a factor variable and clustered standard errors (regress)
res = ws.OLS.from_formula("mpg ~ weight + C(foreign)", df).fit(
    cov_type="cluster", cov_kwds={"groups": df["foreign"]},
)
print(res.summary(style="stata"))          # Stata layout

Coefficient names keep the formula's spelling, so a factor level appears as C(foreign)[T.Foreign]. Use that exact string:

res.test("weight")                                  # single restriction
res.test("weight", "C(foreign)[T.Foreign]")         # joint test of both
res.test("weight = 0")                              # same, as an equation
res.lincom("weight + 1000*C(foreign)[T.Foreign]")   # linear combination

test takes one specification per argument; a comma inside a single string is not a separator.

Other models, on the same data:

# Quantile regression (qreg) — median
ws.QuantileRegression.from_formula("mpg ~ weight + C(foreign)", df).fit()

# Robust regression (rreg)
ws.RobustRegression.from_formula("mpg ~ weight", df).fit()

# Systems of equations (sureg): one list entry per equation
ws.Sureg.from_formula(["price ~ weight + mpg", "mpg ~ weight + foreign"], df).fit()
ws.Sureg.from_formula(
    ["price ~ weight + mpg", "mpg ~ weight + foreign"], df
).fit(cov_type="robust")                     # robust VCE

# Censored regression (tobit). `lower` is a from_formula argument, not a fit one.
# `rep78` has missing values, which are dropped.
ws.Tobit.from_formula("rep78 ~ weight + mpg", df, lower=0).fit()

A publication-style table:

base = ws.OLS.from_formula("mpg ~ weight", df).fit()
full = ws.OLS.from_formula("mpg ~ weight + C(foreign)", df).fit()
print(ws.etable([base, full], model_names=["Base", "Full"]))

Stata's estat vce, for the coefficient covariance or correlation matrix:

print(full.estat_vce())                      # covariance matrix
print(full.estat_vce(correlation=True))      # estat vce, correlation
full.estat_vce().to_frame()                  # full symmetric matrix as a DataFrame

Command coverage

Family	westata	Stata
Linear	OLS	regress (4 weight types, i. factors)
IV	IV2SLS / IVLIML / IVGMM	ivregress 2sls/liml/gmm
HDFE	FEOLS	reghdfe
HDFE	Areg	areg (requires absorb=, reports _cons)
Panel	PanelOLS / RandomEffects / BetweenOLS / FirstDifferenceOLS	xtreg fe/re/be/fd
Binary	Logit / Probit	logit / probit
Count	Poisson / NegativeBinomial	poisson / nbreg
Ordered	OrderedLogit / OrderedProbit	ologit / oprobit
Censored	Tobit (self-implemented)	tobit
Time series	ARIMA / VAR / VECM / GARCH / unitroot	arima / var / vec / arch / dfuller
Postestimation	test / lincom / margins / predict	same names
Reporting	etable	esttab

Weight semantics

Following Stata [U] 20.24:

Type	Point estimate	Variance	Effective N
fweight	literal replication	conventional	sum(w)
aweight	normalized to sum(w) = N	conventional	N
pweight	same as aweight	robust (HC1)	N
iweight	not normalized	conventional	sum(w)

weight_type is mandatory whenever weights are supplied — westata never guesses, because guessing silently changes the VCE.

Evidence levels

Level	Meaning	Count
E3	Agrees with full-precision e(b)/e(V) from real Stata, on every invocation the independent sweep exercised	52
E2	Approximate, or exact only on a subset of invocations: real but bounded differences	26
E1	Materially differs on at least one measured invocation, or not covered by the sweep	9

These levels come from the independent 386-point sweep, not from the library's own tests. That distinction matters: the library's tests passed at 1247 while the sweep found 24 defects, because a test written from the same understanding as the code cannot detect a misunderstanding shared by both.

Measured relative agreement with StataNow 19.5:

Command	Parameters	Standard errors
regress	5.1e-15	6.8e-16
regress, vce(robust)	5.1e-15	4.6e-15
regress, cluster	5.1e-15	7.1e-15
regress, i.factor	3.7e-14	5.1e-16
ivregress 2sls	7.8e-15	7.5e-15
ivregress liml	3.7e-16	3.7e-16
xtreg, fe	1.8e-16	1.9e-16
xtreg, re	6.8e-16	1.0e-15
xtreg, be	3.6e-16	1.7e-16
xtreg, fd	not comparable (command removed)	—
areg	4.6e-16	1.8e-16
reghdfe (1 FE)	4.5e-16	3.3e-16
reghdfe (2 FE)	4.4e-16	1.9e-16
reghdfe (2 FE, cluster)	4.4e-16	2.0e-16
var	2.9e-15	—
poisson	3.9e-13	5.8e-14
logit	~1e-7	~1e-7
probit	1.5e-7	—
tobit	8.5e-8	6.4e-7
margins, dydx	1.6e-7	5.2e-6
ologit / oprobit	~2e-5	—
nbreg	1.7e-5	9.1e-7
arima	4.8e-5	8.2e-5
predict (OLS xb)	1.6e-15	—
predict (logit pr)	3.1e-11	—
predict (tobit e/ystar)	5.8e-10	—
etable vs esttab	0.0	—
vecrank (all lag orders)	0.0 (4 dp)	—
dfuller	2.3e-4	—
ivregress gmm	4.9e-16	6.9e-18
arch	1.4e-04	z values not comparable

The linear, panel-linear, multi-equation and complex-survey families agree to machine precision (1e-15). Non-linear maximum-likelihood models are where the differences live, and almost always in the second-order quantities rather than the coefficients: an exact point estimate with a standard error that is 5% off still gives a wrong z, a wrong p, and a wrong confidence interval.

What was measured wrong, and then fixed:

• ologit/oprobit e(V): a numpy view alias let exp() act on the cut points, and the slope-by-cut block was never transformed — e(V) off by 20x, all cut-point z and CI wrong.
• nl covariance: pinv truncation made standard errors wrong by 6–7 orders of magnitude, so every t statistic was inflated.
• xtreg, fe vce(robust) was mapped to HC1; Stata defines it as vce(cluster panelvar) (verified: the two are identical in Stata, max|V_robust − V_cluster| = 0). The constant's SE was 18x off.
• tobit se(/sigma) lacked the delta-method Jacobian; margins used a fixed numerical step; podiscrete VCEs silently accepted and ignored vce(robust).

What remains approximate is listed per command in docs/USING_WITH_STATA.md, which is organised around the question a user actually has: can I quote this number as Stata's?

Two commands could not be compared at all because StataNow 19 removed them: xtreg, fd (r(198)) and kpss (r(199)). Their westata counterparts are validated by other means.

Behaviours worth knowing before you rely on a result

Six properties that are easy to assume wrongly. Each is deliberate, and each would change a number if you got it backwards:

• df_resid and e(df_r) are different quantities. Under clustering, t tests use G − 1 while e(df_r) stays N − K, matching Stata. Conflating them silently gives wrong p-values.
• Frequency weights reproduce the expanded-data fit. That is Stata's stated invariant; treating fweight as an analytic weight inflates the standard errors substantially.
• margins differentiates the response, not the linear index. For a logit the marginal effect is mean(p(1-p))·b, not the coefficient itself.
• PanelOLS.from_formula ignores entity_effects. Passing it to fit() gives pooled OLS; construct the model directly to get fixed effects.
• ivregress reports z and Wald chi2, while OLS and the panel estimators use t and F.
• e() keys differ per command. There is no uniform set: regress stores no e(aic), tobit stores none either, while the ML families do. Assuming uniformity is how a key gets invented.

Layout

src/westata/
  data/            weight normalization, panel/time-series metadata
  formula/         patsy design matrices, C(x, base=)
  vce/             HC0-HC3, CR1, ssc
  estimators/      linear, iv, hdfe, panel, discrete, tobit, tseries
  results/         StataResults, E(RETURN), Stata/statsmodels summaries
  postestimation/  test, lincom, margins, predict
  report/          etable, LaTeX/Markdown renderers
  conformance/     evidence registry, tolerances, invariants
tests/conformance/ registry.yml + published-output fixtures
docs/              design spec and implementation plans

Documentation

• Complete API reference: docs/API.md (English) and docs/API.zh-CN.md (中文) — all 172 public names, each with the Stata command it corresponds to, a description, and a runnable example whose output is captured from a real run. Both files are generated from the code (docs/render_api_reference.py), so signatures and example outputs cannot drift from it.
• Comparing westata with Stata: docs/USING_WITH_STATA.md — which quantities are interchangeable and which are not, for the four E2 commands.
• Evidence registry: tests/conformance/registry.yml — per-command level, measured gap, and reasoning.
• Known gaps: docs/GAP_REGISTER.md — what is deliberately not implemented, and why.
• Independent sweep: retest/README.md — the 386-point check against live Stata, and how to re-run it. This is where the evidence levels come from.
• Design spec: docs/superpowers/specs/2026-10-22-westata-design.md

Known deviations

Nine commands match Stata less closely than the rest. Each figure is the worst case across an independent 386-point sweep against live Stata; 250 of those 386 points are at machine precision, so these are the outliers rather than the norm.

Command	Worst measured relative difference	Why
xtgls	3.4e+00	Only when the cross-equation Sigma is rank-deficient (more equations than periods). No unique generalised inverse exists; Stata warns here too. On well-posed systems this is machine precision.
teffects nnmatch	7.1e-02	Standard error only — the point estimate is exact (3.6e-15). Uses compiled Mata internally.
frontier	4.3e-02	distribution(exponential) only, where sigma_u is nearly unidentified on the test data. The likelihood at Stata's estimates matches exactly.
bsqreg	3.0e-02	Bootstrap VCE. Coefficients are exact (3.4e-16); a resampling estimator has not converged at any finite replicate count.
permute	1.2e-02	Stata's shuffle order for the permutation draws is not reproducible from its output.
arima	5.5e-03	Optimiser stopping point; our likelihood is higher than Stata's.
rreg	1.0e-03	Final biweight case weights.
xtlogit, pa	1.1e-03	GEE working-correlation handling.
heckman	1.9e-05	Two-step standard errors. The ML path added later is a separate, exact estimate.

One capability is absent rather than approximate: teffects nnmatch, biasadj() raises UnsupportedError.

Known gaps

Not implemented, and they raise instead of approximating: Mata, the graphics system, data management, SEM/GSEM, survival analysis, Bayesian estimation, mi, xtabond2, csdid, rdrobust, synth, general gmm, nlsur, multi-way clustering (CGM), vecrank's fixed-N exact critical values, and Stata's arbitrary r() statistics outside the [U] exp_list grammar.

docs/GAP_REGISTER.md lists every one with its reason. Those are boundary statements, not a to-do list.