from _api_doc_utils import *OLS
Ordinary least squares with robust covariance options
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 Implementation walkthrough
The linear-estimator stack shares a small set of native array and covariance primitives.
fit()converts NumPy inputs to owned dense arrays, validates dimensions and finiteness, and prepends an explicit column of ones.fit_weighted()additionally validates nonnegative weights and multiplies every design row and outcome by \(\sqrt{w_i}\). The original unweighted arrays and weights are retained after fitting.- The point estimate is obtained with
ndarray-linalgleast squares on the design itself. The code does not form \(X'X\) or invert it to obtain coefficients. The returned parameter vector is split into its first intercept entry and remaining slope entries only after the solve succeeds. - CountSketch assigns each original row independently to one of \(s\) buckets and gives it a Rademacher sign \(\pm1\). Signed rows of \(X\) and signed outcomes are accumulated into the same buckets, using the supplied seed or fixed fallback seed
0xBAD5EED. There is no \(1/\sqrt{s}\) normalization; a common rescaling would not change the least-squares minimizer. fit_sketch()solves only the \(s\times p\) sketched system for the point estimate, then stores the untouched full sample. Consequently,predict(), residuals, covariance, and all reported fit statistics are evaluated on the original observations at the approximate coefficients rather than on sketch residuals.summary()rebuilds the design, computes residuals, and dispatches to a shared covariance routine. That routine forms and inverts \(X'X\), explicitly loops over leverage values for HC2/HC3, or first maps observation moments \(x_i e_i\) through the bread to parameter-score rows for HAC and cluster aggregation.- The pairs bootstrap draws integer row indices with replacement, copies matching \(X\) and \(y\) rows, and reruns the exact least-squares path. It does not reproduce a previous sketch. A weighted fit likewise resamples the corresponding weight entries.
The split between QR-style least squares for coefficients and explicit cross-product inversion for inference is intentional but important: a rank-deficient fit can return a least-squares solution while covariance construction still fails or is unstable. CountSketch accelerates only the former.
4 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.
5 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.
6 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) |
7 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]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.
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 |
() |