OLS Example

Weighted fits use analytic weights. Zero-weight rows are excluded from residual degrees of freedom, cluster counts, and the observation sequence used by HAC inference. Positive weights are not interpreted as replicated observations. Covariance uses the design’s QR factorization without explicitly forming its normal equations. A failed refit clears fitted state. Sketched fits report sketch_method, sketch_size, sketch_seed, original_nobs, and inference_approximate=True; their covariance does not include sketch randomness.

1 Fit A Basic Linear Model

import numpy as np
from pprint import pprint

from crabbymetrics import OLS

np.set_printoptions(precision=4, suppress=True)
rng = np.random.default_rng(0)
n = 500
k = 3
beta = np.array([1.5, -2.0, 0.5])
intercept = 0.7

x = rng.normal(size=(n, k))
y = intercept + x @ beta + rng.normal(scale=0.5, size=n)

model = OLS()
model.fit(x, y)

print("true intercept:", intercept)
print("true coef:", beta)
pprint(model.summary())
true intercept: 0.7
true coef: [ 1.5 -2.   0.5]
{'anytime_valid': False,
 'coef': array([ 1.5206, -1.9933,  0.5314]),
 'coef_se': array([0.0224, 0.0249, 0.0223]),
 'inference_approximate': False,
 'intercept': 0.6686882817624615,
 'intercept_se': 0.023367266952062427,
 'nobs': 500,
 'original_nobs': 500,
 'sketch_method': None,
 'sketch_seed': None,
 'sketch_size': None,
 'vcov': array([[ 0.0005,  0.    , -0.0001, -0.    ],
       [ 0.    ,  0.0005, -0.    , -0.    ],
       [-0.0001, -0.    ,  0.0006, -0.    ],
       [-0.    , -0.    , -0.    ,  0.0005]]),
 'vcov_type': 'hc1',
 'weight_type': 'analytic'}

2 Robust Covariance Options

OLS.summary() now shares the same covariance interface used by the other linear estimators:

  • vcov="vanilla" for homoskedastic standard errors
  • vcov="hc1" for heteroskedasticity-robust Eicker-Huber-White standard errors
  • vcov="newey_west" with a lag choice for HAC inference
  • vcov="cluster" with one-way cluster labels
clusters = np.repeat(np.arange(25, dtype=np.int64), n // 25)

vanilla = model.summary(vcov="vanilla")
hac = model.summary(vcov="newey_west", lags=4)
cluster = model.summary(vcov="cluster", clusters=clusters)

print("vanilla SE:", np.round(vanilla["coef_se"], 4))
print("Newey-West SE:", np.round(hac["coef_se"], 4))
print("cluster SE:", np.round(cluster["coef_se"], 4))
vanilla SE: [0.0236 0.0242 0.0227]
Newey-West SE: [0.0229 0.0245 0.02  ]
cluster SE: [0.0266 0.0214 0.0177]