Backtesting ML signals with costs and Deflated Sharpe.

From IC to money

The walk-forward of lesson 5.4.3 ended with a ridge with an IC of 0.0185 (t = 4.9) over 574 weeks: a real, stable and statistically clear signal. The question of this lesson is the one that matters to whoever has to trade it: how much does a portfolio built on that signal earn, after paying to change it every week? The path has four steps, and a fifth that tends to be forgotten:

flowchart TB
  S[Score] --> C[Portfolio rule]
  C --> W[Weights per week]
  W --> T[Turnover and costs]
  T --> N[Net return]
  N --> R[Log]
  R --> D[Deflated Sharpe]

The fifth is the log: every model and portfolio variant you try is a trial, and the Sharpe of the best one only means something after discounting how many trials there were (lesson 2.5.4).

The scores: the walk-forward of 5.4.3

Each lesson runs separately, so this one redoes in a single cell what it needs from 5.4.3: the panel as arrays, the ridge from sufficient statistics and the walk-forward with the time rule. It adds two reference models: the oracle and the ridge with revision_analistas, the leaky feature.

import numpy as np
import pandas as pd
from scipy.stats import norm, rankdata

feat = pd.read_parquet("panel5_features.parquet")
feat = feat.sort_values(["fecha", "ticker"])
etiq = pd.read_parquet("panel5_etiquetas.parquet")
etiq = etiq.sort_values(["fecha", "ticker"])
verdad = pd.read_parquet("panel5_verdad.parquet")
verdad = verdad.sort_values(["fecha", "ticker"])
TODAS = [c for c in feat.columns if c not in
         ("fecha", "ticker", "sector")]
fechas = pd.DatetimeIndex(feat["fecha"].unique())
nW, nA = len(fechas), feat["ticker"].nunique()

F = feat[TODAS].to_numpy(dtype=float).reshape(nW, nA, -1)
with np.errstate(invalid="ignore", divide="ignore"):
    n = np.isfinite(F).sum(axis=1, keepdims=True)
    G = np.nan_to_num(norm.ppf(
        (rankdata(F, axis=1, nan_policy="omit") - 0.5) / n))
R = etiq["ret_sig_semana"].to_numpy().reshape(nW, nA)
fin = etiq["fin_sig_semana"].to_numpy().reshape(nW, nA)[:, 0]
con_feat = np.isfinite(F[:, :, 0]).all(axis=1)
valida = con_feat & np.isfinite(R).all(axis=1)
Yc = np.zeros_like(R)
Yc[valida] = R[valida] - R[valida].mean(axis=1, keepdims=True)
anyo = fechas.year.to_numpy()
semanas = np.flatnonzero(valida & (anyo >= 2015))  # backtest

def estadisticos(cols):
    X1 = np.concatenate([np.ones((nW, nA, 1)), G[:, :, cols]], 2)
    pv = valida.astype(float)
    return (np.einsum("wni,wnj->wij", X1, X1) * pv[:, None, None],
            np.einsum("wni,wn->wi", X1, Yc) * pv[:, None])

def walk_forward(cols, ventana=None, cada=None, anual=False,
                 semivida=None, alfa=10.0):
    S_xx, S_xy = estadisticos(cols)
    k = len(cols)
    pred, beta = np.full((nW, nA), np.nan), None
    for i, w in enumerate(semanas):
        if (beta is None or (cada and i % cada == 0) or
                (anual and anyo[w] != anyo[semanas[i - 1]])):
            u = int(np.searchsorted(fin, fechas[w])) - 1
            p = valida.astype(float)
            p[u + 1:] = 0.0
            if ventana:
                p[:max(0, u + 1 - ventana)] = 0.0
            if semivida:
                p *= 0.5 ** ((u - np.arange(nW)) / semivida)
            A, b = np.tensordot(p, S_xx, 1), np.tensordot(p, S_xy, 1)
            mx, my = A[0, 1:] / A[0, 0], b[0] / A[0, 0]
            C = A[1:, 1:] - A[0, 0] * np.outer(mx, mx)
            c = b[1:] - A[0, 0] * mx * my
            beta = np.linalg.solve(C + alfa * np.eye(k), c)
        pred[w] = G[w][:, cols] @ beta
    return pred

sin_fuga = [i for i, c in enumerate(TODAS)
            if c != "revision_analistas"]
modelos = {
    "exp. annual": walk_forward(sin_fuga, anual=True),
    "exp. weekly": walk_forward(sin_fuga, cada=1),
    "rolling 3 years": walk_forward(sin_fuga, ventana=156, cada=1),
    "half-life 2 years": walk_forward(sin_fuga, semivida=104, cada=1),
}
fuga = walk_forward(list(range(len(TODAS))), anual=True)
oraculo = verdad["alfa_verdadero"].to_numpy().reshape(nW, nA)
print(f"{len(semanas)} backtest weeks: "
      f"{fechas[semanas[0]].date()} to {fechas[semanas[-1]].date()}")

574 backtest weeks: 2015-01-02 to 2025-12-26

From scores to weights, with costs

The portfolio rule is the README one: each week, long the top-score decile and short the bottom one, equally weighted, one dollar per leg (two of gross notional). Between one week and the next, weights drift with returns; turnover is the sum of the absolute weight changes needed to get back to the target portfolio, and the cost is 10 bp per unit of turnover:

with $\tilde w_{i,t-1} = w_{i,t-1}(1 + r_{i,t})$ the previous week's weight after moving with the market and $c$ = 0.001. The function supports three variants that reduce turnover and are used further down:

• cada: rebalance only every n weeks (in between, the weights drift).
• suavizado: trade on an exponential average of the score's rank, $s_t = \lambda\,\text{rank}_t + (1 - \lambda)\,s_{t-1}$, instead of on that week's score. It is the idea of Gârleanu and Pedersen (2013): move only part of the way toward the target portfolio.
• banda: enter in the extreme 10% and do not exit while still in the 30% (a band between buying and holding, like those studied by Novy-Marx and Velikov, 2016).

def cartera(score, q=0.1, cada=1, suavizado=None, banda=None,
            coste=0.001):
    """Gross and net weekly returns and mean turnover."""
    w_ant = np.zeros(nA)
    s_prev = None
    bruto, neto, rot = [], [], []
    for i, w in enumerate(semanas):
        rk = (rankdata(score[w]) - 0.5) / nA
        if suavizado:
            s_prev = rk if s_prev is None else (
                suavizado * rk + (1 - suavizado) * s_prev)
            rk = (rankdata(s_prev) - 0.5) / nA
        if i % cada == 0:
            largo, corto = rk >= 1 - q, rk <= q
            if banda:
                largo |= (w_ant > 0) & (rk >= 1 - banda)
                corto |= (w_ant < 0) & (rk <= banda)
            peso = np.zeros(nA)
            peso[largo] = 1 / largo.sum()
            peso[corto] = -1 / corto.sum()
        else:
            peso = w_ant
        r = float(peso @ R[w])
        t = np.abs(peso - w_ant).sum()
        bruto.append(r)
        neto.append(r - coste * t)
        rot.append(t)
        w_ant = peso * (1 + R[w])
    return np.array(bruto), np.array(neto), float(np.mean(rot[1:]))

def sharpe(x):
    return x.mean() / x.std(ddof=1) * np.sqrt(52)

filas = {}
for nombre, sc in (("oracle", oraculo),
                   ("ridge (exp. annual)", modelos["exp. annual"]),
                   ("ridge + leak", fuga)):
    b, nt, rot = cartera(sc)
    filas[nombre] = {"gross ret.": f"{b.mean() * 52:.1%}",
                     "vol.": f"{b.std(ddof=1) * np.sqrt(52):.1%}",
                     "gross Sharpe": round(sharpe(b), 2),
                     "turnover": f"{rot:.0%}",
                     "net ret.": f"{nt.mean() * 52:.1%}",
                     "net Sharpe": round(sharpe(nt), 2)}
pd.DataFrame(filas).T

                    gross ret.  vol. gross Sharpe turnover net ret. net Sharpe
oracle                   17.1%  9.6%         1.78     259%     3.6%       0.38
ridge (exp. annual)      12.1%  9.7%         1.25     265%    -1.7%      -0.18
ridge + leak            136.5%  6.9%        19.89     287%   121.6%      17.74

The three rows are the README ones. The oracle, which knows the true alpha, has a gross Sharpe of 1.78 and a net one of 0.38: even the best possible information barely survives a weekly turnover of 259%. The honest ridge goes from 1.25 gross to −0.18 net. The leaky one, from 19.9 to 17.7: an impossible number, which in a real backtest should set off the alarm before any test.

The figure at the top shows the difference over time for the ridge: the gross curve rises steadily; the net one, with the same shape, ends below where it started.

The cost arithmetic

The arithmetic is simple and worth doing before any backtest. With a weekly turnover $\tau$ and a cost $c$ per unit of turnover, the annual cost is $52\,\tau\,c$. The signal can bear the cost as long as its annual gross return exceeds it:

b, nt, rot = cartera(modelos["exp. annual"])
coste_anual = 52 * rot * 0.001
equilibrio = b.mean() * 52 / (52 * rot)
print(f"annual gross return {b.mean() * 52:.1%}; "
      f"annual cost {coste_anual:.1%}")
print(f"break-even cost: {equilibrio * 1e4:.1f} bp per side")

annual gross return 12.1%; annual cost 13.8%
break-even cost: 8.8 bp per side

The ridge earns 12.1% gross a year and pays 13.8% in costs. It could bear about 8.8 bp per side: with the README's 10, it loses. There are two levers: more gross return per unit of turnover (a better signal) or less turnover for the same return (a better portfolio rule). The second is the one you can touch without modeling again.

An experiment log

Before trying variants, the log. Each experiment stores its full configuration, its metrics and its series of net returns: without the series you cannot compute the dispersion of the Sharpes or the PBO at the end (lesson 2.2.4). The log also stores the variants that do not work, which are the ones that tell you most about luck.

The grid is small and reasonable, the one anyone would try: the four walk-forward models of 5.4.3 (expanding annual and weekly, three-year rolling, weights with a two-year half-life) times five portfolio rules (deciles, quintiles, deciles every 4 weeks, smoothed ranks with λ = 0.25 and a 30% band). Twenty experiments.

reglas = {"deciles": {}, "quintiles": {"q": 0.2},
          "every 4 weeks": {"cada": 4},
          "smoothed 0.25": {"suavizado": 0.25},
          "band 30%": {"banda": 0.3}}
registro, series = [], {}
for m, sc in modelos.items():
    for r, kw in reglas.items():
        b, nt, rot = cartera(sc, **kw)
        clave = f"{m} | {r}"
        registro.append({"experiment": clave, "model": m,
                         "rule": r, "params": kw,
                         "gross Sharpe": sharpe(b),
                         "net Sharpe": sharpe(nt), "turnover": rot})
        series[clave] = nt
registro = pd.DataFrame(registro).set_index("experiment")
netos = pd.DataFrame(series, index=fechas[semanas])
print(f"logged experiments: {len(registro)}")
(registro[["gross Sharpe", "net Sharpe", "turnover"]]
 .sort_values("net Sharpe", ascending=False).head(8).round(2))

logged experiments: 20
                                   gross Sharpe  net Sharpe  turnover
experiment
exp. annual | smoothed 0.25                0.88        0.37      0.86
exp. weekly | smoothed 0.25                0.79        0.25      0.88
exp. annual | every 4 weeks                0.59        0.18      0.73
exp. annual | band 30%                     1.09        0.12      1.44
exp. weekly | every 4 weeks                0.46        0.05      0.73
half-life 2 years | smoothed 0.25          0.56        0.05      0.86
exp. weekly | band 30%                     1.00        0.03      1.49
exp. annual | deciles                      1.25       -0.18      2.65

The figure places the twenty experiments on the turnover-net Sharpe plane. The rule is visible: the high-turnover variants (weekly deciles and quintiles, on the right) lose money with any model; the low-turnover ones (rebalancing every 4 weeks, smoothing, band) stay close to zero or somewhat above it with the expanding models, and below it with the rolling window, which lesson 5.4.3 already ruled out because of its IC. The best is the annual expanding ridge with smoothed ranks: a net Sharpe of 0.37, with a turnover of 86%. Almost the same as the oracle with weekly deciles.

The Deflated Sharpe of the best one

Is 0.37 a result? Three questions, from least to most demanding. The functions are those of lesson 2.5.4, with the Sharpes at the frequency of the data (weekly):

def psr(sr, umbral, t, asim, curt):
    den = np.sqrt(1 - asim * sr + (curt - 1) / 4 * sr ** 2)
    return float(norm.cdf((sr - umbral) * np.sqrt(t - 1) / den))

def sr_max_esperado(n, var_sr):
    g = np.euler_gamma
    return np.sqrt(var_sr) * ((1 - g) * norm.ppf(1 - 1 / n)
                              + g * norm.ppf(1 - 1 / (n * np.e)))

sr_sem = netos.mean() / netos.std()       # weekly Sharpes
mejor = sr_sem.idxmax()
x = netos[mejor]
T, asim, curt = len(x), x.skew(), x.kurt() + 3
sr = sr_sem[mejor]
et_anual = np.sqrt((1 + sr ** 2 / 2) / T) * np.sqrt(52)
print(f"best: {mejor}")
print(f"net Sharpe {sr * np.sqrt(52):.2f} ± {et_anual:.2f} "
      f"(standard error, {T} weeks)")
print(f"PSR(0) = {psr(sr, 0, T, asim, curt):.3f}")

autoval = np.linalg.eigvalsh(netos.corr().to_numpy())
n_ef = autoval.sum() ** 2 / (autoval ** 2).sum()
for etiqueta, n_int in (("effective trials", n_ef),
                        ("20 logged", 20),
                        ("with the 13 from 5.4.3 x 5 rules", 65)):
    sr0 = sr_max_esperado(n_int, sr_sem.var())
    print(f"N = {n_int:5.1f} ({etiqueta}): SR_0 = "
          f"{sr0 * np.sqrt(52):.2f} annual, "
          f"DSR = {psr(sr, sr0, T, asim, curt):.3f}")

best: exp. annual | smoothed 0.25
net Sharpe 0.37 ± 0.30 (standard error, 574 weeks)
PSR(0) = 0.890
N =   2.0 (effective trials): SR_0 = 0.20 annual, DSR = 0.713
N =  20.0 (20 logged): SR_0 = 0.71 annual, DSR = 0.128
N =  65.0 (with the 13 from 5.4.3 x 5 rules): SR_0 = 0.89 annual, DSR = 0.042

First question: can it be told apart from zero, as if it were the only trial? The standard error of an annual Sharpe with eleven years of weekly data is ±0.30, almost as large as the Sharpe itself, and the PSR(0) is 0.89: even without correction it does not reach 95%.

Second: does it beat the best one expected among the trials? The 20 experiments are highly correlated (they share models and rules), and the participation ratio of the eigenvalues of their correlation matrix gives only 2 effective trials: with such a favorable N, the threshold is a net Sharpe of 0.20 and the DSR, 0.71. With the 20 declared, the threshold rises to 0.71 and the DSR drops to 0.13. And the honest count includes the trials of the previous lesson: the walk-forward of 5.4.3 tested thirteen schemes before choosing the four used here, and each one could have been combined with five rules; with 65, the DSR drops to 0.04. No way of counting reaches 0.95.

Third: what do you present? What the numbers allow: the honest ridge gives a net Sharpe compatible with zero, between −0.4 and 0.4 depending on the portfolio rule with the expanding models, with a standard error of ±0.3. It is the README's answer (“net Sharpe ≈ 0”) and it is the correct conclusion, even though the gross IC is clearly positive.

The figure sums up the second question: the DSR of the best experiment depending on how many trials are declared. No reasonable way of counting brings it close to the 0.95 line.

The DSR does not see leaks

The ridge with revision_analistas needs no correction to look good, and no correction uncovers it:

_, nt_fuga, rot_fuga = cartera(fuga)
sr_f = nt_fuga.mean() / nt_fuga.std(ddof=1)
var_con = pd.concat([sr_sem, pd.Series({"leak": sr_f})]).var()
sr0 = sr_max_esperado(21, var_con)
xf = pd.Series(nt_fuga)
dsr_f = psr(sr_f, sr0, len(xf), xf.skew(), xf.kurt() + 3)
print(f"net Sharpe {sr_f * np.sqrt(52):.1f}; SR_0 with N = 21: "
      f"{sr0 * np.sqrt(52):.2f}; DSR = {dsr_f:.3f}")

net Sharpe 17.7; SR_0 with N = 21: 7.58; DSR = 1.000

A net Sharpe of 17.7 leaves any threshold behind, even the 7.6 that comes out when the leak itself is included in the dispersion of the trials: the DSR is 1.000. The DSR discounts selection (having chosen the best of N trials with no edge), not information from the future: with a leak, the edge is real inside the backtest, because the backtest sees it. The only things that give it away are the order of magnitude (a weekly net Sharpe of 17 in stocks does not exist) and, above all, the audit of when each piece of data is known (lesson 5.2.2). Validation (5.4.1-5.4.3) and the correction for selection (this lesson) are applied after that audit, not instead of it.

What to present

The last row is not a failure: it is information. It says that the next effort is not in another model (the ridge already captures a sizable part of what the oracle captures) but in execution: lower costs (lesson 4.3.1), optimized portfolios with a turnover penalty (lesson 4.1.5) or a signal combined with other, less correlated ones (lesson 5.10.3).

Exercises

1. Apply the best portfolio rule (smoothed ranks with λ = 0.25) to the oracle. What is the best net Sharpe this panel allows with that rule, and how far is the ridge from it?

for nombre, sc in (("oracle", oraculo),
                   ("ridge exp. annual", modelos["exp. annual"])):
    b, nt, rot = cartera(sc, suavizado=0.25)
    print(f"{nombre:<18} gross {sharpe(b):.2f}  net {sharpe(nt):.2f}"
          f"  turnover {rot:.0%}")

oracle             gross 1.14  net 0.64  turnover 91%
ridge exp. annual  gross 0.88  net 0.37  turnover 86%

2. Compute the break-even cost per side of the best variant (annual expanding ridge, smoothing 0.25). At what real costs would it be tradable?

b, nt, rot = cartera(modelos["exp. annual"], suavizado=0.25)
print(f"turnover {rot:.0%} weekly; annual gross {b.mean() * 52:.1%};"
      f" break-even {b.mean() / rot * 1e4:.1f} bp per side")

turnover 86% weekly; annual gross 7.7%; break-even 17.2 bp per side

3. Write down what fields the log of an experiment would have in a real project, beyond those in the cell.