OLS from first principles · print(Purple)

Ordinary least squares is the model to reach for when a correlation is not enough: when the question needs a slope in units, or a coefficient on one variable holding another fixed. It is also short enough to write from scratch, which matters when the only libraries available are numpy and pandas.

The problem

Given $n$ observations of an outcome $y$ and $k$ regressors (including a column of ones for the intercept) stacked into an $n \times k$ matrix $X$ , find the coefficient vector $\beta$ that makes the fitted values $X\beta$ as close as possible to $y$ in the squared sense:

\hat\beta = \arg\min_\beta \; (y - X\beta)^\top (y - X\beta).

The one-line derivation

Expand the objective, differentiate with respect to $\beta$ , and set the gradient to zero:

\frac{\partial}{\partial \beta}\left[ y^\top y - 2\beta^\top X^\top y + \beta^\top X^\top X \beta \right] = -2X^\top y + 2X^\top X\beta = 0

which gives the normal equations and their solution:

X^\top X \, \hat\beta = X^\top y \qquad\Longrightarrow\qquad \hat\beta = (X^\top X)^{-1} X^\top y .

Geometrically, $\hat\beta$ is the vector for which the residual $e = y - X\hat\beta$ is orthogonal to every column of $X$ : $X^\top e = 0$ . The fitted values are the projection of $y$ onto the column space of $X$ .

For a single regressor with an intercept, the same algebra collapses to

\hat\beta_1 = \frac{\operatorname{Cov}(x, y)}{\operatorname{Var}(x)}, \qquad \hat\beta_0 = \bar y - \hat\beta_1 \bar x .

Standard errors

Under the classical assumptions (errors with mean zero, constant variance $\sigma^2$ , uncorrelated with each other and with $X$ ):

\operatorname{Var}(\hat\beta) = \sigma^2 (X^\top X)^{-1}, \qquad \hat\sigma^2 = \frac{e^\top e}{n - k}.

The standard error of each coefficient is the square root of the corresponding diagonal element, and $t_j = \hat\beta_j / \operatorname{se}(\hat\beta_j)$ . $R^2 = 1 - e^\top e / \sum (y_i - \bar y)^2$ .

In numpy, with standard errors

import numpy as np

def ols(X, y, add_intercept=True):
    """OLS with classical standard errors. X: (n, k) array, y: (n,) array."""
    X = np.asarray(X, dtype=float)
    y = np.asarray(y, dtype=float)
    if X.ndim == 1:
        X = X[:, None]
    if add_intercept:
        X = np.column_stack([np.ones(len(X)), X])
    n, k = X.shape
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)          # solves the normal equations, numerically stable
    resid = y - X @ beta
    sigma2 = resid @ resid / (n - k)
    cov = sigma2 * np.linalg.inv(X.T @ X)
    se = np.sqrt(np.diag(cov))
    r2 = 1 - resid @ resid / ((y - y.mean()) @ (y - y.mean()))
    return {"beta": beta, "se": se, "t": beta / se, "r2": r2, "n": n, "resid": resid}

np.linalg.lstsq is used rather than forming $(X^\top X)^{-1}$ directly for the coefficients because it is better conditioned; the inverse is still needed for the covariance, and for a handful of regressors that is fine.

On real data

Does today’s return predict tomorrow’s, pooled across every ticker and day? A slope with a t-statistic answers a slightly different question from the daily rank IC, and the two should agree in sign.

import pandas as pd

df = pd.read_csv("all_stocks_5yr.csv")
df["date"] = pd.to_datetime(df["date"])
df = df.rename(columns={"Name": "ticker"}).sort_values(["ticker", "date"]).reset_index(drop=True)
g = df.groupby("ticker")
df["ret"] = g["close"].pct_change()
df["fwd_ret"] = g["ret"].shift(-1)
d = df.dropna(subset=["ret", "fwd_ret"])
d = d[d["ret"].abs() < 0.25]                 # the pooled regression is not rank-based, so remove the unadjusted jumps

fit = ols(d["ret"].values, d["fwd_ret"].values)
print("n = %d   slope = %+.4f   se = %.4f   t = %+.2f   R2 = %.5f" % (fit["n"], fit["beta"][1], fit["se"][1], fit["t"][1], fit["r2"]))

n = 617984   slope = -0.0016   se = 0.0013   t = -1.26   R2 = 0.00000

The sign is reversal: a 1% move today is followed, on average, by 0.0016% in the other direction. The R² is zero to five decimals, as it must be for a one-day-ahead return, and the pooled t of −1.3 does not reach significance even with 618,000 rows. The daily rank IC on the same data does (t ≈ 2 on the hold-out), because ranking within each date removes the market-wide co-movement and the fat tails that a pooled regression on raw returns has to fight through. That is the usual lesson: on a panel, the pooled regression is the weaker and the less well-specified test.

What heteroskedasticity does

If the error variance is not constant (it never is for returns: volatile stocks and volatile periods have larger errors), then:

$\hat\beta$ is still unbiased and consistent. The slope is fine.
The classical standard errors are wrong, usually too small, so the t-statistics are too large.

The fix is a heteroskedasticity-robust (White, “HC”) covariance:

\widehat{\operatorname{Var}}_{HC0}(\hat\beta) = (X^\top X)^{-1} \Big( \sum_i e_i^2 \, x_i x_i^\top \Big) (X^\top X)^{-1}.

def ols_robust_se(X, y, add_intercept=True):
    """White (HC0) heteroskedasticity-robust standard errors."""
    X = np.asarray(X, dtype=float)
    if X.ndim == 1:
        X = X[:, None]
    if add_intercept:
        X = np.column_stack([np.ones(len(X)), X])
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    e = y - X @ beta
    XtX_inv = np.linalg.inv(X.T @ X)
    meat = (X * e[:, None] ** 2).T @ X
    cov = XtX_inv @ meat @ XtX_inv
    return beta, np.sqrt(np.diag(cov))

beta, se_hc = ols_robust_se(d["ret"].values, d["fwd_ret"].values)
print("slope = %+.4f   classical se = %.4f   robust se = %.4f   robust t = %+.2f" % (beta[1], fit["se"][1], se_hc[1], beta[1] / se_hc[1]))

slope = -0.0016   classical se = 0.0013   robust se = 0.0022   robust t = -0.74

The robust standard error is 70% larger than the classical one and the t-statistic falls from −1.26 to −0.74. Returns are heteroskedastic, the classical errors were too small, and this is the size of the effect on real data.

On a panel there is a second problem the robust errors do not fix: rows on the same date are correlated with each other through the market. The clean answer is to compute the slope per date and take the mean and its standard error across dates, which is the Fama–MacBeth construction used in does a signal add anything beyond a baseline. For a single-series regression, the concern is autocorrelation of the residuals rather than cross-sectional dependence.

With statsmodels, if installed

statsmodels is not part of the base stack this site assumes, but where it is available it gives all of the above in three lines, plus the robust errors with cov_type.

try:
    import statsmodels.api as sm
    X = sm.add_constant(d["ret"].values)
    res = sm.OLS(d["fwd_ret"].values, X).fit()
    res_hc = sm.OLS(d["fwd_ret"].values, X).fit(cov_type="HC0")
    print("statsmodels: slope %+.4f  classical t %+.2f  HC0 t %+.2f" % (res.params[1], res.tvalues[1], res_hc.tvalues[1]))
except ImportError:
    print("statsmodels not installed; the numpy version above is sufficient")

statsmodels: slope -0.0016  classical t -1.26  HC0 t -0.74

The numbers match the numpy implementation to the printed precision. If they do not, the numpy version has a bug, which is the point of checking.

sklearn, for completeness

sklearn.linear_model.LinearRegression().fit(X, y) returns the same coefficients and no standard errors. It is fine for prediction, and the wrong tool when the question is “is this slope distinguishable from zero”.

Pearson and Spearman in one sentence

The OLS slope on standardised variables is the Pearson correlation; Spearman is the Pearson correlation of the ranks, which is what to use when a few extreme rows would otherwise dominate. See Pearson vs Spearman.