"""Residual books. JKP Factors CTF submission. Each month a factor model is estimated from the past 252 daily returns. For scores s, the dollar- and factor-neutral mean-variance portfolio earns s'z, where z is the standardized GLS residual of next month's returns; z is the training target. Two books: "open" hedges the return factors, "neutral" also hedges the principal components of the characteristics. One LightGBM model per book, refit yearly. The books are combined with long-only max-Sharpe weights estimated from their past returns. CTF Admin Modifications (2026-10-01): -------------------------------------- 1. Lowered the minimum number of stocks per month from 100 to 50 when the CTF_EXECUTION_MODE environment variable is "validation". Reason: the validation dataset has about 96 stocks a month, so no test month reached 100 and the validation run returned no rows. Full runs still use 100, so the full-run weights are unchanged. """ from __future__ import annotations import os import time import lightgbm as lgb import numpy as np import pandas as pd SEED = 42 N_THREADS = min(32, os.cpu_count() or 1) TARGET = "ret_exc_lead1m" WINDOW_DAYS = 252 # one year of daily returns MIN_OBS = 200 # returns required in the window # CTF-FIX: the validation dataset has about 96 stocks a month, below the full-run minimum VALIDATION_RUN = os.environ.get("CTF_EXECUTION_MODE", "full").strip().lower() == "validation" MIN_NAMES = 50 if VALIDATION_RUN else 100 DAYS_PER_MONTH = 21.0 PSI_FLOOR = 0.05 LABEL_CLIP = 4.0 REFIT_EVERY = 12 MIN_LABELLED_MONTHS = 36 VALID_MONTHS = 60 LEARNING_RATE = 0.05 NUM_LEAVES = 31 MAX_ROUNDS = 2000 PATIENCE = 200 BOOKS = ("open", "neutral") BLEND_HALFLIFE = 60.0 BLEND_MIN_MONTHS = 36 BLEND_SHRINK = 0.5 SCALE = "risk" # "risk": unit predicted volatility; "gross": unit gross exposure def log(msg: str) -> None: print(msg, flush=True) def to_days(values) -> np.ndarray: return pd.to_datetime(values).to_numpy().astype("datetime64[D]").astype(np.int64) def mp_edge(p: int, n: int) -> float: """Upper edge of the Marchenko-Pastur law for a p-variable correlation matrix from n observations.""" return (1.0 + np.sqrt(p / n)) ** 2 # --------------------------------------------------------------------------------------- # Risk model and portfolio # --------------------------------------------------------------------------------------- def daily_matrix(daily_ret: pd.DataFrame, ids: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """Dense (days x ids) daily excess returns, NaN where missing; ids must be sorted and unique.""" did = daily_ret["id"].to_numpy(np.int64) col = np.minimum(np.searchsorted(ids, did), len(ids) - 1) keep = ids[col] == did day = to_days(daily_ret["date"].to_numpy()[keep]) days, row = np.unique(day, return_inverse=True) M = np.full((len(days), len(ids)), np.nan, dtype=np.float32) M[row, col[keep]] = daily_ret["ret_exc"].to_numpy(np.float32)[keep] return days, M def risk_model(R: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """Factor model of a (T, n) window of daily returns (NaN = missing). Returns the eligibility mask and, for eligible stocks, loadings B, specific variances psi and monthly volatilities sigma. """ obs = ~np.isnan(R) n_obs = obs.sum(axis=0) X = np.where(obs, R, 0.0).astype(np.float64) mean = X.sum(axis=0) / np.maximum(n_obs, 1) X = np.where(obs, X - mean, 0.0) ss = np.square(X).sum(axis=0) ok = (n_obs >= MIN_OBS) & (ss > 0) X, ss, n_obs = X[:, ok], ss[ok], n_obs[ok] sigma = np.sqrt(ss / (n_obs - 1) * DAYS_PER_MONTH) X /= np.sqrt(ss) # unit-norm columns: X'X is a correlation matrix T, n = X.shape _, s, Vt = np.linalg.svd(X, full_matrices=False) k = int(np.sum(np.square(s) > mp_edge(n, T))) B = Vt[:k].T * s[:k] psi = np.maximum(1.0 - np.square(B).sum(axis=1), PSI_FLOOR) return ok, B, psi, sigma def char_factors(Xm: np.ndarray) -> np.ndarray: """Principal components of one month's ranked characteristics above the Marchenko-Pastur edge.""" Z = pd.DataFrame(Xm).rank(pct=True).to_numpy() - 0.5 Z = np.nan_to_num(Z, nan=0.0) Z -= Z.mean(axis=0) sd = Z.std(axis=0) Z = Z[:, sd > 0] / sd[sd > 0] n, p = Z.shape U, s, _ = np.linalg.svd(Z / np.sqrt(n), full_matrices=False) k = int(np.sum(np.square(s) > mp_edge(p, n))) return U[:, :k] * s[:k] def residualize(a: np.ndarray, G: np.ndarray, psi: np.ndarray) -> np.ndarray: """Q a: the GLS residual of a on G with weights 1/psi.""" Gw = G / psi[:, None] return a - G @ np.linalg.solve(G.T @ Gw, Gw.T @ a) def unit_risk_book(s: np.ndarray, G: np.ndarray, psi: np.ndarray, sigma: np.ndarray) -> np.ndarray: """Weights for scores s, scaled to unit predicted volatility.""" v = residualize(np.sqrt(psi) * s, G, psi) / psi return v / np.sqrt(np.sum(np.square(v) * psi)) / sigma def blend_coefficients(past: np.ndarray) -> np.ndarray: """Long-only max-Sharpe weights on the books from their past returns (one row per month).""" k = past.shape[1] if len(past) < BLEND_MIN_MONTHS: return np.full(k, 1.0 / k) age = np.arange(len(past))[::-1] omega = 0.5 ** (age / BLEND_HALFLIFE) omega /= omega.sum() mu = omega @ past dev = past - mu cov = (dev * omega[:, None]).T @ dev / (1.0 - np.square(omega).sum()) vol = np.sqrt(np.diag(cov)) if not np.all(vol > 0): return np.full(k, 1.0 / k) corr = cov / np.outer(vol, vol) x = np.linalg.solve((1.0 - BLEND_SHRINK) * corr + BLEND_SHRINK * np.eye(k), mu / vol) x = np.clip(x, 0.0, None) if x.sum() <= 0: x = np.ones(k) return x / x.sum() / vol # --------------------------------------------------------------------------------------- # Learner # --------------------------------------------------------------------------------------- def book_sharpe_metric(month_codes: np.ndarray, z: np.ndarray, n_months: int): """Validation metric: Sharpe ratio of the monthly returns s'z / ||s||.""" def feval(preds, _data): num = np.bincount(month_codes, weights=preds * z, minlength=n_months) den = np.sqrt(np.bincount(month_codes, weights=preds * preds, minlength=n_months)) R = np.where(den > 0, num / np.where(den > 0, den, 1.0), 0.0) sd = R.std(ddof=1) return "book_sharpe", (R.mean() / sd if sd > 0 else 0.0), True return feval def lgb_params(min_leaf: int) -> dict: return dict(objective="regression", learning_rate=LEARNING_RATE, num_leaves=NUM_LEAVES, min_data_in_leaf=min_leaf, seed=SEED, num_threads=N_THREADS, deterministic=True, force_row_wise=True, verbose=-1, metric="None") def fit(X: np.ndarray, label: np.ndarray, z: np.ndarray, month_ix: np.ndarray, rows: np.ndarray): """Choose the number of rounds on the most recent months, then refit on all rows.""" months = np.unique(month_ix[rows]) n_valid = min(VALID_MONTHS, len(months) // 2) cut = months[-n_valid] tr, va = rows[month_ix[rows] < cut], rows[month_ix[rows] >= cut] min_leaf = int(np.median(np.bincount(month_ix[tr])[months[months < cut]])) codes = month_ix[va] - cut booster = lgb.train( lgb_params(min_leaf), lgb.Dataset(X[tr], label=label[tr]), num_boost_round=MAX_ROUNDS, valid_sets=[lgb.Dataset(X[va], label=label[va])], feval=book_sharpe_metric(codes, z[va], n_valid), callbacks=[lgb.early_stopping(PATIENCE, first_metric_only=True, verbose=False)]) rounds = max(int(booster.best_iteration), 1) min_leaf_all = int(np.median(np.bincount(month_ix[rows])[months])) final = lgb.train(lgb_params(min_leaf_all), lgb.Dataset(X[rows], label=label[rows]), num_boost_round=rounds) return final, rounds, float(booster.best_score["valid_0"]["book_sharpe"]) # --------------------------------------------------------------------------------------- # Model # --------------------------------------------------------------------------------------- def build_weights(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame, return_books: bool = False): """Weights for every month with a fitted model.""" np.random.seed(SEED) t0 = time.perf_counter() feats = [str(f) for f in features["features"].tolist()] panel = chars[["id", "eom", TARGET] + feats].copy() panel["id"] = panel["id"].astype(np.int64) panel["_day"] = to_days(panel["eom"]) panel = panel.sort_values(["_day", "id"], kind="stable").reset_index(drop=True) months, month_ix = np.unique(panel["_day"].to_numpy(), return_inverse=True) starts = np.searchsorted(month_ix, np.arange(len(months) + 1)) ids_all = panel["id"].to_numpy(np.int64) X = panel[feats].to_numpy(np.float32) ret = panel[TARGET].to_numpy(np.float64) del panel log(f"{len(X):,} rows, {len(feats)} characteristics, {len(months)} months, {N_THREADS} threads, " f"books {BOOKS}, scale {SCALE}") uid = np.unique(ids_all) days, M = daily_matrix(daily_ret, uid) col_all = np.searchsorted(uid, ids_all) log(f"daily matrix {M.shape[0]:,} days x {M.shape[1]:,} ids ({time.perf_counter() - t0:.0f}s)") # Risk model and targets for every month. The target of month t uses r_{t+1} and is only used by # refits after t. risk = [None] * len(months) z = {b: np.full(len(X), np.nan) for b in BOOKS} label = {b: np.full(len(X), np.nan) for b in BOOKS} n_ret, n_char = [], [] for t in range(len(months)): lo, hi = starts[t], starts[t + 1] end = int(np.searchsorted(days, months[t], side="right")) if end < MIN_OBS: continue ok, B, psi, sigma = risk_model(M[max(0, end - WINDOW_DAYS):end, col_all[lo:hi]]) if ok.sum() < MIN_NAMES: continue rows = lo + np.flatnonzero(ok) G = {"open": np.column_stack([B, 1.0 / sigma])} if "neutral" in BOOKS: C = char_factors(X[rows]) G["neutral"] = np.column_stack([G["open"], C]) n_char.append(C.shape[1]) n_ret.append(B.shape[1]) risk[t] = (rows, psi, sigma, G) x = ret[rows] / sigma for b in BOOKS: zt = residualize(x, G[b], psi) / np.sqrt(psi) z[b][rows] = zt label[b][rows] = np.clip(zt / zt.std(), -LABEL_CLIP, LABEL_CLIP) log(f"risk pass: {sum(r is not None for r in risk)} months, return factors median {np.median(n_ret):.0f} " f"(range {min(n_ret)}-{max(n_ret)})" + (f", characteristic factors median {np.median(n_char):.0f} (range {min(n_char)}-{max(n_char)})" if n_char else "") + f" ({time.perf_counter() - t0:.0f}s)") # Refits, books and blend. frames, book_frames, past_returns = [], [], [] models, last_refit = {}, None labelled = [t for t in range(len(months)) if risk[t] is not None] for t in range(len(months)): past = [m for m in labelled if m <= t - 1] if len(past) >= MIN_LABELLED_MONTHS and (last_refit is None or t - last_refit >= REFIT_EVERY): tf = time.perf_counter() rows = np.concatenate([risk[m][0] for m in past]) notes = [] for b in BOOKS: models[b], rounds, val_sr = fit(X, label[b], z[b], month_ix, rows) notes.append(f"{b} rounds {rounds} valid SR {val_sr:.2f}") last_refit = t log(f"refit {pd.Timestamp(months[t], unit='D').date()}: {len(rows):,} rows; " + "; ".join(notes) + f"; {time.perf_counter() - tf:.0f}s (elapsed {(time.perf_counter() - t0) / 60:.1f} min)") if not models or risk[t] is None: continue rows, psi, sigma, G = risk[t] W = np.vstack([unit_risk_book(models[b].predict(X[rows]), G[b], psi, sigma) for b in BOOKS]) lam = blend_coefficients(np.array(past_returns) if past_returns else np.empty((0, len(BOOKS)))) w = lam @ W w = w / np.sqrt(np.sum(np.square(w * sigma) * psi)) # unit predicted volatility if SCALE == "gross": w = w / np.abs(w).sum() past_returns.append(W @ ret[rows]) # known at t+1 frames.append(pd.DataFrame({"id": ids_all[rows], "eom": months[t], "w": w})) if return_books: book_frames.append(pd.DataFrame(W.T, columns=list(BOOKS)).assign(id=ids_all[rows], eom=months[t])) if not frames: raise ValueError(f"no tradable month: {len(months)} months, need > {MIN_LABELLED_MONTHS} labelled months " f"with >= {MIN_NAMES} names having {MIN_OBS} daily returns") out = pd.concat(frames, ignore_index=True) out["id"] = out["id"].astype(np.int64) out["eom"] = pd.to_datetime(out["eom"], unit="D").dt.strftime("%Y-%m-%d") assert np.isfinite(out["w"]).all(), "non-finite weights" log(f"done: {len(out):,} rows in {(time.perf_counter() - t0) / 60:.1f} min") if return_books: books = pd.concat(book_frames, ignore_index=True) books["id"] = books["id"].astype(np.int64) books["eom"] = pd.to_datetime(books["eom"], unit="D").dt.strftime("%Y-%m-%d") return out[["id", "eom", "w"]], books return out[["id", "eom", "w"]] def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame: np.random.seed(SEED) weights = build_weights(chars, features, daily_ret) test = chars["ctff_test"].astype(str).str.strip().str.lower().isin(["1", "1.0", "true"]) if not test.any(): # validation panel return weights keys = pd.DataFrame({ "id": chars.loc[test, "id"].astype(np.int64), "eom": pd.to_datetime(chars.loc[test, "eom"]).dt.strftime("%Y-%m-%d"), }) return weights.merge(keys, on=["id", "eom"], how="inner")