Fixed Effects OLS Example

This page shows how to estimate slope coefficients after partialling out one-way or multi-way categorical fixed effects with within.

1 Simulate A Two-Way Fixed Effects Problem

import numpy as np
from pprint import pprint

from crabbymetrics import FixedEffectsOLS, OLS

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(11)
n = 800
beta_true = np.array([1.2, -0.9])

worker = rng.integers(0, 40, size=n, dtype=np.uint32)
firm = rng.integers(0, 18, size=n, dtype=np.uint32)
fe = np.column_stack([worker, firm]).astype(np.uint32)

x = rng.normal(size=(n, 2))
worker_effect = rng.normal(scale=0.8, size=40)
firm_effect = rng.normal(scale=0.5, size=18)
y = (
    x @ beta_true
    + worker_effect[worker]
    + firm_effect[firm]
    + rng.normal(scale=0.1, size=n)
)

model = FixedEffectsOLS()
model.fit(x, fe, y)

print("true coef:", beta_true)
pprint(model.summary())
true coef: [ 1.2 -0.9]
{'absorbed_df': 57,
 'absorbed_df_method': 'exact_two_way',
 'absorption_iterations': [18, 18, 18],
 'absorption_residual_norms': [7.019116151938364e-09,
                               9.263346110384339e-09,
                               8.198622018204918e-09],
 'cluster_correction': True,
 'coef': array([ 1.2034, -0.9   ]),
 'coef_se': array([0.0037, 0.0035]),
 'converged': True,
 'counted_fe_df': 57,
 'fe_df': 'full',
 'residual_df': 741.0,
 'ssc': 'cluster',
 'vcov': array([[ 0., -0.],
       [-0.,  0.]]),
 'vcov_type': 'hc1'}

2 Compare To Explicit Dummy OLS

FixedEffectsOLS is doing the Frisch-Waugh-Lovell version of the same regression. A useful check is to compare it against an OLS fit with worker and firm dummies added directly to the design matrix.

def dummy_block(levels):
    n_levels = int(levels.max()) + 1
    block = np.zeros((levels.shape[0], max(n_levels - 1, 0)))
    for level in range(1, n_levels):
        block[:, level - 1] = (levels == level).astype(float)
    return block


dummy_design = np.column_stack([x, dummy_block(worker), dummy_block(firm)])

baseline = OLS()
baseline.fit(dummy_design, y)
baseline_summary = baseline.summary()

print("FixedEffectsOLS slopes:", model.summary()["coef"])
print("Dummy OLS slopes:      ", baseline_summary["coef"][: x.shape[1]])
print("Absolute difference:   ", np.abs(model.summary()["coef"] - baseline_summary["coef"][: x.shape[1]]))
FixedEffectsOLS slopes: [ 1.2034 -0.9   ]
Dummy OLS slopes:       [ 1.2034 -0.9   ]
Absolute difference:    [0. 0.]

3 Robust Covariance Options

Once the fixed effects are partialled out, summary() exposes the same covariance choices as the other linear estimators.

worker_cluster = worker.astype(np.int64)

hc1 = model.summary()
cluster = model.summary(vcov="cluster", clusters=worker_cluster)
hac = model.summary(vcov="newey_west", lags=3)

print("HC1 SE:", np.round(hc1["coef_se"], 4))
print("cluster SE:", np.round(cluster["coef_se"], 4))
print("Newey-West SE:", np.round(hac["coef_se"], 4))
HC1 SE: [0.0037 0.0035]
cluster SE: [0.0035 0.0041]
Newey-West SE: [0.0036 0.0036]

4 Notes

  • fe must be a 2D uint32 matrix with one column per factor.
  • Category codes should be zero-based and contiguous within each factor.
  • The estimator absorbs the constant. predict(x, fe) includes identified one/two-way FE levels and rejects unseen levels or combinations crossing disconnected components. Three-or-more-way level prediction is not supported.
  • Cluster labels passed to summary(vcov="cluster", clusters=...) should be a separate 1D int64 array with one label per observation.

5 Solver controls and covariance conventions

FixedEffectsOLS(tolerance=1e-8, max_iterations=1000) exposes the absorption controls. An (N, 0) design fits FE only. Summary reports per-right-hand-side iterations and residual norms. Invalid settings and failed absorption raise errors.

For clustered inference, summary() and wald_test() accept:

  • fe_df="full" (unchanged default), "non_nested" (exclude cluster-nested FE, counting an intercept), or "none" (slope count only). Multiway FE rank remains conservative beyond two FE dimensions.
  • ssc="cluster": (N-1)/(N-K); "hc1": N/(N-K); "none": no observation count correction. Here K includes slopes and the selected FE count.
  • cluster_correction=True: multiply by G/(G-1). Set false to omit this factor.

ssc and cluster_correction apply only to cluster covariance. These are explicit conventions, not claims that packages share a universal default. On the fatalities lecture, full FE gives SE 0.385787; fe_df="non_nested" reproduces fixest’s 0.357078; fe_df="none", ssc="hc1", cluster_correction=False reproduces plm HC1’s 0.350150. The slope is unchanged at -0.639980.

summary() retains absorbed_df (full model rank) and separately reports counted_fe_df, residual_df, and the selected correction policies.