Why I Built fereg: reghdfe and fixest Results in a Fraction of the Time

I run reghdfe almost every day. It is the standard tool for regressions with high-dimensional fixed effects in Stata, and together with ivreghdfe and ppmlhdfe it has shaped how a generation of empirical papers in economics and finance is written.

But my datasets kept growing, and reghdfe did not keep up.

On 10 million observations with individual, firm and year effects, one regression takes half a minute. A Poisson model with ppmlhdfe takes four. On some panels, reghdfe does not finish within fifteen minutes. When a paper needs dozens of specifications, robustness checks and subsamples, that adds up to days of waiting.

fixest in R and pyfixest in Python are much faster. But moving to them means switching languages, and the numbers do not quite match. Each package has its own rules for dropping singleton observations, counting the degrees of freedom absorbed by the fixed effects, and correcting the standard errors. The same model can give different samples and standard errors in Stata, R and Python, and a project that uses all three has to reconcile them.

So I built fereg: one estimation engine written in Rust, with Stata, R and Python front ends. Each front end reproduces the package its users already know, and in my benchmarks each one is faster than the package it replaces.

1. One Engine, Three Front Ends

Each front end follows the syntax, defaults, small-sample corrections and output of the package it mirrors:

Existing code needs at most a different command name:

* Stata
reghdfe y x1 x2, absorb(firm year) cluster(firm)
fereg   y x1 x2, absorb(firm year) cluster(firm)

# R
fixest::feols(y ~ x1 + x2 | firm + year, df, cluster = ~firm)
fereg::feols(y ~ x1 + x2 | firm + year, df, cluster = ~firm)

# Python
pyfixest.feols("y ~ x1 + x2 | firm + year", data=df, vcov={"CRV1": "firm"})
fereg.feols("y ~ x1 + x2 | firm + year", data=df, vcov={"CRV1": "firm"})

All three cover OLS, IV/2SLS, Poisson and logit with any number of fixed effects, interactions of fixed effects, varying slopes (such as firm-specific trends), weights, and iid, robust or multi-way clustered standard errors, plus residuals, predictions and the estimated fixed effects. The estimation sample, the degrees of freedom, the variance estimator and the stored results follow the reference package. The test suites check each front end against reghdfe, ivreghdfe, ppmlhdfe, Stata’s logit, fixest and pyfixest.

2. Why Fixed-Effects Regressions Are Slow

With millions of fixed-effect levels, the dummy variables cannot be formed. Instead, the outcome and each regressor are residualized on the fixed effects (the Frisch–Waugh–Lovell theorem), and the coefficients are estimated from the residuals. Each residualization is a linear system with one unknown per fixed-effect level: millions of unknowns, one system per variable. Poisson and logit models solve these systems again at every iteration.

reghdfe, fixest and pyfixest solve them by alternating projections: demean by individual, then by firm, then by year, and repeat until nothing changes. This converges in a few sweeps when the fixed effects are well connected, for example when many workers move between firms. It slows to a crawl when they are not, as in matched employer–employee data with little mobility. Every sweep is a full pass over the data.

3. How fereg Solves It

fereg starts from the linear algebra of the problem instead of from the algorithm. Four ideas do most of the work.

A. Work on cells, not observations

The system depends on the data only through the cells: the distinct combinations of fixed-effect levels that occur, their weights, and the sums of each variable within them. fereg sorts the observations by cell once. After that, every solver iteration is one pass over the cells, with all variables processed together. The observations are read again only at the end, to subtract the fitted fixed effects. With state and year effects on individual-level data, there are far fewer cells than observations.

B. Eliminate the largest fixed effect exactly

The part of the system that belongs to one set of fixed effects is diagonal, so it can be solved in closed form. fereg eliminates the fixed effect with the most levels, usually individuals or workers, and iterates only on the rest. With two fixed effects, what remains is a graph Laplacian, a matrix whose structure is well understood.

C. A solver that watches itself

The remaining system is solved by preconditioned conjugate gradient, for all variables at once. After 6 and 12 iterations, fereg measures how fast each variable converges and projects the remaining work. When iterating would be slow, a cost model decides whether to factor the system exactly (sparse Cholesky) or approximately, with the approximate Cholesky methods developed for graph Laplacians (Kyng and Sachdeva). The solve then needs only a few more steps. Run time therefore depends little on how well connected the data are, which is exactly where alternating projections struggle most.

D. Parallel, and deterministic

Every pass is parallel over observations or cells, not over variables, so a regression with one regressor still uses every core. Sums are taken over fixed chunks and combined in a fixed order, so the estimates are bit-for-bit identical whatever the number of threads: your results do not change when you move from a laptop to a server.

For Poisson and logit, only the weights change between iterations. The cell structure and orderings are built once, each iteration starts from the previous fixed effects, and a factorization is kept until refactoring it is cheaper than iterating.

4. Benchmarks

Seconds per regression at 10 million observations on 16 threads (AMD Ryzen 9 9950X), data in memory: the median of three timed runs after a warm-up (a single run when the warm-up takes over a minute). The data are the benchmark panels of fixest and pyfixest: individuals over 10 years, with firms assigned at random (“simple”) or cyclically (“difficult”, a fixed-effect structure that is hard for iterative solvers).

Green multiples: how many times faster fereg is, against the fereg version for the same language (reghdfe family and xhdfe vs fereg Stata, pyfixest vs fereg Python, fixest vs fereg R). “reghdfe family”: reghdfe for OLS, ivreghdfe for IV, ppmlhdfe for Poisson. pyfixest: the faster of its two solvers (MAP and LSMR). >900: did not finish within 15 minutes. —: no such estimator.

That is about 70 to 110 times faster than reghdfe, ivreghdfe and ppmlhdfe, 8 to 76 times faster than pyfixest, and 3 to 17 times faster than fixest on the simple panel. The gap is widest where it hurts most:

On the difficult panel with three fixed effects, reghdfe did not finish in 15 minutes and fixest took 36 seconds. fereg takes about a third of a second.

The table is a slice of a larger grid of 273 configurations: OLS with iid and clustered errors, IV, Poisson, logit, varying slopes, and worker–firm (AKM) designs with high and low mobility, from 10,000 to 10 million observations, on 16 threads and, up to 1 million observations, on a single thread. In every configuration, the fereg Stata commands are faster than every competitor that finished (reghdfe, ivreghdfe, ppmlhdfe, xhdfe, fixest and pyfixest), and the R and Python packages are faster than fixest and pyfixest. The Stata times include the transfer of the data from Stata to the plugin, which Stata does one value at a time. Full tables and the benchmark code are in the repository.

5. Installing fereg

fereg is open source under the MIT license. The Stata package includes binaries for Linux, macOS (Apple Silicon and Intel) and Windows:

net install fereg, from("https://raw.githubusercontent.com/leoliu0/fereg/gh-pages/stata") replace

fereg y x1 (x2 = z), absorb(firm year) vce(robust)
feppml trade dist, absorb(exporter#year importer#year) cluster(pair)
felogit default x1, absorb(firm year)

Python wheels cover CPython 3.10–3.14 on Linux, macOS and Windows; data can be pandas or polars:

pip install fereg --find-links https://github.com/leoliu0/fereg/releases/expanded_assets/v0.3.0

In R, the source package builds with Rust installed (cargo and rustc); binaries for macOS and Windows are on the release page:

install.packages("https://github.com/leoliu0/fereg/releases/download/v0.3.0/fereg_0.3.0.tar.gz",
                 repos = NULL, type = "source")

6. What Is Not There Yet

fereg covers the models that most applied work runs every day, not everything the reference packages offer. Not implemented yet: probit; HC2/HC3, Newey–West and Driscoll–Kraay standard errors; pyfixest’s wild bootstrap, randomization inference and difference-in-differences helpers; and fixest’s panel lag operators. The full list is in the README.

The solver also stops on its own criterion: when the estimated relative error of each residualized variable falls below 10⁻⁸. Results therefore agree with reghdfe, fixest and pyfixest to numerical precision, not bit for bit.

fereg would not exist without the packages it reproduces. reghdfe and ppmlhdfe (Sergio Correia and coauthors), fixest (Laurent Bergé) and pyfixest (the pyfixest team) define what the right answer is; fereg only tries to reach it faster.

The project is on GitHub: github.com/leoliu0/fereg. If it saves you time, starring the repo helps others find it.