Combinatorial purged cross-validation for strategies.

What k-fold assumes

K-fold cross-validation is the standard way to estimate how a model will behave on data it has not seen: the data is split into k chunks, the model is trained on k − 1 and tested on the remaining one, k times. It works because it assumes the observations are independent and identically distributed: it does not matter which rows fall in each chunk. The scikit-learn documentation says so explicitly and warns that in time series that assumption does not hold.

In a trading problem the assumption fails on two fronts at once:

• Overlapping labels. The label of an observation usually spans several sessions: “will the price rise over the next 60 days?” Monday's label and Tuesday's share 59 of their 60 days.
• Overlapping features. Features are rolling windows: Monday's 120-day return and Tuesday's share 119 days. Two consecutive days have almost identical features.

Together, the two things mean that a test observation has almost identical neighbors in the training set with almost the same answer. A flexible model does not need to predict: recognizing them is enough. The result is an accuracy that does not exist.

flowchart TB
  F[Overlapping features] --> V[Identical neighbors]
  E[Long labels] --> R[Shared answer]
  V --> M[The model recognizes]
  R --> M
  M --> P[Fictitious accuracy]

Measuring the leak on a series with nothing to predict

To measure the leak without a real effect hiding it, you need a series in which you know there is nothing to predict. The ruptura series of regimenes2.parquet is exactly that: independent normal returns with a single structural break on May 22, 2013. From that day on, daily mean 0 and volatility 0.013: pure noise, with no memory. Any out-of-sample AUC other than 0.5 is luck or leakage.

The label is the sign of the return of the next 60 sessions; the features, four returns and two rolling volatilities. Each observation also carries the interval of its label, t0 and t1 (positions of the first and last session it uses), which is what purging needs:

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import KFold

reg = pd.read_parquet("regimenes2.parquet")
r = (reg[reg["serie"] == "ruptura"].set_index("fecha")
     ["rendimiento_log"].loc["2013-05-22":])

H = 60  # label horizon, in sessions
rasgos = pd.DataFrame({f"r{n}": r.rolling(n).sum()
                       for n in (5, 20, 60, 120)})
rasgos["v20"] = r.rolling(20).std()
rasgos["v60"] = r.rolling(60).std()
futuro = r[::-1].rolling(H).sum()[::-1].shift(-1)  # r from t+1 to t+H
datos = rasgos.assign(futuro=futuro).dropna()

X = datos[rasgos.columns].to_numpy()
y = (datos["futuro"] > 0).astype(int).to_numpy()
pos = np.arange(len(datos))
t0, t1 = pos + 1, pos + H  # sessions used by each label
print(len(datos), "observations; share of 1 labels:",
      round(y.mean(), 3))

2994 observations; share of 1 labels: 0.586

The model is a random forest, flexible on purpose: it is the kind of model that best “remembers” neighbors. The auc function trains on some positions and measures the AUC (area under the ROC curve: 0.5 is chance, 1 is perfect) on others:

def modelo():
    return RandomForestClassifier(n_estimators=100, max_features=2,
                                  random_state=0, n_jobs=-1)

def auc(entreno, prueba):
    m = modelo().fit(X[entreno], y[entreno])
    return roc_auc_score(y[prueba], m.predict_proba(X[prueba])[:, 1])

for barajar in (True, False):
    kf = KFold(10, shuffle=barajar,
               random_state=0 if barajar else None)
    res = [auc(a, b) for a, b in kf.split(X)]
    print(f"k-fold {'shuffled' if barajar else 'unshuffled'}: "
          f"mean AUC {np.mean(res):.3f}")

k-fold shuffled: mean AUC 0.928
k-fold unshuffled: mean AUC 0.532

Shuffled, the k-fold “predicts” the sign of the next 60 days of a noise series with an AUC close to 0.9. It is leakage in its purest form: each test row has its neighbor from yesterday and from tomorrow in the training set. Without shuffling, the folds are contiguous blocks of time and most of the leak disappears: the AUC stays at 0.53, close to chance. That 0.53 has a lot of dispersion: with 3,000 overlapping observations, the number of independent labels is about 50 (3,000 / 60), and with so few the AUC of a model with no information can come out several points above or below 0.5.

But without shuffling there is still leakage, at the boundaries of each block. The last training observation before a test block has a label that spans the first 60 days of the block; the first one after it has features computed with the last days of the block. Purging and embargo are there to remove it.

Purging and embargo

• Purging. Remove from the training set every observation whose label $[t_0, t_1]$ overlaps the interval covered by the labels of the test block. Before the block, these are the last H − 1 observations; after it, the first H − 1.
• Embargo. Also remove the observations that start right after the test block, for a fixed number of sessions. Their labels no longer overlap, but their features (windows of up to 120 sessions) were computed with test data. López de Prado proposes a small embargo, on the order of 1% of the sample; with long features it is worth covering the longest window.

The full scheme of a split, with the test in the middle of the series, is that of one row of the figure at the top: training, purge, test, purge, embargo, training.

def purgar(t0, t1, bloques, embargo=0):
    """Training positions: out go the observations whose label
    [t0, t1] overlaps that of any test block, and those that start
    in the `embargo` days after each block."""
    fuera = np.zeros(len(t0), bool)
    for b in bloques:
        ini, fin = t0[b].min(), t1[b].max()
        fuera |= (t0 <= fin) & (t1 >= ini)              # purging
        fuera |= (t0 > fin) & (t0 <= fin + embargo)     # embargo
    return np.flatnonzero(~fuera)

bloque = np.arange(1000, 1300)
entreno = purgar(t0, t1, [bloque], embargo=120)
quitadas = np.setdiff1d(np.setdiff1d(pos, entreno), bloque)
print(f"Test: {len(bloque)}; removed from training: "
      f"{len(quitadas)} (from {quitadas.min()} to {quitadas.max()})")

Test: 300; removed from training: 238 (from 941 to 1478)

A test block of 300 observations forces removing 59 before (purging), 59 after (purging) and 120 more (embargo, the window of the longest feature). Now the comparison: the same unshuffled k-fold, with and without purging and embargo, measuring the AUC separately on the boundary observations (the first 60 and the last 60 of each block, which are the ones with leaked neighbors) and on the center ones:

def k_fold(k, embargo, purga):
    p_borde, y_borde, p_centro, y_centro = [], [], [], []
    for b in np.array_split(pos, k):
        if purga:
            entreno = purgar(t0, t1, [b], embargo)
        else:
            entreno = np.setdiff1d(pos, b)
        m = modelo().fit(X[entreno], y[entreno])
        p = m.predict_proba(X[b])[:, 1]
        borde = (b < b[0] + H) | (b > b[-1] - H)
        p_borde += list(p[borde]); y_borde += list(y[b][borde])
        p_centro += list(p[~borde]); y_centro += list(y[b][~borde])
    return (roc_auc_score(y_borde, p_borde),
            roc_auc_score(y_centro, p_centro))

filas = {"no purging": k_fold(10, 0, False),
         "purging and embargo": k_fold(10, 120, True)}
pd.DataFrame(filas, index=["boundary AUC", "center AUC"]).T.round(3)

                     boundary AUC  center AUC
no purging                  0.558       0.499
purging and embargo         0.497       0.459

Without purging, the boundary observations have an AUC of 0.56, above the 0.50 of the center ones: that is the leak from the overlaps. With purging and embargo, the boundaries drop to 0.50. The difference is of the order of the chance dispersion of this data (about 50 independent labels), so the conclusive test is the shuffled k-fold; the boundary test shows where the leak is when you do not shuffle. It looks small averaged over the whole fold, but in a problem with a weak signal, which is the usual case, a leak of a few AUC points is enough to turn a useless model into one that looks profitable. And it grows with the label horizon and with the number of folds, because there are more boundaries.

From k-fold to combinatorial validation

A purged k-fold gives k estimates, one per fold, which together form one out-of-sample path: each observation is tested once. The walk-forward of lesson 2.5.2 also gives a single path. Combinatorial purged cross-validation (CPCV, López de Prado, 2018) breaks out of that limitation:

1. Split the data into N contiguous groups.
2. Test every combination of k groups as the test set: $\binom{N}{k}$ splits. In each one, train on the remaining N − k groups, purged and embargoed around each test group.
3. Each group appears in the test in $\binom{N-1}{k-1}$ splits. Taking, for each group, its j-th appearance, you can reconstruct complete paths that cover the whole out-of-sample history.

flowchart TB
  D[Data] --> G[N groups]
  G --> C[C N k splits]
  C --> P[Purging and embargo]
  P --> E[Choose on the rest]
  E --> T[Test k groups]
  T --> R[Rebuild paths]

The number of paths is

because there are $\binom{N}{k}$ splits with k test groups each, $k\binom{N}{k}$ tested groups in total, and each path needs N of them. The figure at the top draws the case N = 6, k = 2: 15 splits and 5 paths, with the number of the path each test group goes to.

from itertools import combinations
from math import comb

tabla = {(N, k): {"splits": comb(N, k),
                  "paths": comb(N - 1, k - 1),
                  "train (fraction)": round((N - k) / N, 2)}
         for N in (6, 8, 10) for k in (1, 2, 3)}
pd.DataFrame(tabla).T.rename_axis(["N", "k"])

      splits  paths  train (fraction)
N  k
6  1     6.0    1.0              0.83
   2    15.0    5.0              0.67
   3    20.0   10.0              0.50
8  1     8.0    1.0              0.88
   2    28.0    7.0              0.75
   3    56.0   21.0              0.62
10 1    10.0    1.0              0.90
   2    45.0    9.0              0.80
   3   120.0   36.0              0.70

With k = 1, CPCV is a purged k-fold: N splits, one path. Raising k gives more paths, but at the cost of training on less data (with N = 6 and k = 3 you train on half). N = 10 and k = 2 is a reasonable compromise: 45 splits, 9 paths and 80% of the data for training.

The two CPCV functions, generic: one generates the splits with their purged training set; the other reconstructs the paths from the out-of-sample returns of each (split, group).

def cpcv(n, N, k, t0, t1, embargo):
    """For each combination of k test groups: (index of the
    split, test groups, training positions)."""
    grupos = np.array_split(np.arange(n), N)
    for d, prueba in enumerate(combinations(range(N), k)):
        bloques = [grupos[g] for g in prueba]
        yield d, prueba, grupos, purgar(t0, t1, bloques, embargo)

def caminos(N, k, oos):
    """oos[(split, group)] = returns of that group when it was
    tested in that split. Returns the complete paths."""
    divisiones = list(combinations(range(N), k))
    salida = []
    for j in range(comb(N - 1, k - 1)):
        partes = []
        for g in range(N):
            con_g = [d for d, p in enumerate(divisiones) if g in p]
            partes.append(oos[(con_g[j], g)])
        salida.append(pd.concat(partes))
    return salida

CPCV to choose a parameter

The natural application for a strategy is not a classifier, but the procedure of lesson 2.5.2: choose the lookback L of the futures momentum with the Sharpe of the training data and apply it to the test groups. The label of each daily observation is the next day's return ($t_0 = t_1 = t + 1$), so purging removes almost nothing; the embargo, of 63 sessions, covers part of the feature overlap (the longest lookback is 378).

fut = pd.read_parquet("futuros2_continuos.parquet")
precio = fut.pivot(index="fecha", columns="mercado",
                   values="cierre_ajustado")
rf = precio.pct_change(fill_method=None)
vol = rf.ewm(span=60, min_periods=60).std() * np.sqrt(252)

def sharpe(x):
    return x.mean() / x.std() * np.sqrt(252)

variantes = {}
for n in (10, 21, 42, 63, 126, 189, 252, 378):
    senal = np.sign(precio.pct_change(n, fill_method=None))
    peso = (senal * 0.40 / vol / precio.shape[1]).shift(1)
    variantes[n] = (peso * rf).sum(axis=1, min_count=1)
variantes = pd.DataFrame(variantes).loc["2002":]

def cpcv_eleccion(var, N=10, k=2, embargo=63):
    n = len(var)
    p = np.arange(n)
    oos, elegidos = {}, []
    for d, prueba, grupos, entreno in cpcv(n, N, k, p + 1, p + 1,
                                           embargo):
        mejor = sharpe(var.iloc[entreno]).idxmax()
        elegidos.append(mejor)
        for g in prueba:
            oos[(d, g)] = var.iloc[grupos[g]][mejor]
    sh = pd.Series([sharpe(c) for c in caminos(N, k, oos)])
    return sh, pd.Series(elegidos).value_counts()

sh_fut, elegidos = cpcv_eleccion(variantes)
print("Sharpe of the 9 paths:", sh_fut.round(2).tolist())
print(f"median {sh_fut.median():.2f}, min {sh_fut.min():.2f}, "
      f"max {sh_fut.max():.2f}")
print("Choices in the 45 splits:", elegidos.to_dict())

Sharpe of the 9 paths: [1.06, 1.06, 0.97, 0.84, 1.06, 1.06, 1.06, 1.06, 1.06]
median 1.06, min 0.84, max 1.06
Choices in the 45 splits: {252: 43, 189: 2}

All nine paths are positive and the choice is very stable: 252, the planted value, wins in 43 of the 45 splits, and 189, its neighbor, in the other two. That is why seven paths give exactly 1.06, the Sharpe of the 252 variant over the whole sample (lesson 2.5.2): if every split chooses the same, each path is the same series, cut into groups and glued back together. Only the two paths that go through the splits that chose 189 deviate, down to 0.97 and 0.84. The dispersion between paths is the information a walk-forward does not give, but it has to be read properly: in a CPCV it measures how much the choice changes from one split to another, not which years fall together in the test. Here it is small because the effect is real and stable; the cost of choosing comes down to those two splits.

The contrast with noise, again on ruido2_estrategias.parquet: choose the best of 1,000 in the training set and test it on the test groups. With 1,003 sessions, each group has about 100.

est = pd.read_parquet("ruido2_estrategias.parquet").pivot(
    index="fecha", columns="estrategia", values="rendimiento")
sh_ruido, eleg_ruido = cpcv_eleccion(est, embargo=5)
print("Sharpe of the 9 paths:", sh_ruido.round(2).tolist())
print(f"median {sh_ruido.median():.2f}; "
      f"{len(eleg_ruido)} distinct strategies in 45 splits")

Sharpe of the 9 paths: [-1.3, 0.14, -0.75, -0.35, 0.56, -0.46, -0.2, 0.31, -0.27]
median -0.27; 11 distinct strategies in 45 splits

On the noise, the median of the paths is negative and the range goes from heavy losses to moderate gains: a distribution centered close to zero, as befits a procedure with no edge. The choice changes from one split to another (11 different strategies), and that is why the paths spread out so much. A single path, the one a walk-forward would have given, could have landed anywhere in the range; had it landed on the best one, the procedure would have looked promising.

What CPCV is and what it is not

CPCV does not simulate history. In most splits it also trains on data after the test: for the 2005 group the parameter is chosen with data from 2010 or 2020. That is legitimate for estimating the performance of a procedure if the process generating the data is stable, and it is optimistic if it is not (an effect that is discovered in 2015 and arbitraged away afterwards helps to “predict” 2005). That is why the two tools complement each other:

CPCV paths are not independent (they share training and test data), so their dispersion is not a confidence interval in the strict sense. It is useful for comparing procedures and for feeding measures such as the probability of backtest overfitting, which lesson 2.5.4 computes with a related technique (CSCV).

Exercises

1. The shuffled k-fold leak needs overlapping labels. Check it: repeat the shuffled k-fold with labels of 1, 5, 20 and 60 sessions, with the same features. What happens with a single-session label?

def auc_barajado(h):
    fut_h = r[::-1].rolling(h).sum()[::-1].shift(-1)
    d = rasgos.assign(futuro=fut_h).dropna()
    Xh = d[rasgos.columns].to_numpy()
    yh = (d["futuro"] > 0).astype(int).to_numpy()
    res = []
    for a, b in KFold(10, shuffle=True, random_state=0).split(Xh):
        m = modelo().fit(Xh[a], yh[a])
        res.append(roc_auc_score(yh[b], m.predict_proba(Xh[b])[:, 1]))
    return np.mean(res)

pd.Series({h: auc_barajado(h) for h in (1, 5, 20, 60)},
          name="shuffled AUC").rename_axis("H").round(3)

H
1     0.507
5     0.790
20    0.872
60    0.928
Name: shuffled AUC, dtype: float64

2. Compute how many paths and how many splits CPCV gives with N = 12 and k = 3, and what fraction of the data is left for training before purging. Does it seem a good configuration for 1,003 sessions?

N, k = 12, 3
print(f"splits {comb(N, k)}, paths {comb(N - 1, k - 1)}, "
      f"train {(N - k) / N:.0%}, sessions per group "
      f"{1003 // N}")

splits 220, paths 55, train 75%, sessions per group 83

3. In the futures CPCV, the 63-session embargo does not cover the longest lookback (378). Explain what leak remains and how you would check whether it matters.