“Honest backtest” checklist before moving to paper trading.

import itertools
import numpy as np
import pandas as pd
from dataclasses import dataclass, field
from scipy import stats

aj = pd.read_parquet("universo2_ajustados.parquet")
px = aj.pivot(index="fecha", columns="ticker", values="cierre_ajustado")
fin_mes = px.index.to_series().groupby(px.index.to_period("M")).max()
mensual = px.loc[fin_mes.to_numpy()]
R = mensual.pct_change(fill_method=None)

# the month of a delisting: up to the last price and then the payout
eventos = pd.read_csv("universo2_eventos.csv", parse_dates=["fecha"])
deslist = eventos.dropna(subset=["rentabilidad_deslistado"])
for _, ev in deslist.iterrows():
    tk = ev["ticker"]
    ultima = px[tk].last_valid_index()
    mes = R.index[R.index >= ultima][0]
    antes = R.index[R.index < ultima][-1]
    R.loc[mes, tk] = (px.loc[ultima, tk] / mensual.loc[antes, tk]
                      * (1 + ev["rentabilidad_deslistado"]) - 1)
print(R.shape, f"{len(deslist)} delisted stocks with their final return")
R.loc["2008-10", ["QXQT"]]

car = pd.read_parquet("universo2_caracteristicas.parquet").set_index(
    ["fecha", "ticker"])
def panel(col):
    return car[col].unstack().reindex(R.index)

pn, cap = panel("patrimonio_neto"), panel("capitalizacion")
datos = {
    "R": R,
    "value": np.log(pn.where(pn > 0) / cap),
    "esg": panel("puntuacion_esg").astype(float),
    "horquilla": panel("horquilla_pb"),
    "adv_usd": panel("volumen_medio_21d") * panel("cierre"),
    "cotiza": panel("cierre").notna(),
}
datos["value"].iloc[-1].describe().round(2)

def z_rango(df):
    """Rank-based normal score on each date."""
    rk = df.rank(axis=1)
    n = df.notna().sum(axis=1)
    u = (rk - 0.5).div(n, axis=0).to_numpy()
    return pd.DataFrame(stats.norm.ppf(u), index=df.index,
                        columns=df.columns)

def momentum(R, L):
    """Return of months m-L+1 to m-1 (skips the last one)."""
    with np.errstate(divide="ignore"):    # log1p(-1) = -inf: bankruptcy
        lr = np.log1p(R)
    return np.expm1(lr.shift(1).rolling(L - 1).sum())

def largo_corto(puntos, puede, q):
    """Long the top q and short the bottom q, equally weighted."""
    rk = puntos.where(puede).rank(axis=1, pct=True)
    lar, cor = (rk > 1 - q).astype(float), (rk <= q).astype(float)
    return (lar.div(lar.sum(axis=1), axis=0)
            - cor.div(cor.sum(axis=1), axis=0)).fillna(0)

INICIO = "2007-08"   # the longest window (18 months) already exists

def backtest(d, puntos, q, universo, coste=1.0, retardo=1):
    """Weights at the close of month t; return of month t + retardo."""
    puede = d["cotiza"] & d["cotiza"].columns.isin(list(universo))
    w = largo_corto(puntos, puede, q)
    bruto = (w.shift(retardo) * d["R"]).sum(axis=1)
    giro = (w - w.shift(1)).abs()
    pb = d["horquilla"].fillna(20) / 2 / 1e4      # half spread
    gasto = (giro * pb).sum(axis=1).shift(retardo) * coste
    return (bruto - gasto).loc[INICIO:], w.loc[INICIO:]

def sharpe(x):
    x = np.asarray(x, dtype=float)
    return x.mean() / x.std(ddof=1) * np.sqrt(12)

def senal_vm(d, L):
    return z_rango(momentum(d["R"], L)) + z_rango(d["value"])

TODOS = set(R.columns)                       # point-in-time universe
rejilla_a = pd.DataFrame({
    (L, q): backtest(datos, senal_vm(datos, L), q, TODOS)[0]
    for L in (3, 6, 9, 12, 15, 18) for q in (0.1, 0.2, 0.3, 0.4, 0.5)})
CORTE = "2016-07"                            # in sample: 2007-08 to 2016-07
sr_in = rejilla_a.loc[:CORTE].apply(sharpe)
elegido_a = sr_in.idxmax()
print(f"Chosen in sample: L = {elegido_a[0]} months, "
      f"q = {elegido_a[1]} (Sharpe {sr_in.max():.2f})")
sr_in.unstack().round(2)

def senal_esg(d, k, signo):
    return signo * d["esg"].rolling(k).mean()

sect = pd.read_csv("universo2_sectores.csv")
SUPERVIVIENTES = set(sect.loc[sect["ultima_sesion"].isna(), "ticker"])
rejilla_b = pd.DataFrame({
    (k, s, q): backtest(datos, senal_esg(datos, k, s), q,
                        SUPERVIVIENTES)[0]
    for k in (1, 2, 3, 6, 9, 12) for s in (1, -1)
    for q in (0.1, 0.2, 0.3, 0.4, 0.5)})
sr_in_b = rejilla_b.loc[:CORTE].apply(sharpe)
elegido_b = sr_in_b.idxmax()
print(f"{len(SUPERVIVIENTES)} survivors out of {len(TODOS)}")
print(f"Chosen: {elegido_b[0]}-month average, sign {elegido_b[1]}, "
      f"q = {elegido_b[2]} (Sharpe {sr_in_b.max():.2f})")
sr_in_b.unstack([1, 2]).round(2).iloc[:, :5]

@dataclass
class ResultadoBacktest:
    nombre: str
    rend: pd.Series             # monthly, net of costs
    pesos: pd.DataFrame         # at each month's close
    ejecutar: object            # (datos, coste, retardo) -> rend
    senal: object               # datos -> scores
    rejilla: pd.DataFrame       # ALL the combinations tried
    elegido: tuple
    corte: str                  # last in-sample month
    n_pruebas: int              # including those not written down
    universo: set
    fijos: tuple = ()           # grid axes that have no order
    capital: float = 50e6       # planned capital (USD)

def resultado(nombre, senal, q, universo, rejilla, elegido, n,
              fijos=()):
    rend, w = backtest(datos, senal(datos), q, universo)
    return ResultadoBacktest(
        nombre, rend, w,
        lambda d, c=1.0, t=1: backtest(d, senal(d), q, universo,
                                       c, t)[0],
        senal, rejilla, elegido, CORTE, n, universo, fijos)

L, q = elegido_a
A = resultado("value + momentum", lambda d: senal_vm(d, L), q,
              TODOS, rejilla_a, elegido_a, rejilla_a.shape[1])
k, s, qb = elegido_b
B = resultado("placebo ESG", lambda d: senal_esg(d, k, s), qb,
              SUPERVIVIENTES, rejilla_b, elegido_b, rejilla_b.shape[1],
              fijos=(1,))                    # the sign has no neighbors
for res in (A, B):
    print(f"{res.nombre}: {len(res.rend)} months, Sharpe "
          f"{sharpe(res.rend):.2f}, {res.n_pruebas} tests")

def prueba(valor, umbral, pasa, texto):
    return {"valor": valor, "umbral": umbral, "pasa": bool(pasa),
            "detalle": texto}

def p_ejecucion(res, fechas=("2010-06", "2015-12", "2020-03")):
    """Truncation: the signal for d does not change if what comes after d
    is deleted. Lag: with one more month of delay, the edge does not vanish."""
    completa = res.senal(datos)
    iguales = all(
        np.allclose(res.senal({c: v.loc[:f] for c, v in datos.items()})
                    .loc[f].to_numpy(), completa.loc[f].to_numpy(),
                    equal_nan=True) for f in fechas)
    s1, s2 = sharpe(res.rend), sharpe(res.ejecutar(datos, 1.0, 2))
    return prueba(s2 / s1, ">= 0.5 and truncation identical",
                  iguales and s2 >= 0.5 * s1,
                  f"truncation {'identical' if iguales else 'DIFFERENT'}; "
                  f"Sharpe t+1 {s1:.2f}, t+2 {s2:.2f}")

def p_costes(res):
    s = sharpe(res.ejecutar(datos, 2.0, 1))
    return prueba(s, ">= 0.5", s >= 0.5,
                  f"Sharpe with costs x2 {s:.2f}")

def p_fuera(res):
    r = res.rejilla[res.elegido]
    s_in, s_out = sharpe(r.loc[:res.corte]), sharpe(r.loc[res.corte:]
                                                     .iloc[1:])
    return prueba(s_out / s_in, ">= 0.5", s_out >= 0.5 * s_in > 0,
                  f"in sample {s_in:.2f}, out of sample {s_out:.2f}")

for res in (A, B):
    for f in (p_ejecucion, p_costes, p_fuera):
        print(f"{res.nombre:>16} | {f.__name__:<11} | "
              f"{f(res)['detalle']}")

EULER = np.euler_gamma    # Euler-Mascheroni constant

def p_dsr(res):
    """Deflated Sharpe Ratio (Bailey and López de Prado, 2014)."""
    r = res.rend.to_numpy()
    T, sr = len(r), r.mean() / r.std(ddof=1)        # per month
    g3, g4 = stats.skew(r), stats.kurtosis(r, fisher=False)
    pruebas = res.rejilla.apply(lambda c: c.mean() / c.std(ddof=1))
    N = res.n_pruebas
    sr0 = np.sqrt(pruebas.var(ddof=1)) * (
        (1 - EULER) * stats.norm.ppf(1 - 1 / N)
        + EULER * stats.norm.ppf(1 - 1 / (N * np.e)))
    z = (sr - sr0) * np.sqrt(T - 1) / np.sqrt(
        1 - g3 * sr + (g4 - 1) / 4 * sr ** 2)
    dsr = stats.norm.cdf(z)
    return prueba(dsr, ">= 0.95", dsr >= 0.95,
                  f"Sharpe {sr * np.sqrt(12):.2f} against a maximum "
                  f"expected by chance of {sr0 * np.sqrt(12):.2f}")

def p_pbo(res, S=10):
    """CSCV: half of the blocks in sample, the other half out of sample."""
    M = res.rejilla.to_numpy()
    T, N = M.shape
    bloques = np.array_split(np.arange(T), S)
    logits, pierde = [], []
    for dentro in itertools.combinations(range(S), S // 2):
        i = np.concatenate([bloques[j] for j in dentro])
        o = np.concatenate([bloques[j] for j in range(S)
                            if j not in dentro])
        sr_i = M[i].mean(0) / M[i].std(0, ddof=1)
        sr_o = M[o].mean(0) / M[o].std(0, ddof=1)
        n = sr_i.argmax()
        w = stats.rankdata(sr_o)[n] / (N + 1)
        logits.append(np.log(w / (1 - w)))
        pierde.append(sr_o[n] <= 0)
    pbo, p_perd = np.mean(np.array(logits) <= 0), np.mean(pierde)
    return prueba(p_perd, "loss <= 10%", p_perd <= 0.10,
                  f"PBO {pbo:.2f}; the chosen one loses out of sample in "
                  f"{p_perd:.0%} of the {len(logits)} partitions")

for res in (A, B):
    for f in (p_dsr, p_pbo):
        print(f"{res.nombre:>16} | {f.__name__:<5} | "
              f"{f(res)['detalle']}")

def p_meseta(res):
    """Neighbors of the chosen one in the in-sample grid."""
    sr_in = res.rejilla.loc[:res.corte].apply(sharpe)
    ejes = [sorted({k[i] for k in sr_in.index})
            for i in range(len(res.elegido))]
    pos = {k: [ejes[i].index(v) for i, v in enumerate(k)]
           for k in sr_in.index}
    yo = pos[res.elegido]
    vec = [sr_in[k] for k, p in pos.items() if k != res.elegido
           and max(abs(a - b) for a, b in zip(p, yo)) <= 1
           and all(k[i] == res.elegido[i] for i in res.fijos)]
    ratio = np.mean(vec) / sr_in[res.elegido]
    return prueba(ratio, ">= 0.7", ratio >= 0.7,
                  f"{len(vec)} neighbors: mean {np.mean(vec):.2f} "
                  f"against {sr_in[res.elegido]:.2f}")

def p_capacidad(res, limite=0.10):
    """Capital at which the 95th percentile of the orders is 10%
    of the stock's average daily volume."""
    giro = (res.pesos - res.pesos.shift(1)).abs()
    part = (giro / datos["adv_usd"].loc[INICIO:]).stack()
    p95 = part[part > 0].quantile(0.95)
    cap = limite / p95
    return prueba(cap / 1e6, f">= {res.capital / 1e6:.0f} M USD",
                  cap >= res.capital,
                  f"capacity {cap / 1e6:.0f} M USD")

def p_supervivencia(res):
    """Does the universe include the stocks that stopped trading?"""
    fuera = set(deslist["ticker"]) - res.universo
    return prueba(len(fuera), "0 excluded", not fuera,
                  f"{len(fuera)} delisted stocks outside the universe")

for res in (A, B):
    for f in (p_meseta, p_capacidad, p_supervivencia):
        print(f"{res.nombre:>16} | {f.__name__:<15} | "
              f"{f(res)['detalle']}")

PRUEBAS = {"ejecución t+1": p_ejecucion, "costes x2": p_costes,
           "fuera de muestra": p_fuera, "Deflated Sharpe": p_dsr,
           "PBO y pérdida": p_pbo, "meseta": p_meseta,
           "capacidad": p_capacidad, "supervivencia": p_supervivencia}

def auditar(res):
    filas = {n: f(res) for n, f in PRUEBAS.items()}
    tabla = pd.DataFrame(filas).T[["valor", "umbral", "pasa"]]
    veredicto = "APPROVED" if tabla["pasa"].all() else "NOT APPROVED"
    return tabla, veredicto

informe = {}
for res in (A, B):
    tabla, veredicto = auditar(res)
    informe[res.nombre] = tabla
    print(f"{res.nombre}: {veredicto} "
          f"({tabla['pasa'].sum()} of {len(tabla)} tests)")
pd.concat({n: t["pasa"].map({True: "pass", False: "FAIL"})
           for n, t in informe.items()}, axis=1)

t = informe["placebo ESG"].copy()
t["valor"] = t["valor"].astype(float).round(2)
t

con_todos, _ = backtest(datos, senal_esg(datos, k, s), qb, TODOS)
print(f"Placebo on survivors: {sharpe(B.rend):.2f}; "
      f"with the full universe: {sharpe(con_todos):.2f}")