Purged k-fold con embargo y CPCV para modelos de ML.

Qué cambia respecto a 2.5.3

La lección 2.5.3 trabajó con una serie: una fila por sesión, intervalos de etiqueta por posición (t0, t1) y una CPCV que elegía el mejor de ocho lookbacks. El ML sobre un panel complica las tres piezas:

El panel es el de la lección 5.4.1 (80 acciones, 725 semanas, etiqueta «por encima de la mediana de la semana a 21 sesiones») y el modelo, el mismo random forest. Datos sintéticos con la respuesta conocida: sirven para medir el método.

import zlib
from itertools import combinations
from math import comb
import numpy as np
import pandas as pd
from scipy.stats import norm
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import (TimeSeriesSplit,
                                     cross_val_score)

feat = pd.read_parquet("panel5_features.parquet")
etiq = pd.read_parquet("panel5_etiquetas.parquet")
d = feat.merge(etiq, on=["fecha", "ticker"])
completa = (d.groupby("fecha")[["ret_1s", "ret_sig_21d"]]
            .transform(lambda s: s.notna().all()).all(axis=1))
d = d[completa].copy()
d["med_21d"] = d.groupby("fecha")["ret_sig_21d"].transform("median")
tickers = sorted(d["ticker"].unique())
rng = np.random.default_rng([2029, zlib.crc32(b"n5_panel_cv")])
elegidas = [tickers[i] for i in
            np.sort(rng.choice(400, 80, replace=False))]
d = (d[d["ticker"].isin(elegidas)]
     .sort_values(["fecha", "ticker"]).reset_index(drop=True))

FEATS = [c for c in feat.columns if c not in
         ("fecha", "ticker", "sector", "revision_analistas")]
X = np.nan_to_num(d[FEATS].to_numpy(dtype=float))
X_fuga = np.nan_to_num(
    d[FEATS + ["revision_analistas"]].to_numpy(dtype=float))
y = (d["ret_sig_21d"] > d["med_21d"]).astype(int).to_numpy()

t0 = d["fecha"].to_numpy()        # se conoce la fila
t1 = d["fin_sig_21d"].to_numpy()  # termina su etiqueta
fechas = np.unique(t0)            # las 725 fechas de features
sesiones = pd.DatetimeIndex(pd.read_parquet(
    "panel5_precios.parquet", columns=["fecha"])["fecha"].unique())
sesiones = sesiones.sort_values()

n_ses = sesiones.searchsorted(t1) - sesiones.searchsorted(t0)
n_fechas = (np.searchsorted(fechas, t1, side="right")
            - np.searchsorted(fechas, t0))
print("sesiones por etiqueta:", np.unique(n_ses))
print("fechas de features dentro de cada etiqueta:")
print(pd.Series(n_fechas).value_counts().sort_index().to_string())
seis = np.unique(t0[n_fechas == 6])
print("etiquetas que cubren 6 fechas:", seis.astype("datetime64[D]"))

sesiones por etiqueta: [21]
fechas de features dentro de cada etiqueta:
1       80
2       80
3       80
4       80
5    57600
6       80
etiquetas que cubren 6 fechas: ['2024-12-20']

Todas las etiquetas miden 21 sesiones exactas. En fechas de features, casi todas (57.600) cubren la suya y las cuatro siguientes: la regla de «4 semanas» del README es correcta para ellas. Pero 80 filas, las del 20 de diciembre de 2024 (con Navidad y Año Nuevo en las semanas siguientes), cubren 6; y las últimas filas de la muestra cubren menos porque no hay más fechas de features detrás. Una purga escrita en semanas fallaría en silencio en el primer caso. Escrita con fin_sig_21d, no puede fallar.

Purga y embargo con fechas

La regla es la de 2.5.3, en fechas. Para cada tramo de prueba (un conjunto de fechas de features contiguas), su intervalo va del primer fecha al último fin_sig_21d de sus filas. Fuera del entrenamiento:

• Purga: toda fila cuyo intervalo [fecha, fin_sig_21d] se solape con el del tramo.
• Embargo: además, las filas que empiezan dentro de las embargo sesiones posteriores al fin del tramo. Su etiqueta ya no se solapa, pero sus features de ventana larga (ret_1m, ret_3m, mom_12_1) se calcularon con rendimientos que forman parte de las etiquetas de la prueba.

con $e$ el embargo en sesiones. Un conjunto de prueba puede tener varios tramos (en la CPCV, dos grupos separados), así que la función los separa primero.

def tramos(prueba):
    """(inicio, fin) de cada tramo de fechas contiguas de `prueba`."""
    f = np.unique(t0[prueba])
    pos = np.searchsorted(fechas, f)
    for bloque in np.split(f, np.flatnonzero(np.diff(pos) > 1) + 1):
        en = (t0[prueba] >= bloque[0]) & (t0[prueba] <= bloque[-1])
        yield t0[prueba][en].min(), t1[prueba][en].max()

def purgar(prueba, embargo=0):
    """Posiciones de entrenamiento: fuera la prueba, las filas que
    se solapan con un tramo y las del embargo tras cada tramo."""
    fuera = np.zeros(len(t0), bool)
    fuera[prueba] = True
    for ini, fin in tramos(prueba):
        j = min(sesiones.searchsorted(fin) + embargo,
                len(sesiones) - 1)
        fuera |= (t0 <= sesiones[j]) & (t1 >= ini)
    return np.flatnonzero(~fuera)

class PurgedKFold:
    """k-fold por bloques de fechas con purga y embargo. Tiene
    split() y get_n_splits(): vale como cv= en scikit-learn."""
    def __init__(self, n_splits=5, embargo=0):
        self.n_splits, self.embargo = n_splits, embargo

    def get_n_splits(self, X=None, y=None, groups=None):
        return self.n_splits

    def split(self, X=None, y=None, groups=None):
        for b in np.array_split(fechas, self.n_splits):
            prueba = np.flatnonzero(np.isin(t0, b))
            yield purgar(prueba, self.embargo), prueba

La clase no guarda los datos: lee t0, t1, fechas y sesiones del entorno, lo justo para la lección. En un proyecto, esos cuatro arrays irían en el constructor (como en el PurgedKFold del capítulo 7 de López de Prado, que recibe la serie de fines de etiqueta).

Comprobar la purga contra los datos

Una función que purga y otra que comprueba, escrita de otra forma, se vigilan entre sí. La comprobación no usa tramos: para cada fecha de prueba toma el fin de etiqueta más tardío de sus filas y busca filas de entrenamiento que se solapen con ese intervalo. Mide también el hueco del embargo: cuántas sesiones hay entre el último fin de etiqueta de la prueba y el inicio de la primera fila de entrenamiento posterior.

def comprobar(entreno, prueba):
    """(filas de entreno solapadas con la prueba, sesiones entre el
    último fin de etiqueta de la prueba y el siguiente entreno)."""
    fin_por_fecha = pd.Series(t1[prueba]).groupby(t0[prueba]).max()
    solapes = np.zeros(len(entreno), bool)
    for f, fin in fin_por_fecha.items():
        solapes |= (t0[entreno] <= fin) & (t1[entreno] >= f)
    fin_prueba = t1[prueba].max()
    despues = t0[entreno][t0[entreno] > fin_prueba]
    hueco = (sesiones.searchsorted(despues.min())
             - sesiones.searchsorted(fin_prueba)
             if len(despues) else None)
    return int(solapes.sum()), hueco

# la regla del README: 4 semanas a cada lado del bloque
semana = np.searchsorted(fechas, t0)
regla_4s = []
for b in np.array_split(np.arange(len(fechas)), 5):
    lejos = (semana < b[0] - 4) | (semana > b[-1] + 4)
    regla_4s.append((np.flatnonzero(lejos),
                     np.flatnonzero(np.isin(semana, b))))

filas = []
for nombre, div in [
        ("regla de 4 semanas", regla_4s),
        ("purga por fechas", list(PurgedKFold(5).split())),
        ("purga + embargo 21", list(PurgedKFold(5, 21).split()))]:
    for i, (ent, pru) in enumerate(div):
        sol, hueco = comprobar(ent, pru)
        filas.append((nombre, i, len(ent), sol, hueco))
pd.DataFrame(filas, columns=["esquema", "pliegue", "entreno",
                             "solapes", "hueco (ses.)"]
             ).set_index(["esquema", "pliegue"])

                            entreno  solapes  hueco (ses.)
esquema            pliegue
regla de 4 semanas 0          46080        0           4.0
                   1          45760        0           4.0
                   2          45760        0           3.0
                   3          45760        0           3.0
                   4          46080        0           NaN
purga por fechas   0          46080        0           4.0
                   1          45760        0           4.0
                   2          45760        0           3.0
                   3          45760        0           3.0
                   4          46080        0           NaN
purga + embargo 21 0          45760        0          23.0
                   1          45440        0          23.0
                   2          45440        0          22.0
                   3          45440        0          22.0
                   4          46080        0           NaN

La purga por fechas deja cero solapes en los cinco pliegues, y el embargo de 21 sesiones aparece en la última columna como el hueco que debe ser. La regla de 4 semanas del README también da cero aquí, porque ninguna frontera de los cinco bloques cae en la semana de los festivos: es correcta por suerte, no por construcción. Fíjate además en el pliegue 4, el último: no tiene entrenamiento después, así que el hueco no existe.

La figura muestra, para una acción, el inicio del tercer bloque de prueba (el pliegue 2, que empieza el 4 de agosto de 2017) con la purga por fechas: las cuatro semanas anteriores al bloque salen del entrenamiento porque su etiqueta llega a la primera fecha de prueba, y las cuatro posteriores, porque empiezan antes de que termine la última etiqueta de la prueba. Lo que se aparta es lo que se solapa, ni más ni menos.

Cuánto cuesta el embargo

El embargo no sale gratis: cada sesión de embargo es una semana de 80 filas menos para entrenar en cada frontera posterior. La ventana más larga de las features es la del momentum, 252 sesiones. Cubrirla entera sería lo más seguro y lo más caro:

sin_purga = sum(len(t0) - len(p) for _, p in PurgedKFold(5).split())
coste = {}
for e in (0, 5, 21, 63, 252):
    n = sum(len(ent) for ent, _ in PurgedKFold(5, e).split())
    coste[e] = {"filas de entreno": n,
                "perdidas": f"{1 - n / sin_purga:.1%}"}
pd.DataFrame(coste).T.rename_axis("embargo (sesiones)")

                   filas de entreno perdidas
embargo (sesiones)
0                            229440     1.1%
5                            229120     1.2%
21                           228160     1.7%
63                           225280     2.9%
252                          212800     8.3%

La purga sola cuesta un 1,1 % del entrenamiento. Un embargo de 21 sesiones (un mes) sube la pérdida al 1,7 %; uno de 63 (tres meses, la ventana de ret_3m), al 2,9 %; cubrir el año del momentum, al 8,3 %. Con más pliegues hay más fronteras y el coste se multiplica: en la CPCV de abajo, cada división tiene hasta dos tramos de prueba con su embargo. La recomendación de López de Prado, un embargo pequeño del orden del 1 % de la muestra, es un compromiso: con 725 semanas, unas 7 semanas. Aquí se usarán 63 sesiones, que cubren todas las ventanas salvo las de 252.

Lo que cambian la purga y el embargo en el acierto

Con los pliegues por bloques, la fuga gruesa del k-fold barajado ya no está (lección 5.4.1). Queda la de las fronteras: las primeras y las últimas semanas de cada bloque de prueba tienen vecinas solapadas en el entrenamiento si no se purga. Por eso, además del acierto de todo el pliegue, se mide el de las filas de frontera (las 4 primeras y las 4 últimas semanas de cada bloque).

def bosque():
    return RandomForestClassifier(
        n_estimators=60, min_samples_leaf=5, max_features="sqrt",
        max_samples=0.3, n_jobs=4, random_state=0)

bloques = [(np.flatnonzero(~np.isin(t0, b)),
            np.flatnonzero(np.isin(t0, b)))
           for b in np.array_split(fechas, 5)]
res = {}
for nombre, div in [("bloques sin purga", bloques),
                    ("purga", list(PurgedKFold(5).split())),
                    ("purga + embargo 63",
                     list(PurgedKFold(5, 63).split()))]:
    total, borde = [], []
    for ent, pru in div:
        ok = bosque().fit(X[ent], y[ent]).predict(X[pru]) == y[pru]
        s = semana[pru]
        en_borde = (s < s.min() + 4) | (s > s.max() - 4)
        total.append(ok.mean())
        borde.append(ok[en_borde].mean())
    res[nombre] = {"acierto": np.mean(total),
                   "acierto en fronteras": np.mean(borde)}
pd.DataFrame(res).T.map(lambda a: f"{a:.1%}")

                   acierto acierto en fronteras
bloques sin purga    50.0%                50.8%
purga                50.1%                49.2%
purga + embargo 63   50.1%                49.3%

En el conjunto de cada pliegue, las tres versiones dan lo mismo: un 50,0-50,1 %. En las fronteras, sin purga el bosque acierta el 50,8 % y con purga el 49,2 %: es la huella de la fuga que queda en los bordes, que el promedio diluye porque las fronteras son solo 8 de las 145 semanas de cada bloque. La diferencia es del orden del ruido de estas cifras (unas 12.800 filas de frontera, muy solapadas entre sí), así que no se puede leer como una medida precisa; lo que sí es seguro es la dirección, y que crece con el número de pliegues, porque hay más fronteras.

Un splitter propio dentro de scikit-learn

Como PurgedKFold tiene split y get_n_splits, cualquier función de scikit-learn que acepte cv= lo usa sin más: cross_val_score, cross_validate, GridSearchCV. Sirve para presentar la baseline: una regresión logística sobre las features en rango normal dentro de cada semana (el mismo tratamiento que la lección 5.2.4), con mucha regularización. Es rápida y lineal: si el bosque no la supera, no aporta nada.

rango = d.groupby("fecha")[FEATS].rank()
n_sem = d.groupby("fecha")[FEATS].transform("count")
G = np.nan_to_num(norm.ppf((rango - 0.5) / n_sem))

def logistica():
    return LogisticRegression(C=0.01, max_iter=500)

cv = PurgedKFold(5, embargo=63)
for nombre, modelo, Xm in [("logística (rangos)", logistica(), G),
                           ("random forest", bosque(), X)]:
    a = cross_val_score(modelo, Xm, y, cv=cv, scoring="accuracy")
    print(f"{nombre:<20}acierto {a.mean():.2%} "
          f"(pliegues {np.round(a, 3).tolist()})")

logística (rangos)  acierto 50.55% (pliegues [0.51, 0.503, 0.503, 0.505, 0.507])
random forest       acierto 50.07% (pliegues [0.502, 0.505, 0.503, 0.491, 0.503])

La logística acierta algo más que el bosque, con la misma validación. Las dos cifras están cerca del 50 %, como debe ser con un techo del 51,1 % (lección 5.4.1), y la diferencia, de unas décimas, es del orden del ruido. Para saber si es consistente hace falta más de un número: es lo que da la CPCV.

CPCV para un modelo de machine learning

La CPCV de 2.5.3 elegía un parámetro con los datos de entrenamiento y aplicaba la variante elegida a los grupos de prueba. Con un modelo, el procedimiento es:

1. Partir las fechas en N grupos contiguos.
2. Para cada una de las $\binom{N}{k}$ combinaciones de k grupos de prueba: purgar y embargar alrededor de cada grupo (la función purgar ya separa los tramos), entrenar un modelo nuevo con el resto y guardar sus predicciones para cada grupo de prueba por separado.
3. Recomponer $\binom{N-1}{k-1}$ caminos: para cada grupo, la predicción de su j-ésima aparición en la prueba.
4. Medir cada camino con las métricas del problema: aquí el acierto y el IC (correlación de rangos entre la probabilidad y el rendimiento a 21 sesiones, semana a semana, y su media).

flowchart TB
  F[Fechas en N grupos] --> C[C N k divisiones]
  C --> P[Purga y embargo]
  P --> E[Entrenar un modelo]
  E --> G[Predecir grupos]
  G --> R[Recomponer caminos]
  R --> M[IC por camino]

Con N = 6 y k = 2 salen 15 divisiones, 15 modelos y 5 caminos; cada modelo entrena con unas 4/6 partes de las fechas, menos la purga y el embargo alrededor de los dos grupos de prueba.

N, K = 6, 2
grupos = np.array_split(fechas, N)
divisiones = list(combinations(range(N), K))
fila_grupo = [np.flatnonzero(np.isin(t0, g)) for g in grupos]

def cpcv(modelo, Xm, embargo=63):
    """{(división, grupo): probabilidades de las filas del grupo}."""
    pred = {}
    for i, prueba_g in enumerate(divisiones):
        prueba = np.concatenate([fila_grupo[g] for g in prueba_g])
        ent = purgar(prueba, embargo)
        m = modelo().fit(Xm[ent], y[ent])
        for g in prueba_g:
            pred[(i, g)] = m.predict_proba(Xm[fila_grupo[g]])[:, 1]
    return pred

def ic_semanal(p, filas):
    """Media de la correlación de rangos semanal entre p y el
    rendimiento a 21 sesiones."""
    df = pd.DataFrame({"s": semana[filas], "p": p,
                       "r": d["ret_sig_21d"].to_numpy()[filas]})
    rk = df.groupby("s")[["p", "r"]].rank()
    return rk.groupby(df["s"]).corr().xs("p", level=1)["r"].mean()

def caminos(pred):
    filas = []
    for j in range(comb(N - 1, K - 1)):
        p = np.concatenate([
            pred[([i for i, dv in enumerate(divisiones)
                   if g in dv][j], g)] for g in range(N)])
        todas = np.concatenate(fila_grupo)
        filas.append({"acierto": ((p > 0.5) == y[todas]).mean(),
                      "IC": ic_semanal(p, todas)})
    return pd.DataFrame(filas).rename_axis("camino")

tabla = pd.concat({
    "logística": caminos(cpcv(logistica, G)),
    "bosque": caminos(cpcv(bosque, X)),
    "bosque + fuga": caminos(cpcv(bosque, X_fuga))}, axis=1)
tabla.round(4)

       logística          bosque         bosque + fuga
         acierto      IC acierto      IC       acierto      IC
camino
0         0.5041  0.0123  0.5009 -0.0011        0.5126  0.0442
1         0.5023  0.0105  0.5000  0.0008        0.5145  0.0451
2         0.5069  0.0190  0.5026  0.0040        0.5167  0.0503
3         0.5048  0.0148  0.5025  0.0137        0.5162  0.0499
4         0.5037  0.0125  0.5005 -0.0008        0.5123  0.0432

Para leer la tabla falta la referencia: el IC del oráculo (el alfa verdadero frente al mismo rendimiento a 21 sesiones) en las mismas 80 acciones. Como todos los caminos cubren todas las fechas, es el mismo para los cinco.

verdad = pd.read_parquet("panel5_verdad.parquet",
                         columns=["fecha", "ticker", "alfa_verdadero"])
alfa = d[["fecha", "ticker"]].merge(verdad)["alfa_verdadero"]
ic_oraculo = ic_semanal(alfa.to_numpy(), np.arange(len(d)))
print(f"IC del oráculo: {ic_oraculo:.4f}")
dif = tabla["bosque"]["IC"] - tabla["logística"]["IC"]
print("bosque menos logística, por camino:",
      dif.round(4).tolist())
print(f"caminos en los que gana la logística: {(dif < 0).sum()}/5")

IC del oráculo: 0.0270
bosque menos logística, por camino: [-0.0135, -0.0096, -0.015, -0.001, -0.0134]
caminos en los que gana la logística: 5/5

La figura de la cabecera pone los quince números en una escala. Tres lecturas:

• La baseline gana. La logística sobre rangos tiene un IC de 0,011 a 0,019 por camino y el bosque, de −0,001 a 0,014. La comparación es pareada (mismas divisiones, mismas filas de prueba) y la logística queda por encima en los cinco caminos. Con esta relación señal/ruido, la flexibilidad del bosque se gasta en ruido; la lección 5.5.2 vuelve sobre cómo regularizarlo.
• Nadie pasa del techo. El oráculo da el IC máximo alcanzable; los dos modelos honestos quedan por debajo, como debe ser. La logística recoge la mitad, más o menos.
• La fuga gana en todos los caminos. Con revision_analistas, el bosque tiene un IC de 0,043 a 0,050 y un acierto del 51,2-51,7 %, por encima del oráculo en las cinco. La CPCV, con purga y embargo, no lo puede ver: la información del futuro está dentro de la feature (lección 5.2.2).

Hiperparámetros dentro de la CPCV

Si los hiperparámetros del bosque se eligieran mirando los caminos de la CPCV, los caminos dejarían de ser fuera de muestra: serían el conjunto con el que se eligió. Hay dos salidas correctas:

• Validación anidada: dentro de cada división de la CPCV, un PurgedKFold sobre el entrenamiento elige los hiperparámetros (por ejemplo, con GridSearchCV(cv=PurgedKFold(…))), y el modelo elegido predice los grupos de prueba. Cuesta (número de divisiones) × (combinaciones) × (pliegues internos) entrenamientos.
• Contar los intentos: si se prueban varias configuraciones y se elige la mejor por los caminos, se registran todas y el resultado se descuenta con el Deflated Sharpe o la PBO (lecciones 2.5.4 y 5.4.4).

Y los pesos por unicidad de la lección 5.3.4 entran en el mismo bucle: fit(X[ent], y[ent], sample_weight=w[ent]), con los pesos calculados solo con las etiquetas de entrenamiento.

Ejercicios

1. Comprueba TimeSeriesSplit con comprobar: cinco pliegues sin hueco y con gap=4 * 80. ¿Cuántas filas de entrenamiento se solapan con la prueba en cada caso?

for gap in (0, 4 * 80):
    tss = TimeSeriesSplit(5, gap=gap)
    sol = [comprobar(a, b)[0] for a, b in tss.split(X)]
    print(f"gap={gap:>3}: solapes por pliegue {sol}")

gap=  0: solapes por pliegue [390, 376, 362, 348, 334]
gap=320: solapes por pliegue [70, 56, 42, 28, 14]

2. Calcula, sin entrenar nada, cuántos modelos y caminos salen con N = 8 y k = 2, y qué fracción media de las filas queda para entrenar tras purgar y embargar 63 sesiones alrededor de los dos grupos de prueba.

N8, K8 = 8, 2
grupos8 = np.array_split(fechas, N8)
fracs = []
for dv in combinations(range(N8), K8):
    prueba = np.flatnonzero(np.isin(t0, np.concatenate(
        [grupos8[g] for g in dv])))
    fracs.append(len(purgar(prueba, 63)) / len(t0))
print(f"modelos {comb(N8, K8)}, caminos {comb(N8 - 1, K8 - 1)}, "
      f"entreno medio {np.mean(fracs):.1%} "
      f"(sin purga {(N8 - K8) / N8:.1%})")

modelos 28, caminos 7, entreno medio 70.7% (sin purga 75.0%)

3. Un compañero elige el número de árboles y el tamaño de hoja del bosque con los cinco caminos de la CPCV (prueba doce combinaciones y se queda con la de mejor IC medio) y presenta ese IC como resultado fuera de muestra. ¿Qué está mal y cómo lo arreglarías?