Backtest de señales de ML con costes y Deflated Sharpe.

Del IC al dinero

El walk-forward de la lección 5.4.3 terminó con un ridge de IC 0,0185 (t = 4,9) en 574 semanas: una señal real, estable y estadísticamente clara. La pregunta de esta lección es la que importa a quien tiene que operarla: ¿cuánto gana una cartera construida con esa señal, después de pagar por cambiarla cada semana? El camino tiene cuatro pasos, y un quinto que se suele olvidar:

flowchart TB
  S[Puntuación] --> C[Regla de cartera]
  C --> W[Pesos por semana]
  W --> T[Rotación y costes]
  T --> N[Rendimiento neto]
  N --> R[Registro]
  R --> D[Deflated Sharpe]

El quinto es el registro: cada variante de modelo y de cartera que se prueba es un intento, y el Sharpe de la mejor solo significa algo descontando cuántos intentos hubo (lección 2.5.4).

Las puntuaciones: el walk-forward de 5.4.3

Cada lección se ejecuta por separado, así que esta rehace en una celda lo necesario de la 5.4.3: el panel en arrays, el ridge desde estadísticos suficientes y el walk-forward con la regla del tiempo. Añade dos modelos de referencia: el oráculo y el ridge con revision_analistas, la feature con fuga.

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. anual": walk_forward(sin_fuga, anual=True),
    "exp. semanal": walk_forward(sin_fuga, cada=1),
    "móvil 3 años": walk_forward(sin_fuga, ventana=156, cada=1),
    "semivida 2 años": 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)} semanas de backtest: "
      f"{fechas[semanas[0]].date()} a {fechas[semanas[-1]].date()}")

574 semanas de backtest: 2015-01-02 a 2025-12-26

De puntuaciones a pesos, con costes

La regla de cartera es la del README: cada semana, larga en el decil de puntuación más alta y corta en el más baja, equiponderada, un dólar por pata (dos de nominal bruto). Entre una semana y la siguiente, los pesos derivan con los rendimientos; la rotación es la suma de los cambios absolutos de peso necesarios para volver a la cartera objetivo, y el coste, 10 pb por cada unidad de rotación:

con $\tilde w_{i,t-1} = w_{i,t-1}(1 + r_{i,t})$ el peso de la semana anterior después de moverse con el mercado y $c$ = 0,001. La función admite tres variantes que reducen la rotación y se usan más abajo:

• cada: rebalancear solo cada n semanas (entre medias, los pesos derivan).
• suavizado: operar con una media exponencial del rango de la puntuación, $s_t = \lambda\,\text{rango}_t + (1 - \lambda)\,s_{t-1}$, en lugar de con la puntuación de esa semana. Es la idea de Gârleanu y Pedersen (2013): moverse solo una parte del camino hacia la cartera objetivo.
• banda: entrar en el 10 % de los extremos y no salir mientras se siga en el 30 % (una banda entre comprar y mantener, como las que estudian Novy-Marx y Velikov, 2016).

def cartera(score, q=0.1, cada=1, suavizado=None, banda=None,
            coste=0.001):
    """Rendimientos semanales brutos y netos y rotación media."""
    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 (("oráculo", oraculo),
                   ("ridge (exp. anual)", modelos["exp. anual"]),
                   ("ridge + fuga", fuga)):
    b, nt, rot = cartera(sc)
    filas[nombre] = {"rend. bruto": f"{b.mean() * 52:.1%}",
                     "vol.": f"{b.std(ddof=1) * np.sqrt(52):.1%}",
                     "Sharpe bruto": round(sharpe(b), 2),
                     "rotación": f"{rot:.0%}",
                     "rend. neto": f"{nt.mean() * 52:.1%}",
                     "Sharpe neto": round(sharpe(nt), 2)}
pd.DataFrame(filas).T

                   rend. bruto  vol. Sharpe bruto rotación rend. neto Sharpe neto
oráculo                  17.1%  9.6%         1.78     259%       3.6%        0.38
ridge (exp. anual)       12.1%  9.7%         1.25     265%      -1.7%       -0.18
ridge + fuga            136.5%  6.9%        19.89     287%     121.6%       17.74

Las tres filas son las del README. El oráculo, que conoce el alfa verdadero, tiene un Sharpe bruto de 1,78 y neto de 0,38: incluso la mejor información posible apenas sobrevive a una rotación del 259 % semanal. El ridge honesto pasa de 1,25 bruto a −0,18 neto. El de la fuga, de 19,9 a 17,7: un número imposible, que en un backtest real tendría que disparar la alarma antes que cualquier contraste.

La figura de la cabecera muestra la diferencia en el tiempo para el ridge: la curva bruta sube con regularidad; la neta, con la misma forma, termina por debajo de donde empezó.

La cuenta del coste

La aritmética es sencilla y conviene hacerla antes de ningún backtest. Con una rotación semanal $\tau$ y un coste $c$ por unidad de rotación, el coste anual es $52\,\tau\,c$. La señal soporta el coste mientras su rendimiento bruto anual lo supere:

b, nt, rot = cartera(modelos["exp. anual"])
coste_anual = 52 * rot * 0.001
equilibrio = b.mean() * 52 / (52 * rot)
print(f"rendimiento bruto anual {b.mean() * 52:.1%}; "
      f"coste anual {coste_anual:.1%}")
print(f"coste de equilibrio: {equilibrio * 1e4:.1f} pb por lado")

rendimiento bruto anual 12.1%; coste anual 13.8%
coste de equilibrio: 8.8 pb por lado

El ridge gana un 12,1 % bruto al año y paga un 13,8 % en costes. Soportaría unos 8,8 pb por lado: con los 10 del README, pierde. Hay dos palancas: más rendimiento bruto por unidad de rotación (mejor señal) o menos rotación para el mismo rendimiento (mejor regla de cartera). La segunda es la que se puede tocar sin volver a modelizar.

Un registro de experimentos

Antes de probar variantes, el registro. Cada experimento guarda su configuración completa, sus métricas y su serie de rendimientos netos: sin la serie no se puede calcular la dispersión de los Sharpes ni la PBO al final (lección 2.2.4). El registro guarda también las variantes que no funcionan, que son las que más información dan sobre la suerte.

La malla es pequeña y razonable, la que cualquiera probaría: los cuatro modelos del walk-forward de 5.4.3 (expansivo anual y semanal, móvil de tres años, pesos con semivida de dos años) por cinco reglas de cartera (deciles, quintiles, deciles cada 4 semanas, rangos suavizados con λ = 0,25 y banda del 30 %). Veinte experimentos.

reglas = {"deciles": {}, "quintiles": {"q": 0.2},
          "cada 4 semanas": {"cada": 4},
          "suavizado 0,25": {"suavizado": 0.25},
          "banda 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({"experimento": clave, "modelo": m,
                         "regla": r, "parámetros": kw,
                         "Sharpe bruto": sharpe(b),
                         "Sharpe neto": sharpe(nt), "rotación": rot})
        series[clave] = nt
registro = pd.DataFrame(registro).set_index("experimento")
netos = pd.DataFrame(series, index=fechas[semanas])
print(f"experimentos registrados: {len(registro)}")
(registro[["Sharpe bruto", "Sharpe neto", "rotación"]]
 .sort_values("Sharpe neto", ascending=False).head(8).round(2))

experimentos registrados: 20
                                  Sharpe bruto  Sharpe neto  rotación
experimento
exp. anual | suavizado 0,25               0.88         0.37      0.86
exp. semanal | suavizado 0,25             0.79         0.25      0.88
exp. anual | cada 4 semanas               0.59         0.18      0.73
exp. anual | banda 30 %                   1.09         0.12      1.44
exp. semanal | cada 4 semanas             0.46         0.05      0.73
semivida 2 años | suavizado 0,25          0.56         0.05      0.86
exp. semanal | banda 30 %                 1.00         0.03      1.49
exp. anual | deciles                      1.25        -0.18      2.65

La figura pone los veinte experimentos en el plano rotación-Sharpe neto. Se ve la regla: las variantes de mucha rotación (deciles y quintiles semanales, a la derecha) pierden dinero con cualquier modelo; las de poca rotación (rebalanceo cada 4 semanas, suavizado, banda) quedan cerca de cero o algo por encima con los modelos expansivos, y por debajo con la ventana móvil, que la lección 5.4.3 ya descartó por su IC. La mejor es el ridge expansivo anual con rangos suavizados: Sharpe neto de 0,37, con una rotación del 86 %. Casi lo mismo que el oráculo con deciles semanales.

El Deflated Sharpe del mejor

¿Es 0,37 un resultado? Tres preguntas, de menos a más exigentes. Las funciones son las de la lección 2.5.4, con los Sharpes en la frecuencia de los datos (semanal):

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()       # Sharpes semanales
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"mejor: {mejor}")
print(f"Sharpe neto {sr * np.sqrt(52):.2f} ± {et_anual:.2f} "
      f"(error típico, {T} semanas)")
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 (("intentos efectivos", n_ef),
                        ("20 registrados", 20),
                        ("con los 13 de 5.4.3 x 5 reglas", 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} anual, "
          f"DSR = {psr(sr, sr0, T, asim, curt):.3f}")

mejor: exp. anual | suavizado 0,25
Sharpe neto 0.37 ± 0.30 (error típico, 574 semanas)
PSR(0) = 0.890
N =   2.0 (intentos efectivos): SR_0 = 0.20 anual, DSR = 0.713
N =  20.0 (20 registrados): SR_0 = 0.71 anual, DSR = 0.128
N =  65.0 (con los 13 de 5.4.3 x 5 reglas): SR_0 = 0.89 anual, DSR = 0.042

Primera pregunta: ¿se distingue de cero, como si fuera el único intento? El error típico de un Sharpe anual con once años de datos semanales es de ±0,30, casi tanto como el propio Sharpe, y el PSR(0) es 0,89: ni siquiera sin corregir llega al 95 %.

Segunda: ¿supera al mejor esperado entre los intentos? Los 20 experimentos están muy correlacionados (comparten modelos y reglas), y la razón de participación de los autovalores de su matriz de correlación da solo 2 intentos efectivos: con ese N tan favorable, el umbral es un Sharpe neto de 0,20 y el DSR, 0,71. Con los 20 declarados, el umbral sube a 0,71 y el DSR baja a 0,13. Y la cuenta honesta incluye los intentos de la lección anterior: el walk-forward de 5.4.3 probó trece esquemas antes de elegir los cuatro de aquí, y cada uno podría haberse combinado con cinco reglas; con 65, el DSR baja a 0,04. Ninguna forma de contar llega al 0,95.

Tercera: ¿qué se presenta? Lo que los números permiten: el ridge honesto da un Sharpe neto compatible con cero, entre −0,4 y 0,4 según la regla de cartera con los modelos expansivos, con un error típico de ±0,3. Es la respuesta del README («Sharpe neto ≈ 0») y es la conclusión correcta, aunque el IC bruto sea claramente positivo.

La figura resume la segunda pregunta: el DSR del mejor experimento según cuántos intentos se declaren. Ninguna forma razonable de contar lo acerca a la línea de 0,95.

El DSR no ve las fugas

El ridge con revision_analistas no necesita ninguna corrección para parecer bueno, y ninguna corrección lo descubre:

_, nt_fuga, rot_fuga = cartera(fuga)
sr_f = nt_fuga.mean() / nt_fuga.std(ddof=1)
var_con = pd.concat([sr_sem, pd.Series({"fuga": 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"Sharpe neto {sr_f * np.sqrt(52):.1f}; SR_0 con N = 21: "
      f"{sr0 * np.sqrt(52):.2f}; DSR = {dsr_f:.3f}")

Sharpe neto 17.7; SR_0 con N = 21: 7.58; DSR = 1.000

Un Sharpe neto de 17,7 deja atrás cualquier umbral, incluso el de 7,6 que sale al meter la propia fuga en la dispersión de los intentos: el DSR es 1,000. El DSR descuenta la selección (haber elegido el mejor de N intentos sin ventaja), no la información del futuro: con una fuga, la ventaja es real dentro del backtest, porque el backtest la ve. Lo único que la delata es el orden de magnitud (un Sharpe neto semanal de 17 en acciones no existe) y, sobre todo, la auditoría de cuándo se conoce cada dato (lección 5.2.2). La validación (5.4.1-5.4.3) y la corrección por selección (esta lección) se aplican después de esa auditoría, no en su lugar.

Qué presentar

La última fila no es un fracaso: es información. Dice que el siguiente esfuerzo no está en otro modelo (el ridge ya recoge una parte apreciable de lo que recoge el oráculo) sino en la ejecución: costes más bajos (lección 4.3.1), carteras optimizadas con penalización de la rotación (lección 4.1.5) o una señal combinada con otras menos correlacionadas (lección 5.10.3).

Ejercicios

1. Aplica la mejor regla de cartera (rangos suavizados con λ = 0,25) al oráculo. ¿Cuál es el mejor Sharpe neto que permite este panel con esa regla, y cuánto le falta al ridge?

for nombre, sc in (("oráculo", oraculo),
                   ("ridge exp. anual", modelos["exp. anual"])):
    b, nt, rot = cartera(sc, suavizado=0.25)
    print(f"{nombre:<18} bruto {sharpe(b):.2f}  neto {sharpe(nt):.2f}"
          f"  rotación {rot:.0%}")

oráculo            bruto 1.14  neto 0.64  rotación 91%
ridge exp. anual   bruto 0.88  neto 0.37  rotación 86%

2. Calcula el coste por lado de equilibrio de la mejor variante (ridge expansivo anual, suavizado 0,25). ¿Con qué costes reales sería operable?

b, nt, rot = cartera(modelos["exp. anual"], suavizado=0.25)
print(f"rotación {rot:.0%} semanal; bruto anual {b.mean() * 52:.1%};"
      f" equilibrio {b.mean() / rot * 1e4:.1f} pb por lado")

rotación 86% semanal; bruto anual 7.7%; equilibrio 17.2 pb por lado

3. Escribe qué campos tendría el registro de un experimento en un proyecto real, más allá de los de la celda.