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 zero-based category codes, one column per fixed-effect dimension.

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. Fixed-effect coefficients are not recovered or stored.

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 regressor and fixed-effect dimension, and optional nonnegative weights. It passes the fixed-effect code matrix and a contiguous weight slice to within::Solver using that crate’s default convergence parameters.
  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; nuisance fixed-effect coefficients are never reconstructed. This is why the class can report slope inference but deliberately cannot predict level outcomes for new observations.
  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.

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. Prediction is intentionally unavailable because the nuisance effects are discarded.

6 Python API

Constructor: cm.FixedEffectsOLS

Call fit(x, fe, y) or fit_weighted(x, fe, y, sample_weight). There is no predict() because the class is estimation-first and does not materialize fixed-effect coefficients. 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))
()
Code 
cls = cm.FixedEffectsOLS
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
FixedEffectsOLS()
bootstrap(self, /, n_bootstrap, seed=None)
fit(self, /, x, fe, y)
fit_weighted(self, /, x, fe, y, sample_weight)
summary(self, /, vcov='hc1', lags=None, clusters=None)
wald_test(self, /, r, q=None, vcov=None, lags=None, clusters=None)

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, 'residual_df': 257.0, 'absorbed_df_method': 'exact_two_way'}

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 ()
residual_df ()
absorbed_df_method ()