Research agents: hypotheses, alpha mining and review.

import ast
import re
import warnings

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

precios = pd.read_parquet("acciones6_precios.parquet")
C = precios.pivot(index="fecha", columns="ticker", values="cierre")
V = precios.pivot(index="fecha", columns="ticker",
                  values="volumen").astype(float)
R = C.pct_change()
sesiones = C.index

noticias = pd.read_parquet("noticias6_corpus.parquet")
historia = pd.read_csv("acciones6_tickers.csv")
empresas = pd.read_csv("acciones6_empresas.csv")
actual = empresas.set_index("id_empresa")["ticker"]
a_actual = dict(zip(historia["ticker"],
                    historia["id_empresa"].map(actual)))
noticias["tk"] = noticias["ticker"].map(a_actual)

# Close of each session: 16:00 New York time, except the 15 half
# sessions of 2019-2025 (13:00), according to exchange_calendars 4.13.2.
MEDIAS = pd.to_datetime([
    "2019-07-03", "2019-11-29", "2019-12-24", "2020-11-27",
    "2020-12-24", "2021-11-26", "2022-11-25", "2023-07-03",
    "2023-11-24", "2024-07-03", "2024-11-29", "2024-12-24",
    "2025-07-03", "2025-11-28", "2025-12-24"])
hora = np.where(sesiones.isin(MEDIAS), 13, 16)
cierres = (sesiones + pd.to_timedelta(hora, unit="h")).tz_localize(
    "America/New_York").tz_convert("UTC")

lex = pd.read_csv("texto6_lexico.csv").query("lista == 'financiera'")
POS = set(lex.loc[lex["categoria"] == "positiva", "palabra"])
NEG = set(lex.loc[lex["categoria"] == "negativa", "palabra"])

def tono(texto):
    w = re.findall(r"[a-záéíóúüñ]+", texto.lower())
    p, n = sum(x in POS for x in w), sum(x in NEG for x in w)
    return (p - n) / (p + n + 1)

noticias["tono"] = [tono(t + " " + c) for t, c in
                    zip(noticias["titular"], noticias["cuerpo"])]
print(f"{C.shape[1]} stocks, {len(sesiones)} sessions; "
      f"{len(noticias)} news items, "
      f"{noticias['tk'].isna().sum()} without a company")

def tono_por_sesion(noticias):
    """Each news item counts in the first session that closes AFTER
    it is published: what is published after the close, in the next."""
    i = np.searchsorted(cierres.values, noticias["ts"].values,
                        side="right")
    ok = i < len(sesiones)
    n = noticias[ok].assign(sesion=sesiones[i[ok]])
    return (n.groupby(["sesion", "tk"])["tono"].sum().unstack()
            .reindex(index=sesiones, columns=C.columns).fillna(0.0))

T = tono_por_sesion(noticias)
i = np.searchsorted(cierres.values, noticias["ts"].values,
                    side="right")
fecha_ny = (noticias["ts"].dt.tz_convert("America/New_York")
            .dt.tz_localize(None).dt.normalize())
distinta = (sesiones[np.minimum(i, len(sesiones) - 1)]
            != fecha_ny.values).mean()
print(f"news items whose session is not their New York date: "
      f"{distinta:.0%}")

semanas = pd.Series(sesiones, index=sesiones).groupby(
    sesiones.to_period("W")).max()
W = pd.DatetimeIndex(semanas.values)
W = W[(W >= "2019-01-04") & (W <= "2025-12-26")]
fin = C.reindex(semanas.values)        # close of each week
Y = (fin.shift(-1) / fin - 1).reindex(W).to_numpy()  # the next one

def g(df):
    """Weekly cross-sectional normal rank."""
    x = df.reindex(W).to_numpy(float)
    n = np.isfinite(x).sum(axis=1, keepdims=True)
    return np.nan_to_num(norm.ppf(
        (rankdata(x, axis=1, nan_policy="omit") - 0.5) / n))

BASE = {
    "ret": lambda n: C / C.shift(n) - 1,
    "vol": lambda n: R.rolling(n).std(),
    "volu": lambda n: np.log(V.rolling(n).mean()
                             / V.rolling(63).mean()),
    "tono": lambda n: T.rolling(n).sum(),
}
_cache = {}

def evalua(expr):
    """Evaluates a signal of the language; error if it is not in it."""
    def nodo(x):
        if isinstance(x, ast.Expression):
            return nodo(x.body)
        if isinstance(x, ast.UnaryOp) and isinstance(x.op, ast.USub):
            return -nodo(x.operand)
        if isinstance(x, ast.BinOp) and isinstance(
                x.op, (ast.Sub, ast.Mult, ast.Div)):
            a, b = nodo(x.left), nodo(x.right)
            if isinstance(x.op, ast.Sub):
                return a - b
            if isinstance(x.op, ast.Mult):
                return a * b
            return a / np.where(np.abs(b) < 0.1, np.nan, b)
        if (isinstance(x, ast.Call) and isinstance(x.func, ast.Name)
                and x.func.id in BASE and len(x.args) == 1
                and isinstance(x.args[0], ast.Constant)
                and type(x.args[0].value) is int
                and 1 <= x.args[0].value <= 252 and not x.keywords):
            clave = (x.func.id, x.args[0].value)
            if clave not in _cache:
                _cache[clave] = g(BASE[clave[0]](clave[1]))
            return _cache[clave]
        raise ValueError(f"outside the language: {ast.dump(x)[:40]}")
    return nodo(ast.parse(expr.strip(), mode="eval"))

def valida(linea):
    try:
        evalua(linea)
        return True
    except (SyntaxError, ValueError):
        return False

print(len(W), "weeks;", valida("-ret(5)"), valida("ret(5) + 1"),
      valida("**Momentum**: ret(252)"))

import torch
from transformers import AutoModelForCausalLM, AutoTokenizer

torch.manual_seed(0)
torch.use_deterministic_algorithms(True)
MODELO = "Qwen/Qwen2.5-0.5B-Instruct"
tok = AutoTokenizer.from_pretrained(MODELO)
llm = AutoModelForCausalLM.from_pretrained(MODELO,
                                           dtype=torch.float32)

def genera(mensajes, max_tokens):
    x = tok.apply_chat_template(mensajes, add_generation_prompt=True,
                                return_tensors="pt", return_dict=True)
    with torch.inference_mode():
        y = llm.generate(**x, do_sample=False,
                         max_new_tokens=max_tokens,
                         pad_token_id=tok.eos_token_id)
    n = x["input_ids"].shape[1]
    return tok.decode(y[0, n:], skip_special_tokens=True)

INSTRUCCIONES = (
    "Propones señales para ordenar acciones. Usa SOLO este lenguaje: "
    "ret(n) rentabilidad de n sesiones; vol(n) volatilidad de n "
    "sesiones; volu(n) volumen medio de n sesiones; tono(n) tono "
    "medio de las noticias de n sesiones; operadores - * / y "
    "numeros. Escribe 5 señales, una por linea, sin explicaciones.")
texto = genera([{"role": "system", "content": INSTRUCCIONES},
                {"role": "user", "content":
                 "Hipotesis para un horizonte de una semana."}], 80)
lineas = [l for l in texto.splitlines() if l.strip()]
for l in lineas:
    print(l[:86] + ("…" if len(l) > 86 else ""))
print(f"valid: {sum(map(valida, lineas))} of {len(lineas)}")

EJEMPLOS = "\n".join(["-ret(5)", "ret(252)-ret(21)", "tono(5)",
                      "-vol(63)", "ret(5)*volu(5)"])
texto = genera([
    {"role": "system", "content":
     "Escribe señales en este lenguaje: ret(n), vol(n), volu(n), "
     "tono(n), con - * /. Solo expresiones, una por linea."},
    {"role": "user", "content": "5 señales"},
    {"role": "assistant", "content": EJEMPLOS},
    {"role": "user", "content": "5 señales distintas"}], 60)
lineas = [l.strip() for l in texto.splitlines() if l.strip()]
unicas = list(dict.fromkeys(l for l in lineas if valida(l)))
copias = [l for l in unicas if l in EJEMPLOS.split("\n")]
print(lineas)
print(f"valid {sum(map(valida, lineas))}, distinct {len(unicas)},"
      f" copied from the examples {len(copias)}: {copias}")

FUNCIONES = ([f"ret({n})" for n in (1, 5, 10, 21, 63, 126, 252)]
             + [f"vol({n})" for n in (5, 21, 63, 126)]
             + ["volu(5)", "volu(21)", "tono(5)", "tono(21)"])
espacio = [f for f in FUNCIONES] + [f"-{f}" for f in FUNCIONES]
for i, a in enumerate(FUNCIONES):
    for b in FUNCIONES[i + 1:]:
        espacio += [f"{a}*{b}", f"-{a}*{b}", f"{a}-{b}", f"{b}-{a}"]
rng = np.random.default_rng(646)
candidatas = [str(c) for c in
              rng.choice(espacio, size=200, replace=False)]
print(f"{len(espacio)} possible signals; {len(candidatas)} chosen")
print(candidatas[:6])

Yr = rankdata(Y, axis=1)
Yc = Yr - Yr.mean(axis=1, keepdims=True)

def ic_semanal(score):
    s = rankdata(np.nan_to_num(score, nan=0.0), axis=1)
    s = s - s.mean(axis=1, keepdims=True)
    den = np.sqrt((s * s).sum(1) * (Yc * Yc).sum(1))
    with np.errstate(invalid="ignore", divide="ignore"):
        return (s * Yc).sum(1) / den

def t_ic(x):
    x = x[np.isfinite(x)]
    return x.mean() / x.std(ddof=1) * np.sqrt(len(x))

IC = pd.DataFrame({c: ic_semanal(evalua(c)) for c in candidatas},
                  index=W)
corte = pd.Timestamp("2022-06-30")
descubre = IC.index < corte
embargo = W[descubre][-1:]   # its label ends in the 1st validation week
descubre &= ~IC.index.isin(embargo)
valida_ = IC.index > corte
res = pd.DataFrame({
    "IC desc.": IC[descubre].mean(), "t desc.": IC[descubre].apply(
        lambda x: t_ic(x.to_numpy())),
    "IC valid.": IC[valida_].mean(), "t valid.": IC[valida_].apply(
        lambda x: t_ic(x.to_numpy()))})
print(f"discovery: {descubre.sum()} weeks up to "
      f"{W[descubre][-1].date()}; validation: {valida_.sum()} from "
      f"{W[valida_][0].date()}")
res.sort_values("t desc.", ascending=False).head(10).round(3)

pd.DataFrame({c: {"t desc.": t_ic(ic_semanal(evalua(c))[descubre]),
                  "t valid.": t_ic(ic_semanal(evalua(c))[valida_])}
              for c in unicas}).T.round(2)

from scipy.stats import t as t_dist

m, n_d = len(res), int(descubre.sum())
p = 2 * t_dist.sf(res["t desc."].abs(), n_d - 1)
bh = np.flatnonzero(np.sort(p) <= 0.10 * np.arange(1, m + 1) / m)
print(f"|t| > 1.96: {(res['t desc.'].abs() > 1.96).sum()} of {m}"
      f" (by chance {0.05 * m:.0f} would be expected)")
print(f"Bonferroni at 5%: {(p < 0.05 / m).sum()}; "
      f"Benjamini-Hochberg at 10%: {bh.max() + 1 if len(bh) else 0}")

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

autoval = np.linalg.eigvalsh(IC[descubre].fillna(0).corr())
n_ef = autoval.sum() ** 2 / (autoval ** 2).sum()
print(f"effective trials: {n_ef:.0f}; maximum t expected with no "
      f"signal at all: {max_esperado(m):.2f} (N = {m}), "
      f"{max_esperado(n_ef):.2f} (N = {n_ef:.0f})")

verdad = pd.read_parquet("panel5_verdad.parquet").pivot(
    index="fecha", columns="ticker", values="alfa_verdadero")
ic_or = pd.Series(ic_semanal(
    verdad.reindex(index=W, columns=C.columns).to_numpy()), index=W)
print(f"level 5 oracle: IC {ic_or[descubre].mean():.3f} "
      f"(t {t_ic(ic_or[descubre].to_numpy()):.2f}) in discovery,"
      f" {ic_or[valida_].mean():.3f} in validation")

def larga_corta(score, q=0.2):
    """Top quintile minus bottom quintile, equally weighted, weekly."""
    rk = rankdata(np.nan_to_num(score), axis=1) / score.shape[1]
    n = score.shape[1]
    w = ((rk > 1 - q).astype(float) - (rk <= q)) / (q * n)
    return pd.Series((w * np.nan_to_num(Y)).sum(1), index=W)

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

carteras = pd.DataFrame({c: larga_corta(evalua(c))
                         for c in candidatas})
sr = carteras[descubre].mean() / carteras[descubre].std()
mejor = res["t desc."].idxmax()
x = carteras.loc[descubre, mejor]
for n in (1, n_ef, m):
    sr0 = np.sqrt(sr.var()) * max_esperado(n) if n > 1 else 0.0
    print(f"N = {n:5.0f}: threshold {sr0 * np.sqrt(52):.2f} annual,"
          f" DSR = {psr(sr[mejor], sr0, x):.3f}")
for nombre, z in (("discovery", x),
                  ("validation", carteras.loc[valida_, mejor])):
    print(f"{mejor}, {nombre}: annual Sharpe "
          f"{z.mean() / z.std() * np.sqrt(52):.2f}")

bloques = np.array_split(np.arange(len(W)), 6)
filas = []
for k, prueba in enumerate(bloques):
    entrena = np.ones(len(W), bool)
    entrena[max(prueba[0] - 1, 0):prueba[-1] + 2] = False  # embargo
    t_ent = IC[entrena].apply(lambda x: t_ic(x.to_numpy()))
    elegida = t_ent.idxmax()
    filas.append({"bloque": f"{W[prueba[0]]:%Y-%m} to "
                            f"{W[prueba[-1]]:%Y-%m}",
                  "elegida": elegida,
                  "t entrena": round(t_ent.max(), 2),
                  "IC prueba": round(IC.iloc[prueba][elegida].mean(),
                                     4),
                  "IC tono(5)": round(IC.iloc[prueba]["tono(5)"].mean(),
                                      4)})
cv = pd.DataFrame(filas).set_index("bloque")
print(cv)
print(f"mean: chosen {cv['IC prueba'].mean():.4f}; "
      f"tono(5) {cv['IC tono(5)'].mean():.4f}")

propuesta_a = """
def tono_suave(T):
    # smooths the tone to remove noise
    return T.rolling(5, center=True).sum()
"""
propuesta_b = """
def tono_diario(noticias):
    # tone of each New York day, per company
    fecha = (noticias["ts"].dt.tz_convert("America/New_York")
             .dt.tz_localize(None).dt.normalize())
    return (noticias.groupby([fecha, "tk"])["tono"].sum().unstack()
            .reindex(index=sesiones, columns=C.columns).fillna(0.0))
"""

def revision_estatica(codigo):
    """Look-ahead patterns that can be seen in the code's text."""
    avisos = []
    for n in ast.walk(ast.parse(codigo)):
        if isinstance(n, ast.keyword) and n.arg == "center":
            avisos.append(f"line {n.value.lineno}: center=")
        if (isinstance(n, ast.Call)
                and isinstance(n.func, ast.Attribute)
                and n.func.attr == "shift" and n.args
                and isinstance(n.args[0], ast.UnaryOp)):
            avisos.append(f"line {n.lineno}: negative shift")
        if isinstance(n, ast.Attribute) and n.attr in ("bfill",
                                                       "backfill"):
            avisos.append(f"line {n.lineno}: {n.attr}")
    return avisos or ["no warnings"]

for nombre, codigo in (("A", propuesta_a), ("B", propuesta_b)):
    print(nombre, revision_estatica(codigo))

espacio_agente = {"np": np, "pd": pd, "sesiones": sesiones, "C": C}
exec(propuesta_b, espacio_agente)
tono_diario = espacio_agente["tono_diario"]

def prueba_pit(funcion, fechas):
    """What is computed for date D with all the news has to match
    what is computed with the news published before D's close.
    Returns the dates and the cells that change."""
    completo = funcion(noticias)
    malas, celdas = 0, 0
    for d in fechas:
        corte = cierres[sesiones.get_loc(d)]
        parcial = funcion(noticias[noticias["ts"] < corte])
        dif = ~np.isclose(completo.loc[d], parcial.loc[d])
        malas, celdas = malas + dif.any(), celdas + dif.sum()
    return malas, celdas

muestra = W[::18]                       # 21 Fridays spread out
for nombre, f in (("tono_por_sesion", tono_por_sesion),
                  ("tono_diario (B)", tono_diario)):
    malas, celdas = prueba_pit(f, muestra)
    print(f"{nombre:<16} dates with differences: {malas} of "
          f"{len(muestra)}; cells: {celdas}")

for nombre, T_ in (("correct", T),
                   ("proposal B", tono_diario(noticias))):
    ic = ic_semanal(g(T_.rolling(5).sum()))
    print(f"tono(5) {nombre:<12} IC {np.nanmean(ic):.4f}  "
          f"t {t_ic(ic):.2f}")