import matplotlib
matplotlib.use("Agg") # headless, non-interactive backend
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels
import statsmodels.api as sm
rng = np.random.default_rng(42) # one generator, seeded, for the whole notebook
np.set_printoptions(precision=6, suppress=True)
print("numpy", np.__version__, "| pandas", pd.__version__, "| statsmodels", statsmodels.__version__)nb2.1 — OLS from Scratch with NumPy vs. statsmodels
Companion to Ch 2.1, “OLS in Matrix Form.”
The chapter rebuilt regression from the floor up and landed on a single line, \hat{\boldsymbol\beta} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}. It also made a string of promises that are easy to state and more convincing to check on a computer: that the residuals sum to zero and sit perpendicular to every regressor, that the hat matrix \mathbf{H} and residual maker \mathbf{M} are symmetric and idempotent, and that a trustworthy library like statsmodels returns exactly the coefficients our hand-built formula does.
This notebook does all of that. We start with Sam’s four-day toy dataset — small enough that the chapter inverted the matrices by hand and got \hat{\boldsymbol\beta} = (1,\,7/6) — reproduce that number from scratch, then scale the identical code up to a larger simulated return series. Along the way we will see one genuinely important numerical lesson: why you should solve the normal equations with np.linalg.solve rather than forming an explicit inverse, and what near-collinearity does to the whole machine.
Setup
One seeded generator for the whole notebook (the reproducibility rule from CONVENTIONS.md), and matplotlib forced onto its headless Agg backend so the notebook runs anywhere, with or without a screen.
Sam’s tiny dataset
Recall the four trading days from the chapter. On each day Sam records the market return x_i and his stock’s return y_i, both in percent:
| Day i | Market x_i | Stock y_i |
|---|---|---|
| 1 | -2 | -1 |
| 2 | 0 | 0 |
| 3 | 1 | 3 |
| 4 | 1 | 2 |
We build the design matrix \mathbf{X} by hand: a column of ones for the intercept, then the market-return column. The column of ones is the intercept’s “regressor” — a variable equal to 1 for every observation, so its coefficient \beta_0 adds the same baseline to every prediction.
x = np.array([-2.0, 0.0, 1.0, 1.0])
y = np.array([-1.0, 0.0, 3.0, 2.0])
N = x.shape[0]
ones = np.ones(N)
X = np.column_stack([ones, x]) # N x K, first column = intercept
K = X.shape[1]
print("X =\n", X)
print("\ny =", y)
print(f"\nN = {N} observations, K = {K} regressors (including intercept)")Building \mathbf{X}'\mathbf{X} and \mathbf{X}'\mathbf{y}
The whole estimator is assembled from two cross-products. In NumPy @ is matrix multiplication and .T is the transpose, so \mathbf{X}'\mathbf{X} is X.T @ X and \mathbf{X}'\mathbf{y} is X.T @ y. The chapter computed these by hand and got a diagonal \mathbf{X}'\mathbf{X} — a small gift from the market returns happening to sum to zero in this sample, which makes the two columns orthogonal. Let us confirm.
XtX = X.T @ X
Xty = X.T @ y
print("X'X =\n", XtX)
print("\nX'y =", Xty)
# Cross-check the hand arithmetic from the chapter:
# X'X = [[4, 0], [0, 6]], X'y = [4, 7]
assert np.allclose(XtX, [[4, 0], [0, 6]])
assert np.allclose(Xty, [4, 7])
print("\nMatches the chapter's hand computation.")Solving the normal equations: solve, not inv
The chapter’s boxed formula is \hat{\boldsymbol\beta} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}. Written that way it looks like you should form the inverse (\mathbf{X}'\mathbf{X})^{-1} and multiply. You can — and for this tiny diagonal matrix it is harmless — but it is the wrong habit. The formula is really shorthand for “solve the linear system \mathbf{X}'\mathbf{X}\,\hat{\boldsymbol\beta} = \mathbf{X}'\mathbf{y} for the unknown \hat{\boldsymbol\beta},” and that is what np.linalg.solve(XtX, Xty) does directly.
Why prefer solve? Computing an explicit inverse and then multiplying does strictly more arithmetic than solving the system once, and every extra floating-point operation is a chance to lose precision. solve uses a factorization (an LU decomposition with pivoting) tuned to the single job of solving \mathbf{A}\mathbf{z} = \mathbf{b}; it is faster and numerically more stable, and the gap widens dramatically when \mathbf{X}'\mathbf{X} is ill-conditioned (Section “Near-collinearity” below). The rule of thumb: if you ever write inv(A) @ b, you almost certainly wanted solve(A, b). We show both here only to prove they agree on well-behaved data.
# Preferred: solve the normal equations directly.
beta_hat = np.linalg.solve(XtX, Xty)
# Educational comparison only: the explicit-inverse route.
beta_via_inv = np.linalg.inv(XtX) @ Xty
print("beta_hat (solve) =", beta_hat)
print("beta_hat (inv) =", beta_via_inv)
print("\nbeta_0 = %.6f, beta_1 = %.6f (chapter: 1, 7/6 = %.6f)" %
(beta_hat[0], beta_hat[1], 7/6))
# Reproduce the chapter's exact answer (1, 7/6) to machine precision.
assert np.allclose(beta_hat, [1.0, 7/6]), "should match Sam's hand result"
assert np.allclose(beta_hat, beta_via_inv)
print("\nSam's beta recovered exactly: hat-beta = (1, 7/6).")Fitted values, residuals, and the orthogonality promises
With \hat{\boldsymbol\beta} in hand the fitted values are \hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol\beta} and the residuals are \hat{\mathbf{e}} = \mathbf{y} - \hat{\mathbf{y}}. The chapter promised two algebraic identities that hold exactly, by construction, whenever the model includes an intercept:
- the residuals sum to zero, \sum_i \hat e_i = 0;
- the residuals are orthogonal to every regressor, \mathbf{X}'\hat{\mathbf{e}} = \mathbf{0}.
These are not lucky outcomes — they are what “least squares” mechanically forces. Let us see them appear, and check the hand-computed residual vector (1/3,\,-1,\,5/6,\,-1/6) from the chapter.
y_hat = X @ beta_hat
e_hat = y - y_hat
print("y_hat =", y_hat)
print("e_hat =", e_hat)
print("\nchapter e_hat = [1/3, -1, 5/6, -1/6] =", np.array([1/3, -1, 5/6, -1/6]))
print("\nsum of residuals :", e_hat.sum()) # expect ~0
print("X' e_hat (both ~ 0) :", X.T @ e_hat) # expect [~0, ~0]
print("y_hat' e_hat (~ 0) :", y_hat @ e_hat) # fitted _|_ residual
assert np.allclose(e_hat, [1/3, -1, 5/6, -1/6])
assert abs(e_hat.sum()) < 1e-12
assert np.allclose(X.T @ e_hat, 0, atol=1e-12)
print("\nAll orthogonality identities hold to machine precision.")The hat matrix \mathbf{H} and the residual maker \mathbf{M}
The chapter peeled the matrix that acts on \mathbf{y} off of the fitted values and named it the hat matrix, \mathbf{H} = \mathbf{X}(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}', \qquad \hat{\mathbf{y}} = \mathbf{H}\mathbf{y}. Its complement is the residual maker \mathbf{M} = \mathbf{I} - \mathbf{H}, with \hat{\mathbf{e}} = \mathbf{M}\mathbf{y}. Here we do need a matrix that looks like (\mathbf{X}'\mathbf{X})^{-1} sitting inside a product, so we build the inner solve as \texttt{solve}(\mathbf{X}'\mathbf{X},\ \mathbf{X}') — i.e. solve \mathbf{X}'\mathbf{X}\,\mathbf{Z} = \mathbf{X}' for the K\times N matrix \mathbf{Z} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}' — then left-multiply by \mathbf{X}. Same solve-over-inv discipline as before.
# Z = (X'X)^{-1} X' via solve (no explicit inverse), then H = X Z.
Z = np.linalg.solve(XtX, X.T) # K x N
H = X @ Z # N x N hat matrix
M = np.eye(N) - H # residual maker
print("H =\n", H)
print("\nM = I - H =\n", M)
# H reproduces the fitted values and M the residuals.
assert np.allclose(H @ y, y_hat)
assert np.allclose(M @ y, e_hat)
print("\nH y == y_hat and M y == e_hat: confirmed.")Symmetric and idempotent — the fingerprints of a projection
A projection matrix has two telltale algebraic properties, and the chapter proved both for \mathbf{H} and \mathbf{M}:
- symmetric: \mathbf{H}' = \mathbf{H} (and \mathbf{M}' = \mathbf{M});
- idempotent: \mathbf{H}\mathbf{H} = \mathbf{H} (and \mathbf{M}\mathbf{M} = \mathbf{M}) — applying the projection twice does nothing new, because once you have dropped onto the plane you are already on it.
Two more identities fall out for free: \mathbf{H}\mathbf{M} = \mathbf{0} (the explained and unexplained pieces meet at a right angle) and \mathbf{M}\mathbf{X} = \mathbf{0} (the residual maker annihilates the regressors — hence its other name, the annihilator). We verify all of it numerically.
checks = {
"H symmetric (H' == H)" : np.allclose(H, H.T),
"M symmetric (M' == M)" : np.allclose(M, M.T),
"H idempotent (HH == H)" : np.allclose(H @ H, H),
"M idempotent (MM == M)" : np.allclose(M @ M, M),
"HM == 0" : np.allclose(H @ M, 0, atol=1e-12),
"MX == 0 (annihilator)" : np.allclose(M @ X, 0, atol=1e-12),
}
for name, ok in checks.items():
print(f" {'PASS' if ok else 'FAIL':4} {name}")
assert all(checks.values()), "a projection identity failed"
print("\nH and M behave exactly like the orthogonal projections they are.")R^2 and the leverages h_{ii}
Two scalars worth pulling out. The coefficient of determination R^2 is the fraction of the variation in \mathbf{y} the regression explains. Because the model has an intercept we measure variation around the mean: R^2 = 1 - \frac{\sum_i \hat e_i^2}{\sum_i (y_i - \bar y)^2} = 1 - \frac{\text{SSR}}{\text{TSS}}. The diagonal entries of the hat matrix, h_{ii} = H_{ii}, are the leverages: h_{ii} measures how much observation i’s own outcome y_i pulls on its own fitted value \hat y_i. A point with an unusual x has high leverage. They run between 0 and 1 and sum to K (the number of regressors), a fact we check. The leverages drive the robust standard errors of Chapter 2.4, so we meet them now.
SSR = e_hat @ e_hat
TSS = ((y - y.mean()) ** 2).sum()
R2 = 1.0 - SSR / TSS
h = np.diag(H) # leverages h_ii
print(f"SSR = {SSR:.6f}, TSS = {TSS:.6f}, R^2 = {R2:.6f}")
print("\nleverages h_ii =", h)
print(f"sum of leverages = {h.sum():.6f} (should equal K = {K})")
print(f"most influential day (highest leverage): day {int(np.argmax(h)) + 1}")
assert np.isclose(h.sum(), K) # trace(H) = K
assert np.all((h >= -1e-12) & (h <= 1 + 1e-12))
print("\nLeverages sum to K and lie in [0, 1], as projection theory requires.")Validation against statsmodels
The whole point of building it from scratch is to trust it — and the cleanest trust test is to hand the same data to a battle-tested library and demand agreement. statsmodels wants the design matrix with the intercept column already present (which we have), so we feed it X directly rather than calling sm.add_constant. We compare coefficients, fitted values, residuals, R^2, and leverages.
model = sm.OLS(y, X) # X already has the intercept column
result = model.fit()
print("from scratch beta_hat :", beta_hat)
print("statsmodels params :", result.params)
print("\nmax abs coefficient difference :", np.max(np.abs(beta_hat - result.params)))
# statsmodels exposes leverage via get_influence().hat_matrix_diag
infl = result.get_influence()
assert np.allclose(beta_hat, result.params, atol=1e-10), "coefficients must match to ~1e-10"
assert np.allclose(y_hat, result.fittedvalues, atol=1e-10)
assert np.allclose(e_hat, result.resid, atol=1e-10)
assert np.isclose(R2, result.rsquared, atol=1e-10)
assert np.allclose(h, infl.hat_matrix_diag, atol=1e-10)
print("\nEverything matches statsmodels to ~1e-10: coefficients, fitted values,")
print("residuals, R^2, and leverages. The from-scratch machine is correct.")Same code, bigger data: Sam’s beta on a simulated return series
Nothing above was special to four observations. To prove it, we simulate N = 500 days of a market return and a stock return generated with a true beta of 1.2 and an intercept (alpha) of 0.05, plus idiosyncratic noise. We reuse the identical solve-based estimator and once more check it against statsmodels. This is the chapter’s promise that “one derivation covers every linear regression” made concrete.
N_big = 500
alpha_true, beta_true = 0.05, 1.2
mkt = rng.normal(0.0, 1.0, size=N_big) # market return, percent
eps = rng.normal(0.0, 0.8, size=N_big) # idiosyncratic noise
ret = alpha_true + beta_true * mkt + eps # stock return
Xb = np.column_stack([np.ones(N_big), mkt])
yb = ret
beta_big = np.linalg.solve(Xb.T @ Xb, Xb.T @ yb) # same estimator, scaled up
res_big = sm.OLS(yb, Xb).fit()
print(f"true (alpha, beta) = ({alpha_true}, {beta_true})")
print(f"from-scratch (alpha, beta)= ({beta_big[0]:.4f}, {beta_big[1]:.4f})")
print(f"statsmodels (alpha, beta)= ({res_big.params[0]:.4f}, {res_big.params[1]:.4f})")
print(f"\nmax abs difference vs statsmodels: {np.max(np.abs(beta_big - res_big.params)):.2e}")
assert np.allclose(beta_big, res_big.params, atol=1e-10)
print("\nIdentical code, 500 rows, still matches statsmodels to ~1e-10.")A preview of trouble: near-collinearity blows up (\mathbf{X}'\mathbf{X})^{-1}
The chapter warned that the dangerous case in practice is not perfect collinearity (a coding mistake software catches) but near-collinearity: two regressors that are almost — but not quite — redundant. Then \mathbf{X}'\mathbf{X} is still invertible, but barely, and \hat{\boldsymbol\beta} becomes wildly sensitive to tiny changes in the data.
The diagnostic is the condition number of \mathbf{X}'\mathbf{X}: the ratio of its largest to smallest eigenvalue. A small condition number (near 1) means a well-behaved, roughly spherical bowl; a huge one means a long, nearly flat trough where the location of the minimum is barely pinned down. We build a third regressor that is the market column plus a tiny perturbation and watch the condition number explode — and with it the gap between solve and the explicit inverse.
def condition_and_fit(eps_scale):
# Add a near-duplicate of the market column at the given noise scale; report condition #.
x3 = mkt + rng.normal(0.0, eps_scale, size=N_big) # nearly == mkt when eps_scale tiny
Xc = np.column_stack([np.ones(N_big), mkt, x3])
XtXc = Xc.T @ Xc
cond = np.linalg.cond(XtXc)
try:
b_solve = np.linalg.solve(XtXc, Xc.T @ yb)
b_inv = np.linalg.inv(XtXc) @ (Xc.T @ yb)
gap = np.max(np.abs(b_solve - b_inv))
except np.linalg.LinAlgError:
# So ill-conditioned the factorization gives up: numerically singular.
gap = np.nan
return cond, gap
print(f"{'noise scale':>12} | {'cond(X^T X)':>14} | {'|solve - inv|':>14}")
print("-" * 48)
for s in [1.0, 1e-2, 1e-4, 1e-6, 1e-8]:
cond, gap = condition_and_fit(s)
gap_str = "singular!" if np.isnan(gap) else f"{gap:.3e}"
print(f"{s:>12.0e} | {cond:>14.3e} | {gap_str:>14}")
print("\nAs the third column approaches a copy of the market column, cond(X'X) explodes")
print("and solve vs. inv start to disagree -- the precision problem Chapter 2.4 quantifies.")Visualizing the condition number’s blow-up
A picture of the same table: as the near-duplicate column’s noise shrinks toward zero, the condition number of \mathbf{X}'\mathbf{X} climbs without bound on a log scale. This is the near-singular wall the chapter promised — the regression cousin of last week’s nearly singular covariance matrix.
scales = np.logspace(0, -10, 25) # 1 down to 1e-10
conds = np.array([np.linalg.cond((lambda s: (
np.column_stack([np.ones(N_big), mkt, mkt + rng.normal(0, s, N_big)]).T @
np.column_stack([np.ones(N_big), mkt, mkt + rng.normal(0, s, N_big)])
))(s)) for s in scales])
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.loglog(scales, conds, marker="o")
ax.invert_xaxis() # shrinking noise -> moving right
ax.set_xlabel("noise scale of the near-duplicate column (smaller = more collinear)")
ax.set_ylabel(r"condition number of $\mathbf{X}'\mathbf{X}$")
ax.set_title("Near-collinearity blows up the condition number")
ax.grid(True, which="both", alpha=0.3)
fig.tight_layout()
fig.savefig("nb21_condition_number.png", dpi=110)
plt.close(fig)
print("Saved figure: nb21_condition_number.png")
print(f"condition number ranges from {conds.min():.2e} to {conds.max():.2e}")Your Turn
You now have a complete, verified OLS engine in a handful of NumPy lines. Extend it:
Add a real regressor. Append a third genuine column to Sam’s bigger design
Xb— say a lagged market return,np.concatenate([[0.0], mkt[:-1]])— refit with both yoursolve-based estimator andstatsmodels, and confirm they still agree to1e-10. How do the leverages change once you have three regressors (remember they must still sum to K)?Break it on purpose. Build a design with a perfectly collinear column — for example the chapter’s suggestion x_{i3} = 1 + x_i, which is the intercept column plus the market column. Pass
Xbad = np.column_stack([ones, x, 1 + x])tonp.linalg.solve(Xbad.T @ Xbad, Xbad.T @ y)and watch it raiseLinAlgError: Singular matrix. Then find the nonzero vector \mathbf{v} with \mathbf{X}\mathbf{v} = \mathbf{0} (hint: it has entries 1, 1, -1) and verifyXbad @ vis the zero vector. This is perfect multicollinearity — the machine cannot choose how to split the explained variation, so it refuses.Leverage by hand. Using Sam’s original four-day
X, compute H_{11} directly as \mathbf{x}_1'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_1 with \mathbf{x}_1' = (1, -2) and confirm it equalsnp.diag(H)[0]from this notebook. Which day has the highest leverage, and does that match your intuition about which market return is most unusual?