---
title: "5.5.5 PCA, clustering and autoencoders"
description: "PCA, stock and regime clustering and a Keras autoencoder on 400 synthetic stocks with known factors: what they recover and what they do not."
url: "https://pythonparatrading.com/en/python-trading-course/level-5/pca-clustering-autoencoders/"
lang: "en"
image: "https://pythonparatrading.com/og/en/python-trading-course/level-5/pca-clustering-autoencoders.png"
---
# 5.5.5 · Unsupervised learning: PCA, clustering and autoencoders
> Without a label there is no return to predict, but there is structure to discover: common factors, groups of assets and market regimes. This lesson applies PCA, clustering and a Keras autoencoder to 400 synthetic stocks whose factor model and sectors are known, and measures what each method recovers, with how much data and with what delay.
- Course: [Level 5: Machine learning and deep learning (Keras/PyTorch)](https://pythonparatrading.com/en/python-trading-course/level-5/)
- Module: 5.5 · Classical models
- Duration: 65 min
- Prerequisites: [5.1.1 · Why financial machine learning is different](https://pythonparatrading.com/en/python-trading-course/level-5/why-financial-machine-learning-is-different/), [4.1.2 · Covariance estimation: shrinkage, factor models and denoising](https://pythonparatrading.com/en/python-trading-course/level-4/covariance-estimation-shrinkage-and-denoising/), [4.1.6 · Factor models: CAPM, Fama-French, PCA and attribution](https://pythonparatrading.com/en/python-trading-course/level-4/factor-models-capm-fama-french-pca/), [4.1.4 · Risk parity, minimum variance, HRP and HERC](https://pythonparatrading.com/en/python-trading-course/level-4/risk-parity-minimum-variance-hrp-herc/), [2.1.7 · Regime detection: hidden Markov models and structural breaks](https://pythonparatrading.com/en/python-trading-course/level-2/regime-detection-hmm-structural-breaks/)
## By the end of this lesson
- Compute PCA on the correlation matrix of returns, decide how many components are structure with the Marchenko-Pastur edge and identify the market and the sectors in the loadings.
- Group stocks with k-means on loadings and with hierarchical clustering on correlations, measure the agreement with the true sectors (ARI) and know why the market has to be removed first.
- Detect regimes with a Gaussian mixture fitted only with the past, and measure its agreement with the true regimes and its delay.
- Train an autoencoder in Keras 3 with the PyTorch backend, check that the linear one is equivalent to PCA and compare the nonlinear one out of sample.
## Structure without a label
The previous lessons in the module learn a function from the features to the label. **Unsupervised** learning has no
label: it looks for structure in the data themselves. In trading it has four common uses, and none of them predicts
returns on its own:
```mermaid
flowchart TB
A[N returns] --> B[PCA]
A --> C[Groups of assets]
A --> D[Groups of dates]
A --> E[Autoencoder]
B --> F[Factors, residuals]
C --> G[Diversification]
D --> H[Regimes]
E --> I[Anomalies]
```
- **Statistical factors** (PCA): summarize 400 series in a few common directions; the residuals are the basis of
statistical arbitrage (Avellaneda and Lee, 2010) and the covariance matrix is cleaned with them (lesson 4.1.2).
- **Groups of assets**: “real” sectors, not those of the classification, to diversify or neutralize a portfolio
(HRP, lesson 4.1.4).
- **Regimes**: group dates with similar behavior (lesson 2.1.7 with hidden Markov models).
- **Compression**: represent each day with a few variables and measure what does not fit (anomalies).
The problem with evaluating unsupervised learning is that there is no metric like the IC. Here there is, because the
data are **synthetic** and the model that generated them is published (README of `public/curso/datos/nivel-5/`): the
daily return of each of the 400 stocks is
$$
r_{i,t} = \beta_{i,t}\, f_{m,t} + f_{s(i),t} + \varepsilon_{i,t} + \ldots
$$
a market factor with one beta per stock, one factor per sector (10 sectors of 40 stocks), idiosyncratic noise with
fat tails and small effects (weekly alpha, earnings). Each method is judged by what it recovers of that truth.
## The data: 400 daily returns
We work with the simple daily returns of the closes, standardized with the mean and standard deviation of
**2011-2018**, the fitting period. 2019-2025 is kept to measure out of sample. Standardizing with the whole sample
would be the leakage of lesson 5.2.4 in its unsupervised version.
```python
import numpy as np
import pandas as pd
pr = pd.read_parquet("panel5_precios.parquet",
columns=["fecha", "ticker", "cierre"])
R = (pr.pivot(index="fecha", columns="ticker", values="cierre")
.pct_change().iloc[1:])
activos = pd.read_csv("panel5_activos.csv", index_col="ticker")
verdad = pd.read_csv("panel5_activos_verdad.csv", index_col="ticker")
sector = activos.loc[R.columns, "sector"].to_numpy()
beta = verdad.loc[R.columns, "beta"].to_numpy()
ent = (R.index < "2019-01-01") # 2011-2018
mu, sd = R[ent].mean(), R[ent].std() # training only
Z = ((R - mu) / sd).to_numpy()
Ze, Zp = Z[ent], Z[~ent]
print("sessions x stocks:", R.shape)
print("training:", Ze.shape[0], " test:", Zp.shape[0])
print("sectors:", pd.Series(sector).value_counts().unique())
```
```output
sessions x stocks: (3771, 400)
training: 2011 test: 1760
sectors: [40]
```
## PCA and the Marchenko-Pastur edge
PCA looks for the directions of maximum variance: the eigenvectors of the correlation matrix, sorted by their
eigenvalue. The practical question is how many components are structure and how many are noise. If returns were
independent noise, the eigenvalues of their correlation matrix with $T$ sessions and $N$ assets would fall, for large
$T$ and $N$, below
$$
\lambda_+ = \left(1 + \sqrt{N/T}\right)^2
$$
the upper edge of the Marchenko-Pastur distribution. Lesson 4.1.2 used it to clean covariances; here it serves to
count factors.
```python
C = np.corrcoef(Ze, rowvar=False)
val, vec = np.linalg.eigh(C)
val, vec = val[::-1], vec[:, ::-1] # largest to smallest
T, N = Ze.shape
borde = (1 + np.sqrt(N / T)) ** 2 # λ+ with variance 1
print(f"N/T = {N / T:.3f}, Marchenko-Pastur edge = {borde:.2f}")
print("first 12 eigenvalues:", val[:12].round(2))
print(f"above the edge: {(val > borde).sum()}")
print(f"explained variance: PC1 {val[0] / N:.1%}, "
f"PC1-PC10 {val[:10].sum() / N:.1%}")
```
```output
N/T = 0.199, Marchenko-Pastur edge = 2.09
first 12 eigenvalues: [120.51 9.79 7.4 6.62 5.88 5.74 5.44 4.71 4.28 3.72
1.37 1.35]
above the edge: 10
explained variance: PC1 30.1%, PC1-PC10 43.5%
```
The header figure plots the first twenty. The first, 120.5, takes 30.1% of the variance: it is the market. From the
second to the tenth they fall from 9.8 to 3.7 and stay clearly above the edge (2.09). From the eleventh onward they
pile up between 1.2 and 1.4: noise. **Ten** components, which together explain 43.5% of the variance.
The generator has **eleven** factors. The missing one has not been lost: the market moves all stocks in the same
direction, and so does the average of the ten sector factors, so PCA merges those two directions into the first
component and what remains distinguishable are nine **contrasts** between sectors. A principal component is a
direction of variance, not a factor of the model one by one; naming each one requires looking at its loadings.
## What each component is
```python
pc1 = vec[:, 0] * np.sign(vec[:, 0].sum())
carga = pc1 * sd.to_numpy() # in return units
f1 = Ze @ pc1 # component series
ew = R[ent].mean(axis=1).to_numpy() # equal-weighted portfolio
print(f"corr(PC1 loading, true beta) = "
f"{np.corrcoef(carga, beta)[0, 1]:.3f}")
print(f"corr(PC1 loading not rescaled, beta) = "
f"{np.corrcoef(pc1, beta)[0, 1]:.3f}")
print(f"corr(PC1 series, equal-weighted port.) = "
f"{np.corrcoef(f1, ew)[0, 1]:.3f}")
```
```output
corr(PC1 loading, true beta) = 0.987
corr(PC1 loading not rescaled, beta) = 0.620
corr(PC1 series, equal-weighted port.) = 0.999
```
The PC1 series has a correlation of 0.999 with the equal-weighted portfolio. Its loadings, once multiplied by each
stock's volatility (PCA of correlations works with standardized returns), have a correlation of 0.987 with the **true
beta**; without rescaling, only 0.620, because the loading of a correlation measures $\beta_i
\sigma_m/\sigma_i$ and the most volatile stocks are penalized. It is a detail often forgotten when using PC1 as a
statistical beta.
```python
cargas = pd.DataFrame(vec[:, 1:10] * np.sqrt(val[1:10]),
columns=[f"PC{i}" for i in range(2, 11)])
cargas["sector"] = sector
media = cargas.groupby("sector").mean()
print(media.round(2).to_string())
```
```output
PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
sector
consumo_basico -0.03 -0.04 -0.04 -0.06 -0.01 0.03 -0.25 -0.16 -0.02
consumo_ciclico -0.05 -0.03 -0.03 0.03 -0.01 0.01 -0.06 0.19 -0.20
energia 0.46 0.03 0.00 -0.00 0.02 -0.01 0.01 -0.01 -0.01
finanzas -0.10 0.35 0.10 -0.06 0.05 -0.00 0.04 -0.03 -0.02
industria -0.04 0.00 -0.03 -0.00 0.01 -0.04 -0.09 0.16 0.21
materiales -0.06 -0.02 -0.05 0.30 0.03 -0.13 0.05 -0.08 0.00
salud -0.03 -0.18 0.34 -0.02 0.02 0.00 0.03 -0.02 0.01
servicios_publicos -0.04 -0.05 -0.06 -0.15 -0.25 -0.16 0.09 -0.03 -0.01
tecnologia -0.03 -0.03 -0.06 0.05 -0.12 0.29 0.10 -0.03 0.03
telecomunicaciones -0.06 -0.13 -0.14 -0.14 0.24 -0.01 0.11 -0.03 -0.01
```
Each of components 2 to 10 loads heavily on one or two sectors and little on the rest: PC2 is energy (`energia`,
0.46) against all the others; PC3, financials (`finanzas`, 0.35) against health care and telecommunications; PC4,
health care (`salud`, 0.34). The last ones mix more (PC8 opposes consumer staples, −0.25, to telecommunications and
technology). PCA knows nothing about sectors; it finds them because they are the structure of the covariance. In real
data, loadings mix sectors, styles (value, momentum) and size, and their interpretation changes from one window to
another (lesson 4.1.6).
## Grouping stocks
If each sector occupies its own direction in the space of components 2 to 10, a k-means on those loadings
(normalized so that the direction counts, not the size) should find the groups. To measure it, the **adjusted Rand
index** (ARI): 1 if two partitions coincide (with any group names), around 0 if they are as alike as they would be by
chance.
```python
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score
L = vec[:, 1:10] * np.sqrt(val[1:10])
L = L / np.linalg.norm(L, axis=1, keepdims=True)
km = KMeans(n_clusters=10, n_init=10, random_state=0).fit(L)
print(f"ARI k-means (PC2-PC10 loadings) versus sectors: "
f"{adjusted_rand_score(sector, km.labels_):.3f}")
print(pd.crosstab(km.labels_, sector).max(axis=1).tolist())
```
```output
ARI k-means (PC2-PC10 loadings) versus sectors: 1.000
[40, 40, 40, 40, 40, 40, 40, 40, 40, 40]
```
ARI of 1.000: the ten groups are the ten sectors, with 40 stocks each. The most common alternative in finance is
**hierarchical clustering** on the correlation distance $d_{ij} = \sqrt{2(1-\rho_{ij})}$, the same one used by the
HRP of lesson 4.1.4. The cell tries it with and without the first component:
```python
from sklearn.cluster import AgglomerativeClustering
def jerarquico(Zv, quitar_mercado=True, k=10):
"""Clusters by correlation (distance sqrt(2(1-rho)))."""
if quitar_mercado:
w, v = np.linalg.eigh(np.corrcoef(Zv, rowvar=False))
Zv = Zv - np.outer(Zv @ v[:, -1], v[:, -1])
rho = np.corrcoef(Zv, rowvar=False)
dist = np.sqrt(np.clip(2 * (1 - rho), 0, None))
return AgglomerativeClustering(
n_clusters=k, metric="precomputed",
linkage="average").fit(dist).labels_
for quitar in (False, True):
lab = jerarquico(Ze, quitar)
print(f"remove PC1={quitar!s:5}: ARI "
f"{adjusted_rand_score(sector, lab):.3f}, sizes "
f"{sorted(np.bincount(lab).tolist(), reverse=True)}")
```
```output
remove PC1=False: ARI 0.001, sizes [386, 5, 2, 1, 1, 1, 1, 1, 1, 1]
remove PC1=True : ARI 1.000, sizes [40, 40, 40, 40, 40, 40, 40, 40, 40, 40]
```
Without removing anything, the result is useless: one group of 386 stocks and nine of one to five (ARI 0.001). With an
average correlation of 0.32 because of the market, all stocks look alike and the sector differences stay buried;
average linkage keeps merging everything into a single group. Removing the projection on PC1, the residuals keep only
the sector part and the clustering recovers the ten exact sectors.
That it comes out perfect says more about the data than about the method: the synthetic sectors are clean (40 stocks,
one factor each, no company straddling two) and there are eight years of daily data. The next cell repeats the
clustering with shorter windows, in six different stretches of 2011-2018:
```python
filas = {}
for n in (21, 63, 126, 252, 504):
aris = [adjusted_rand_score(sector, jerarquico(Ze[s:s + n]))
for s in range(0, 1800, 300)]
filas[n] = {"ARI medio": np.mean(aris), "mínimo": np.min(aris),
"máximo": np.max(aris)}
pd.DataFrame(filas).T.rename_axis("sesiones").round(2)
```
```output
ARI medio mínimo máximo
sesiones
21 0.21 0.11 0.29
63 0.39 0.29 0.55
126 0.58 0.46 0.82
252 0.84 0.74 0.92
504 0.97 0.93 0.99
```
*Fig. — Hierarchical clustering of the 400 stocks (correlation without the first component, 10 groups) with windows of different length within 2011-2018: adjusted Rand index against the true sectors, mean and range of six windows.* A line rising with window length, with bars from minimum to maximum: 0.21 with 21 sessions, 0.39 with 63, 0.58 with 126, 0.84 with 252 and 0.97 with 504 (highlighted). With one month of data the groups barely resemble the sectors; with two years they recover them almost entirely.
With one month of data (21 sessions), the average ARI is 0.21; with a quarter, 0.39; with a year, 0.84; with two,
0.97. The groups of a clustering done with little data depend on the window, and an HRP recomputed with 63 sessions
reorganizes its tree every time. In real data, with a less clean structure that also changes, the window length is a
trade-off between noise and recency.
## The first component in crises
The variance taken by PC1 is not constant. Computed quarter by quarter, it measures how much stocks move together:
it is the idea of the *absorption ratio* of Kritzman et al. (2011) as a measure of systemic risk.
```python
Rv = R.to_numpy()
trim = {}
for per, idx in R.groupby(R.index.to_period("Q")).indices.items():
Zq = Rv[idx]
Zq = (Zq - Zq.mean(0)) / Zq.std(0)
ev = np.linalg.eigvalsh(np.corrcoef(Zq, rowvar=False))
trim[str(per)] = ev[-1] / ev.sum()
trim = pd.Series(trim)
print(f"PC1 by quarter: median {trim.median():.2f}, "
f"minimum {trim.min():.2f} ({trim.idxmin()}), "
f"maximum {trim.max():.2f} ({trim.idxmax()})")
print(trim.nlargest(5).round(2).to_string())
```
```output
PC1 by quarter: median 0.31, minimum 0.13 (2019Q2), maximum 0.68 (2020Q1)
2020Q1 0.68
2020Q2 0.68
2015Q4 0.53
2018Q4 0.53
2016Q1 0.52
```
*Fig. — Fraction of the 400 stocks' variance captured by the first component, computed quarter by quarter (about 63 sessions) from 2011 to 2025: it rises in the generator's stress spells.* A quarterly line oscillating between 13% and 53% most of the time, with peaks in 2011, late 2015, late 2018 and 2022, and a maximum of 68% in the first and second quarters of 2020 (highlighted). The median is 31%.
The median is 31%, with a minimum of 13% in the second quarter of 2019 and a maximum of 68% in the first two quarters
of 2020, the generator's hardest stress stretch (market volatility of 80% annualized). The other peaks (late 2015,
late 2018, 2022) coincide with the other stress stretches. It is what you expect of any portfolio in a crisis:
diversification across stocks disappears just when it is needed (lesson 4.2.4, correlation stress testing). Note
that quarters of 63 sessions give a Marchenko-Pastur edge of about 13: with short windows only the market stands out
from the noise (exercise 2 counts it year by year).
## Regimes with clustering
What groups stocks can also group **dates**. Each Friday is described with four features computed only with the
past: the logarithm of market volatility over 21 sessions, the approximate average correlation (market variance
divided by the average variance of the stocks, 63 sessions), the market return over 21 sessions and the logarithm of
the cross-sectional dispersion between stocks. A **mixture of three Gaussians** is fitted with 2012-2018 and
classifies all the weeks; the generator publishes the true regime of each one (`calma` for calm, `estres` for stress,
`rebote` for rebound) in `panel5_coef_verdad.csv`.
```python
from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler
reg = pd.read_csv("panel5_coef_verdad.csv", parse_dates=["fecha"],
index_col="fecha")["regimen"]
mkt = R.mean(axis=1)
var_i = R.rolling(63).var().mean(axis=1)
F = pd.DataFrame({
"log_vol_21": np.log(mkt.rolling(21).std() * np.sqrt(252)),
"correlacion_63": mkt.rolling(63).var() / var_i,
"ret_21": mkt.rolling(21).sum(),
"log_dispersion_21": np.log(R.std(axis=1).rolling(21).mean()),
}).reindex(reg.index, method="ffill") # Friday of each week
e2 = F.index < "2019-01-01"
esc = StandardScaler().fit(F[e2])
gm = GaussianMixture(3, covariance_type="full", n_init=5,
random_state=0).fit(esc.transform(F[e2]))
grupo = pd.Series(gm.predict(esc.transform(F)), index=F.index)
print(pd.DataFrame(esc.inverse_transform(gm.means_),
columns=F.columns).round(3))
print(pd.crosstab(reg[~e2], grupo[~e2],
rownames=["2019-2025"], colnames=["grupo"]))
```
```output
log_vol_21 correlacion_63 ret_21 log_dispersion_21
0 -1.691 0.324 0.035 -4.070
1 -2.224 0.186 0.008 -4.158
2 -1.242 0.425 -0.052 -3.936
grupo 0 1 2
2019-2025
calma 72 190 1
estres 15 4 33
rebote 36 0 15
```
The group numbers are assigned by the algorithm; their means say what each one is. Group 2 has the highest
volatility ($e^{-1.24} \approx$ 29% annualized), the highest correlation (0.43) and a negative return: it is stress.
Group 1 is quiet calm (volatility of 11%, correlation 0.19). Group 0 sits between the two. In 2019-2025, outside the
fitting sample:
- **Stress**: 33 of its 52 weeks fall in group 2 and 15 in the intermediate one.
- **Calm**: 190 of 263 in group 1 and 72 in the intermediate one; a single one in the stress group.
- **Rebound**: none in calm, but spread between the intermediate group (36) and stress (15). A rebound is volatile
like stress and is only distinguished by the sign of the return, which over three months of data is noisy.
Grouping dates also has its version of delay: the features are backward-looking windows, so the regime has to last to
move them.
```python
estres = grupo == int(np.argmax(gm.means_[:, 0]))
inicios = reg.index[(reg == "estres") & (reg.shift() != "estres")]
for ini in inicios[inicios >= "2019-01-01"]:
tras = estres[ini:]
primera = tras.index[tras.to_numpy()][0]
print(f"stress from {ini:%Y-%m-%d}: the group flags it on "
f"{primera:%Y-%m-%d} ({(primera - ini).days // 7} weeks)")
```
```output
stress from 2020-02-21: the group flags it on 2020-03-06 (2 weeks)
stress from 2022-01-07: the group flags it on 2022-02-11 (5 weeks)
stress from 2025-02-21: the group flags it on 2025-03-14 (3 weeks)
```
The February 2020 stress appears in group 2 two weeks after it starts; that of January 2022, five (it was a mild
stress, 25% volatility); that of February 2025, three. A detector without future information always arrives late:
the useful question is how late and how much that delay costs the strategy that uses it. Lesson 2.1.7 did the same
with an HMM, which adds the persistence of the states, and lesson 5.10.2 returns to regimes as data drift.
## Autoencoders
An **autoencoder** is a network that learns to reproduce its input by passing through a narrow bottleneck: an encoder
compresses the day's 400 returns into $k$ numbers and a decoder reconstructs them. It is trained with squared error,
without a label.
```mermaid
flowchart TB
A[400 returns] --> B[tanh layer 64]
B --> C[Bottleneck of 10]
C --> D[tanh layer 64]
D --> E[400 reconstructed]
E --> F[Squared error]
```
If all the layers are linear, the best autoencoder with a bottleneck of $k$ reconstructs from the same subspace as
the first $k$ principal components: it is PCA with another algorithm. The promise of nonlinear layers is to capture
structure that PCA does not see. The cell trains both in **Keras 3 with the PyTorch backend** (the `KERAS_BACKEND`
variable is set by the runner; lesson 5.6.1 explains the backends), with a seed (`keras.utils.set_random_seed`),
deterministic torch algorithms and 60 epochs on 2011-2018. In production you would try more architectures,
regularization and early stopping with a purged validation.
```python
import os
import keras
import torch
torch.use_deterministic_algorithms(True)
print("backend:", keras.backend.backend(), os.environ["KERAS_BACKEND"])
def r2(A, B):
"""Fraction of the (standardized) variance reconstructed."""
return 1 - ((A - B) ** 2).sum() / (A ** 2).sum()
def autoencoder(k, oculta=None, epocas=60):
keras.utils.set_random_seed(0)
x = entrada = keras.Input((N,))
if oculta:
x = keras.layers.Dense(oculta, activation="tanh")(x)
x = keras.layers.Dense(k, name="cuello")(x)
if oculta:
x = keras.layers.Dense(oculta, activation="tanh")(x)
salida = keras.layers.Dense(N)(x)
m = keras.Model(entrada, salida)
m.compile(optimizer=keras.optimizers.Adam(1e-3), loss="mse")
m.fit(Ze.astype("float32"), Ze.astype("float32"), epochs=epocas,
batch_size=64, shuffle=True, verbose=0)
return m
k = 10
V = vec[:, :k]
res = {"PCA": (r2(Ze, Ze @ V @ V.T), r2(Zp, Zp @ V @ V.T))}
for nombre, oculta in [("autoencoder lineal", None),
("autoencoder tanh 64", 64)]:
m = autoencoder(k, oculta)
pe = m.predict(Ze.astype("float32"), verbose=0)
pp = m.predict(Zp.astype("float32"), verbose=0)
res[nombre] = (r2(Ze, pe), r2(Zp, pp))
pd.DataFrame(res, index=["R2 2011-2018", "R2 2019-2025"]).T.round(4)
```
```output
backend: torch torch
R2 2011-2018 R2 2019-2025
PCA 0.4352 0.4811
autoencoder lineal 0.4350 0.4811
autoencoder tanh 64 0.4301 0.4408
```
The measure is the fraction of the standardized variance that is reconstructed ($R^2$), in sample (2011-2018) and out
of sample (2019-2025):
- The **linear autoencoder** matches PCA to the fourth decimal out of sample (0.4811 both): it learns the same
subspace, more slowly.
- The **nonlinear** one (tanh layers of 64) loses little in sample (0.4301 versus 0.4352) and quite a lot out of
sample (0.4408 versus 0.4811). It has nothing nonlinear to find, because the true model is a linear factor model,
and its extra parameters fit peculiarities of 2011-2018 that do not repeat.
```python
curva = {}
for kk in (1, 2, 3, 5, 10, 15, 20, 30, 50):
Vk = vec[:, :kk]
curva[kk] = {"dentro": r2(Ze, Ze @ Vk @ Vk.T),
"fuera": r2(Zp, Zp @ Vk @ Vk.T)}
pd.DataFrame(curva).T.rename_axis("componentes").round(4)
```
```output
dentro fuera
componentes
1 0.3013 0.3678
2 0.3257 0.3894
3 0.3442 0.4046
5 0.3755 0.4281
10 0.4352 0.4811
15 0.4517 0.4889
20 0.4672 0.4963
30 0.4963 0.5116
50 0.5500 0.5412
```
*Fig. — Fraction of the 400 stocks' standardised variance that PCA reconstructs by number of components, fitted on 2011-2018 and measured in and out of sample (2019-2025), against two Keras autoencoders with a bottleneck of 10.* Two rising PCA curves: out of sample, from 0.368 with one component to 0.541 with 50; in sample, from 0.301 to 0.550. With 10 components, PCA out of sample gives 0.481 and the linear autoencoder matches it (0.481); the non-linear autoencoder with tanh layers of 64 (highlighted) falls below, at 0.441.
The figure places the two autoencoders on the PCA curve. With a bottleneck of 10 (the components above the edge), the
nonlinear one reconstructs out of sample little more than PCA with 5 components (0.428). And the PCA curve has its own
lesson: going from 10 to 50 components raises the in-sample $R^2$ (from 0.435 to 0.550) more than the out-of-sample one
(from 0.481 to 0.541), because those 40 extra components are the noise that the Marchenko-Pastur edge had already
discarded.
That the $R^2$ is higher out of sample than in sample is not leakage: 2019-2025 contains 2020, when the market took up
to 68% of the variance, and a period with more correlation is reconstructed better with the same components. That is
why methods are compared over the same period.
## Common mistakes
- **PCA on non-standardized returns or with the whole sample.** Without standardizing, volatile stocks dominate the
components; standardizing with the whole sample, the test data enter the fit.
- **Counting factors by eye.** The Marchenko-Pastur edge gives a reference; the components below it are noise and
adding them raises the in-sample fit, not the out-of-sample one.
- **Clustering correlations without removing the market.** You get one giant group and a few stragglers.
- **Reading group numbers as fixed labels.** k-means or a mixture may number the groups differently when refitted;
a group is identified by its means.
- **Regimes fitted with the whole sample.** A detector that “gets” past regimes right with later information does not
say when it would have detected them.
- **Using an autoencoder because it is more modern.** If it does not improve on PCA out of sample with the same
dimension, it is PCA that is more expensive and less stable.
## Exercises
1. Repeat the hierarchical clustering without removing the market but with 11 groups (in case the market takes one).
Does it improve? Then try Ward linkage on the residuals without PC1.
**Solution.**
```python
lab = jerarquico(Ze, quitar_mercado=False, k=11)
print(sorted(np.bincount(lab).tolist(), reverse=True))
print(f"ARI: {adjusted_rand_score(sector, lab):.3f}")
lab_w = AgglomerativeClustering(n_clusters=10, linkage="ward").fit(
(Ze - np.outer(Ze @ vec[:, 0], vec[:, 0])).T).labels_
print(f"Ward on the residuals: ARI "
f"{adjusted_rand_score(sector, lab_w):.3f}")
```
```output
[385, 5, 2, 1, 1, 1, 1, 1, 1, 1, 1]
ARI: 0.001
Ward on the residuals: ARI 0.977
```
With 11 groups the result is the same disaster (385 stocks in one group, ARI 0.001): the problem is not the number of
groups, it is that the correlation distance with the market left in does not separate sectors. Ward on the residuals
gives an ARI of 0.977, almost perfect: Ward merges the groups that least increase the within-group variance and, with
clean residuals, works as well as average linkage on the correlation.
2. Compute year by year how many eigenvalues exceed the Marchenko-Pastur edge and what fraction of the variance PC1
takes. What happens in 2020?
**Solution.**
```python
anios = {}
for a, idx in R.groupby(R.index.year).indices.items():
Za = Rv[idx]
Za = (Za - Za.mean(0)) / Za.std(0)
ev = np.linalg.eigvalsh(np.corrcoef(Za, rowvar=False))
b = (1 + np.sqrt(N / len(idx))) ** 2
anios[a] = {"sesiones": len(idx), "borde": round(b, 2),
"encima": int((ev > b).sum()),
"PC1": round(ev[-1] / N, 2)}
pd.DataFrame(anios).T.astype({"sesiones": int, "encima": int})
```
```output
sesiones borde encima PC1
2011 251 5.12 8 0.31
2012 250 5.13 7 0.32
2013 252 5.11 8 0.23
2014 252 5.11 7 0.28
2015 252 5.11 8 0.38
2016 252 5.11 7 0.37
2017 251 5.12 8 0.20
2018 251 5.12 8 0.34
2019 252 5.11 8 0.26
2020 253 5.10 2 0.58
2021 252 5.11 7 0.28
2022 251 5.12 6 0.44
2023 250 5.13 9 0.23
2024 252 5.11 8 0.29
2025 250 5.13 7 0.35
```
With about 252 sessions per year, $N/T \approx 1.6$ and the edge rises to 5.1: with one year of data only 6 to 9
components can be distinguished, not the 10 of the long sample. In 2020, only **2**: the market takes 58% of the
variance and crushes the sectors. In crises the structure simplifies to one factor; just when diversification across
sectors matters most is when it is least visible in the data.
3. A colleague computes the first five principal components with the returns of 2011-2025, uses each stock's exposure
to each component as a feature of a 5.5.3 model and validates it with walk-forward from 2019. What is wrong?
**Solution.**
The components were computed with returns from 2019-2025, the test period: each stock's loadings incorporate how it
moved afterward, and the walk-forward inherits that information even if the supervised model is retrained properly.
The standardization has the same problem. You have to recompute the standardization and the PCA in each training
window (with the data up to the retraining date) and also align the components across windows, because their sign and
order can change from one to the next (for example, fixing the sign of PC1 so that its sum is positive, as in the
lesson, and matching the others by correlation of loadings).
## Self-assessment
1. The generator has 11 factors (market and 10 sectors), but only 10 eigenvalues exceed the Marchenko-Pastur edge. Why?
- A. Because the market and the average of the ten sector factors move all stocks in the same direction: PCA merges them into the first component and nine contrasts between sectors remain
- B. Because the edge is miscalculated
- C. Because one sector has no factor of its own
- D. Because PCA can only find 10 components with 400 stocks
*Answer A: Because the market and the average of the ten sector factors move all stocks in the same direction: PCA merges them into the first component and nine contrasts between sectors remain.* The “all stocks at once” directions of the market and of the sum of sectors almost coincide. What remains distinguishable are the differences between sectors: nine components. Counting factors with PCA gives directions, not the generator's factors one by one.
2. A hierarchical clustering of the correlation of the 400 stocks gives one group of 386 and nine almost empty ones (ARI 0.001). Removing the first component, it recovers the 10 sectors (ARI 1.000). What happened?
- A. Clustering without removing anything overfits
- B. Removing PC1 introduces leakage
- C. The market dominates the correlations (0.32 on average between stocks): all stocks look alike because of it and the sector differences stay buried until it is removed
- D. Average linkage does not work for stocks
*Answer C: The market dominates the correlations (0.32 on average between stocks): all stocks look alike because of it and the sector differences stay buried until it is removed.* With a strong common factor, the distance between two stocks from different sectors is almost the same as between two from the same sector. The residuals without the market keep only the sector part, and there the groups appear.
3. The Gaussian mixture fitted with 2012-2018 flags the start of the February 2020 stress two weeks late and that of January 2022, five. What explains it?
- A. That the mixture is badly trained
- B. That the true regimes are wrongly dated
- C. That more groups are needed
- D. That the features are past windows (21-session volatility, 63-session correlation): a regime has to last long enough to move their averages, and a mild stress takes longer
*Answer D: That the features are past windows (21-session volatility, 63-session correlation): a regime has to last long enough to move their averages, and a mild stress takes longer.* A detector that only uses past information always arrives late: the question is how late. The 2022 stress is mild (25% vol.), and the window takes longer to separate from the calm.
4. With a bottleneck of 10, the linear autoencoder reconstructs 48.11% of the variance out of sample (the same as PCA) and the nonlinear one with tanh layers, 44.08%. What do you conclude?
- A. That autoencoders do not work with financial data
- B. That, when the true structure is linear (a factor model), an autoencoder cannot improve on PCA and its flexibility is spent on noise; the nonlinear one has to be justified with data that have nonlinear structure
- C. That more epochs are needed to match PCA
- D. That PCA overfits
*Answer B: That, when the true structure is linear (a factor model), an autoencoder cannot improve on PCA and its flexibility is spent on noise; the nonlinear one has to be justified with data that have nonlinear structure.* A linear autoencoder with squared error learns the same subspace as PCA. The nonlinear one adds parameters that fit peculiarities of 2011-2018 that do not repeat: in sample it loses little (0.430 versus 0.435) and out of sample quite a lot more.
5. Why is the PCA reconstruction R² higher in 2019-2025 (0.481) than in 2011-2018 (0.435), where it was fitted?
- A. Because 2019-2025 includes 2020, when the market explains much more variance (the first component reaches 68% in one quarter): a period with more correlation is easier to reconstruct with the same components
- B. Because of leakage from the test period
- C. Because PCA improves over time
- D. Because the standardization uses the test data
*Answer A: Because 2019-2025 includes 2020, when the market explains much more variance (the first component reaches 68% in one quarter): a period with more correlation is easier to reconstruct with the same components.* The R² depends on how much variance is common in the test period. Comparing in and out of sample only makes sense between methods measured over the same period, not between periods.
## Key points
- PCA of the correlation matrix of 400 stocks: 10 eigenvalues above the Marchenko-Pastur edge (2.09). PC1 is the market (loading × volatility versus beta: correlation 0.987) and PC2-PC10, contrasts between sectors.
- The true model has 11 factors and PCA sees 10: the components are directions of variance, not the generator's factors one by one.
- To group stocks you have to remove the market first: with it, hierarchical clustering gives one group of 386 (ARI 0.001); without it, the exact 10 sectors (ARI 1.000). With one month of data, the average ARI drops to 0.21.
- A Gaussian mixture on volatility, correlation, return and dispersion, fitted with the past, separates stress from calm, does not distinguish the rebound well and arrives 2 to 5 weeks late.
- A linear autoencoder is PCA (out-of-sample R² 0.4811 in both). The nonlinear one falls below (0.4408): with a linear structure, more flexibility is more noise.
- Everything unsupervised that feeds a model is fitted only with training data: standardization, components, groups and regimes.
## Resources
### Research
- [Price Discovery & Arbitrage Price Theory](https://pythonparatrading.com/en/arbitrage-pricing-theory-and-price-discovery-2/) · Research-0006 · 6 min
### Documentation
- [scikit-learn: decomposition (PCA) and clustering (k-means, hierarchical)](https://scikit-learn.org/stable/modules/clustering.html) · scikit-learn.org
- [scikit-learn: Gaussian mixture models](https://scikit-learn.org/stable/modules/mixture.html) · scikit-learn.org
- [Keras 3: getting started and choosing the backend](https://keras.io/getting_started/) · keras.io
### Books
- [Avellaneda, M. and Lee, J.-H. (2010): Statistical Arbitrage in the US Equities Market. Quantitative Finance 10(7), 761-782](https://doi.org/10.1080/14697680903124632) · doi.org
- [Kritzman, M., Li, Y., Page, S. and Rigobon, R. (2011): Principal Components as a Measure of Systemic Risk. Journal of Portfolio Management 37(4), 112-126](https://doi.org/10.3905/jpm.2011.37.4.112) · doi.org
- [Hinton, G. and Salakhutdinov, R. (2006): Reducing the Dimensionality of Data with Neural Networks. Science 313, 504-507](https://doi.org/10.1126/science.1127647) · doi.org
- [Gu, S., Kelly, B. and Xiu, D. (2021): Autoencoder Asset Pricing Models. Journal of Econometrics 222(1), 429-450](https://doi.org/10.1016/j.jeconom.2020.07.009) · doi.org
### Data
- [Level 5 data: panel prices (Parquet)](https://pythonparatrading.com/curso/datos/nivel-5/panel5_precios.parquet) · Download · panel5_precios.parquet
Previous lesson: [5.5.4 · Optuna and the risk of overfitting the validation](https://pythonparatrading.com/en/python-trading-course/level-5/optuna-validation-overfitting/)
Next lesson: [5.5.6 · Learning to rank for stock selection](https://pythonparatrading.com/en/python-trading-course/level-5/learning-to-rank-stock-selection/)