OLS

Ordinary least squares with robust covariance options

from _api_doc_utils import *

1 Where it fits

Group: Regression

OLS estimates the linear projection

\[ y_i = \alpha + x_i'\beta + u_i \]

by least squares. The class stores the fitted intercept and coefficient vector and exposes the common linear-model covariance surface: classical, HC1, Newey-West, and cluster-robust standard errors.

2 Criterion and computation

Write \(\tilde x_i=(1,x_i')'\) and let \(\theta=(\alpha,\beta')'\). The ordinary fit solves

\[ \hat\theta = \arg\min_\theta \sum_{i=1}^n (y_i-\tilde x_i'\theta)^2. \]

The weighted method instead minimizes \(\sum_i w_i(y_i-\tilde x_i'\theta)^2\) by applying least squares to \(\sqrt{w_i}\tilde x_i\) and \(\sqrt{w_i}y_i\). Both paths use a rank-revealing least-squares solve rather than explicitly forming \((X'X)^{-1}\) for the point estimate.

The sketched method draws a CountSketch matrix \(S\in\mathbb R^{s\times n}\) and solves

\[ \hat\theta_S=\arg\min_\theta\|S(y-\tilde X\theta)\|_2^2. \]

This is an approximate point estimator. The fitted object still stores the full original design and computes residuals and covariance estimates on all \(n\) observations at \(\hat\theta_S\).

3 Inference

Let \(e_i=y_i-\tilde x_i'\hat\theta\), \(H=\tilde X(\tilde X'\tilde X)^{-1}\tilde X'\), and \(d_f=n-p\), where \(p\) includes the intercept. The implemented covariance choices are:

\[ \widehat V_{\mathrm{vanilla}} = \frac{e'e}{d_f}(\tilde X'\tilde X)^{-1}, \]

and

\[ \widehat V_{\mathrm{HC}k} = (\tilde X'\tilde X)^{-1} \left(\sum_i a_{ki}e_i^2\tilde x_i\tilde x_i'\right) (\tilde X'\tilde X)^{-1}, \]

with \(a_{0i}=1\), \(a_{1i}=n/d_f\), \(a_{2i}=(1-h_{ii})^{-1}\), and \(a_{3i}=(1-h_{ii})^{-2}\). Newey-West uses the Bartlett-weighted autocovariances of the parameter scores \((\tilde X'\tilde X)^{-1}\tilde x_i e_i\) and defaults to \(\max\{1,\lfloor4(n/100)^{2/9}\rfloor\}\) lags. Cluster covariance sums those scores within clusters and applies

\[ \frac{G}{G-1}\frac{n-1}{d_f}. \]

Wald tests use the selected covariance. The bootstrap is a pairs bootstrap over rows and refits the model in every draw.

With anytime-valid output enabled, the class converts each squared \(t\) statistic to the implemented mixture evidence value

\[ \log G_t = \frac12\log r+ \frac{\nu+1}{2} \left[ \log(1+t^2/\nu)-\log(1+rt^2/\nu) \right], \qquad r=\frac{g}{g+n},\quad \nu=n-p, \]

and reports \(\min\{1,\exp(-\log G_t)\}\). The confidence-sequence half-width is \(c_{n,\alpha,g}\,\mathrm{SE}\), where

\[ c_{n,\alpha,g}^2 = \nu\frac{1-q}{q-r}, \qquad q=(r\alpha^2)^{1/(\nu+1)}. \]

This sequential calculation reuses the requested fixed-sample covariance; it does not change the coefficient fit.

4 Performance and numerical behavior

For a dense \(n\times p\) design, fitting costs about \(O(np^2+p^3)\) and stores \(O(np+p^2)\) values. HC2/HC3 explicitly evaluate every leverage and therefore also cost \(O(np^2)\). Newey-West adds \(O(Lnp)\) score work for \(L\) lags; cluster aggregation is linear in \(n\) after the scores are formed. A CountSketch reduces the least-squares row dimension to \(s\), but full-data prediction, residual construction, and inference remain; it is most useful when the point-estimation solve dominates. Ill-conditioned or rank-deficient designs may produce unstable covariance inversions even when the least-squares point estimate succeeds.

5 Python API

Constructor: cm.OLS

Use fit(x, y) for the standard estimator, fit_weighted(x, y, sample_weight) for weighted least squares, and fit_sketch(x, y, sketch_size, seed=None) when a randomized row sketch is acceptable for a large tall design. predict(x) returns fitted values. bootstrap(B, seed=None) returns bootstrap draws with the intercept in the first column.

print(inspect.signature(cm.OLS))

()

cls = cm.OLS
display(HTML(html_table(["Public method"], public_methods(cls))))

Public method
bootstrap(self, /, n_bootstrap, seed=None)
fit(self, /, x, y)
fit_sketch(self, /, x, y, sketch_size, seed=None)
fit_weighted(self, /, x, y, sample_weight)
predict(self, /, x)
summary(self, /, vcov='hc1', lags=None, clusters=None, anytime_valid=False, g=1.0, level=0.95)
wald_test(self, /, r, q=None, vcov=None, lags=None, clusters=None)

6 Minimal example

rng = np.random.default_rng(1)
x = rng.normal(size=(200, 3))
y = 0.5 + x @ np.array([1.0, -0.7, 0.25]) + rng.normal(scale=0.3, size=200)
model = cm.OLS()
model.fit(x, y)
clusters = np.repeat(np.arange(20, dtype=np.int64), 10)
print(model.summary(vcov='hc1')['coef'])
print(model.summary(vcov='cluster', clusters=clusters)['coef_se'])
print(model.predict(x[:3]))

[ 1.00664118 -0.70425283  0.224367  ]
[0.02367873 0.0222344  0.0228217 ]
[ 0.29427264 -1.39837835 -0.41709421]

7 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.

rng = np.random.default_rng(101)
x = rng.normal(size=(80, 3))
y = 0.5 + x @ np.array([1.0, -0.7, 0.25]) + rng.normal(scale=0.3, size=80)
model = cm.OLS()
model.fit(x, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))

summary() key	shape
intercept	()
coef	(3,)
intercept_se	()
coef_se	(3,)
vcov	(4, 4)
vcov_type	()
anytime_valid	()