OLS Standard Errors and t-Statistics From Scratch in NumPy

Given a design matrix X (n × p, full column rank, n > p; include a column of ones yourself if you want an intercept) and targets y, fit ordinary least squares and return [beta, se, t]: the p coefficients, their standard errors under homoskedastic errors, and the t-statistics beta / se.

Estimate the error variance as RSS / (n − p). The fit is never exact in the tests.

Examples

Input X = [[1, 0], [1, 1], [1, 2], [1, 3]], y = [1, 3, 2, 5]

Output [[1.1, 1.1], [0.9721, 0.5196], [1.1316, 2.117]]

An intercept and a slope from four points.

Input X = [[1, 1], [1, 2], [1, 3], [1, 4], [1, 5]], y = [2.1, 3.9, 6.2, 7.8, 10.1]

Output [[0.05, 1.99], [0.1981, 0.0597], [0.2524, 33.3213]]

A strong slope: its t-statistic is large.

Input X = [[1], [1], [1]], y = [2, 4, 9]

Output [[5], [2.0817], [2.4019]]

Intercept only: beta is the mean, and se is the standard error of the mean.

Constraints

2 ≤ n ≤ 1000, 1 ≤ p ≤ 20, n > p
X has full column rank
Checked to a relative tolerance of 10⁻⁷

Hints

Hint 1

Get beta first. Solving the normal equations works; a QR factorisation of X is more accurate.

Hint 2

The residual variance is the residual sum of squares divided by n − p, not n.

Hint 3

The covariance of beta is sigma² (XᵀX)⁻¹. With X = QR, (XᵀX)⁻¹ = R⁻¹R⁻ᵀ, whose diagonal is the squared row norms of R⁻¹.

Where candidates go wrong

Dividing the residual sum of squares by n: the standard errors come out too small, worst when p is large relative to n.
Forming np.linalg.inv(X.T @ X): it squares the condition number of X, and with nearly collinear columns the standard errors lose most of their digits.
Taking the square root of the whole covariance matrix instead of its diagonal.

The idea

Factor X = QR with Q having orthonormal columns and R upper triangular. Then XᵀX = RᵀR, the normal equations become Rβ = Qᵀy, solved by back substitution, and the condition number involved is that of X, not its square.

The residuals y − Xβ̂ have n − p degrees of freedom, so σ̂² = RSS / (n − p) is unbiased for σ². The covariance of β̂ is σ²(XᵀX)⁻¹ = σ² R⁻¹R⁻ᵀ, so the variance of β̂ᵢ is σ̂² times the squared norm of row i of R⁻¹. The t-statistic divides each coefficient by its standard error.

Time O(np²) · Space O(np)

Reference solution

import numpy as np
def ols_inference(X, y):
    X = np.asarray(X, dtype=float)
    y = np.asarray(y, dtype=float)
    n, p = X.shape
    Q, Rm = np.linalg.qr(X)                       # X = QR, Rm is p x p upper triangular
    beta = np.linalg.solve(Rm, Q.T @ y)           # R beta = Q^T y
    resid = y - X @ beta
    sigma2 = resid @ resid / (n - p)              # p degrees of freedom used by the fit
    Rinv = np.linalg.solve(Rm, np.eye(p))
    se = np.sqrt(sigma2 * np.sum(Rinv ** 2, axis=1))   # diag of (X^T X)^-1 = row norms of R^-1, squared
    return [beta.tolist(), se.tolist(), (beta / se).tolist()]

Console⌘/Ctrl + Enter runs

Python runs in your browser with NumPy. The first run downloads the runtime, which takes a few seconds.

Why interviewers ask it

Every signal in quant research is judged by a regression and its t-statistic. Interviewers ask for the standard errors because that is where candidates who have only called a library run out of road.

The theory behind it is in the lesson Statistical inference in code: standard errors, permutation tests, the Kalman filter.

Keep going

Easy10 min

Regression coefficients with standard errors and t-statistics

Examples

Constraints

Hints

Where candidates go wrong

The idea

Reference solution

Follow-ups an interviewer asks

Why interviewers ask it

Keep going

The autocorrelation function of a series

Regression coefficients with standard errors and t-statistics

An exact permutation test for a difference in means

A one-dimensional Kalman filter