Capstone C3: stock and ETF portfolio with ML, HRP and Spanish taxation.

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

feat = pd.read_parquet("panel5_features.parquet")
etiq = pd.read_parquet("panel5_etiquetas.parquet",
                       columns=["fecha", "ticker", "ret_sig_21d",
                                "fin_sig_21d"])
d = feat.merge(etiq, on=["fecha", "ticker"])
# the last panel date of each calendar month: rebalancing day
fin_mes = d.groupby(d["fecha"].dt.to_period("M"))["fecha"].max()
d = d[d["fecha"].isin(fin_mes)].sort_values(["fecha", "ticker"])
TODAS = list(feat.columns[3:])
fechas = pd.DatetimeIndex(d["fecha"].unique())
nM, nA = len(fechas), d["ticker"].nunique()
tickers = d["ticker"].unique()[:nA]
F = d[TODAS].to_numpy(float).reshape(nM, nA, -1)
Y = d["ret_sig_21d"].to_numpy().reshape(nM, nA)
fin = d["fin_sig_21d"].to_numpy().reshape(nM, nA)[:, 0]


def rango_normal(X):
    """Cross-sectional normal rank of each month (5.2.4); no data, 0."""
    n = np.isfinite(X).sum(axis=1, keepdims=True)
    with np.errstate(invalid="ignore", divide="ignore"):
        u = (rankdata(X, axis=1, nan_policy="omit") - 0.5) / n
        return np.nan_to_num(norm.ppf(u))


G = np.stack([rango_normal(F[:, :, j]) for j in range(len(TODAS))], 2)
# audit: monthly IC of each feature in 2012-2014 (only the first
# training window) and its t; a huge t is an alarm
tramo = (fechas.year < 2015) & np.isfinite(Y).all(axis=1)
RY = rango_normal(Y)
ic = np.array([[np.corrcoef(G[m, :, j], RY[m])[0, 1]
                for m in np.flatnonzero(tramo)]
               for j in range(len(TODAS))])
t = ic.mean(1) / ic.std(1, ddof=1) * np.sqrt(ic.shape[1])
auditoria = pd.DataFrame({"ic": ic.mean(1), "t": t}, index=TODAS)
print(f"{nM} months x {nA} stocks; audited window: "
      f"{tramo.sum()} months")
print(auditoria.sort_values("t", key=abs, ascending=False)
      .head(4).round(3).to_string())
SOSPECHOSAS = auditoria.index[auditoria["t"].abs() > 5].tolist()
FE = [j for j, c in enumerate(TODAS) if c not in SOSPECHOSAS]
print("out:", SOSPECHOSAS, "| remaining", len(FE), "features")

import lightgbm as lgb
from sklearn.linear_model import Ridge

ok = np.isfinite(Y).all(axis=1)
RY = np.where(ok[:, None], rango_normal(Y), np.nan)  # label
quintil = np.where(ok[:, None], (rankdata(Y, axis=1) - 1) * 5 // nA,
                   -1).astype(int)
ARB = dict(n_estimators=200, learning_rate=0.03, num_leaves=15,
           min_child_samples=200, subsample=0.7, subsample_freq=1,
           colsample_bytree=0.7, random_state=0, n_jobs=4,
           deterministic=True, force_row_wise=True, verbose=-1)
anyo = fechas.year.to_numpy()
punt = {k: np.full((nM, nA), np.nan)
        for k in ("ridge", "lgbm_reg", "lgbm_rank")}
for y in range(2015, 2026):
    hoy = np.flatnonzero(anyo == y)
    # time rule: labels closed before the first rebalance of the
    # year (purging the overlapping months)
    tr = np.flatnonzero(ok & (fin < fechas[hoy[0]]))
    X = G[tr][:, :, FE].reshape(-1, len(FE))
    mod = {"ridge": Ridge(alpha=10.0).fit(X, RY[tr].ravel()),
           "lgbm_reg": lgb.LGBMRegressor(**ARB).fit(X, RY[tr].ravel()),
           "lgbm_rank": lgb.LGBMRanker(
               objective="lambdarank", label_gain=[0, 1, 2, 3, 4],
               lambdarank_truncation_level=40, **ARB).fit(
               X, quintil[tr].ravel(), group=[nA] * len(tr))}
    for k, m in mod.items():
        punt[k][hoy] = m.predict(G[hoy][:, :, FE].reshape(
            -1, len(FE))).reshape(len(hoy), nA)
oos = np.flatnonzero((anyo >= 2015) & ok)
verdad = pd.read_parquet("panel5_verdad.parquet",
                         columns=["fecha", "ticker", "alfa_verdadero"])
verdad = verdad[verdad["fecha"].isin(fechas)].sort_values(
    ["fecha", "ticker"])
punt["oraculo"] = verdad["alfa_verdadero"].to_numpy().reshape(nM, nA)
vol = F[:, :, TODAS.index("vol_63d")]


def corr_r(a, b):
    return np.mean([np.corrcoef(rankdata(a[i]), rankdata(b[i]))[0, 1]
                    for i in oos])


filas = {}
for k, p in punt.items():
    x = np.array([np.corrcoef(rankdata(p[i]), rankdata(Y[i]))[0, 1]
                  for i in oos])
    filas[k] = {"IC": x.mean(),
                "t": x.mean() / x.std(ddof=1) * np.sqrt(len(x)),
                "corr_vol": corr_r(p, vol)}
print(f"{len(oos)} out-of-sample months (2015-2025)")
pd.DataFrame(filas).T.round(4)

from scipy.cluster.hierarchy import leaves_list, linkage
from scipy.spatial.distance import squareform

px = pd.read_parquet("panel5_precios.parquet",
                     columns=["fecha", "ticker", "cierre"])
px = px.pivot(index="fecha", columns="ticker", values="cierre")
px = px[tickers].astype(float)
ret = px.pct_change()
TOP, COSTE = 40, 0.001            # 40 stocks; 10 bp per side


def var_ivp(S, ix):
    c = S[np.ix_(ix, ix)]
    w = 1 / np.diag(c)
    w /= w.sum()
    return w @ c @ w


def hrp(S):
    """HRP from lesson 4.1.4 (single linkage, bisection)."""
    sd = np.sqrt(np.diag(S))
    dist = np.sqrt(np.clip((1 - S / np.outer(sd, sd)) / 2, 0, 1))
    orden = leaves_list(linkage(squareform(dist, checks=False),
                                "single"))
    w, grupos = np.ones(len(S)), [orden]
    while grupos:
        sig = []
        for g in grupos:
            if len(g) > 1:
                a, b = g[: len(g) // 2], g[len(g) // 2:]
                va, vb = var_ivp(S, a), var_ivp(S, b)
                w[a] *= vb / (va + vb)
                w[b] *= va / (va + vb)
                sig += [a, b]
        grupos = sig
    return w


REB = [i for i in range(nM) if anyo[i] >= 2015]
BANDA = 120                       # held while within the top 120


def seleccion(p, i, tenidas, vetadas=()):
    """The TOP best, keeping the ones already held while they stay
    within the top BANDA (5.4.4); `vetadas` do not come in."""
    orden = [j for j in np.argsort(-p[i]) if tickers[j] not in vetadas]
    rango = {j: n for n, j in enumerate(orden)}
    sel = sorted((j for j in tenidas if rango.get(j, nA) < BANDA),
                 key=rango.get)[:TOP]
    sel += [j for j in orden if j not in sel][:TOP - len(sel)]
    return np.array(sel)


def pesos(sel, i, metodo):
    w = np.zeros(nA)
    if metodo == "igual":
        w[sel] = 1 / len(sel)
    else:
        hist = ret.loc[:fechas[i]].iloc[-252:].to_numpy()[:, sel]
        w[sel] = hrp(np.cov(hist, rowvar=False))
    return w


def cartera(p, metodo):
    """Monthly return net of costs, and turnover."""
    w_ant, filas = np.zeros(nA), []
    for a, b in zip(REB[:-1], REB[1:]):
        w = pesos(seleccion(p, a, np.flatnonzero(w_ant)), a, metodo)
        rot = np.abs(w - w_ant).sum()
        r = (px.loc[fechas[b]] / px.loc[fechas[a]] - 1).to_numpy()
        bruto = float(w @ r)
        filas.append((fechas[b], bruto - COSTE * rot, rot))
        w_ant = w * (1 + r) / (1 + bruto)
    return pd.DataFrame(filas, columns=["fecha", "neto", "rot"]
                        ).set_index("fecha")


mercado = pd.Series(
    [(px.loc[fechas[b]] / px.loc[fechas[a]] - 1).mean()
     for a, b in zip(REB[:-1], REB[1:])],
    index=fechas[REB[1:]])                 # equal-weighted universe
pruebas = {f"{k} + {m}": cartera(punt[k], m)
           for k in ("ridge", "lgbm_reg", "lgbm_rank")
           for m in ("igual", "hrp")}
activo = pd.DataFrame({k: v["neto"] - mercado
                       for k, v in pruebas.items()})
registro = pd.DataFrame({
    "activo_anual_%": 100 * 12 * activo.mean(),
    "ir": activo.mean() / activo.std() * np.sqrt(12),
    "sharpe": [v["neto"].mean() / v["neto"].std() * np.sqrt(12)
               for v in pruebas.values()],
    "rotacion_mes_%": [100 * v["rot"].iloc[1:].mean()
                       for v in pruebas.values()]})
registro.to_csv("registro_pruebas.csv")
techo = {m: 12 * (cartera(punt["oraculo"], m)["neto"]
                  - mercado).mean() for m in ("igual", "hrp")}
print(f"{len(mercado)} months; universe: "
      f"{100 * ((1 + mercado).prod() ** (12 / len(mercado)) - 1):.1f}"
      f" % a year, Sharpe "
      f"{mercado.mean() / mercado.std() * np.sqrt(12):.2f}")
print("ceiling (oracle), annual active return:",
      {k: f"{100 * v:.1f} %" for k, v in techo.items()})
registro.round(2)

from collections import defaultdict, deque
from dateutil.relativedelta import relativedelta

DOS_MESES = relativedelta(months=2)


class Lotes:
    """FIFO lots by value (art. 37.2) with the two-month rule
    (art. 33.5 f, lesson 7.4.6), simplified: see the text."""

    def __init__(self):
        # lot: [date, units, cost, deferred, free units]
        self.cola = defaultdict(deque)
        self.pend = defaultdict(list)    # loss sales still to match
        self.lineas = []                 # (year, amount, type)
        self.encadenadas = 0

    def comprar(self, f, t, uds, coste):
        dif, libres = 0.0, uds
        for v in self.pend[t]:           # purchase after a sale
            if v["fecha"] + DOS_MESES >= f and v["uds"] > 0 and libres:
                u = min(v["uds"], libres)
                dif += u * v["perdida_ud"]
                self.encadenadas += v["encadenada"]
                self.lineas.append((v["fecha"].year,
                                    u * v["perdida_ud"], "difiere"))
                v["uds"] -= u
                libres -= u
        self.cola[t].append([f, uds, coste, dif, libres])

    def vender(self, f, t, uds, importe):
        cola, falta, coste, liberada = self.cola[t], uds, 0.0, 0.0
        while falta:
            lote = cola[0]
            u = min(falta, lote[1])
            x = u / lote[1]
            coste += lote[2] * x
            liberada += lote[3] * x
            lote[1] -= u
            lote[2] -= lote[2] * x
            lote[3] -= lote[3] * x
            lote[4] = min(lote[4], lote[1])
            falta -= u
            if lote[1] == 0:
                cola.popleft()
        res = importe - coste
        self.lineas.append((f.year, res, "venta"))
        if liberada:
            self.lineas.append((f.year, -liberada, "integra"))
            res -= liberada
        if res >= 0:
            return
        quedan = sum(l[1] for l in cola)
        n = 0                            # earlier purchases, up to
        for lote in reversed(cola):      # the remaining balance
            if lote[0] + DOS_MESES >= f and quedan:
                u = min(lote[4], quedan - n)
                lote[3] += -res / uds * u
                lote[4] -= u
                n += u
        if n:
            self.lineas.append((f.year, -res / uds * n, "difiere"))
            self.encadenadas += liberada > 0
        if uds > n:
            self.pend[t].append({"fecha": f, "uds": uds - n,
                                 "perdida_ud": -res / uds,
                                 "encadenada": liberada > 0})


# test: buy 100, sell 100 at a loss, buy back 60 a month later
L = Lotes()
L.comprar(pd.Timestamp("2025-01-10"), "ABCD", 100, 10_000.0)
L.vender(pd.Timestamp("2025-03-10"), "ABCD", 100, 8_000.0)
L.comprar(pd.Timestamp("2025-04-08"), "ABCD", 60, 4_900.0)
L.vender(pd.Timestamp("2025-09-15"), "ABCD", 60, 5_400.0)
for linea in L.lineas:
    print(linea)
print("year 2025:", round(sum(i for a, i, _ in L.lineas
                              if a == 2025), 2))

TRAMOS = [(6000, 0.19), (50000, 0.21), (200000, 0.23),
          (300000, 0.27), (np.inf, 0.30)]      # 2026 schedule


def cuota(base):
    c, desde = 0.0, 0.0
    for hasta, t in TRAMOS:
        c += max(0.0, min(base, hasta) - desde) * t
        desde = hasta
    return c


def impuesto(lineas, y, pend):
    """Tax due for the year on gains and losses (no dividends: prices
    are total return), with the losses of the previous four years
    (art. 49). Modifies `pend`."""
    saldo = sum(i for a, i, _ in lineas if a == y)
    for a in sorted(pend):
        if a < y - 4:
            del pend[a]
        elif saldo > 0:
            x = min(pend[a], saldo)
            pend[a] -= x
            saldo -= x
    if saldo < 0:
        pend[y] = -saldo
    return cuota(max(saldo, 0.0))


def simular(p, metodo, fiscal=True, bloqueo=False, capital=100_000.0):
    """Portfolio in whole shares: sells, buys, pays 10 bp per side
    and, at the first July rebalance, the previous year's IRPF. At
    the end it sells everything and settles the last year.
    `bloqueo`: what was sold at a loss is not bought or increased for
    two months (repurchase warning of 7.4.6)."""
    caja, uds, L, pend = capital, defaultdict(int), Lotes(), {}
    vetos, pagado, valor = {}, {}, {}
    for i in REB:
        f, precio = fechas[i], px.loc[fechas[i]]
        y = f.year - 1
        if fiscal and f.month >= 7 and y >= 2015 and y not in pagado:
            pagado[y] = impuesto(L.lineas, y, pend)
            caja -= pagado[y]
        valor[f] = caja + sum(n * precio[t] for t, n in uds.items())
        vet = {t for t, h in vetos.items() if h > f and bloqueo}
        if i == REB[-1]:
            obj = {}                     # final liquidation
        else:
            tenidas = [j for j, t in enumerate(tickers) if uds[t] > 0]
            fuera = [t for t in vet if uds[t] == 0]
            sel = seleccion(p, i, tenidas, fuera)
            w = pesos(sel, i, metodo)
            obj = {tickers[j]: int(w[j] * valor[f] * 0.995
                                   / precio.iloc[j]) for j in sel}
        for t in [t for t in uds if uds[t] > obj.get(t, 0)]:
            n = uds[t] - obj.get(t, 0)
            imp = n * precio[t] * (1 - COSTE)
            k = len(L.lineas)
            L.vender(f, t, n, imp)
            if any(x[2] == "difiere" for x in L.lineas[k:]) or (
                    L.pend[t] and L.pend[t][-1]["fecha"] == f):
                vetos[t] = f + DOS_MESES + pd.Timedelta(days=1)
            caja += imp
            uds[t] -= n
        for t, n_obj in obj.items():
            n = n_obj - uds[t]
            if n > 0 and t not in vet:
                coste = n * precio[t] * (1 + COSTE)
                L.comprar(f, t, n, coste)
                caja -= coste
                uds[t] += n
    if fiscal:
        pagado[f.year] = impuesto(L.lineas, f.year, pend)
        caja -= pagado[f.year]
    return pd.Series(valor), caja, L, pagado


def lineas_anuales(L, regla=True):
    """Balance of gains and losses per year; without the rule, like
    an engine that defers nothing."""
    tipos = ("venta", "difiere", "integra") if regla else ("venta",)
    return pd.Series([i for _, i, k in L.lineas if k in tipos],
                     index=[a for a, _, k in L.lineas if k in tipos]
                     ).groupby(level=0).sum()


t_res = {}
for nombre, kw in (("no IRPF", dict(fiscal=False)),
                   ("with IRPF", {}),
                   ("with IRPF and lockout", dict(bloqueo=True))):
    v, neto, L, pagado = simular(punt["lgbm_rank"], "hrp", **kw)
    anos = (v.index[-1] - v.index[0]).days / 365.25
    perd = -sum(i for _, i, k in L.lineas if k == "venta" and i < 0)
    t_res[nombre] = {
        "final_neto": neto,
        "tae_%": 100 * ((neto / v.iloc[0]) ** (1 / anos) - 1),
        "irpf": sum(pagado.values()),
        "perdidas": perd,
        "diferidas_%": 100 * sum(i for _, i, k in L.lineas
                                 if k == "difiere") / perd}
    if nombre == "with IRPF":
        L_ref, v_ref = L, v
print(pd.DataFrame(t_res).T.round(1).to_string())
con, sin = lineas_anuales(L_ref), lineas_anuales(L_ref, False)
print("Year's balance with the rule minus without it (€):")
print(pd.DataFrame({"": (con - sin).round(0).astype(int)}).T
      .to_string())

import hashlib
import json


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


def max_esperado(n):
    gm = np.euler_gamma
    return ((1 - gm) * norm.ppf(1 - 1 / n)
            + gm * norm.ppf(1 - 1 / (n * np.e)))


ir = activo.mean() / activo.std()          # per month
mejor = ir.idxmax()
x = activo[mejor]
sr0 = ir.std() * max_esperado(len(ir))
mom = (len(x), x.skew(), x.kurt() + 3)
dsr = psr(ir[mejor], sr0, *mom)
print(f"best of {len(ir)}: {mejor}, IR {ir[mejor] * np.sqrt(12):.2f};"
      f" threshold {sr0 * np.sqrt(12):.2f}; DSR {dsr:.3f}")
pre = "lgbm_rank + hrp"                     # the one of the brief
xp = activo[pre]
print(f"brief ({pre}): IR {ir[pre] * np.sqrt(12):.2f}, PSR(0) "
      f"{psr(ir[pre], 0, len(xp), xp.skew(), xp.kurt() + 3):.3f}")


def sha(ruta):
    return hashlib.sha256(open(ruta, "rb").read()).hexdigest()[:12]


informe = {
    "proyecto": "C3 ML portfolio, HRP and IRPF", "config": pre,
    "datos": {a: sha(a) for a in ("panel5_features.parquet",
                                  "panel5_precios.parquet")},
    "fuera": SOSPECHOSAS,
    "validacion": {"n_pruebas": len(ir), "mejor": mejor,
                   "dsr_mejor": round(float(dsr), 3),
                   "ir_encargo": round(float(ir[pre]) * 12**0.5, 3)},
    "fiscal": {k: {c: round(float(v), 1) for c, v in d.items()}
               for k, d in t_res.items()}}
with open("informe_c3.json", "w") as fh:
    json.dump(informe, fh, ensure_ascii=False, indent=1)
print(json.dumps(informe["validacion"], ensure_ascii=False))

from ib_async import IB, Order, Stock

ib = IB()
ib.connect("127.0.0.1", 4002, clientId=34)       # IB Gateway, paper
cartera = {p.contract.symbol: p.position for p in ib.positions()}
objetivo = calcular_objetivo(cartera)              # signal + HRP + band
ventas = {t: n - objetivo.get(t, 0) for t, n in cartera.items()
          if n > objetivo.get(t, 0)}
compras = {t: n - cartera.get(t, 0) for t, n in objetivo.items()
           if n > cartera.get(t, 0)}
for lado, lista in (("SELL", ventas), ("BUY", compras)):
    for t, n in sorted(lista.items()):
        c = ib.qualifyContracts(Stock(t, "SMART", "USD"))[0]
        orden = Order(action=lado, totalQuantity=n, orderType="MOC",
                      orderRef=f"C3-{t}-{lado}-20261030")
        ib.placeOrder(c, orden)                    # at the close
ib.sleep(5)
print([(t.contract.symbol, t.orderStatus.status) for t in ib.trades()])