reghdfe in Python: regression with high-dimensional fixed effects

Many applied regressions need several sets of fixed effects at once: worker and firm effects in a wage equation, exporter and importer effects in trade data, firm and year effects in a corporate panel. Each set can have thousands of levels, too many to add as dummy variables. A high-dimensional fixed-effects estimator removes all of them by iterative demeaning and reports the slopes that the full dummy-variable regression would give, without ever building the dummies.

What reghdfe estimates

In a matched worker–firm panel, the wage of worker i in year t depends on a worker effect, the effect of the firm j(i, t) that employs the worker, a year effect and the regressors:

yᵢₜ = xᵢₜ′β + αᵢ + ψⱼ₍ᵢ,ₜ₎ + δₜ + eᵢₜ

By the Frisch–Waugh–Lovell theorem, β equals the slope from regressing y on x after both have been projected off all the dummy sets together. With one set, that projection is the familiar within transformation. With two or more there is no closed form, so the variables are demeaned by each set in turn until nothing changes (Guimarães and Portugal 2010; Gaure 2013).

The key assumption is that eᵢₜ is uncorrelated with the regressors once all the effects are removed. The effects absorb anything that is constant within a worker, within a firm or within a year, whether or not it is correlated with x. They do not absorb a confounder that changes over time within a worker–firm match.

Three details matter in practice.

• Identification from movers. Firm effects are separated from worker effects only through workers who change firms (Abowd, Kramarz and Margolis 1999). Firms linked by movers form connected sets, and within each set one effect is redundant. The number of redundant effects enters the degrees of freedom.
• Singletons. An observation that is alone in its level, such as a worker seen only once, is fitted exactly by its own effect. It carries no information about β, and keeping it overstates the degrees of freedom and can bias clustered standard errors downward. Such observations are dropped, repeatedly, until none remain.
• Clustering. Errors are often correlated within firms, or within both firms and years. Cluster-robust standard errors, one-way or two-way (Cameron, Gelbach and Miller 2011), allow for that. An absorbed effect nested within a cluster variable is not counted in the degrees of freedom.

In OpenEconometrics

oe.reghdfe absorbs any number of fixed-effect columns. The call below corresponds to the community Stata command reghdfe log_wage log_hours, absorb(worker firm year) vce(cluster firm):

oe.reghdfe(data=df, y="log_wage", x=["log_hours"], absorb=["worker", "firm", "year"], cluster="firm")

• No Intercept is reported; the constant is absorbed by the fixed effects.
• covariance="nonrobust" (the default), "robust" (HC1) or "cluster" with cluster= naming one or two columns. Passing cluster= alone selects the cluster covariance. Coefficient tests are Student t with G − 1 degrees of freedom under clustering (the smallest G with two columns).
• drop_singletons=True is the default; drop_singletons=False keeps singletons, like reghdfe's keepsingletons. tolerance and max_iterations control the demeaning.
• Weights may be "aweight", "fweight" or "pweight"; sampling weights need a robust or cluster covariance.
• extra["absorbed"] lists each dimension's levels, redundant effects and whether it is nested in a cluster. metrics include r_squared_within, df_absorbed and n_singletons_dropped.

Two related functions use the same absorption. oe.ivreghdfe fits instrumental-variables regressions (2SLS, LIML or GMM) with absorbed effects; see the instrumental variables page. oe.ppmlhdfe fits Poisson pseudo-maximum likelihood with absorbed effects, widely used for gravity equations of trade; see Poisson, negative binomial and PPML. For a single set of panel effects, oe.xtreg also offers random effects and Driscoll–Kraay errors. The linear-family guide gives the degrees-of-freedom rules and every error code.

Supported options

Options of oe.reghdfe

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates a matched panel of 3,000 workers and 150 firms over six years. Each worker enters in a random year, stays for a random number of years and changes firm with probability 0.2 each year, so some workers appear only once. Hours depend on the worker effect and the firm effect, so workers in high-wage firms work more. The true elasticity of the wage with respect to hours is 0.3, and the errors are independent.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(workers=3000, firms=150, years=6, seed=5):
    """Matched worker-firm panel: wages depend on worker, firm and year effects."""
    rng = np.random.default_rng(seed)
    worker_effect = rng.normal(scale=0.3, size=workers)
    firm_effect = rng.normal(scale=0.2, size=firms)
    rows = []
    for w in range(workers):
        start = rng.integers(0, years)
        length = rng.integers(1, years - start + 1)      # some workers appear only once
        firm = rng.integers(firms)
        for t in range(start, start + length):
            if t > start and rng.random() < 0.2:          # 20% of workers change firm each year
                firm = rng.integers(firms)
            rows.append((w, firm, 2015 + t))
    df = pd.DataFrame(rows, columns=["worker", "firm", "year"])
    theta = worker_effect[df["worker"]]
    psi = firm_effect[df["firm"]]
    trend = 0.03 * (df["year"] - 2015)
    df["log_hours"] = 0.8 * theta + 1.0 * psi + trend + rng.normal(scale=0.3, size=len(df))
    df["log_wage"] = (1 + 0.3 * df["log_hours"] + theta + psi + trend
                      + rng.normal(scale=0.15, size=len(df)))
    return df


if __name__ == "__main__":
    df = simulate()
    fit = oe.reghdfe(data=df, y="log_wage", x=["log_hours"],
                     absorb=["worker", "firm", "year"], cluster="firm")
    print(fit.summary())
    print()
    for d in fit.extra["absorbed"]:
        print(f"{d['column']:<7} levels {d['levels']:>5}  redundant {d['redundant']:>4}  "
              f"nested in cluster {d['nested']}")
    print()
    pooled = oe.ols(data=df, y="log_wage", x=["log_hours"], covariance="cluster", cluster="firm")
    print(f"{'absorbed effects':<22}{'slope':>8}{'s.e.':>8}{'N':>7}")
    print(f"{'none (OLS)':<22}{pooled.coefficients[1].estimate:>8.4f}"
          f"{pooled.coefficients[1].std_error:>8.4f}{pooled.nobs:>7}")
    for absorb in (["worker"], ["worker", "year"], ["worker", "firm", "year"]):
        r = oe.reghdfe(data=df, y="log_wage", x=["log_hours"], absorb=absorb, cluster="firm")
        print(f"{' + '.join(absorb):<22}{r.coefficients[0].estimate:>8.4f}"
              f"{r.coefficients[0].std_error:>8.4f}{r.nobs:>7}")

With openecon 0.3.20a1 from PyPI the script prints:

Linear regression with high-dimensional fixed effects — log_wage
Observations: 5485  |  Covariance: cluster  |  Confidence: 95%

Term       Estimate  Std. error       t     P>|stat|  CI lower  CI upper
---------  --------  ----------  ------  -----------  --------  --------
log_hours  0.327251  0.00912299  35.871  3.45621e-75  0.309223  0.345278

Excluded observations: 1239

r_squared: 0.940345  |  adjusted_r_squared: 0.91201  |  r_squared_within: 0.320204  |  adjusted_r_squared_within: 0.320021
rmse: 0.143709  |  df_model: 1  |  df_resid: 3718  |  df_absorbed: 1766
n_singletons_dropped: 1239
Model F test (slopes): F(1, 149) = 1286.73, p = 3.456e-75
Warning: Dropped 1239 singleton observation(s) alone in a level of an absorbed dimension (reghdfe default; drop_singletons=False keeps them).

worker  levels  1761  redundant    0  nested in cluster False
firm    levels   150  redundant  150  nested in cluster True
year    levels     6  redundant    1  nested in cluster False

absorbed effects         slope    s.e.      N
none (OLS)              0.8958  0.0137   6724
worker                  0.4390  0.0117   5485
worker + year           0.4242  0.0117   5485
worker + firm + year    0.3273  0.0091   5485

Reading the output

• Singletons are dropped first. 1,239 of the 6,724 rows were alone in a level: they are the 1,239 workers seen only once (3,000 workers less the 1,761 that remain). That leaves 5,485 observations; the OLS row uses all 6,724.
• The degrees of freedom follow the absorbed structure. The 1,761 remaining worker effects are all counted. The 150 firm effects are nested in the firm clusters, so none is counted. One of the six year effects is redundant. That gives df_absorbed = 1,766 and df_resid = 3,718, and the clustered tests use 149 degrees of freedom, as in F(1, 149).
• Each set of effects removes part of the bias. OLS gives 0.8958 for an elasticity of 0.3. Worker effects bring it to 0.4390, and year effects to 0.4242. The estimate stays too high until the firm effects are absorbed: a worker who moves to a higher-wage firm also works more hours there, so the change in hours within a worker carries the firm effect with it.
• This sample misses the true value. With all three sets the slope is 0.3273, with a 95% interval from 0.3092 to 0.3453, which does not contain 0.3; the estimate is about three standard errors too high. The model is correctly specified, and a dummy-variable regression on the same data gives the same slope, so this is sampling error in this draw, not bias or a convergence failure. Draws this far off are rare: in 40 further samples from the same design (seeds 100 to 139), the estimates averaged 0.3000 and every 95% interval contained 0.3.
• The fit is mostly the effects. The overall R² of 0.940 includes the absorbed effects; the within R² of 0.320 is the share of the remaining variation that hours explain.

Limitations

• Cluster-robust covariance for oe.reghdfe allows at most 2 clustering dimensions.
• Covariances are limited to nonrobust, HC1 and one- or two-way cluster; HC0, HC2, HC3 and HAC covariances are not offered.
• iweights are not accepted.
• From the third absorbed dimension on, redundant effects are counted with a pairwise rule that gives a lower bound, so the degrees of freedom and the inference are conservative.
• With fewer clusters than slopes, the model F test is reduced to the available rank and flagged in the result.
• Demeaning is iterative: a fit that does not reach tolerance within max_iterations sweeps raises absorption_nonconvergence.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Abowd, J. M., Kramarz, F. and Margolis, D. N. (1999). High wage workers and high wage firms. Econometrica, 67(2), 251–333. doi:10.1111/1468-0262.00020
• Cameron, A. C., Gelbach, J. B. and Miller, D. L. (2011). Robust inference with multiway clustering. Journal of Business & Economic Statistics, 29(2), 238–249. doi:10.1198/jbes.2010.07136
• Gaure, S. (2013). OLS with multiple high dimensional category variables. Computational Statistics & Data Analysis, 66, 8–18. doi:10.1016/j.csda.2013.03.024
• Guimarães, P. and Portugal, P. (2010). A simple feasible procedure to fit models with high-dimensional fixed effects. Stata Journal, 10(4), 628–649. doi:10.1177/1536867X1101000406