Clustering — k-means, Gaussian Mixtures & Market Regimes

Two Algorithms, One Family

Unsupervised learning has no target and no loss to hold it to account — which is precisely what makes it worth doing carefully. k-means assigns each point to its nearest of $k$ centroids and recomputes the centroids until nothing moves (Lloyd's algorithm), making hard assignments and implicitly assuming spherical, equal-size clusters. A Gaussian mixture generalises it: every cluster is a full Gaussian with its own mean and covariance, fit by Expectation–Maximization, giving elliptical clusters and soft assignments. k-means is exactly the limit of a Gaussian mixture whose clusters are forced to be perfectly round, all the same size, and then shrunk until they have no width at all: once a cluster has no width, the probability of belonging to it collapses to yes-or-no and only the distance to its centre is left, which is k-means. The two are therefore one family rather than two methods, and are built here as one module.

The application is regime discovery: 3,459 trading days of S&P 500 data from 2000 to 2013, each described by just its return and its log realized volatility. No day carries a label. Whether the algorithm can recover the crises is a question with an answer we happen to know, which is what makes it a fair test of a method that normally has none.

Matching a reference solver, and what that does not prove

Both implementations are checked against scikit-learn, and the check is worth reading closely because the two things it can mean come apart. On fit they agree: k-means inertia — the total squared distance from every point to the centre of its own cluster, which is exactly the quantity k-means is trying to make small — comes to 2982.9 against 2983.0, GMM log-likelihood −8728.2 against −8732.1. On the partition they do not quite — the two 3-cluster labelings agree at ARI 0.896, not 1.000 — the adjusted Rand index counts the pairs of days that two partitions agree to place together or apart, corrected so that chance agreement scores 0 and an identical partition scores 1. EM is a local optimiser and this likelihood surface has several modes, so matching a reference solver on a scalar summary is weaker evidence than it looks. Here the from-scratch run reaches the higher likelihood of the two, so it is not the poorer solution; it is a different one.

How Many Clusters?

The number of clusters is the real decision, and the two standard tools disagree. The elbow points to three — inertia falls 36% going to $k=3$ and only 23% going to $k=4$. BIC — the Bayesian information criterion, which scores a model as its likelihood minus a penalty for the number of parameters it used, so extra clusters have to earn their cost — keeps falling to $k=7$ before turning up. Both are in-sample fit measures, so neither can arbitrate against the other, and running the BIC grid only as far as $k=7$ would make it look like it never stops improving at all.

What settles it is a criterion neither of them uses: stability. A clustering algorithm cannot decline to return clusters, so the question is whether a partition is a property of the data or of one particular fit. Refit on two overlapping 80% subsamples and compare the labels they give the days both of them saw, scored by the adjusted Rand index — 1 for identical partitions, 0 for chance agreement.

Stability does not name a single $k$ by itself: three, four and five all reproduce well. What it does decisively is eliminate the top of the range — from $k=6$ the GMM falls to 0.45, barely above chance, and BIC's preferred $k=7$ scores 0.57. Those extra components land somewhere different every time they are fitted, which is what an over-split looks like from the outside. The criteria compose rather than compete: stability rules out everything above five, and the elbow picks three from what survives. BIC on its own would have been wrong.

A control for why BIC overshoots

Why does BIC overshoot? Because a Gaussian mixture will spend components on non-Gaussian shape, not only on genuinely distinct groups. A control makes this concrete: simulate from a single fat-tailed distribution — one population, no regimes in it whatsoever — and score it identically. BIC asks for three components. The Gaussian control, correctly, asks for one. Daily equity returns are famously fat-tailed, so some of BIC's appetite for regimes here is kurtosis being fitted rather than structure being found.

The Regimes

The three regimes are economically legible, and not in the way the volatility axis alone would give. The two turbulent clusters sit at essentially the same volatility (+0.31 against +0.47 log-vol); what separates them is direction — 76% of the first are down days against 1% of the second. A sell-off regime and a rebound regime, at comparable turbulence. That is the return dimension earning its place, and it is worth confirming directly: agreement between this partition and a plain three-way split of log-volatility is only ARI 0.270, so these are not volatility terciles wearing a disguise.

Colouring the history by regime dates the crises with no dates supplied. Against an unconditional turbulent share of 34%, 2008 runs at 68% — but so does 2002 at 64% and 2009 at 62%. The dot-com unwind is fully comparable to the financial crisis rather than a footnote to it, and the mid-decade calm (2005 at 4%, 2006 at 6%) is as sharply separated as the crises are. This is the static, distributional shadow of a Markov-switching model, which adds what clustering cannot see: the dynamics of how one regime becomes another.

Where this sits

The same realized-volatility series is used in Realized Volatility, Recurrent Networks & LSTMs — Forecasting Volatility and Financial Returns Predictability, so what an unlabelled clustering finds in it can be set beside what the supervised and Bayesian treatments do. Volatility Persistence: Regimes vs Long Memory asks the sharper version of the question raised by the fat-tailed control — whether apparent regimes in volatility are regimes at all. The dynamics this notebook lacks are supplied by Bayesian Markov-Switching AR. On choosing $k$, Dirichlet-Process Mixtures — How Many Components? and Dirichlet-Process & Pitman–Yor Mixtures replace the grid search with inference, putting a prior over the cluster count and reporting a posterior; Latent Class Analysis — Choosing the Number of Classes is the same problem with categorical indicators, where a mixture model is called a latent class model.

Notebook

Downloads

Clustering Module — Source Code

"""From-scratch clustering: k-means (Lloyd's algorithm) and a Gaussian mixture
model fit by Expectation-Maximization.

k-means makes HARD assignments to the nearest centroid and implicitly assumes
spherical, equal-size clusters (it minimizes within-cluster squared distance).
A GMM generalizes it: each cluster is a full Gaussian (its own mean AND
covariance -- so clusters can be elliptical and overlapping), and EM produces
SOFT assignments (a posterior probability of belonging to each cluster).  In
fact k-means is the limit of GMM-EM with shared spherical covariance shrunk to
zero, so this file builds the two as one family.

Both are validated against scikit-learn in the notebook.  The Bayesian relative
-- a Dirichlet-process mixture that INFERS the number of clusters instead of
fixing k -- lives in the nonparametrics arc; here k is chosen by the elbow (SSE)
or BIC.
"""
import numpy as np
from scipy.stats import multivariate_normal


def kmeans(X, k, n_iter=100, seed=0):
    """Lloyd's algorithm. Returns (labels, centroids, inertia)."""
    X = np.asarray(X, float); rng = np.random.default_rng(seed)
    C = X[rng.choice(len(X), k, replace=False)].copy()          # random data points as seeds
    for _ in range(n_iter):
        D = ((X[:, None, :] - C[None, :, :]) ** 2).sum(2)       # squared distances to each centroid
        lab = D.argmin(1)
        newC = np.array([X[lab == j].mean(0) if np.any(lab == j) else C[j] for j in range(k)])
        if np.allclose(newC, C):
            break
        C = newC
    inertia = ((X - C[lab]) ** 2).sum()                         # within-cluster SSE
    return lab, C, inertia


class GMM:
    """Gaussian mixture model fit by Expectation-Maximization (full covariance).

    E-step: responsibility r_ij = pi_j N(x_i | mu_j, Sig_j) / sum_l (...).
    M-step: pi_j, mu_j, Sig_j <- responsibility-weighted proportions/means/covs.
    Initialized from k-means; `reg` adds a ridge to each covariance for stability.
    """

    def __init__(self, k, n_iter=200, tol=1e-6, reg=1e-6, seed=0):
        self.k = k; self.n_iter = n_iter; self.tol = tol; self.reg = reg; self.seed = seed

    def fit(self, X):
        X = np.asarray(X, float); n, d = X.shape
        lab, C, _ = kmeans(X, self.k, seed=self.seed)           # k-means initialization
        self.mu = C.copy()
        self.Sig = np.array([np.cov(X[lab == j].T) + self.reg * np.eye(d)
                             if np.sum(lab == j) > 1 else np.eye(d) for j in range(self.k)])
        self.pi = np.array([max((lab == j).mean(), 1e-3) for j in range(self.k)]); self.pi /= self.pi.sum()
        self.loglik_ = []
        ll_old = -np.inf
        for _ in range(self.n_iter):
            R = np.column_stack([self.pi[j] * multivariate_normal.pdf(X, self.mu[j], self.Sig[j], allow_singular=True)
                                 for j in range(self.k)])       # E-step (unnormalized)
            ll = np.log(R.sum(1) + 1e-300).sum(); self.loglik_.append(ll)
            R = R / (R.sum(1, keepdims=True) + 1e-300)          # responsibilities
            Nk = R.sum(0) + 1e-10                               # M-step
            self.pi = Nk / n
            self.mu = (R.T @ X) / Nk[:, None]
            for j in range(self.k):
                Xc = X - self.mu[j]
                self.Sig[j] = (R[:, j, None] * Xc).T @ Xc / Nk[j] + self.reg * np.eye(d)
            if abs(ll - ll_old) < self.tol:
                break
            ll_old = ll
        self.n_params = self.k * (d + d * (d + 1) / 2) + (self.k - 1)   # means + covs + weights
        return self

    def predict_proba(self, X):
        X = np.asarray(X, float)
        R = np.column_stack([self.pi[j] * multivariate_normal.pdf(X, self.mu[j], self.Sig[j], allow_singular=True)
                             for j in range(self.k)])
        return R / (R.sum(1, keepdims=True) + 1e-300)

    def predict(self, X):
        return self.predict_proba(X).argmax(1)

    def bic(self, X):
        ll = self.loglik_[-1]
        return -2 * ll + self.n_params * np.log(len(X))

References

• Lloyd, S. P. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory 28(2), 129–137. — the k-means algorithm
• Dempster, A. P., Laird, N. M. & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society B 39(1), 1–38. — EM
• Hubert, L. & Arabie, P. (1985). Comparing partitions. Journal of Classification 2(1), 193–218. — the adjusted Rand index
• Ben-Hur, A., Elisseeff, A. & Guyon, I. (2002). A stability based method for discovering structure in clustered data. Pacific Symposium on Biocomputing, 6–17. — resampling stability as a criterion for k
• Schwarz, G. (1978). Estimating the dimension of a model. Annals of Statistics 6(2), 461–464. — BIC
• Hamilton, J. D. (1989). A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica 57(2), 357–384. — regimes with dynamics