FixedEffectsOLS

Within-estimator for high-dimensional fixed effects

1 Where it fits

Group: Regression

FixedEffectsOLS partials out one or more categorical fixed effects, then runs least squares on residualized variables:

M_F y = M_F X\beta + M_F u.

The fixed-effect matrix fe is a 2D uint32 array of category IDs, one column per fixed-effect dimension. IDs are compacted internally; they need not be contiguous.

2 Criterion and absorption

For factor incidence matrix D, the conceptual weighted problem is

\min_{\beta,\alpha} \sum_{i=1}^n w_i(y_i-x_i'\beta-d_i'\alpha)^2.

The class uses the Frisch-Waugh-Lovell representation. It iteratively demeans y and every column of X over each fixed-effect dimension using the within solver, producing M_Dy and M_DX, and then solves

\hat\beta = \arg\min_\beta \|W^{1/2}M_D(y-X\beta)\|_2^2.

The iterative absorption must converge for the outcome and every regressor; otherwise the fit raises an error. One- and two-way level effects are recovered for identified prediction. Their normalization is not itself an economic estimand.

3 Implementation walkthrough

Weights have analytic, not frequency, semantics. Zero-weight rows are removed before absorption and excluded from observation counts, fixed-effect level/rank counts, HAC lag ordering, and cluster counts. Summary and Wald inference use the same active sample. A failed refit clears prior fitted state.

This class has a narrow delegation boundary: the within Rust crate performs absorption, while crabbymetrics owns orchestration, the slope solve, degrees of freedom, covariance, and Python behavior.

  1. The wrapper validates aligned rows, at least one fixed-effect dimension (zero regressors are allowed), and optional nonnegative weights. It passes the fixed-effect code matrix and a contiguous weight slice to within::Solver using the constructor’s tolerance and max_iterations controls.
  2. Instead of invoking absorption separately, it creates one batch whose first right-hand side is y and whose remaining right-hand sides are the columns of X. One solver instance therefore uses the same encoded factors and numerical settings for every variable.
  3. The fit inspects the convergence flag for every right-hand side. If even one regressor has not converged, it raises an error containing the complete convergence and final-residual vectors. It never combines partially converged columns into a regression.
  4. The returned demeaned vectors are copied into y^*=M_Dy and X^*=M_DX. The package then applies optional square-root weights again in the slope least-squares solve. This matches the within solver’s weighted projection convention and the weighted Frisch-Waugh-Lovell objective.
  5. The slope vector and residualized arrays are retained. The fit recovers one- or two-way level effects and their connected-component information. Prediction rejects unseen levels and two-way combinations crossing disconnected components; three-or-more-way level prediction is unsupported.
  6. For degrees of freedom, one factor contributes its observed level count. Two factors are represented as a bipartite graph and a union-find pass computes L_1+L_2-C exactly. For three or more factors the code counts observed levels and subtracts J-1, labeling the result as conservative rather than claiming exact graph rank.
  7. Covariance uses X^* and residuals y^*-X^*\hat\beta in the shared linear sandwich code. The bootstrap resamples X, y, every fixed-effect code, and optional weights together, then repeats absorption from scratch.

The implementation is optimized for many levels and few reported slopes: it avoids dummy expansion but keeps a dense residualized copy of every regressor. Convergence quality belongs to the delegated absorber; inferential conventions and failure policy belong to this wrapper.

4 Inference

All OLS covariance options are applied to the residualized design and outcome: homoskedastic, HC0-HC3, Newey-West, and cluster robust. The residual degrees of freedom are

d_f=n-p-r_D,

where p is the number of reported slopes and r_D is the absorbed fixed-effect rank. The implementation computes r_D exactly for one factor as its observed level count. For two factors it uses

r_D=L_1+L_2-C,

where C is the number of connected components in the bipartite level graph. With three or more factors it uses the conservative approximation \sum_jL_j-(J-1). Standard errors can therefore be conservative in disconnected multiway designs. Wald tests use the selected covariance, and the pairs bootstrap resamples rows, including all fixed-effect identifiers, before re-absorption and refitting.

The default fe_df="full" preserves the previous full-FE correction. fe_df="none" excludes absorbed effects from the inference count; fe_df="non_nested" excludes effects nested in the supplied cluster partition and requires cluster covariance. Cluster inference also exposes ssc="cluster", "hc1", or "none", and cluster_correction=True/False. Use the same settings in summaries and Wald tests. The worked example and replication report explain reference-library conventions; changing them changes the statistical correction, not the fitted slopes.

5 Performance and numerical behavior

Absorption avoids an n\times\sum_jL_j dummy matrix. Each demeaning sweep is approximately O(n(p+1)J) for J fixed-effect dimensions, with total cost multiplied by the number of sweeps required for convergence. Memory is dominated by the residualized n\times p design. Poorly connected factor graphs can converge slowly. Covariance construction remains dense in the slope dimension, and bootstrap cost is the cost of full absorption and regression times the number of draws. Level prediction adds the recovered effects to the slope prediction when the requested combination is identified.

6 Python API

Constructor: cm.FixedEffectsOLS

Call fit(x, fe, y) or fit_weighted(x, fe, y, sample_weight). predict(x, fe) returns identified level predictions for one- and two-way effects. For an FE-only fit, pass an (n, 0) regressor array. summary() supports the same covariance options as the other linear estimators and reports absorbed_df, residual_df, and absorbed_df_method. Rank is exact for one- and two-way fixed effects; three or more dimensions use a labeled conservative count.

Code 
print(inspect.signature(cm.FixedEffectsOLS))
(tolerance=1e-08, max_iterations=1000)
Code 
cls = cm.FixedEffectsOLS
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
FixedEffectsOLS(tolerance=1e-08, max_iterations=1000)
bootstrap(self, /, n_bootstrap, seed=None)
fit(self, /, x, fe, y)
fit_weighted(self, /, x, fe, y, sample_weight)
predict(self, /, x, fe)
summary(self, /, vcov='hc1', lags=None, clusters=None, *, fe_df='full', ssc='cluster', cluster_correction=True)
wald_test(self, /, r, q=None, vcov=None, lags=None, clusters=None, *, fe_df='full', ssc='cluster', cluster_correction=True)

7 Minimal example

rng = np.random.default_rng(3)
n = 300
x = rng.normal(size=(n, 2))
worker = rng.integers(0, 30, size=n, dtype=np.uint32)
firm = rng.integers(0, 12, size=n, dtype=np.uint32)
fe = np.column_stack([worker, firm]).astype(np.uint32)
y = x @ np.array([0.8, -0.5]) + rng.normal(size=30)[worker] + rng.normal(size=12)[firm] + rng.normal(scale=0.2, size=n)
model = cm.FixedEffectsOLS()
model.fit(x, fe, y)
print(model.summary(vcov='cluster', clusters=worker.astype(np.int64)))
{'coef': array([ 0.78646351, -0.52522162]), 'coef_se': array([0.01354886, 0.01305264]), 'vcov': array([[1.83571547e-04, 4.64880268e-05],
       [4.64880268e-05, 1.70371338e-04]]), 'vcov_type': 'cluster', 'absorbed_df': 41, 'counted_fe_df': 41, 'fe_df': 'full', 'ssc': 'cluster', 'cluster_correction': True, 'residual_df': 257.0, 'absorbed_df_method': 'exact_two_way', 'converged': True, 'absorption_iterations': [22, 22, 21], 'absorption_residual_norms': [9.652144187311186e-09, 6.984514961255232e-09, 8.894950479535544e-09]}

8 summary() contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

Code 
rng = np.random.default_rng(103)
n = 100
x = rng.normal(size=(n, 2))
g = rng.integers(0, 10, size=n, dtype=np.uint32)
h = rng.integers(0, 5, size=n, dtype=np.uint32)
fe = np.column_stack([g, h]).astype(np.uint32)
y = x @ np.array([0.8, -0.5]) + rng.normal(size=10)[g] + rng.normal(size=n) * 0.2
model = cm.FixedEffectsOLS()
model.fit(x, fe, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))
summary() key shape
coef (2,)
coef_se (2,)
vcov (2, 2)
vcov_type ()
absorbed_df ()
counted_fe_df ()
fe_df ()
ssc ()
cluster_correction ()
residual_df ()
absorbed_df_method ()
converged ()
absorption_iterations (3,)
absorption_residual_norms (3,)