"""CTF model: cross-sectional rank features -> ridge + elastic net ensemble -> long-short weights. Pipeline (everything uses only information available at the time it is used): 1. Each month, every characteristic is turned into its cross-sectional ECDF rank, centred at zero; missing values are set to 0 (the cross-sectional median). The target (next-month excess return) is rank-transformed the same way, per month. 2. Models are trained on months with ctff_test == False only. Because the models are linear, training only needs the per-segment cross-products X'X, X'y, y'y, so the full panel never has to be held in memory as one big matrix. 3. Hyper-parameters are chosen by expanding-window time-series validation inside the training sample (5 folds; each fold trains on the past and validates on the next block). 4. The final predictor is the average of the standardised ridge and elastic-net predictions. The map from predictions to weights (rank / z-score / decile) is chosen by Sharpe ratio on the two most recent validation blocks. 5. Estimation uses an expanding window with equal month weights. During the output (test) period the model is re-estimated every 12 months on an expanding window that contains only months strictly before the month being predicted (their realised returns are known by then). Hyper-parameters stay fixed at the values chosen in step 3. 6. Mean-variance weights: with the single-factor risk model Sigma = s_mkt^2 * b b' + diag(s_u^2) (b = market beta, s_u = idiosyncratic volatility, both characteristics in the data) and alpha = volatility * score, the weights are Sigma^{-1} alpha (Woodbury identity), capped at 4x the mean absolute weight and scaled to gross exposure 2. The weight map (rank / z-score / decile / mean-variance) is reported on the validation blocks. (Without these columns: within the long and the short leg, weights are divided by the stock's recent volatility (a characteristic already in the data), so that no single volatile stock dominates portfolio risk; each leg is re-normalised to sum to +1 / -1.) 7. Volatility management: the portfolio is scaled each month by (long-run strategy volatility) / (volatility of its last L realised monthly returns), bounded to [0.25, 2]. Whether to use it and which L (6/12/24) is chosen on the validation blocks, never on the output period. Realised returns used at month t are those of earlier months only. 8. Unscaled weights are dollar-neutral with gross exposure 2 (long 1, short 1); rebalanced monthly. """ import numpy as np import pandas as pd SEED = 42 N_MIN = 30 # months with fewer stocks are skipped for training CUTS = [0.5, 0.6, 0.7, 0.8, 0.9, 1.0] # expanding-window fold boundaries (share of months) RIDGE_REL = [0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10.0, 30.0] # ridge penalty relative to mean(diag(Q)) EN_REL = [0.5, 0.2, 0.1, 0.05, 0.02, 0.01] # elastic-net penalty relative to its lambda_max EN_ALPHA = [0.5, 1.0] # l1 share (1.0 = lasso) BETA_PROXIES = ["betabab_1260d", "beta_60m", "betadown_252d", "beta_dimson_21d"] # market beta, first found is used IVOL_PROXIES = ["ivol_capm_21d", "ivol_ff3_21d", "ivol_capm_252d", "ivol_hxz4_21d"] # idiosyncratic volatility USE_MVO = True # use the mean-variance weights (single-factor risk model) if the columns exist HALF_LIFE = None # months; exponential down-weighting of old months (None = equal weights; 120 was tested: worse in validation) HALF_LIVES_INFO = [None, 120, 60, 36] # reported for information only, not used for selection REFIT = 12 # re-estimate the model every 12 output months SCHEMES = ["rank", "zscore", "decile"] VOL_PROXIES = ["rvol_21d", "ivol_capm_21d", "ivol_ff3_21d", "ivol_capm_252d", "rvol_252d"] # stock risk proxy, first found is used VOL_WINDOWS = [None, 6, 12, 24] # volatility-management windows (months); None = off VOL_LO, VOL_HI = 0.25, 2.0 # bounds of the volatility scale factor DIAG = {} # diagnostics for the write-up (ignored by the pipeline) # ----------------------------------------------------------------------------- helpers def _to_bool(col): """ctff_test may arrive as bool, 0/1 numbers or strings ('0'/'1'/'True'/'False').""" if col.dtype == bool: return col.to_numpy() v = col.astype(str).str.strip().str.lower() return v.isin(["1", "1.0", "true", "t", "yes", "y"]).to_numpy() # ----------------------------------------------------------------------------- transforms def _rank_features(X): """Cross-sectional ECDF rank of each column, centred at 0; NaN -> 0.""" r = pd.DataFrame(X).rank(axis=0, pct=True, method="average").to_numpy(dtype=np.float32) return np.nan_to_num(r - np.float32(0.5), nan=0.0) def _rank_target(y): return pd.Series(y).rank(pct=True, method="average").to_numpy(dtype=np.float64) - 0.5 # ----------------------------------------------------------------------------- models def _ridge(Q, c, lam): return np.linalg.solve(Q + lam * np.eye(len(c)), c) def _enet(Q, c, lam, a, beta0=None, tol=1e-8, max_sweeps=300): """Cyclic coordinate descent for 0.5 b'Qb - c'b + lam*(a|b|_1 + (1-a)/2 |b|_2^2).""" p = len(c) beta = np.zeros(p) if beta0 is None else beta0.copy() g = Q @ beta d = np.diag(Q) + lam * (1.0 - a) l1 = lam * a for _ in range(max_sweeps): biggest = 0.0 for j in range(p): if d[j] <= 0: continue bj = beta[j] rho = c[j] - g[j] + Q[j, j] * bj new = np.sign(rho) * max(abs(rho) - l1, 0.0) / d[j] if new != bj: g += Q[j] * (new - bj) beta[j] = new biggest = max(biggest, abs(new - bj)) if biggest < tol: break return beta def _sse(G, b, yy, beta): return float(yy - 2.0 * beta @ b + beta @ G @ beta) # ----------------------------------------------------------------------------- weights def _zs(p): s = p.std() return (p - p.mean()) / s if s > 0 else np.zeros_like(p) def _clean_sigma(v): """Stock risk proxy for one month: invalid -> median, winsorised at the 5th/95th percentile.""" v = np.asarray(v, dtype=float).copy() v[~np.isfinite(v) | (v <= 0)] = np.nan if np.isnan(v).all(): return None v[np.isnan(v)] = np.nanmedian(v) lo, hi = np.percentile(v, [5, 95]) return np.clip(v, lo, hi) def _risk_inputs(rv, iv, be): """One month of risk inputs: total vol, idiosyncratic vol, beta (cleaned) and a market-variance estimate.""" rv, iv = _clean_sigma(rv), _clean_sigma(iv) be = np.asarray(be, dtype=float).copy() if rv is None or iv is None: return None be[~np.isfinite(be)] = np.nan if np.isnan(be).all(): return None be[np.isnan(be)] = np.nanmedian(be) be = np.clip(be, -0.5, 3.0) est = max(float(np.sum(rv ** 2 - iv ** 2)) / float(np.sum(be ** 2)), 1e-8) # daily market variance return {"rv": rv, "iv": iv, "be": be, "est": est} def _mvo_weights(p, rk, sm): """Sigma^{-1} alpha for Sigma = sm * b b' + diag(iv^2); alpha = total vol * score.""" n = len(p) if n < 5 or rk is None: return np.zeros(n) mu = rk["rv"] * _zs(p) d_inv = 1.0 / rk["iv"] ** 2 be = rk["be"] db = d_inv * be w = d_inv * mu - db * (sm / (1.0 + sm * float(be @ db))) * float(be @ (d_inv * mu)) s = np.abs(w).sum() if not np.isfinite(s) or s <= 0: return np.zeros(n) w = 2.0 * w / s c = 4.0 * np.abs(w).mean() w = np.clip(w, -c, c) return 2.0 * w / np.abs(w).sum() def _weights(p, scheme, sigma=None, rk=None, sm=None): """Long-short weights: dollar-neutral, long leg sums to +1, short leg to -1 (gross 2).""" n = len(p) if n < 2: return np.zeros(n) if scheme == "mvo": return _mvo_weights(p, rk, sm) if scheme == "rank": v = pd.Series(p).rank(method="average").to_numpy() - (n + 1) / 2.0 elif scheme == "zscore": v = np.clip(_zs(p), -3.0, 3.0) v = v - v.mean() else: # decile: top and bottom 10% k = max(1, n // 10) order = np.argsort(p, kind="stable") v = np.zeros(n) v[order[:k]] = -1.0 v[order[-k:]] = 1.0 if sigma is not None: v = v / sigma pos, neg = v > 0, v < 0 if not pos.any() or not neg.any(): return np.zeros(n) w = np.zeros(n) w[pos] = v[pos] / v[pos].sum() w[neg] = v[neg] / (-v[neg].sum()) return w def _vol_scale(hist, L): h = np.asarray(hist, dtype=float) if len(h) < max(L, 12): return 1.0 s_recent, s_all = h[-L:].std(ddof=1), h.std(ddof=1) if s_recent <= 0 or s_all <= 0: return 1.0 return float(np.clip(s_all / s_recent, VOL_LO, VOL_HI)) def _sharpe(r): r = np.asarray(r, dtype=float) return float(r.mean() / r.std(ddof=1) * np.sqrt(12)) if len(r) > 2 and r.std(ddof=1) > 0 else -np.inf # ----------------------------------------------------------------------------- main def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame: np.random.seed(SEED) banned = {"id", "eom", "ret_exc_lead1m", "ctff_test"} feats = list(dict.fromkeys(f for f in features["features"].tolist() if f in chars.columns and f not in banned)) cols = [chars.columns.get_loc(f) for f in feats] ycol = chars.columns.get_loc("ret_exc_lead1m") p = len(feats) risk_name = next((c for c in VOL_PROXIES if c in chars.columns), None) vcol = chars.columns.get_loc(risk_name) if risk_name else None beta_name = next((c for c in BETA_PROXIES if c in chars.columns), None) ivol_name = next((c for c in IVOL_PROXIES if c in chars.columns), None) risk_ok = bool(risk_name and beta_name and ivol_name) rcols = [vcol, chars.columns.get_loc(ivol_name), chars.columns.get_loc(beta_name)] if risk_ok else None schemes = SCHEMES + (["mvo"] if risk_ok else []) DIAG["risk_proxy"], DIAG["beta_proxy"], DIAG["ivol_proxy"] = risk_name, beta_name, ivol_name print(f"[model] risk columns: total vol {risk_name}, idiosyncratic vol {ivol_name}, beta {beta_name}", flush=True) groups = chars.groupby("eom", sort=True).indices months = sorted(groups) flag = _to_bool(chars["ctff_test"]) is_test = {m: bool(flag[groups[m]].all()) for m in months} test_m = [m for m in months if is_test[m]] or months[-12:] n_min = int(min(N_MIN, max(5, np.median([len(groups[m]) for m in months]) // 2))) train_m = [m for m in months if not is_test[m] and len(groups[m]) >= n_min] if not train_m: # e.g. a small validation sample with only test months print("[warn] no months flagged ctff_test == False; training on all months", flush=True) train_m = [m for m in months if len(groups[m]) >= n_min] ids_all, eom_all = chars["id"].to_numpy(), chars["eom"].to_numpy() def zero_output(): parts = [pd.DataFrame({"id": ids_all[groups[m]], "eom": eom_all[groups[m]], "w": np.zeros(len(groups[m]))}) for m in test_m] return pd.concat(parts, ignore_index=True) if p == 0 or not train_m: return zero_output() M = len(train_m) print(f"[model] {len(months)} months, {M} training months, {len(test_m)} output months, {p} features", flush=True) # ---- fold layout: expanding window, validation block k+1 follows training blocks 0..k use_cv = M >= 60 if use_cv: cut = [int(round(M * c)) for c in CUTS] seg_of = np.searchsorted(np.array(cut), np.arange(M), side="right") nseg = len(cut) else: seg_of, nseg = np.zeros(M, dtype=int), 1 Gm = np.zeros((M, p, p), dtype=np.float32) # per-month cross-products X'X, X'y, y'y, n bm = np.zeros((M, p)) ym = np.zeros(M) nm = np.zeros(M) lam = 0.5 ** (1.0 / HALF_LIFE) if HALF_LIFE else 1.0 def agg(lo, hi, end=None, lam_=None): """Sum of the month statistics in [lo, hi); weight lam**(end-1-s) if `end` is given, else 1.""" lm = lam if lam_ is None else lam_ Gs, bs, ys, ns = np.zeros((p, p)), np.zeros(p), 0.0, 0.0 for s_ in range(lo, hi): wt = 1.0 if end is None else lm ** (end - 1 - s_) Gs += wt * Gm[s_] bs += wt * bm[s_] ys += wt * ym[s_] ns += wt * nm[s_] return Gs, bs, ys, ns cache = {} # months of the last two validation blocks sm_hist = [] # monthly market-variance estimates (past only) def month_risk(ix): if not risk_ok: return None V = chars.iloc[ix, rcols].to_numpy(dtype=np.float64) return _risk_inputs(V[:, 0], V[:, 1], V[:, 2]) # ---- pass 1: cross-products per segment for i, m in enumerate(train_m): ix = groups[m] X = _rank_features(chars.iloc[ix, cols].to_numpy(dtype=np.float32)) y = chars.iloc[ix, ycol].to_numpy(dtype=np.float64) ok = ~np.isnan(y) rk = month_risk(ix) if rk is not None: sm_hist.append(rk["est"]) if ok.sum() < n_min: continue Xv, yv = X[ok].astype(np.float64), _rank_target(y[ok]) s = seg_of[i] Gm[i] = Xv.T @ Xv bm[i] = Xv.T @ yv ym[i] = yv @ yv nm[i] = len(yv) if use_cv and s >= nseg - 2: # raw returns kept for validation performance sg = _clean_sigma(chars.iloc[ix, vcol].to_numpy()) if vcol is not None else None rko = None if rk is None else {k_: rk[k_][ok] for k_ in ("rv", "iv", "be")} sm_now = float(np.mean(sm_hist[-12:])) if sm_hist else 1e-4 cache[m] = (X[ok], s, y[ok], None if sg is None else sg[ok], rko, sm_now) if nm.sum() == 0: return zero_output() if sm_hist: DIAG["implied_market_vol_annual"] = float(np.sqrt(np.median(sm_hist) * 252)) print(f"[model] implied market volatility (annualised, median): " f"{DIAG['implied_market_vol_annual']:.1%}", flush=True) # ---- time-series validation of hyper-parameters (uses only cross-products) ridge_i, en_ia, en_ir = 3, 1, 3 # defaults if validation is not possible scheme, vol_L, r_val = "rank", None, np.array([]) if use_cv: r_sc = np.zeros(len(RIDGE_REL)) e_sc = np.zeros((len(EN_ALPHA), len(EN_REL))) zero_sc = 0.0 for k in range(nseg - 1): Gt, bt, _, nt = agg(0, cut[k], end=cut[k]) Q, c = Gt / nt, bt / nt d0 = np.mean(np.diag(Q)) Gv, bv, yvv, _ = agg(cut[k], cut[k + 1]) zero_sc += yvv for i, rel in enumerate(RIDGE_REL): r_sc[i] += _sse(Gv, bv, yvv, _ridge(Q, c, rel * d0)) for ia, a in enumerate(EN_ALPHA): beta = None for ir, rel in enumerate(EN_REL): beta = _enet(Q, c, rel * np.max(np.abs(c)) / a, a, beta) e_sc[ia, ir] += _sse(Gv, bv, yvv, beta) ridge_i = int(np.argmin(r_sc)) en_ia, en_ir = np.unravel_index(int(np.argmin(e_sc)), e_sc.shape) DIAG["oos_r2_ridge"] = 1.0 - r_sc[ridge_i] / zero_sc DIAG["oos_r2_enet"] = 1.0 - e_sc[en_ia, en_ir] / zero_sc DIAG["ridge_rel_grid"], DIAG["ridge_val_sse"] = RIDGE_REL, (r_sc / zero_sc).tolist() print(f"[cv] ridge rel={RIDGE_REL[ridge_i]} (R2 {DIAG['oos_r2_ridge']:.5f}); " f"enet alpha={EN_ALPHA[en_ia]} rel={EN_REL[en_ir]} (R2 {DIAG['oos_r2_enet']:.5f})", flush=True) info = {} for H in HALF_LIVES_INFO: # information only: validation R2 of ridge by half-life lm = 0.5 ** (1.0 / H) if H else 1.0 tot, sse = 0.0, 0.0 for k in range(nseg - 1): Gt, bt, _, nt = agg(0, cut[k], end=cut[k], lam_=lm) Q, c = Gt / nt, bt / nt Gv, bv, yvv, _ = agg(cut[k], cut[k + 1]) sse += _sse(Gv, bv, yvv, _ridge(Q, c, RIDGE_REL[ridge_i] * np.mean(np.diag(Q)))) tot += yvv info[str(H)] = 1.0 - sse / tot DIAG["val_r2_by_half_life"] = info print(f"[cv] validation R2 of ridge by half-life (information only) {info}; half-life used {HALF_LIFE}", flush=True) # ---- weight scheme: Sharpe of RAW monthly returns on the last two validation blocks betas = {} for k in (nseg - 3, nseg - 2): # folds whose validation blocks are cached Gt, bt, _, nt = agg(0, cut[k], end=cut[k]) Q, c = Gt / nt, bt / nt a = EN_ALPHA[en_ia] betas[k + 1] = (_ridge(Q, c, RIDGE_REL[ridge_i] * np.mean(np.diag(Q))), _enet(Q, c, EN_REL[en_ir] * np.max(np.abs(c)) / a, a)) rets = {s_: [] for s_ in schemes} dec = [] for m, (X, s, yr, sg, rko, sm_now) in cache.items(): br, be = betas[s] pred = 0.5 * (_zs(X @ br) + _zs(X @ be)) for sc in schemes: rets[sc].append(float(_weights(pred, sc, sg, rko, sm_now) @ yr)) order = np.argsort(pred, kind="stable") dec.append([float(yr[g].mean()) for g in np.array_split(order, 10)]) shp = {sc: _sharpe(rets[sc]) for sc in schemes} scheme = "mvo" if (USE_MVO and "mvo" in shp) else max(shp, key=shp.get) r_val = np.array(rets[scheme]) DIAG["val_sharpe_by_scheme"] = shp DIAG["val_months"] = len(dec) DIAG["val_decile_mean_return"] = np.mean(dec, axis=0).tolist() DIAG["val_returns_by_scheme"] = rets print(f"[cv] validation Sharpe (raw returns) by scheme {shp} -> using '{scheme}'" + (" (design choice)" if scheme == "mvo" else ""), flush=True) # ---- volatility management: choose window (or none) on the validation months wins = [L for L in VOL_WINDOWS if L] start = max(wins) if wins else 0 if wins and len(r_val) >= start + 30: vs = {} for L in VOL_WINDOWS: sr = [(1.0 if L is None else _vol_scale(r_val[:i], L)) * r_val[i] for i in range(start, len(r_val))] vs[str(L)] = _sharpe(sr) best = max(vs, key=vs.get) # ties -> first entry (None = off) vol_L = None if best == "None" else int(best) DIAG["val_sharpe_by_vol_window"] = vs print(f"[cv] validation Sharpe by volatility window {vs} -> using {vol_L}", flush=True) DIAG.update(scheme=scheme, vol_window=vol_L, ridge_rel=RIDGE_REL[ridge_i], en_alpha=EN_ALPHA[en_ia], en_rel=EN_REL[en_ir]) # ---- pass 2: predict month by month; re-estimate every REFIT months on past months only a = EN_ALPHA[en_ia] train_set = set(train_m) run_G, run_b, _, run_n = agg(0, M, end=M) # decayed sums over all training months beta_r = beta_e = None hist = list(r_val) if vol_L else [] # realised (unscaled) strategy returns known so far out, scales, n_refit = [], [], 0 for idx, m in enumerate(test_m): ix = groups[m] if idx % REFIT == 0: # refit: run_* hold months strictly before m Q, c = run_G / run_n, run_b / run_n beta_r = _ridge(Q, c, RIDGE_REL[ridge_i] * np.mean(np.diag(Q))) beta_e = _enet(Q, c, EN_REL[en_ir] * np.max(np.abs(c)) / a, a, beta_e) n_refit += 1 if n_refit == 1: DIAG["enet_nonzero_first"] = int((beta_e != 0).sum()) top = np.argsort(-np.abs(beta_e))[:10] DIAG["top_enet_features_first"] = [(feats[j], float(beta_e[j])) for j in top] X = _rank_features(chars.iloc[ix, cols].to_numpy(dtype=np.float32)) pred = 0.5 * (_zs(X @ beta_r) + _zs(X @ beta_e)) sg = _clean_sigma(chars.iloc[ix, vcol].to_numpy()) if vcol is not None else None rk = month_risk(ix) if rk is not None: sm_hist.append(rk["est"]) sm_t = float(np.mean(sm_hist[-12:])) if sm_hist else 1e-4 w0 = _weights(pred, scheme, sg, rk, sm_t) sc_t = _vol_scale(hist, vol_L) if vol_L else 1.0 scales.append((str(m), sc_t)) out.append(pd.DataFrame({"id": ids_all[ix], "eom": eom_all[ix], "w": sc_t * w0})) # month m's realised returns become usable from the next month on y = chars.iloc[ix, ycol].to_numpy(dtype=np.float64) ok = ~np.isnan(y) if m not in train_set: # one month older: decay the running sums run_G *= lam run_b *= lam run_n *= lam if ok.sum() >= n_min: if vol_L: hist.append(float(w0[ok] @ y[ok])) if m not in train_set: Xv, yv = X[ok].astype(np.float64), _rank_target(y[ok]) run_G += Xv.T @ Xv run_b += Xv.T @ yv run_n += len(yv) DIAG["half_life"] = HALF_LIFE DIAG["n_refits"] = n_refit DIAG["scale_series"] = scales DIAG["enet_nonzero_last"] = int((beta_e != 0).sum()) top = np.argsort(-np.abs(beta_e))[:10] DIAG["top_enet_features_last"] = [(feats[j], float(beta_e[j])) for j in top] print(f"[fit] {n_refit} refits; elastic net keeps {DIAG['enet_nonzero_first']} -> " f"{DIAG['enet_nonzero_last']} of {p} features; volatility window {vol_L}", flush=True) return pd.concat(out, ignore_index=True)