Purged k-fold with embargo and CPCV for ML models.

What changes with respect to 2.5.3

Lesson 2.5.3 worked with a single series: one row per session, label intervals by position (t0, t1) and a CPCV that chose the best of eight lookbacks. ML on a panel complicates all three pieces:

The panel is the one of lesson 5.4.1 (80 stocks, 725 weeks, label “above the week's median at 21 sessions”) and the model, the same random forest. Synthetic data with a known answer: they are for measuring the method.

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()        # the row becomes known
t1 = d["fin_sig_21d"].to_numpy()  # its label ends
fechas = np.unique(t0)            # the 725 feature dates
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("sessions per label:", np.unique(n_ses))
print("feature dates inside each label:")
print(pd.Series(n_fechas).value_counts().sort_index().to_string())
seis = np.unique(t0[n_fechas == 6])
print("labels covering 6 dates:", seis.astype("datetime64[D]"))

sessions per label: [21]
feature dates inside each label:
1       80
2       80
3       80
4       80
5    57600
6       80
labels covering 6 dates: ['2024-12-20']

All labels measure exactly 21 sessions. In feature dates, almost all of them (57,600) cover their own and the next four: the README's “4 weeks” rule is correct for them. But 80 rows, those of December 20, 2024 (with Christmas and New Year in the following weeks), cover 6; and the last rows of the sample cover fewer because there are no more feature dates after them. A purge written in weeks would fail silently in the first case. Written with fin_sig_21d, it cannot fail.

Purge and embargo with dates

The rule is that of 2.5.3, in dates. For each test segment (a set of contiguous feature dates), its interval runs from the first fecha to the last fin_sig_21d of its rows. Out of training:

• Purge: every row whose interval [fecha, fin_sig_21d] overlaps that of the segment.
• Embargo: in addition, the rows that start within the embargo sessions after the end of the segment. Their label no longer overlaps, but their long-window features (ret_1m, ret_3m, mom_12_1) were computed with returns that are part of the test labels.

with $e$ the embargo in sessions. A test set can have several segments (in the CPCV, two separate groups), so the function separates them first.

def tramos(prueba):
    """(start, end) of each segment of contiguous dates in `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):
    """Training positions: out go the test, the rows that overlap
    a segment and those in the embargo after each segment."""
    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 by date blocks with purge and embargo. It has
    split() and get_n_splits(): it works as cv= in 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

The class does not store the data: it reads t0, t1, fechas and sesiones from the environment, just enough for the lesson. In a project, those four arrays would go in the constructor (as in the PurgedKFold of López de Prado's chapter 7, which receives the series of label ends).

Checking the purge against the data

A function that purges and another that checks, written in a different way, keep watch on each other. The check does not use tramos: for each test date it takes the latest label end of its rows and looks for training rows that overlap that interval. It also measures the embargo gap: how many sessions there are between the last label end of the test and the start of the first later training row.

def comprobar(entreno, prueba):
    """(training rows overlapping the test, sessions between the
    last label end of the test and the next training row)."""
    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

# the README's rule: 4 weeks on each side of the block
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 [
        ("4-week rule", regla_4s),
        ("purge by dates", list(PurgedKFold(5).split())),
        ("purge + 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", "gap (ses.)"]
             ).set_index(["esquema", "pliegue"])

                            entreno  solapes  gap (ses.)
esquema            pliegue
4-week rule        0          46080        0         4.0
                   1          45760        0         4.0
                   2          45760        0         3.0
                   3          45760        0         3.0
                   4          46080        0         NaN
purge by dates     0          46080        0         4.0
                   1          45760        0         4.0
                   2          45760        0         3.0
                   3          45760        0         3.0
                   4          46080        0         NaN
purge + 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

The purge by dates leaves zero overlaps in the five folds, and the 21-session embargo shows up in the last column as the gap it should be. The README's 4-week rule also gives zero here, because no boundary of the five blocks falls in the holiday week: it is correct by luck, not by construction. Note also fold 4, the last one: it has no training after it, so the gap does not exist.

The figure shows, for one stock, the start of the third test block (fold 2, which starts on August 4, 2017) with the purge by dates: the four weeks before the block leave training because their label reaches the first test date, and the four after it, because they start before the last test label ends. What is set aside is what overlaps, no more and no less.

How much the embargo costs

The embargo is not free: each session of embargo is one week of 80 rows fewer to train on at each later boundary. The longest feature window is that of momentum, 252 sessions. Covering it in full would be the safest and the most expensive option:

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] = {"training rows": n,
                "perdidas": f"{1 - n / sin_purga:.1%}"}
pd.DataFrame(coste).T.rename_axis("embargo (sessions)")

                   training rows perdidas
embargo (sessions)
0                         229440     1.1%
5                         229120     1.2%
21                        228160     1.7%
63                        225280     2.9%
252                       212800     8.3%

The purge alone costs 1.1% of training. A 21-session embargo (one month) raises the loss to 1.7%; one of 63 (three months, the window of ret_3m), to 2.9%; covering the year of momentum, to 8.3%. With more folds there are more boundaries and the cost multiplies: in the CPCV below, each split has up to two test segments with their embargo. López de Prado's recommendation, a small embargo on the order of 1% of the sample, is a compromise: with 725 weeks, about 7 weeks. Here 63 sessions will be used, which cover all windows except the 252-session ones.

What purge and embargo change in the hit rate

With block folds, the gross leak of shuffled k-fold is gone (lesson 5.4.1). What remains is the leak at the boundaries: the first and last weeks of each test block have overlapping neighbors in training if there is no purge. That is why, in addition to the hit rate of the whole fold, the hit rate of the boundary rows (the first 4 and the last 4 weeks of each block) is measured.

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 [("blocks, no purge", bloques),
                    ("purge", list(PurgedKFold(5).split())),
                    ("purge + 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),
                   "hit rate at boundaries": np.mean(borde)}
pd.DataFrame(res).T.map(lambda a: f"{a:.1%}")

                   acierto hit rate at boundaries
blocks, no purge     50.0%                  50.8%
purge                50.1%                  49.2%
purge + embargo 63   50.1%                  49.3%

Over each fold as a whole, the three versions give the same result: 50.0-50.1%. At the boundaries, without purge the forest's hit rate is 50.8% and with purge 49.2%: it is the trace of the leak that remains at the edges, which the average dilutes because the boundaries are only 8 of the 145 weeks of each block. The difference is on the order of the noise of these figures (about 12,800 boundary rows, heavily overlapping with each other), so it cannot be read as a precise measurement; what is certain is the direction, and that it grows with the number of folds, because there are more boundaries.

A custom splitter inside scikit-learn

Since PurgedKFold has split and get_n_splits, any scikit-learn function that accepts cv= uses it as is: cross_val_score, cross_validate, GridSearchCV. It serves to introduce the baseline: a logistic regression on the features as normal ranks within each week (the same treatment as lesson 5.2.4), with heavy regularization. It is fast and linear: if the forest does not beat it, the forest adds nothing.

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 [("logistic (ranks)", logistica(), G),
                           ("random forest", bosque(), X)]:
    a = cross_val_score(modelo, Xm, y, cv=cv, scoring="accuracy")
    print(f"{nombre:<20}hit rate {a.mean():.2%} "
          f"(folds {np.round(a, 3).tolist()})")

logistic (ranks)    hit rate 50.55% (folds [0.51, 0.503, 0.503, 0.505, 0.507])
random forest       hit rate 50.07% (folds [0.502, 0.505, 0.503, 0.491, 0.503])

The logistic regression gets somewhat more right than the forest, with the same validation. Both figures are close to 50%, as they should be with a ceiling of 51.1% (lesson 5.4.1), and the difference, a few tenths, is on the order of the noise. To know whether it is consistent, more than one number is needed: that is what CPCV provides.

CPCV for a machine learning model

The CPCV of 2.5.3 chose a parameter with the training data and applied the chosen variant to the test groups. With a model, the procedure is:

1. Split the dates into N contiguous groups.
2. For each of the $\binom{N}{k}$ combinations of k test groups: purge and embargo around each group (the purgar function already separates the segments), train a new model on the rest and store its predictions for each test group separately.
3. Reassemble $\binom{N-1}{k-1}$ paths: for each group, the prediction from its j-th appearance in the test.
4. Measure each path with the problem's metrics: here the hit rate and the IC (rank correlation between the probability and the 21-session return, week by week, and its mean).

flowchart TB
  F[Dates in N groups] --> C[C N k splits]
  C --> P[Purge and embargo]
  P --> E[Train a model]
  E --> G[Predict groups]
  G --> R[Reassemble paths]
  R --> M[IC per path]

With N = 6 and k = 2 there are 15 splits, 15 models and 5 paths; each model trains on about 4/6 of the dates, minus the purge and embargo around the two test groups.

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):
    """{(split, group): probabilities of the group's rows}."""
    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):
    """Mean of the weekly rank correlation between p and the
    21-session return."""
    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({
    "logistic": caminos(cpcv(logistica, G)),
    "forest": caminos(cpcv(bosque, X)),
    "forest + leak": caminos(cpcv(bosque, X_fuga))}, axis=1)
tabla.round(4)

       logistic          forest         forest + leak
        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

To read the table the reference is missing: the oracle's IC (the true alpha against the same 21-session return) on the same 80 stocks. Since all paths cover all dates, it is the same for all five.

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"oracle IC: {ic_oraculo:.4f}")
dif = tabla["forest"]["IC"] - tabla["logistic"]["IC"]
print("forest minus logistic, per path:",
      dif.round(4).tolist())
print(f"paths won by the logistic regression: {(dif < 0).sum()}/5")

oracle IC: 0.0270
forest minus logistic, per path: [-0.0135, -0.0096, -0.015, -0.001, -0.0134]
paths won by the logistic regression: 5/5

The header figure puts the fifteen numbers on one scale. Three readings:

• The baseline wins. The logistic regression on ranks has an IC of 0.011 to 0.019 per path and the forest, of −0.001 to 0.014. The comparison is paired (same splits, same test rows) and the logistic regression comes out ahead on all five paths. With this signal-to-noise ratio, the forest's flexibility is spent on noise; lesson 5.5.2 comes back to how to regularize it.
• Nobody goes past the ceiling. The oracle gives the maximum achievable IC; the two honest models stay below it, as they should. The logistic regression captures about half of it.
• The leak wins on every path. With revision_analistas, the forest has an IC of 0.043 to 0.050 and a hit rate of 51.2-51.7%, above the oracle on all five. CPCV, with purge and embargo, cannot see it: the information from the future is inside the feature (lesson 5.2.2).

Hyperparameters inside CPCV

If the forest's hyperparameters were chosen by looking at the CPCV paths, the paths would stop being out of sample: they would be the set used to choose. There are two correct ways out:

• Nested validation: inside each CPCV split, a PurgedKFold on the training data chooses the hyperparameters (for example, with GridSearchCV(cv=PurgedKFold(…))), and the chosen model predicts the test groups. It costs (number of splits) × (combinations) × (inner folds) trainings.
• Counting the trials: if several configurations are tried and the best is chosen by the paths, all of them are recorded and the result is discounted with the Deflated Sharpe or the PBO (lessons 2.5.4 and 5.4.4).

And the uniqueness weights of lesson 5.3.4 go into the same loop: fit(X[ent], y[ent], sample_weight=w[ent]), with the weights computed only with the training labels.

Exercises

1. Check TimeSeriesSplit with comprobar: five folds without a gap and with gap=4 * 80. How many training rows overlap the test in each case?

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}: overlaps per fold {sol}")

gap=  0: overlaps per fold [390, 376, 362, 348, 334]
gap=320: overlaps per fold [70, 56, 42, 28, 14]

2. Compute, without training anything, how many models and paths there are with N = 8 and k = 2, and what average fraction of the rows is left for training after purging and embargoing 63 sessions around the two test groups.

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"models {comb(N8, K8)}, paths {comb(N8 - 1, K8 - 1)}, "
      f"mean training {np.mean(fracs):.1%} "
      f"(no purge {(N8 - K8) / N8:.1%})")

models 28, paths 7, mean training 70.7% (no purge 75.0%)

3. A colleague chooses the number of trees and the leaf size of the forest with the five CPCV paths (tries twelve combinations and keeps the one with the best mean IC) and presents that IC as an out-of-sample result. What is wrong and how would you fix it?