Aprendizaje no supervisado: PCA, clustering y autoencoders.

Estructura sin etiqueta

Las lecciones anteriores del módulo aprenden una función de las features a la etiqueta. El aprendizaje no supervisado no tiene etiqueta: busca estructura en los propios datos. En trading tiene cuatro usos habituales, y ninguno predice rendimientos por sí solo:

flowchart TB
  A[N rendimientos] --> B[PCA]
  A --> C[Grupos de activos]
  A --> D[Grupos de fechas]
  A --> E[Autoencoder]
  B --> F[Factores, residuos]
  C --> G[Diversificación]
  D --> H[Regímenes]
  E --> I[Anomalías]

• Factores estadísticos (PCA): resumir 400 series en unas pocas direcciones comunes; los residuos son la base del arbitraje estadístico (Avellaneda y Lee, 2010) y la matriz de covarianzas se limpia con ellos (lección 4.1.2).
• Grupos de activos: sectores «de verdad», no los de la clasificación, para diversificar o neutralizar una cartera (HRP, lección 4.1.4).
• Regímenes: agrupar fechas con un comportamiento parecido (lección 2.1.7 con modelos ocultos de Markov).
• Compresión: representar cada día con pocas variables y medir lo que no encaja (anomalías).

El problema de evaluar lo no supervisado es que no hay una métrica como el IC. Aquí sí la hay, porque los datos son sintéticos y el modelo que los generó está publicado (README de public/curso/datos/nivel-5/): el rendimiento diario de cada una de las 400 acciones es

un factor de mercado con una beta por acción, un factor por sector (10 sectores de 40 acciones), ruido propio con colas gruesas y efectos pequeños (alfa semanal, resultados). Cada método se juzga por lo que recupera de esa verdad.

Los datos: 400 rendimientos diarios

Se trabaja con los rendimientos diarios simples de los cierres, tipificados con la media y la desviación de 2011-2018, el periodo de ajuste. 2019-2025 queda para medir fuera de muestra. Tipificar con toda la muestra sería la fuga de la lección 5.2.4 en versión no supervisada.

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()       # solo con entrenamiento
Z = ((R - mu) / sd).to_numpy()
Ze, Zp = Z[ent], Z[~ent]
print("sesiones x acciones:", R.shape)
print("entrenamiento:", Ze.shape[0], " prueba:", Zp.shape[0])
print("sectores:", pd.Series(sector).value_counts().unique())

sesiones x acciones: (3771, 400)
entrenamiento: 2011  prueba: 1760
sectores: [40]

PCA y el borde de Marchenko-Pastur

PCA busca las direcciones de máxima varianza: los autovectores de la matriz de correlación, ordenados por su autovalor. La pregunta práctica es cuántos componentes son estructura y cuántos ruido. Si los rendimientos fueran ruido independiente, los autovalores de su matriz de correlación con $T$ sesiones y $N$ activos caerían, para $T$ y $N$ grandes, por debajo de

el borde superior de la distribución de Marchenko-Pastur. La lección 4.1.2 la usó para limpiar covarianzas; aquí sirve para contar factores.

C = np.corrcoef(Ze, rowvar=False)
val, vec = np.linalg.eigh(C)
val, vec = val[::-1], vec[:, ::-1]           # de mayor a menor
T, N = Ze.shape
borde = (1 + np.sqrt(N / T)) ** 2            # λ+ con varianza 1
print(f"N/T = {N / T:.3f}, borde de Marchenko-Pastur = {borde:.2f}")
print("12 primeros autovalores:", val[:12].round(2))
print(f"por encima del borde: {(val > borde).sum()}")
print(f"varianza explicada: PC1 {val[0] / N:.1%}, "
      f"PC1-PC10 {val[:10].sum() / N:.1%}")

N/T = 0.199, borde de Marchenko-Pastur = 2.09
12 primeros autovalores: [120.51   9.79   7.4    6.62   5.88   5.74   5.44   4.71   4.28   3.72
   1.37   1.35]
por encima del borde: 10
varianza explicada: PC1 30.1%, PC1-PC10 43.5%

La figura de la cabecera dibuja los veinte primeros. El primero, 120,5, se lleva el 30,1 % de la varianza: es el mercado. Del segundo al décimo bajan de 9,8 a 3,7 y quedan claramente por encima del borde (2,09). Del undécimo en adelante se amontonan entre 1,2 y 1,4: ruido. Diez componentes, que juntos explican el 43,5 % de la varianza.

El generador tiene once factores. El que falta no se ha perdido: el mercado mueve a todas las acciones en el mismo sentido, y la media de los diez factores sectoriales también, así que PCA funde esas dos direcciones en el primer componente y lo que queda distinguible son nueve contrastes entre sectores. Un componente principal es una dirección de varianza, no un factor del modelo uno a uno; poner nombre a cada uno exige mirar sus cargas.

Qué es cada componente

pc1 = vec[:, 0] * np.sign(vec[:, 0].sum())
carga = pc1 * sd.to_numpy()          # en unidades de rendimiento
f1 = Ze @ pc1                        # serie del componente
ew = R[ent].mean(axis=1).to_numpy()  # cartera equiponderada
print(f"corr(carga PC1, beta verdadera)        = "
      f"{np.corrcoef(carga, beta)[0, 1]:.3f}")
print(f"corr(carga PC1 sin desescalar, beta)   = "
      f"{np.corrcoef(pc1, beta)[0, 1]:.3f}")
print(f"corr(serie PC1, cartera equiponderada) = "
      f"{np.corrcoef(f1, ew)[0, 1]:.3f}")

corr(carga PC1, beta verdadera)        = 0.987
corr(carga PC1 sin desescalar, beta)   = 0.620
corr(serie PC1, cartera equiponderada) = 0.999

La serie del PC1 tiene una correlación de 0,999 con la cartera equiponderada. Sus cargas, una vez multiplicadas por la volatilidad de cada acción (la PCA de correlaciones trabaja con rendimientos tipificados), tienen una correlación de 0,987 con la beta verdadera; sin desescalar, solo 0,620, porque la carga de una correlación mide $\beta_i \sigma_m/\sigma_i$ y las acciones más volátiles quedan penalizadas. Es un detalle que se olvida a menudo al usar el PC1 como beta estadística.

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())

                     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

Cada uno de los componentes 2 a 10 carga fuerte sobre uno o dos sectores y poco sobre el resto: el PC2 es energía (0,46) contra todos; el PC3, finanzas (0,35) contra salud y telecomunicaciones; el PC4, salud (0,34). Los últimos mezclan más (el PC8 opone consumo básico, −0,25, a telecomunicaciones y tecnología). PCA no sabe de sectores; los encuentra porque son la estructura de la covarianza. En datos reales, las cargas mezclan sectores, estilos (value, momentum) y tamaño, y su interpretación cambia de una ventana a otra (lección 4.1.6).

Agrupar acciones

Si cada sector ocupa su propia dirección en el espacio de los componentes 2 a 10, un k-means sobre esas cargas (normalizadas para que cuente la dirección, no el tamaño) debería encontrar los grupos. Para medirlo, el índice de Rand ajustado (ARI): 1 si dos particiones coinciden (con cualquier nombre de grupo), alrededor de 0 si se parecen lo que se parecerían al azar.

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 (cargas PC2-PC10) frente a sectores: "
      f"{adjusted_rand_score(sector, km.labels_):.3f}")
print(pd.crosstab(km.labels_, sector).max(axis=1).tolist())

ARI k-means (cargas PC2-PC10) frente a sectores: 1.000
[40, 40, 40, 40, 40, 40, 40, 40, 40, 40]

ARI de 1,000: los diez grupos son los diez sectores, con 40 acciones cada uno. La alternativa más común en finanzas es el clustering jerárquico sobre la distancia de correlación $d_{ij} = \sqrt{2(1-\rho_{ij})}$, la misma que usa el HRP de la lección 4.1.4. La celda la prueba con y sin el primer componente:

from sklearn.cluster import AgglomerativeClustering

def jerarquico(Zv, quitar_mercado=True, k=10):
    """Clústeres por correlación (distancia 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"quitar PC1={quitar!s:5}: ARI "
          f"{adjusted_rand_score(sector, lab):.3f}, tamaños "
          f"{sorted(np.bincount(lab).tolist(), reverse=True)}")

quitar PC1=False: ARI 0.001, tamaños [386, 5, 2, 1, 1, 1, 1, 1, 1, 1]
quitar PC1=True : ARI 1.000, tamaños [40, 40, 40, 40, 40, 40, 40, 40, 40, 40]

Sin quitar nada, el resultado es inútil: un grupo de 386 acciones y nueve de una a cinco (ARI 0,001). Con una correlación media de 0,32 por el mercado, todas las acciones se parecen y las diferencias sectoriales quedan enterradas; el enlace medio va juntando todo en un único grupo. Quitando la proyección sobre el PC1, los residuos solo conservan la parte sectorial y el clustering recupera los diez sectores exactos.

Que salga perfecto dice más de los datos que del método: los sectores sintéticos son limpios (40 acciones, un factor cada uno, ninguna empresa a caballo entre dos) y hay ocho años de datos diarios. La celda siguiente repite el clustering con ventanas más cortas, en seis tramos distintos de 2011-2018:

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)

          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

Con un mes de datos (21 sesiones), el ARI medio es 0,21; con un trimestre, 0,39; con un año, 0,84; con dos, 0,97. Los grupos de un clustering hecho con pocos datos dependen de la ventana, y un HRP que se recalcula con 63 sesiones reorganiza su árbol cada vez. En datos reales, con estructura menos limpia y que además cambia, la longitud de la ventana es un compromiso entre ruido y actualidad.

El primer componente en las crisis

La varianza que se lleva el PC1 no es constante. Calculada trimestre a trimestre, mide cuánto se mueven las acciones a la vez: es la idea del absorption ratio de Kritzman et al. (2011) como medida de riesgo sistémico.

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 por trimestre: mediana {trim.median():.2f}, "
      f"mínimo {trim.min():.2f} ({trim.idxmin()}), "
      f"máximo {trim.max():.2f} ({trim.idxmax()})")
print(trim.nlargest(5).round(2).to_string())

PC1 por trimestre: mediana 0.31, mínimo 0.13 (2019Q2), máximo 0.68 (2020Q1)
2020Q1    0.68
2020Q2    0.68
2015Q4    0.53
2018Q4    0.53
2016Q1    0.52

La mediana es el 31 %, con un mínimo del 13 % en el segundo trimestre de 2019 y un máximo del 68 % en los dos primeros trimestres de 2020, el tramo de estrés más duro del generador (volatilidad del mercado del 80 % anual). Los demás picos (finales de 2015, finales de 2018, 2022) coinciden con los otros tramos de estrés. Es lo que se espera de cualquier cartera en una crisis: la diversificación entre acciones desaparece justo cuando hace falta (lección 4.2.4, stress testing de correlación). Fíjate en que los trimestres de 63 sesiones dan un borde de Marchenko-Pastur de unos 13: con ventanas cortas solo el mercado sobresale del ruido (el ejercicio 2 lo cuenta año a año).

Regímenes con clustering

Lo mismo que agrupa acciones puede agrupar fechas. Cada viernes se describe con cuatro features calculadas solo con el pasado: el logaritmo de la volatilidad del mercado en 21 sesiones, la correlación media aproximada (varianza del mercado entre la varianza media de las acciones, 63 sesiones), el rendimiento del mercado en 21 sesiones y el logaritmo de la dispersión entre acciones. Una mezcla de tres gaussianas se ajusta con 2012-2018 y clasifica todas las semanas; el generador publica el régimen verdadero de cada una (calma, estres, rebote) en panel5_coef_verdad.csv.

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")    # viernes de cada semana
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"]))

   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

Los números de grupo los pone el algoritmo; sus medias dicen qué es cada uno. El grupo 2 tiene la volatilidad más alta ($e^{-1,24} \approx$ 29 % anual), la correlación más alta (0,43) y rendimiento negativo: es el estrés. El grupo 1 es la calma tranquila (volatilidad del 11 %, correlación 0,19). El grupo 0 queda entre los dos. En 2019-2025, fuera de la muestra de ajuste:

• Estrés: 33 de sus 52 semanas caen en el grupo 2 y 15 en el intermedio.
• Calma: 190 de 263 en el grupo 1 y 72 en el intermedio; una sola en el de estrés.
• Rebote: ninguna en la calma, pero repartido entre el intermedio (36) y el estrés (15). Un rebote es volátil como el estrés y solo se distingue por el signo del rendimiento, que en tres meses de datos es ruidoso.

Agrupar fechas también tiene su versión del retraso: las features son ventanas hacia atrás, así que el régimen tiene que durar para moverlas.

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"estrés desde {ini:%Y-%m-%d}: el grupo lo marca el "
          f"{primera:%Y-%m-%d} ({(primera - ini).days // 7} semanas)")

estrés desde 2020-02-21: el grupo lo marca el 2020-03-06 (2 semanas)
estrés desde 2022-01-07: el grupo lo marca el 2022-02-11 (5 semanas)
estrés desde 2025-02-21: el grupo lo marca el 2025-03-14 (3 semanas)

El estrés de febrero de 2020 aparece en el grupo 2 dos semanas después de empezar; el de enero de 2022, cinco (era un estrés suave, 25 % de volatilidad); el de febrero de 2025, tres. Un detector sin información futura siempre llega tarde: la pregunta útil es cuánto y cuánto cuesta ese retraso en la estrategia que lo usa. La lección 2.1.7 hizo lo mismo con un HMM, que añade la persistencia de los estados, y la lección 5.10.2 vuelve a los regímenes como deriva de los datos.

Autoencoders

Un autoencoder es una red que aprende a reproducir su entrada pasando por un cuello estrecho: un codificador comprime los 400 rendimientos del día a $k$ números y un decodificador los reconstruye. Se entrena con error cuadrático, sin etiqueta.

flowchart TB
  A[400 rendimientos] --> B[Capa tanh 64]
  B --> C[Cuello de 10]
  C --> D[Capa tanh 64]
  D --> E[400 reconstruidos]
  E --> F[Error cuadrático]

Si todas las capas son lineales, el mejor autoencoder con un cuello de $k$ reconstruye desde el mismo subespacio que los $k$ primeros componentes principales: es PCA con otro algoritmo. La promesa de las capas no lineales es capturar estructura que PCA no ve. La celda entrena los dos en Keras 3 con backend PyTorch (la variable KERAS_BACKEND la fija el runner; la lección 5.6.1 explica los backends), con semilla (keras.utils.set_random_seed), algoritmos deterministas de torch y 60 épocas sobre 2011-2018. En producción se probarían más arquitecturas, regularización y parada temprana con una validación purgada.

import os
import keras
import torch

torch.use_deterministic_algorithms(True)
print("backend:", keras.backend.backend(), os.environ["KERAS_BACKEND"])

def r2(A, B):
    """Fracción de la varianza (tipificada) reconstruida."""
    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)

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

La medida es la fracción de la varianza tipificada que se reconstruye ($R^2$), dentro (2011-2018) y fuera (2019-2025):

• El autoencoder lineal iguala a PCA hasta el cuarto decimal fuera de muestra (0,4811 los dos): aprende el mismo subespacio, más despacio.
• El no lineal (capas tanh de 64) pierde poco dentro (0,4301 frente a 0,4352) y bastante fuera (0,4408 frente a 0,4811). No tiene nada no lineal que encontrar, porque el modelo verdadero es un modelo de factores lineal, y sus parámetros extra ajustan peculiaridades de 2011-2018 que no se repiten.

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)

             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

La figura pone los dos autoencoders sobre la curva de PCA. Con un cuello de 10 (los componentes sobre el borde), el no lineal reconstruye fuera poco más que PCA con 5 componentes (0,428). Y la curva de PCA tiene su propia lección: pasar de 10 a 50 componentes sube el $R^2$ dentro (de 0,435 a 0,550) más que fuera (de 0,481 a 0,541), porque esos 40 componentes extra son el ruido que el borde de Marchenko-Pastur ya había descartado.

Que el $R^2$ sea mayor fuera que dentro no es una fuga: 2019-2025 contiene 2020, cuando el mercado se llevó hasta el 68 % de la varianza, y un periodo con más correlación se reconstruye mejor con los mismos componentes. Por eso los métodos se comparan en el mismo periodo.

Errores típicos

• PCA sobre rendimientos sin tipificar o con toda la muestra. Sin tipificar, las acciones volátiles dominan los componentes; tipificando con toda la muestra, los datos de prueba entran en el ajuste.
• Contar factores a ojo. El borde de Marchenko-Pastur da una referencia; los componentes por debajo son ruido y añadirlos sube el ajuste dentro de muestra, no fuera.
• Clustering de correlaciones sin quitar el mercado. Sale un grupo gigante y unos pocos sueltos.
• Leer los números de grupo como etiquetas fijas. k-means o una mezcla pueden numerar los grupos de otra forma al reajustarlos; el grupo se identifica por sus medias.
• Regímenes ajustados con toda la muestra. Un detector que «acierta» los regímenes pasados con información de después no dice cuándo los habría detectado.
• Usar un autoencoder porque es más moderno. Si no mejora a PCA fuera de muestra con la misma dimensión, es PCA más caro y menos estable.

Ejercicios

1. Repite el clustering jerárquico sin quitar el mercado pero con 11 grupos (por si el mercado ocupa uno). ¿Mejora? Prueba después el enlace de Ward sobre los residuos sin el PC1.

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 sobre los residuos: ARI "
      f"{adjusted_rand_score(sector, lab_w):.3f}")

[385, 5, 2, 1, 1, 1, 1, 1, 1, 1, 1]
ARI: 0.001
Ward sobre los residuos: ARI 0.977

2. Calcula año a año cuántos autovalores superan el borde de Marchenko-Pastur y qué fracción de la varianza se lleva el PC1. ¿Qué pasa en 2020?

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})

      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

3. Un compañero calcula los cinco primeros componentes principales con los rendimientos de 2011-2025, usa la exposición de cada acción a cada componente como feature de un modelo de 5.5.3 y lo valida con walk-forward desde 2019. ¿Qué está mal?