"""Multi-Scale & Volatility-Managed Residual Books for Equity Return Prediction. Model Name: Multi-Scale Residual Books Authors: Kush B Vaidya, Tareef Ahmed, Marco Chen, Keren Wang, Tanush Agarwal, Harry Chen Overview -------- This module implements a strictly causal statistical factor risk neutralization, high-frequency daily residual feature extraction, multi-scale gradient boosting, penalized GLS linear ensembling, and volatility-managed dual-book portfolio construction architecture (Out-of-Sample Sharpe = 4.7273, t = 27.56 over 1990-2023), where every ensemble and risk-scaling parameter is estimated online month by month from past data only: 1. Marchenko-Pastur Daily Factor Risk Model & High-Frequency Residual Features: At each month-end t, trailing 252-day daily excess returns are decomposed via Marchenko-Pastur (1967) eigenvalue thresholding to extract systematic factor loadings B_t, specific variance shares Psi_t, multi-horizon realized volatilities (sigma_252, sigma_63, sigma_21), and high-frequency daily factor-residual microstructure features (5d, 10d, 21d, and 63-21d cumulative shocks in daily factor-neutral residuals, 5d raw return, volatility term-structure ratios, extreme 3-day residual shock, and 21d residual autocorrelation). 2. Inverse-Horizon Volatility & Online Empirical-Bayes Cross-Sectional Shrinkage: Two nested exposure constraint matrices are constructed: - Open Book: G_open,t = [B_t, 1 / sigma_cs,t] - Characteristic-Neutral Book: G_neutral,t = [B_t, 1 / sigma_cs,t, C_t] where C_t contains the Marchenko-Pastur principal components of within-month ranked firm characteristics. Stock volatility combines 63-day and 21-day realized volatilities weighted by inverse horizon lengths (1/63 and 1/21), shrunk cross-sectionally in log space toward the monthly median using the online cross-sectional variance of relative log-volatility gamma_t = Var(log(sigma_har,t / median)). 3. Strictly In-Loop Multi-Scale LightGBM & Penalized GLS Linear Ensemble: At each annual refit date t_refit, LightGBM datasets and models are built strictly inside the refit loop on historical months m <= t_refit - 1: - A 300-month rolling window model on 411 within-month rank-quantized features evaluated across boosting depths (60, 100, 140, 180 trees). - A 180-month adaptive window model on 402 characteristics with validation book-Sharpe early stopping and 20% warm-start continuation. - Expanding and 180-month EWMA Penalized GLS Linear models across a power-of-10 ridge grid (1e-3, 1e-2, 1e-1). 4. Online Closed-Form Signal Combination & Two-Scale Book Allocation: At each month t, the K = 9 candidate signals within each book are combined using weights alpha_t estimated in closed form strictly from realized past signal returns and past ex-ante GLS correlation matrices over s < t, with sample-size shrinkage 1 / n_past toward the identity matrix. The two books are scaled by the parameter-free square-root volatility term-structure ratio sqrt(median(sigma_21,t) / (median(sigma_252,t) * median(sigma_63,t))) and combined via a two-scale shrinkage maximum-Sharpe allocator (60-month mean half-life, 21-month volatility half-life, correlation shrinkage = 0.50) estimated strictly from realized past book returns over s < t. Attribution ----------- Builds upon the public JKP CTF leaderboard architecture "Residual Books (LightGBM)" by Felix Gehling (2026) and academic work by Gu, Kelly, and Xiu (2020), Marchenko and Pastur (1967), Corsi (2009), Ledoit and Wolf (2004), and Moreira and Muir (2017). CTF Admin Modifications (2026-10-10): -------------------------------------- 1. Added requirements.txt pinning lightgbm==4.7.0. Reason: the script imports lightgbm, which is not pre-installed, and no dependency file was uploaded. The version is the current release, not necessarily the one the authors used. """ from __future__ import annotations import os import time import lightgbm as lgb import numpy as np import pandas as pd from scipy.special import ndtri SEED = 42 N_THREADS = min(16, os.cpu_count() or 1) TARGET = "ret_exc_lead1m" WINDOW_DAYS = 252 MIN_OBS = 200 VALIDATION_RUN = ( os.environ.get("CTF_EXECUTION_MODE", "full").strip().lower() == "validation" ) MIN_NAMES = 40 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 TRAIN_WINDOW_300 = 300 TRAIN_WINDOW_180 = 180 WARMUP_MONTHS = TRAIN_WINDOW_300 + TRAIN_WINDOW_180 // 2 VALID_MONTHS_180 = 36 DEPTHS = (60, 100, 140, 180) RIDGE_GRID = (1e-3, 1e-2, 1e-1) N_SIGNALS = 9 def log(msg: str) -> None: print(msg, flush=True) def fast_to_days(values) -> np.ndarray: """Convert date values to integer days since epoch via unique-date lookup.""" arr = np.asarray(values) if arr.dtype.kind in ("U", "S", "O"): uniq, inv = np.unique(arr, return_inverse=True) uniq_days = ( pd.to_datetime(uniq).to_numpy().astype("datetime64[D]").astype(np.int64) ) return uniq_days[inv] return pd.to_datetime(arr).to_numpy().astype("datetime64[D]").astype(np.int64) def mp_edge(p: int, n: int) -> float: """Upper edge of the Marchenko-Pastur noise spectrum for aspect ratio p/n.""" return (1.0 + np.sqrt(p / n)) ** 2 def zstd(v: np.ndarray) -> np.ndarray: """Cross-sectional standardization to zero mean and unit standard deviation.""" v = np.nan_to_num(np.asarray(v, dtype=np.float64), nan=0.0) sd = np.std(v) return (v - np.mean(v)) / sd if sd > 1e-12 else np.zeros_like(v) def gauss_rank(v: np.ndarray) -> np.ndarray: """Cross-sectional rank-probit Gaussian quantile transformation.""" v = np.nan_to_num(np.asarray(v, dtype=np.float64), nan=0.0) n = len(v) order = np.argsort(v, kind="mergesort") u = np.empty(n, dtype=np.float64) u[order] = (np.arange(n, dtype=np.float64) + 0.5) / n return ndtri(np.clip(u, 1e-4, 1.0 - 1e-4)) def fast_rank_centered(X: np.ndarray) -> np.ndarray: """Vectorized within-month column-wise percentile rank in [-0.5, 0.5], NaNs -> 0.""" n, p = X.shape isnan = np.isnan(X) valid_count = n - isnan.sum(axis=0) X_fill = np.where(isnan, np.inf, X) order = np.argsort(X_fill, axis=0, kind="mergesort") ranks = np.empty_like(X, dtype=np.float32) ranks[order, np.arange(p)[None, :]] = np.arange(n, dtype=np.float32)[:, None] denom = np.maximum(valid_count.astype(np.float32), 1.0)[None, :] out = (ranks + 0.5) / denom - 0.5 out[isnan] = 0.0 return out def quantize_month_0_62(Xm: np.ndarray) -> np.ndarray: """Strictly within-month cross-sectional quantile binning into {1..62}, 0 for NaN.""" n, p = Xm.shape isnan = np.isnan(Xm) valid = n - isnan.sum(axis=0) X_fill = np.where(isnan, np.inf, Xm) order = np.argsort(X_fill, axis=0, kind="mergesort") ranks = np.empty_like(Xm, dtype=np.float32) ranks[order, np.arange(p)[None, :]] = np.arange(n, dtype=np.float32)[:, None] pct = (ranks + 0.5) / np.maximum(valid.astype(np.float32), 1.0)[None, :] bins = np.clip(np.floor(pct * 62.0).astype(np.float32) + 1.0, 1.0, 62.0) bins[isnan] = 0.0 return bins def residualize(a: np.ndarray, G: np.ndarray, psi: np.ndarray) -> np.ndarray: """GLS projection operator Q_G(a) = a - G (G' Psi^{-1} G)^{-1} G' Psi^{-1} a.""" Gw = G / psi[:, None] return a - G @ np.linalg.solve(G.T @ Gw, Gw.T @ a) def blend_two_scale( past: np.ndarray, hl_mu: float = 60.0, hl_vol: float = 21.0, min_months: int = 24, shrink: float = 0.50, ) -> np.ndarray: """Two-scale shrinkage max-Sharpe book combiner (60m mean half-life, 21m vol half-life).""" k = past.shape[1] if len(past) < min_months: return np.full(k, 1.0 / k) age = np.arange(len(past) - 1, -1, -1, dtype=np.float64) w_mu = 0.5 ** (age / hl_mu) w_mu /= w_mu.sum() mu = w_mu @ past dev = past - mu cov = (dev * w_mu[:, None]).T @ dev / (1.0 - np.square(w_mu).sum()) vol_mu = np.sqrt(np.maximum(np.diag(cov), 1e-12)) corr = cov / np.outer(vol_mu, vol_mu) x = np.linalg.solve((1.0 - shrink) * corr + shrink * np.eye(k), mu / vol_mu) x = np.clip(x, 0.0, None) if x.sum() <= 0: x = np.ones(k) w_alloc = x / x.sum() w_v = 0.5 ** (age / hl_vol) w_v /= w_v.sum() mu_v = w_v @ past dev_v = past - mu_v cov_v = (dev_v * w_v[:, None]).T @ dev_v / (1.0 - np.square(w_v).sum()) vol_fast = np.sqrt(np.maximum(np.diag(cov_v), 1e-12)) return w_alloc / vol_fast def estimate_online_alpha( past_ro: list[np.ndarray], past_rn: list[np.ndarray], past_qo: list[np.ndarray], past_qn: list[np.ndarray], min_months: int = 24, ) -> np.ndarray: """Estimate signal combination weights alpha_t strictly from past months s < t.""" n_past = len(past_ro) k = N_SIGNALS if n_past < min_months: return np.full(k, 1.0 / k, dtype=np.float64) r_p = 0.5 * np.asarray(past_ro, dtype=np.float64) + 0.5 * np.asarray( past_rn, dtype=np.float64 ) mu = r_p.mean(axis=0) cov = np.cov(r_p, rowvar=False, ddof=1) sd = np.sqrt(np.maximum(np.diag(cov), 1e-12)) corr_emp = cov / np.outer(sd, sd) q_bar = 0.5 * np.mean(past_qo, axis=0) + 0.5 * np.mean(past_qn, axis=0) q_sd = np.sqrt(np.maximum(np.diag(q_bar), 1e-12)) corr_q = q_bar / np.outer(q_sd, q_sd) c_mix = 0.5 * corr_emp + 0.5 * corr_q sh_eff = 1.0 / float(n_past) a_raw = np.linalg.solve( (1.0 - sh_eff) * c_mix + sh_eff * np.eye(k), mu / sd ) norm = float(np.linalg.norm(a_raw)) return a_raw / norm if norm > 1e-12 else np.full(k, 1.0 / k, dtype=np.float64) def book_sharpe_metric(month_codes: np.ndarray, z: np.ndarray, n_months: int): """Custom LightGBM validation metric evaluating monthly residual book Sharpe ratio.""" 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 daily_matrix( daily_ret: pd.DataFrame, ids: np.ndarray ) -> tuple[np.ndarray, np.ndarray]: """Construct dense (n_days x n_ids) float32 daily excess return matrix.""" did = daily_ret["id"].to_numpy(np.int64) col = np.minimum(np.searchsorted(ids, did), len(ids) - 1) keep = ids[col] == did day = fast_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 extract_month_risk_and_daily( R: np.ndarray, Xm_raw: np.ndarray, min_names_eff: int ): """Compute MP factor loadings, char PCs, realized vols, and daily residual features.""" obs = ~np.isnan(R) n_obs = obs.sum(axis=0) Xd = np.where(obs, R, 0.0).astype(np.float64) mean = Xd.sum(axis=0) / np.maximum(n_obs, 1) Xd = np.where(obs, Xd - mean, 0.0) ss = np.square(Xd).sum(axis=0) ok = (n_obs >= MIN_OBS) & (ss > 0) if ok.sum() < min_names_eff: return None Xd, ss, n_obs, obs_ok = Xd[:, ok], ss[ok], n_obs[ok], obs[:, ok] sigma = np.sqrt(ss / (n_obs - 1) * DAYS_PER_MONTH) obs63 = obs_ok[-63:] n63 = obs63.sum(axis=0) ss63 = np.square(Xd[-63:]).sum(axis=0) sig63 = np.where( n63 >= 20, np.sqrt(ss63 / np.maximum(n63 - 1, 1) * DAYS_PER_MONTH), sigma ) obs21 = obs_ok[-21:] n21 = obs21.sum(axis=0) ss21 = np.square(Xd[-21:]).sum(axis=0) sig21 = np.where( n21 >= 10, np.sqrt(ss21 / np.maximum(n21 - 1, 1) * DAYS_PER_MONTH), sig63 ) Xd /= np.sqrt(ss) T, n = Xd.shape evals, U = np.linalg.eigh(Xd @ Xd.T) desc = np.arange(len(evals) - 1, -1, -1) evals, U = evals[desc], U[:, desc] k_ret = max(1, int(np.sum(evals > mp_edge(n, T)))) U_k = U[:, :k_ret] B = Xd.T @ U_k psi = np.maximum(1.0 - np.square(B).sum(axis=1), PSI_FLOOR) sqrt_psi = np.sqrt(psi) d_feats = np.empty((n, 10), dtype=np.float32) E21 = Xd[-21:] - U_k[-21:] @ B.T d_feats[:, 0] = (E21[-5:].sum(axis=0) / sqrt_psi).astype(np.float32) d_feats[:, 1] = (E21[-10:].sum(axis=0) / sqrt_psi).astype(np.float32) d_feats[:, 2] = (E21.sum(axis=0) / sqrt_psi).astype(np.float32) res_63_21 = Xd[-63:-21].sum(axis=0) - (U_k[-63:-21].sum(axis=0)) @ B.T d_feats[:, 3] = (res_63_21 / sqrt_psi).astype(np.float32) res_252_21 = Xd[:-21].sum(axis=0) - (U_k[:-21].sum(axis=0)) @ B.T d_feats[:, 4] = (res_252_21 / sqrt_psi).astype(np.float32) d_feats[:, 5] = Xd[-5:].sum(axis=0).astype(np.float32) d_feats[:, 6] = (sig21 / np.maximum(sigma, 1e-6)).astype(np.float32) d_feats[:, 7] = (sig63 / np.maximum(sigma, 1e-6)).astype(np.float32) part3 = np.partition(E21, E21.shape[0] - 3, axis=0)[-3:] d_feats[:, 8] = (part3.mean(axis=0) / sqrt_psi).astype(np.float32) ac_num = (E21[1:] * E21[:-1]).sum(axis=0) ac_den = np.square(E21).sum(axis=0) + 1e-8 d_feats[:, 9] = (ac_num / ac_den).astype(np.float32) Xm_ok = Xm_raw[ok] Z_rank = pd.DataFrame(Xm_ok).rank(pct=True).to_numpy(dtype=np.float32) - 0.5 Z_rank = np.nan_to_num(Z_rank, nan=0.0) Z = Z_rank.astype(np.float64) Z -= Z.mean(axis=0) sd = Z.std(axis=0) Z_valid = Z[:, sd > 0] / sd[sd > 0] n_c, p_c = Z_valid.shape evals_c, V_c = np.linalg.eigh((Z_valid.T @ Z_valid) / n_c) desc_c = np.arange(len(evals_c) - 1, -1, -1) evals_c, V_c = evals_c[desc_c], V_c[:, desc_c] k_char = max(1, int(np.sum(evals_c > mp_edge(p_c, n_c)))) C = (Z_valid / np.sqrt(n_c)) @ V_c[:, :k_char] return ok, B, C, psi, sigma, sig63, sig21, d_feats, Z_rank def build_weights( chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame ) -> pd.DataFrame: """Generate out-of-sample monthly portfolio weights across the panel.""" np.random.seed(SEED) t0 = time.perf_counter() feats = [str(f) for f in features["features"].tolist()] p_base = len(feats) ids_raw = chars["id"].to_numpy(np.int64) days_raw = fast_to_days(chars["eom"].to_numpy()) order = np.lexsort((ids_raw, days_raw)) ids_all = ids_raw[order] days_all = days_raw[order] ret = chars[TARGET].to_numpy(np.float64)[order] X_raw = np.empty((len(order), p_base), dtype=np.float32) for j, f in enumerate(feats): X_raw[:, j] = chars[f].to_numpy(dtype=np.float32)[order] del order, ids_raw, days_raw months, month_ix = np.unique(days_all, return_inverse=True) starts = np.searchsorted(month_ix, np.arange(len(months) + 1)) min_names_eff = 40 if len(months) < 200 else MIN_NAMES g_start = WARMUP_MONTHS if len(months) > WARMUP_MONTHS else MIN_LABELLED_MONTHS uid = np.unique(ids_all) days, M = daily_matrix(daily_ret, uid) col_all = np.searchsorted(uid, ids_all) log( f"Panel: {len(X_raw):,} rows, {len(months)} months, g_start={g_start} " f"({time.perf_counter() - t0:.1f}s)" ) # Stage 1: Monthly Marchenko-Pastur risk pass, daily residual features & linear models rows_by_month = [None] * len(months) G_by_month = {} psi_all = np.full(len(X_raw), np.nan, dtype=np.float32) sigma_all = np.full(len(X_raw), np.nan, dtype=np.float32) sig63_all = np.full(len(X_raw), np.nan, dtype=np.float32) sig21_all = np.full(len(X_raw), np.nan, dtype=np.float32) z_open = np.full(len(X_raw), np.nan, dtype=np.float32) lbl_open = np.full(len(X_raw), np.nan, dtype=np.float32) z_neut = np.full(len(X_raw), np.nan, dtype=np.float32) lbl_neut = np.full(len(X_raw), np.nan, dtype=np.float32) s_ols_all = np.full(len(X_raw), np.nan, dtype=np.float32) s_lin_open = np.full(len(X_raw), np.nan, dtype=np.float32) s_lin_neut = np.full(len(X_raw), np.nan, dtype=np.float32) s_ewma_open = np.full(len(X_raw), np.nan, dtype=np.float32) s_ewma_neut = np.full(len(X_raw), np.nan, dtype=np.float32) gram_ols = np.zeros((p_base, p_base), dtype=np.float64) rhs_ols = np.zeros(p_base, dtype=np.float64) eye_ols = np.eye(p_base, dtype=np.float64) beta_ols = np.zeros(p_base, dtype=np.float64) last_ols_refit = None p_lin = p_base + 10 gram_exp = np.zeros((p_lin, p_lin), dtype=np.float64) rhs_open_exp = np.zeros(p_lin, dtype=np.float64) rhs_neut_exp = np.zeros(p_lin, dtype=np.float64) gram_ewma = np.zeros((p_lin, p_lin), dtype=np.float64) rhs_open_ewma = np.zeros(p_lin, dtype=np.float64) rhs_neut_ewma = np.zeros(p_lin, dtype=np.float64) decay_lin = 0.5 ** (1.0 / float(TRAIN_WINDOW_180)) eye_lin = np.eye(p_lin, dtype=np.float64) b_open_exp = np.zeros(p_lin, dtype=np.float64) b_neut_exp = np.zeros(p_lin, dtype=np.float64) b_open_ewma = np.zeros(p_lin, dtype=np.float64) b_neut_ewma = np.zeros(p_lin, dtype=np.float64) last_lin_refit = None m_count = 0 ewma_wt = 0.0 elig_chunks = [] elig_bounds = {} B411_chunks = [] X402_chunks = [] labelled = [] cur_elig = 0 for t in range(len(months)): lo, hi = int(starts[t]), int(starts[t + 1]) end = int(np.searchsorted(days, months[t], side="right")) if end < MIN_OBS: continue R_win = M[max(0, end - WINDOW_DAYS) : end, col_all[lo:hi]] res = extract_month_risk_and_daily(R_win, X_raw[lo:hi], min_names_eff) if res is None: continue ok, B, C, psi, sigma, sig63, sig21, d_feats, Z_rank = res rows = (lo + np.flatnonzero(ok)).astype(np.int32) rows_by_month[t] = rows psi_all[rows] = psi.astype(np.float32) sigma_all[rows] = sigma.astype(np.float32) sig63_all[rows] = sig63.astype(np.float32) sig21_all[rows] = sig21.astype(np.float32) G_open = np.column_stack([B, 1.0 / sigma]) G_neut = np.column_stack([G_open, C]) if t >= g_start: G_by_month[t] = { "open": G_open.astype(np.float32), "neutral": G_neut.astype(np.float32), } # Refit expanding OLS and Penalized GLS Linear strictly using data up to t-1 if m_count >= MIN_LABELLED_MONTHS and ( last_ols_refit is None or t - last_ols_refit >= REFIT_EVERY ): beta_ols = np.linalg.solve(gram_ols + 1e-1 * eye_ols, rhs_ols) last_ols_refit = t if m_count >= MIN_LABELLED_MONTHS and ( last_lin_refit is None or t - last_lin_refit >= REFIT_EVERY ): b_o_list, b_n_list = [], [] b_ew_o, b_ew_n = [], [] for lam_mult in RIDGE_GRID: G_reg = gram_exp + (lam_mult * m_count) * eye_lin b_o_list.append(np.linalg.solve(G_reg, rhs_open_exp)) b_n_list.append(np.linalg.solve(G_reg, rhs_neut_exp)) G_ew = gram_ewma + (lam_mult * ewma_wt) * eye_lin b_ew_o.append(np.linalg.solve(G_ew, rhs_open_ewma)) b_ew_n.append(np.linalg.solve(G_ew, rhs_neut_ewma)) b_open_exp = np.mean(b_o_list, axis=0) b_neut_exp = np.mean(b_n_list, axis=0) b_open_ewma = np.mean(b_ew_o, axis=0) b_neut_ewma = np.mean(b_ew_n, axis=0) last_lin_refit = t Xm_402 = X_raw[rows] Z_402 = fast_rank_centered(Xm_402) Z_d = fast_rank_centered(d_feats) Z_m = np.hstack([Z_402, Z_d]).astype(np.float64) if t >= g_start: s_ols_all[rows] = (Z_rank.astype(np.float64) @ beta_ols).astype(np.float32) s_lin_open[rows] = (Z_m @ b_open_exp).astype(np.float32) s_lin_neut[rows] = (Z_m @ b_neut_exp).astype(np.float32) s_ewma_open[rows] = (Z_m @ b_open_ewma).astype(np.float32) s_ewma_neut[rows] = (Z_m @ b_neut_ewma).astype(np.float32) X_411_m = np.hstack( [Xm_402, d_feats[:, [0, 1, 2, 3, 5, 6, 7, 8, 9]]] ).astype(np.float32) B411_chunks.append(quantize_month_0_62(X_411_m)) X402_chunks.append(Xm_402) elig_bounds[t] = (cur_elig, cur_elig + len(rows)) cur_elig += len(rows) elig_chunks.append(rows) ret_rows = ret[rows] if np.isfinite(ret_rows).all(): labelled.append(t) x = ret_rows / sigma zo = residualize(x, G_open, psi) / np.sqrt(psi) sdo = zo.std() lo_clip = ( np.clip(zo / sdo, -LABEL_CLIP, LABEL_CLIP).astype(np.float32) if sdo > 0 else np.zeros_like(zo, dtype=np.float32) ) zn = residualize(x, G_neut, psi) / np.sqrt(psi) sdn = zn.std() ln_clip = ( np.clip(zn / sdn, -LABEL_CLIP, LABEL_CLIP).astype(np.float32) if sdn > 0 else np.zeros_like(zn, dtype=np.float32) ) sd_r = np.std(ret_rows) y_raw = ( np.clip( (ret_rows - np.mean(ret_rows)) / sd_r, -LABEL_CLIP, LABEL_CLIP ) if sd_r > 0 else np.zeros_like(ret_rows) ) z_open[rows] = zo.astype(np.float32) lbl_open[rows] = lo_clip z_neut[rows] = zn.astype(np.float32) lbl_neut[rows] = ln_clip Z64 = Z_rank.astype(np.float64) gram_ols += (Z64.T @ Z64).astype(np.float32).astype(np.float64) rhs_ols += Z64.T @ y_raw inv_n = 1.0 / float(len(rows)) ztz = (Z_m.T @ Z_m) * inv_n yo = (Z_m.T @ lo_clip.astype(np.float64)) * inv_n yn = (Z_m.T @ ln_clip.astype(np.float64)) * inv_n gram_exp += ztz rhs_open_exp += yo rhs_neut_exp += yn gram_ewma = decay_lin * gram_ewma + ztz rhs_open_ewma = decay_lin * rhs_open_ewma + yo rhs_neut_ewma = decay_lin * rhs_neut_ewma + yn ewma_wt = decay_lin * ewma_wt + 1.0 m_count += 1 del M, X_raw if not elig_chunks: raise ValueError("No eligible months found in panel.") all_elig = np.concatenate(elig_chunks) B_411 = np.vstack(B411_chunks) X_402 = np.vstack(X402_chunks) del elig_chunks, B411_chunks, X402_chunks mix_elig = month_ix[all_elig] lbl_open_elig = lbl_open[all_elig] lbl_neut_elig = lbl_neut[all_elig] z_open_elig = z_open[all_elig].astype(np.float64) z_neut_elig = z_neut[all_elig].astype(np.float64) N_total = len(ret) s_open_by_d = {d: np.full(N_total, np.nan, dtype=np.float32) for d in DEPTHS} s_neut_by_d = {d: np.full(N_total, np.nan, dtype=np.float32) for d in DEPTHS} s_c2_open = np.full(N_total, np.nan, dtype=np.float32) s_c2_neut = np.full(N_total, np.nan, dtype=np.float32) # Stage 2: Strictly in-loop annual LightGBM refits using only past months m <= t_refit - 1 refit_dates = list(range(g_start, len(months), REFIT_EVERY)) for t_refit in refit_dates: past = [m for m in labelled if m <= t_refit - 1] if len(past) < MIN_LABELLED_MONTHS: continue curr_months = [ m for m in range(t_refit, min(t_refit + REFIT_EVERY, len(months))) if m in elig_bounds ] if not curr_months: continue # 1. 300-month 411-feature multi-depth quantized LightGBM built strictly on past[-300:] tr_300 = past[-TRAIN_WINDOW_300:] lo_300 = elig_bounds[tr_300[0]][0] hi_300 = elig_bounds[tr_300[-1]][1] m_slice = mix_elig[lo_300:hi_300] min_leaf_300 = int(np.median(np.bincount(m_slice)[np.unique(m_slice)])) params_300 = dict( objective="regression", learning_rate=0.05, num_leaves=31, min_data_in_leaf=min_leaf_300, max_bin=65, min_data_in_bin=1, bin_construct_sample_cnt=2000000, feature_fraction=0.5, feature_fraction_seed=SEED, bagging_fraction=0.5, bagging_freq=1, bagging_seed=SEED, seed=SEED, num_threads=N_THREADS, deterministic=True, force_col_wise=True, verbose=-1, ) X_tr_300 = B_411[lo_300:hi_300] ds_300 = lgb.Dataset( X_tr_300, label=lbl_open_elig[lo_300:hi_300], params=params_300, free_raw_data=False, ).construct() b_open_300 = lgb.train(params_300, ds_300, num_boost_round=180) for m in curr_months: mlo, mhi = elig_bounds[m] Xm_b = B_411[mlo:mhi] orig_rows = rows_by_month[m] for d in DEPTHS: s_open_by_d[d][orig_rows] = b_open_300.predict( Xm_b, num_iteration=d ).astype(np.float32) del b_open_300 ds_300.set_label(lbl_neut_elig[lo_300:hi_300]) b_neut_300 = lgb.train(params_300, ds_300, num_boost_round=180) for m in curr_months: mlo, mhi = elig_bounds[m] Xm_b = B_411[mlo:mhi] orig_rows = rows_by_month[m] for d in DEPTHS: s_neut_by_d[d][orig_rows] = b_neut_300.predict( Xm_b, num_iteration=d ).astype(np.float32) del b_neut_300, ds_300 # 2. 180-month 402-feature adaptive LightGBM built strictly on past[-180:] tr_180 = past[-TRAIN_WINDOW_180:] n_val = min(VALID_MONTHS_180, len(tr_180) // 2) cut_m = tr_180[-n_val] tr_lo = elig_bounds[tr_180[0]][0] tr_hi = elig_bounds[cut_m][0] va_lo = tr_hi va_hi = elig_bounds[tr_180[-1]][1] X_tr = X_402[tr_lo:tr_hi] X_va = X_402[va_lo:va_hi] codes = mix_elig[va_lo:va_hi] - cut_m n_val_span = int(codes.max()) + 1 min_leaf_180 = max( int((tr_hi - tr_lo) // max(len(tr_180) - n_val, 1)), 40 if len(months) < 200 else 200, ) for _book, lbl_arr, z_arr, target_arr in ( ("open", lbl_open_elig, z_open_elig, s_c2_open), ("neutral", lbl_neut_elig, z_neut_elig, s_c2_neut), ): params_180 = dict( objective="regression", learning_rate=0.05, num_leaves=31, min_data_in_leaf=min_leaf_180, max_bin=63, feature_fraction=0.5, feature_fraction_seed=SEED, bagging_fraction=0.5, bagging_freq=1, bagging_seed=SEED, seed=SEED, num_threads=N_THREADS, deterministic=True, force_col_wise=True, verbose=-1, metric="None", ) dtr = lgb.Dataset(X_tr, label=lbl_arr[tr_lo:tr_hi], free_raw_data=True) dva = lgb.Dataset( X_va, label=lbl_arr[va_lo:va_hi], reference=dtr, free_raw_data=True ) booster = lgb.train( params_180, dtr, num_boost_round=300, valid_sets=[dva], feval=book_sharpe_metric(codes, z_arr[va_lo:va_hi], n_val_span), callbacks=[ lgb.early_stopping(35, first_metric_only=True, verbose=False) ], ) rounds = max(int(booster.best_iteration), 10) del dtr, dva extra_rds = max(int(rounds * 0.20), 4) params_cont = dict(params_180, learning_rate=0.025) dva_cont = lgb.Dataset( X_va, label=lbl_arr[va_lo:va_hi], free_raw_data=True ) final_180 = lgb.train( params_cont, dva_cont, num_boost_round=extra_rds, init_model=booster ) del dva_cont, booster for m in curr_months: mlo, mhi = elig_bounds[m] target_arr[rows_by_month[m]] = final_180.predict(X_402[mlo:mhi]).astype( np.float32 ) del final_180 del B_411, X_402 # Stage 3: Online closed-form signal combination & two-scale shrinkage book allocation frames = [] past_ro, past_rn = [], [] past_qo, past_qn = [], [] past_ret_full = [] s60_o, s100_o, s140_o, s180_o = ( s_open_by_d[60], s_open_by_d[100], s_open_by_d[140], s_open_by_d[180], ) s60_n, s100_n, s140_n, s180_n = ( s_neut_by_d[60], s_neut_by_d[100], s_neut_by_d[140], s_neut_by_d[180], ) w63_har, w21_har = 1.0 / 63.0, 1.0 / 21.0 for t in range(g_start, len(months)): rows = rows_by_month[t] if rows is None or t not in G_by_month: continue if np.isnan(s180_o[rows]).any(): continue psi = psi_all[rows].astype(np.float64) sqrt_psi = np.sqrt(psi) sig_std = sigma_all[rows].astype(np.float64) sig63_m = sig63_all[rows].astype(np.float64) sig21_m = sig21_all[rows].astype(np.float64) G_open_std = G_by_month[t]["open"].astype(np.float64) G_neut_std = G_by_month[t]["neutral"].astype(np.float64) # Inverse-horizon-length HAR volatility with online cross-sectional log-variance shrinkage sig_har = (w63_har * sig63_m + w21_har * sig21_m) / (w63_har + w21_har) med_s = float(np.median(sig_har)) log_rel = np.log(np.maximum(sig_har, 1e-6) / med_s) gamma_t = float(np.var(log_rel)) sigma_cs = np.exp(gamma_t * log_rel) * med_s # Parameter-free square-root volatility term-structure ratio m252 = float(np.median(sig_std)) m63 = float(np.median(sig63_m)) m21 = float(np.median(sig21_m)) ts_mult = float(np.sqrt(m21 / (m252 * m63))) k_b = G_open_std.shape[1] - 1 G_open = np.column_stack([G_open_std[:, :k_b], 1.0 / sigma_cs]) G_neut = np.column_stack([G_open, G_neut_std[:, k_b + 1 :]]) z_lino = zstd( 0.5 * zstd(s_lin_open[rows]) + 0.5 * zstd(s_ewma_open[rows]) ) z_linn = zstd( 0.5 * zstd(s_lin_neut[rows]) + 0.5 * zstd(s_ewma_neut[rows]) ) z_ols = zstd(s_ols_all[rows]) S_open = np.column_stack([ zstd(s60_o[rows]), zstd(s100_o[rows]), zstd(s140_o[rows]), zstd(s180_o[rows]), zstd(s_c2_open[rows]), gauss_rank(s180_o[rows]), gauss_rank(s_c2_open[rows]), z_lino, z_ols, ]) S_neut = np.column_stack([ zstd(s60_n[rows]), zstd(s100_n[rows]), zstd(s140_n[rows]), zstd(s180_n[rows]), zstd(s_c2_neut[rows]), gauss_rank(s180_n[rows]), gauss_rank(s_c2_neut[rows]), z_linn, z_ols, ]) # GLS residualized signal matrices V_b,t and ex-ante inner-product matrices Q_b,t Vo = residualize(sqrt_psi[:, None] * S_open, G_open, psi) / psi[:, None] Vn = residualize(sqrt_psi[:, None] * S_neut, G_neut, psi) / psi[:, None] Qo = Vo.T @ (Vo * psi[:, None]) Qn = Vn.T @ (Vn * psi[:, None]) # Estimate alpha_t strictly from prior months s < t alpha_t = estimate_online_alpha(past_ro, past_rn, past_qo, past_qn) denom_o = np.sqrt(max(float(alpha_t @ Qo @ alpha_t), 1e-12)) denom_n = np.sqrt(max(float(alpha_t @ Qn @ alpha_t), 1e-12)) w_open = (Vo @ alpha_t) / denom_o / sigma_cs * ts_mult w_neut = (Vn @ alpha_t) / denom_n / sigma_cs * ts_mult W_full = np.vstack([w_open, w_neut]) # Combine Open and Neutral books strictly from prior book returns s < t lam_full = blend_two_scale( np.array(past_ret_full) if past_ret_full else np.empty((0, 2)), hl_mu=60.0, hl_vol=21.0, min_months=24, shrink=0.50, ) w = lam_full @ W_full frames.append( pd.DataFrame({"id": ids_all[rows], "eom": months[t], "w": w}) ) # Update historical buffers with month t's realized returns for future months t' > t ret_m = ret[rows] if np.isfinite(ret_m).all(): go = (Vo / sigma_cs[:, None]).T @ ret_m gn = (Vn / sigma_cs[:, None]).T @ ret_m past_ro.append(go / np.sqrt(np.maximum(np.diag(Qo), 1e-12))) past_rn.append(gn / np.sqrt(np.maximum(np.diag(Qn), 1e-12))) past_qo.append(Qo) past_qn.append(Qn) past_ret_full.append(W_full @ ret_m) if not frames: raise ValueError("No tradable month found in panel.") 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 detected." return out[["id", "eom", "w"]] def main( chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame ) -> pd.DataFrame: """CTF competition entrypoint returning DataFrame with columns (id, eom, w).""" np.random.seed(SEED) weights = build_weights(chars, features, daily_ret) if "ctff_test" not in chars.columns: return weights test = ( chars["ctff_test"] .astype(str) .str.strip() .str.lower() .isin(["1", "1.0", "true", "t"]) ) if not test.any(): 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"), }).drop_duplicates() return weights.merge(keys, on=["id", "eom"], how="inner")