Principal components

Many inputs move together. Ten points on the yield curve, five hundred stock returns, the beds, baths and square footage of a house: most of the variation lives in a few directions. Principal component analysis finds them.

The singular value decomposition

Centre the data matrix $X$ (subtract each column’s mean) and decompose it:

where $U$ and $V$ have orthonormal columns and $D$ is diagonal with $d_1 \ge d_2 \ge \dots \ge 0$. The columns $v_1, v_2, \dots$ of $V$ are the eigenvectors of $X^\top X$ (the covariance matrix up to a factor $N-1$): new axes, each perpendicular to the ones before.

The principal components are the data expressed on those axes, $z_1 = X v_1$, $z_2 = X v_2$, … The variance of $z_j$ is $d_j^2/(N-1)$, so the first component has the largest variance, the second the largest of what’s left, and so on.

Three ways PCA appears

• Unsupervised compression. Transform the inputs to their first $M$ components before fitting anything. No target needed.
• Principal components regression (PCR). Regress $y$ on $z_1, \dots, z_M$ with $M < p$. The components are orthogonal, so this is $M$ separate one-variable regressions. Then map the coefficients back to the original inputs.
• Ridge, seen through PCA. Ridge shrinks the coefficient on component $j$ by the factor $d_j^2/(d_j^2 + \lambda)$: it barely touches the big components and squashes the small ones. PCR is the hard-threshold version: it keeps the first $M$ untouched and drops the rest completely.

In code

import numpy as np
rng = np.random.default_rng(10)

# Fake yield-curve changes: level, slope and curvature factors plus noise
tenors = np.array([0.25, 1, 2, 5, 10, 30])
level, slope, curve = np.ones(6), (tenors - 7) / 10, -((np.log(tenors) - 1.3) ** 2) + 1.5
F = rng.standard_normal((2500, 3)) * [8, 4, 2]
dY = F @ np.vstack([level, slope, curve]) + rng.standard_normal((2500, 6))

dY = dY - dY.mean(axis=0)
U, d, Vt = np.linalg.svd(dY, full_matrices=False)
explained = d**2 / np.sum(d**2)
print(np.round(explained, 3))          # the first three components explain almost everything
print(np.round(Vt[0], 2))              # PC1 ≈ equal weights: a parallel "level" shift