""" Yale SOM - MGT 924 Statistical Foundations, Problem Set 3: JKP Common Task Framework (CTF) model. Model: equal-weight ensemble of three forecasters -> 414-factor Barra risk model -> tangency portfolio. 1. Features (point-in-time, cross-sectional only): - 402 characteristics ranked within month x Fama-French-12 industry (percentile - 0.5, missing -> 0). - "fe" bundle (Ridge only): the 20 characteristics with the largest |pre-1975 IC| x market-wide size rank and x Amihud rank (40), last 1/5/10-day returns from daily data (3), 13 JKP theme composites (13). 2. Targets: - rank_res: next-month return (winsorised 1/99%) net of FF12 dummies + standardised beta_252d and market_equity (monthly OLS), ranked within month x FF12. - spec (NN only): next-month return net of FF12 dummies + 13 JKP theme exposures (monthly OLS), divided by the risk model's predicted specific volatility, ranked within month. 3. Forecasters, refit every year on an expanding window from 1952 (forecasts for return year Y use data through Y-1): - XGBoost on the 402 ranks, target rank_res (hyper-parameters tuned once on return months <= Dec 1974). - Ridge on the 402 ranks + fe bundle, target rank_res (penalty tuned once on return months <= Dec 1974, in code). - Neural net (MLP 256-128-64): inputs = the 402 ranks and their products with trailing 1/6/12-month characteristic payoffs (top-minus-bottom-tercile returns, only payoffs realised by month t); target spec; early stopping on the last two years of the window, then a refit on the full window for the chosen epoch count; average of 5 seeds. Ensemble: equal-weight mean of the three forecasts, each z-scored within month. 4. Risk model and portfolio (Jensen, Kelly, Malamud & Pedersen-style Barra model): - Exposures: 12 industry dummies + 402 standardised characteristic ranks; daily ridge factor regressions. - Factor covariance: last 2,520 days, correlations half-life 504 days, variances half-life 84 days. - Specific risk: EWMA (half-life 84 days) of daily residuals; gaps filled by a ridge model of log vol. - Weights: w = Sigma^-1 mu (Woodbury), scaled to 10% predicted annual volatility each month. No look-ahead: everything for month t uses data dated <= t; all selections/tuning use data <= 1974. CTF Admin Modifications (2026-10-08): -------------------------------------- 1. Added requirements.txt pinning torch==2.13.0 and xgboost==3.4.1. Reason: the script imports both packages, neither is pre-installed, and no dependency file was uploaded. The versions are ones other submissions have run with, not necessarily the authors'. 2. In main(), the test month-ends are now converted to timestamps before they are compared with the daily dates. Reason: the pipeline loads eom as plain dates, and comparing those with timestamps stopped the run with "'<=' not supported between instances of 'datetime.datetime' and 'datetime.date'". The conversion changes nothing when eom already arrives as timestamps. """ from __future__ import annotations import os import random import warnings import numpy as np import pandas as pd warnings.filterwarnings("ignore") ID_COL, DATE_COL, RET_COL, TEST_COL = "id", "eom", "ret_exc_lead1m", "ctff_test" TUNE_LAST_YEAR = 1974 # all tuning / selection uses return months <= Dec 1974 NN_SEEDS = tuple(int(x) for x in os.environ.get("CTF_NN_SEEDS", "42,43,44,45,46").split(",")) # approved: 5 seeds N_THREADS = os.cpu_count() or 4 # ---- frozen hyper-parameters (tuned once on return months <= Dec 1974 by blocked 5-fold CV; see documentation) XGB_PARAMS = dict(tree_method="hist", objective="reg:squarederror", base_score=0.0, eval_metric="rmse", booster="gbtree", learning_rate=0.01, max_depth=4, min_split_loss=0.0, subsample=0.9478722068283343, colsample_bytree=0.9477611940298507, min_child_weight=1.0, reg_lambda=2.211168094478662, seed=1) XGB_TREES = 1103 NN_LR, NN_WD, NN_DROPOUT = 3e-4, 1e-4, 0.1 NN_MAX_EPOCHS, NN_PATIENCE, NN_VAL_YEARS = 30, 5, 2 RIDGE_ALPHAS = np.logspace(0, 8, 17) # ---- risk model settings (as in the Markowitz-ML benchmark) RIDGE_LAMBDA = 1e-4 COV_OBS, HL_COR, HL_VAR, HL_STOCK_VAR, INITIAL_VAR_OBS = 2520, 504, 84, 84, 63 VOL_ANN = 0.10 FF12 = ["BusEq", "Chems", "Durbl", "Enrgy", "Hlth", "Manuf", "Money", "NoDur", "Other", "Shops", "Telcm", "Utils"] RES_FACTORS = ["beta_252d", "market_equity"] THEMES = ["Accruals", "Debt Issuance", "Investment", "Low Leverage", "Low Risk", "Momentum", "Profit Growth", "Profitability", "Quality", "Seasonality", "Short-Term Reversal", "Size", "Value"] # JKP (2023) theme membership and signs THEME_MAP = { 'age': ('Low Leverage', -1), 'aliq_at': ('Investment', -1), 'aliq_mat': ('Low Leverage', -1), 'ami_126d': ('Size', 1), 'at_be': ('Low Leverage', -1), 'at_gr1': ('Investment', -1), 'at_me': ('Value', 1), 'at_turnover': ('Quality', 1), 'be_gr1a': ('Investment', -1), 'be_me': ('Value', 1), 'beta_60m': ('Low Risk', -1), 'beta_dimson_21d': ('Low Risk', -1), 'betabab_1260d': ('Low Risk', -1), 'betadown_252d': ('Low Risk', -1), 'bev_mev': ('Value', 1), 'bidaskhl_21d': ('Low Leverage', 1), 'capex_abn': ('Debt Issuance', -1), 'capx_gr1': ('Investment', -1), 'capx_gr2': ('Investment', -1), 'capx_gr3': ('Investment', -1), 'cash_at': ('Low Leverage', 1), 'chcsho_12m': ('Value', -1), 'coa_gr1a': ('Investment', -1), 'col_gr1a': ('Investment', -1), 'cop_at': ('Quality', 1), 'cop_atl1': ('Quality', 1), 'corr_1260d': ('Seasonality', -1), 'coskew_21d': ('Seasonality', -1), 'cowc_gr1a': ('Accruals', -1), 'dbnetis_at': ('Seasonality', -1), 'debt_gr3': ('Debt Issuance', -1), 'debt_me': ('Value', 1), 'dgp_dsale': ('Quality', 1), 'div12m_me': ('Value', 1), 'dolvol_126d': ('Size', -1), 'dolvol_var_126d': ('Profitability', -1), 'dsale_dinv': ('Profit Growth', 1), 'dsale_drec': ('Profit Growth', -1), 'dsale_dsga': ('Profit Growth', 1), 'earnings_variability': ('Low Risk', -1), 'ebit_bev': ('Profitability', 1), 'ebit_sale': ('Profitability', 1), 'ebitda_mev': ('Value', 1), 'emp_gr1': ('Investment', -1), 'eq_dur': ('Value', -1), 'eqnetis_at': ('Value', -1), 'eqnpo_12m': ('Value', 1), 'eqnpo_me': ('Value', 1), 'eqpo_me': ('Value', 1), 'f_score': ('Profitability', 1), 'fcf_me': ('Value', 1), 'fnl_gr1a': ('Debt Issuance', -1), 'gp_at': ('Quality', 1), 'gp_atl1': ('Quality', 1), 'inv_gr1': ('Investment', -1), 'inv_gr1a': ('Investment', -1), 'iskew_capm_21d': ('Short-Term Reversal', -1), 'iskew_ff3_21d': ('Short-Term Reversal', -1), 'iskew_hxz4_21d': ('Short-Term Reversal', -1), 'ival_me': ('Value', 1), 'ivol_capm_21d': ('Low Risk', -1), 'ivol_capm_252d': ('Low Risk', -1), 'ivol_ff3_21d': ('Low Risk', -1), 'ivol_hxz4_21d': ('Low Risk', -1), 'kz_index': ('Seasonality', 1), 'lnoa_gr1a': ('Investment', -1), 'lti_gr1a': ('Seasonality', -1), 'market_equity': ('Size', -1), 'mispricing_mgmt': ('Investment', 1), 'mispricing_perf': ('Quality', 1), 'ncoa_gr1a': ('Investment', -1), 'ncol_gr1a': ('Debt Issuance', -1), 'netdebt_me': ('Low Leverage', -1), 'netis_at': ('Value', -1), 'nfna_gr1a': ('Debt Issuance', 1), 'ni_ar1': ('Debt Issuance', 1), 'ni_be': ('Profitability', 1), 'ni_inc8q': ('Quality', 1), 'ni_ivol': ('Low Leverage', 1), 'ni_me': ('Value', 1), 'niq_at': ('Quality', 1), 'niq_at_chg1': ('Profit Growth', 1), 'niq_be': ('Profitability', 1), 'niq_be_chg1': ('Profit Growth', 1), 'niq_su': ('Profit Growth', 1), 'nncoa_gr1a': ('Investment', -1), 'noa_at': ('Debt Issuance', -1), 'noa_gr1a': ('Investment', -1), 'o_score': ('Profitability', -1), 'oaccruals_at': ('Accruals', -1), 'oaccruals_ni': ('Accruals', -1), 'ocf_at': ('Profitability', 1), 'ocf_at_chg1': ('Profit Growth', 1), 'ocf_me': ('Value', 1), 'ocfq_saleq_std': ('Low Risk', -1), 'op_at': ('Quality', 1), 'op_atl1': ('Quality', 1), 'ope_be': ('Profitability', 1), 'ope_bel1': ('Profitability', 1), 'opex_at': ('Quality', 1), 'pi_nix': ('Seasonality', 1), 'ppeinv_gr1a': ('Investment', -1), 'prc': ('Size', -1), 'prc_highprc_252d': ('Momentum', 1), 'qmj': ('Quality', 1), 'qmj_growth': ('Quality', 1), 'qmj_prof': ('Quality', 1), 'qmj_safety': ('Quality', 1), 'rd_me': ('Size', 1), 'rd_sale': ('Low Leverage', 1), 'rd5_at': ('Low Leverage', 1), 'resff3_12_1': ('Momentum', 1), 'resff3_6_1': ('Momentum', 1), 'ret_1_0': ('Short-Term Reversal', -1), 'ret_12_1': ('Momentum', 1), 'ret_12_7': ('Profit Growth', 1), 'ret_3_1': ('Momentum', 1), 'ret_6_1': ('Momentum', 1), 'ret_60_12': ('Investment', -1), 'ret_9_1': ('Momentum', 1), 'rmax1_21d': ('Low Risk', -1), 'rmax5_21d': ('Low Risk', -1), 'rmax5_rvol_21d': ('Short-Term Reversal', -1), 'rskew_21d': ('Short-Term Reversal', -1), 'rvol_21d': ('Low Risk', -1), 'sale_bev': ('Quality', 1), 'sale_emp_gr1': ('Profit Growth', 1), 'sale_gr1': ('Investment', -1), 'sale_gr3': ('Investment', -1), 'sale_me': ('Value', 1), 'saleq_gr1': ('Investment', -1), 'saleq_su': ('Profit Growth', 1), 'seas_1_1an': ('Profit Growth', 1), 'seas_1_1na': ('Momentum', 1), 'seas_11_15an': ('Seasonality', 1), 'seas_11_15na': ('Seasonality', -1), 'seas_16_20an': ('Seasonality', 1), 'seas_16_20na': ('Accruals', -1), 'seas_2_5an': ('Seasonality', 1), 'seas_2_5na': ('Investment', -1), 'seas_6_10an': ('Seasonality', 1), 'seas_6_10na': ('Low Risk', -1), 'sti_gr1a': ('Seasonality', 1), 'taccruals_at': ('Accruals', -1), 'taccruals_ni': ('Accruals', -1), 'tangibility': ('Low Leverage', 1), 'tax_gr1a': ('Profit Growth', 1), 'turnover_126d': ('Low Risk', -1), 'turnover_var_126d': ('Profitability', -1), 'z_score': ('Low Leverage', 1), 'zero_trades_126d': ('Low Risk', 1), 'zero_trades_21d': ('Low Risk', 1), 'zero_trades_252d': ('Low Risk', 1), } def log(msg: str): print(f"[model] {msg}", flush=True) # ===================================================================================================================== # 1. Data preparation # ===================================================================================================================== def ff12_class(sic: pd.Series) -> np.ndarray: """Fama-French 12 industries from SIC codes (benchmark ranges and precedence).""" s = pd.to_numeric(sic, errors="coerce").to_numpy(dtype=float) def r(a, b): return (s >= a) & (s <= b) rules = [ ("NoDur", r(100, 999) | r(2000, 2399) | r(2700, 2749) | r(2770, 2799) | r(3100, 3199) | r(3940, 3989)), ("Durbl", r(2500, 2519) | r(3630, 3659) | r(3710, 3711) | (s == 3714) | (s == 3716) | r(3750, 3751) | (s == 3792) | r(3900, 3939) | r(3990, 3999)), ("Manuf", r(2520, 2589) | r(2600, 2699) | r(2750, 2769) | r(3000, 3099) | r(3200, 3569) | r(3580, 3629) | r(3700, 3709) | r(3712, 3713) | (s == 3715) | r(3717, 3749) | r(3752, 3791) | r(3793, 3799) | r(3830, 3839) | r(3860, 3899)), ("Enrgy", r(1200, 1399) | r(2900, 2999)), ("Chems", r(2800, 2829) | r(2840, 2899)), ("BusEq", r(3570, 3579) | r(3660, 3692) | r(3694, 3699) | r(3810, 3829) | r(7370, 7379)), ("Telcm", r(4800, 4899)), ("Utils", r(4900, 4949)), ("Shops", r(5000, 5999) | r(7200, 7299) | r(7600, 7699)), ("Hlth", r(2830, 2839) | (s == 3693) | r(3840, 3859) | r(8000, 8099)), ("Money", r(6000, 6999)), ] out = np.full(len(s), "Other", dtype=object) done = np.zeros(len(s), dtype=bool) for name, m in rules: m = m & ~done out[m] = name done |= m return out def pct_rank(block: pd.DataFrame, by, min_obs: int | None) -> np.ndarray: """rank(ties='max') / n_non_missing within groups, exact zeros -> 0, minus 0.5; NaN kept (benchmark ranking).""" g = block.groupby(by) n = g.transform("count") r = g.rank(method="max") / n r = r.where(block != 0, 0.0).where(block.notna()) if min_obs is not None: r = r.where(n >= min_obs) return (r - 0.5).to_numpy(dtype=np.float32) def rank_within(v: np.ndarray, groups) -> np.ndarray: r = pd.Series(v).groupby(groups).rank(pct=True) - 0.5 return r.fillna(0.0).to_numpy(np.float32) def prepare(chars: pd.DataFrame, features: list[str]): """keys + within-industry ranks (X_ind), market-wide ranks (X_mkt), risk-model exposures (X_fm).""" df = chars[[ID_COL, DATE_COL, "eom_ret", "sic", RET_COL] + features].copy() df[DATE_COL] = pd.to_datetime(df[DATE_COL]) df["eom_ret"] = pd.to_datetime(df["eom_ret"]) df = df.sort_values([DATE_COL, ID_COL], kind="mergesort").reset_index(drop=True) keys = pd.DataFrame({"id": df[ID_COL].astype(np.int64), "eom": df[DATE_COL], "eom_ret": df["eom_ret"], "ret": pd.to_numeric(df[RET_COL], errors="coerce").astype(float), "ff12": ff12_class(df["sic"])}) keys["month"] = (keys["eom_ret"].dt.year * 12 + keys["eom_ret"].dt.month - 1).to_numpy() keys["ret_year"] = keys["eom_ret"].dt.year.to_numpy() n, p = len(df), len(features) X_mkt = np.empty((n, p), np.float32) X_ind = np.empty((n, p), np.float32) R_fm = np.empty((n, p), np.float32) by_ind = [keys["eom"], keys["ff12"]] for j0 in range(0, p, 25): c = features[j0:j0 + 25] blk = df[c].astype(np.float64) X_mkt[:, j0:j0 + len(c)] = pct_rank(blk, keys["eom"], None) X_ind[:, j0:j0 + len(c)] = pct_rank(blk, by_ind, None) R_fm[:, j0:j0 + len(c)] = pct_rank(blk, keys["eom"], 10) del df for A in (X_mkt, X_ind, R_fm): np.nan_to_num(A, copy=False, nan=0.0) # risk-model exposures: standardise within eom; constant columns get N(0,1) noise (benchmark convention) rng = np.random.default_rng(1) eom = keys["eom"].to_numpy() starts = np.flatnonzero(np.r_[True, eom[1:] != eom[:-1]]) ends = np.r_[starts[1:], n] for a, b in zip(starts, ends): blk = R_fm[a:b].astype(np.float64) sd = blk.std(axis=0, ddof=1) const = np.round(sd, 6) == 0 if const.any(): blk[:, const] += rng.standard_normal((b - a, const.sum())) sd = blk.std(axis=0, ddof=1) R_fm[a:b] = ((blk - blk.mean(axis=0)) / sd).astype(np.float32) inds = [i for i in FF12 if (keys["ff12"] == i).any()] D = np.stack([(keys["ff12"] == i).to_numpy(np.float32) for i in inds], axis=1) X_fm = np.concatenate([D, R_fm], axis=1) del R_fm return keys, X_ind, X_mkt, X_fm, inds + features, len(inds) # ===================================================================================================================== # 2. Risk model (daily factor regressions, factor covariances, specific risk) # ===================================================================================================================== def ridge_glmnet(X, y): n, k = X.shape return np.linalg.solve(X.T @ X + n * RIDGE_LAMBDA * np.eye(k), X.T @ y) def ewma_vol(x, lam, start): from scipy.signal import lfilter out = np.full(len(x), np.nan) if len(x) <= start: return out v0 = np.sum(x[:start] ** 2) / (start - 1) out[start] = v0 if len(x) > start + 1: out[start + 1:] = lfilter([1 - lam], [1, -lam], x[start:-1] ** 2, zi=[lam * v0])[0] return np.sqrt(out) def weighted_cov(Z, hl): t = len(Z) w = (0.5 ** (1 / hl)) ** np.arange(t, 0, -1) w = w / w.sum() Zc = Z - w @ Z return (Zc * w[:, None]).T @ Zc / (1 - np.sum(w ** 2)) def risk_model(daily_ret: pd.DataFrame, keys: pd.DataFrame, X_fm: np.ndarray, n_ind: int, eval_eoms: np.ndarray): dt = pd.to_datetime(daily_ret["date"]) d = daily_ret.loc[(dt >= keys["eom"].min()) & (dt <= keys["eom"].max()) & np.isfinite(daily_ret["ret_exc"]), [ID_COL, "date", "ret_exc"]].copy() d["date"] = pd.to_datetime(d["date"]) mk = lambda ids, dates: ids.astype(np.int64) * 100_000 + dates.dt.year * 12 + dates.dt.month - 1 row = pd.Index(mk(keys["id"], keys["eom_ret"])).get_indexer(mk(d[ID_COL], d["date"])) d = d.assign(row=row)[row >= 0].sort_values(["date", ID_COL], kind="mergesort").reset_index(drop=True) # daily cross-sectional ridge regressions on the exposures of the row whose return month contains the day dates = d["date"].to_numpy() starts = np.flatnonzero(np.r_[True, dates[1:] != dates[:-1]]) ends = np.r_[starts[1:], len(d)] rows, r = d["row"].to_numpy(), d["ret_exc"].to_numpy() F = np.empty((len(starts), X_fm.shape[1])) resid = np.empty(len(d)) for i, (a, b) in enumerate(zip(starts, ends)): X = X_fm[rows[a:b]].astype(np.float64) f = ridge_glmnet(X, r[a:b]) resid[a:b] = r[a:b] - X @ f absent = ~X[:, :n_ind].any(axis=0) if absent.any(): f[:n_ind][absent] = np.median(f[:n_ind][~absent]) F[i] = f fct_dates = dates[starts] log(f"risk model: {len(fct_dates):,} daily factor regressions, K = {F.shape[1]}") # factor covariances at the portfolio eoms (data <= eom) covs = {} for e in eval_eoms: Z = F[fct_dates <= e][-COV_OBS:] C = weighted_cov(Z, HL_COR) sdc = np.sqrt(np.diag(C)) cor = C / np.outer(sdc, sdc) sd = np.sqrt(np.diag(weighted_cov(Z, HL_VAR))) covs[pd.Timestamp(e)] = cor * np.outer(sd, sd) # specific risk: EWMA of each stock's residuals, >= 200 own observations in the last 252 trading days sr = pd.DataFrame({"id": d[ID_COL].to_numpy(), "date": d["date"].to_numpy(), "res": resid}) del d sr = sr.sort_values(["id", "date"], kind="mergesort").reset_index(drop=True) lam = 0.5 ** (1 / HL_STOCK_VAR) ids = sr["id"].to_numpy() st = np.flatnonzero(np.r_[True, ids[1:] != ids[:-1]]) en = np.r_[st[1:], len(sr)] res, vol = sr["res"].to_numpy(), np.empty(len(sr)) for a, b in zip(st, en): vol[a:b] = ewma_vol(res[a:b], lam, INITIAL_VAR_OBS) sr["res_vol"] = vol pos = np.searchsorted(fct_dates, sr["date"].to_numpy()) td_252 = np.where(pos >= 252, fct_dates[np.maximum(pos - 252, 0)], np.datetime64("NaT")) sr["date_200"] = sr.groupby("id")["date"].shift(200) sr = sr[(sr["date_200"].to_numpy() >= td_252) & sr["res_vol"].notna()] sr["eom"] = sr["date"] + pd.offsets.MonthEnd(0) sr = sr.groupby(["id", "eom"], sort=False).tail(1) loc = pd.Index(mk(keys["id"], keys["eom"])).get_indexer(mk(sr["id"], sr["eom"])) res_vol = np.full(len(keys), np.nan) res_vol[loc[loc >= 0]] = sr["res_vol"].to_numpy()[loc >= 0] # fill missing specific risk with a monthly ridge of log vol on the same-eom exposures (all eoms) for rows_e in keys.groupby("eom").indices.values(): have = np.isfinite(res_vol[rows_e]) if have.sum() < 30 or have.all(): continue X = X_fm[rows_e].astype(np.float64) b = ridge_glmnet(X[have], np.log(res_vol[rows_e][have])) res_vol[rows_e[~have]] = np.exp(X[~have] @ b) return covs, res_vol # ===================================================================================================================== # 3. Targets, engineered features, characteristic payoffs # ===================================================================================================================== def winsorised_residual(keys, E) -> np.ndarray: r, out = keys["ret"].to_numpy(), np.full(len(keys), np.nan) for rows in keys.groupby("eom").indices.values(): rows = rows[np.isfinite(r[rows])] if len(rows) < 10: continue y = np.clip(r[rows], *np.quantile(r[rows], [0.01, 0.99])) b = np.linalg.lstsq(E[rows], y, rcond=None)[0] out[rows] = y - E[rows] @ b return out def theme_exposures(X_fm, keys, cols, n_ind): pos = {c: j for j, c in enumerate(cols)} T = np.empty((len(keys), len(THEMES))) for t, th in enumerate(THEMES): mem = [c for c in THEME_MAP if THEME_MAP[c][0] == th and c in pos] T[:, t] = np.mean([THEME_MAP[c][1] * X_fm[:, pos[c]] for c in mem], axis=0) T = pd.DataFrame(T).groupby(keys["eom"].to_numpy()).transform(lambda x: (x - x.mean()) / x.std()) return np.concatenate([X_fm[:, :n_ind], T.to_numpy()], axis=1) def targets(keys, X_fm, cols, n_ind, res_vol): E = np.asarray(X_fm[:, [cols.index(c) for c in cols[:n_ind] + RES_FACTORS]], np.float64) res = winsorised_residual(keys, E) rank_res = (pd.Series(res).groupby([keys["eom"].to_numpy(), keys["ff12"].to_numpy()]).rank(pct=True) - 0.5).to_numpy() Xt = theme_exposures(X_fm, keys, cols, n_ind).astype(np.float64) e = winsorised_residual(keys, Xt) v = res_vol * np.sqrt(21) s = np.where(np.isfinite(v) & (v > 0), e / v, np.nan) spec = (pd.Series(s).groupby(keys["eom"].to_numpy()).rank(pct=True) - 0.5).to_numpy() return rank_res.astype(np.float32), spec.astype(np.float32) def fe_bundle(keys, X_ind, X_mkt, features, daily_ret) -> np.ndarray: """Ridge's engineered features; the top-20 characteristics are selected on return months <= Dec 1974.""" pos = {c: j for j, c in enumerate(features)} groups = [keys["eom"].to_numpy(), keys["ff12"].to_numpy()] y = (keys.groupby(["eom", "ff12"])["ret"].rank(pct=True) - 0.5).to_numpy() tr = (keys["ret_year"].to_numpy() <= TUNE_LAST_YEAR) & np.isfinite(y) eom = keys["eom"].to_numpy()[tr] Xt, yt = X_ind[tr].astype(np.float64), y[tr] ics = [] for m in np.unique(eom): r = eom == m xc, yc = Xt[r] - Xt[r].mean(axis=0), yt[r] - yt[r].mean() sd = np.sqrt((xc ** 2).sum(axis=0) * (yc ** 2).sum()) ics.append(np.divide(xc.T @ yc, sd, out=np.zeros(len(sd)), where=sd > 0)) top = np.argsort(-np.abs(np.mean(ics, axis=0)))[:20] blocks = [] for cond in ("market_equity", "ami_126d"): z = X_mkt[:, pos[cond]] for j in top: blocks.append(X_ind[:, j] * z * (2 * np.sqrt(3))) # last 1/5/10-day returns up to eom dt = pd.to_datetime(daily_ret["date"]) d = daily_ret.loc[(dt <= keys["eom"].max()) & daily_ret[ID_COL].isin(keys["id"].unique()), [ID_COL, "date", "ret_exc"]].copy() d["date"] = pd.to_datetime(d["date"]) d = d.sort_values([ID_COL, "date"], kind="mergesort").reset_index(drop=True) lr = np.log1p(d["ret_exc"].clip(lower=-0.99).fillna(0.0).to_numpy()) cs = d.assign(cs=lr).groupby(ID_COL)["cs"].cumsum().to_numpy() p_ = d.groupby(ID_COL).cumcount().to_numpy() out = {"id": d[ID_COL].to_numpy(), "eom": (d["date"] + pd.offsets.MonthEnd(0)).to_numpy()} for h in (1, 5, 10): lag = np.r_[np.full(h, np.nan), cs[:-h]] lag = np.where(p_ >= h, lag, np.where(p_ == h - 1, 0.0, np.nan)) out[f"r{h}"] = np.expm1(cs - lag) dd = pd.DataFrame(out).groupby(["id", "eom"], sort=False).tail(1) k = keys[["id", "eom"]].merge(dd, on=["id", "eom"], how="left") for h in (1, 5, 10): blocks.append(rank_within(k[f"r{h}"].to_numpy(), groups)) for th in THEMES: mem = [(pos[c], s) for c, (t, s) in THEME_MAP.items() if t == th and c in pos] blocks.append(rank_within(np.mean([s * X_ind[:, j] for j, s in mem], axis=0), groups)) return np.column_stack(blocks).astype(np.float32) def payoff_tables(keys, X_ind): """Per eom: P1, P6, P12 trailing characteristic payoffs (formation months <= t-1), normalised across chars.""" eoms = np.sort(keys["eom"].unique()) rows = keys.groupby("eom").indices ret = keys["ret"].to_numpy() pay = np.zeros((len(eoms), X_ind.shape[1])) for i, e in enumerate(eoms): r = rows[pd.Timestamp(e)] r = r[np.isfinite(ret[r])] if len(r) < 20: continue y = np.clip(ret[r], *np.quantile(ret[r], [0.01, 0.99])) x = X_ind[r] top, bot = x > 1 / 6, x < -1 / 6 pay[i] = (np.where(top, y[:, None], 0).sum(0) / np.maximum(top.sum(0), 1) - np.where(bot, y[:, None], 0).sum(0) / np.maximum(bot.sum(0), 1)) pay[i][(top.sum(0) < 5) | (bot.sum(0) < 5)] = 0.0 P = np.zeros((len(eoms), 3, X_ind.shape[1])) for i in range(1, len(eoms)): P[i, 0] = pay[i - 1] P[i, 1] = pay[max(0, i - 6):i].mean(0) P[i, 2] = pay[max(0, i - 12):i].mean(0) P = P / (P.std(2, keepdims=True) + 1e-12) return np.searchsorted(eoms, keys["eom"].to_numpy()), P.astype(np.float32) # ===================================================================================================================== # 4. Forecasters # ===================================================================================================================== def ridge_forecasts(keys, X, y, test_years): """Ridge with unpenalised intercept; alpha tuned by blocked 5-fold CV on return months <= Dec 1974.""" month = keys["month"].to_numpy() ok = np.isfinite(y) months = np.unique(month[ok]) p = X.shape[1] + 1 G, c, yy, nn = np.empty((len(months), p, p)), np.empty((len(months), p)), np.empty(len(months)), np.empty(len(months)) for i, m in enumerate(months): r = np.flatnonzero((month == m) & ok) Xa = np.column_stack([X[r].astype(np.float64), np.ones(len(r))]) G[i], c[i], yy[i], nn[i] = Xa.T @ Xa, Xa.T @ y[r], y[r] @ y[r], len(r) def solve(Gs, cs, a): pen = np.full(p, a) pen[-1] = 0.0 return np.linalg.solve(Gs + np.diag(pen), cs) keep = np.flatnonzero(months <= TUNE_LAST_YEAR * 12 + 11) Gt, ct = G[keep].sum(0), c[keep].sum(0) mse = np.zeros(len(RIDGE_ALPHAS)) for blk in np.array_split(keep, 5): Gv, cv, yv, nv = G[blk].sum(0), c[blk].sum(0), yy[blk].sum(), nn[blk].sum() for j, a in enumerate(RIDGE_ALPHAS): b = solve(Gt - Gv, ct - cv, a) mse[j] += (yv - 2 * b @ cv + b @ Gv @ b) / nv / 5 alpha = float(RIDGE_ALPHAS[np.argmin(mse)]) log(f"ridge: alpha = {alpha:g} (pre-1975 CV)") ry = keys["ret_year"].to_numpy() pred = np.full(len(keys), np.nan) for Y in test_years: tr = months <= (Y - 1) * 12 + 11 b = solve(G[tr].sum(0), c[tr].sum(0), alpha) te = np.flatnonzero(ry == Y) pred[te] = X[te].astype(np.float64) @ b[:-1] + b[-1] return pred def xgb_forecasts(keys, X, y, test_years): """XGBoost with the frozen settings, refit every year on an expanding window.""" import xgboost as xgb X = np.ascontiguousarray(X, dtype=np.float32) y = y.astype(np.float64) ry = keys["ret_year"].to_numpy() params = dict(XGB_PARAMS, nthread=N_THREADS) pred = np.full(len(y), np.nan) for Y in test_years: idx = np.flatnonzero((ry <= Y - 1) & np.isfinite(y)) bst = xgb.train(params, xgb.QuantileDMatrix(X[idx], y[idx], nthread=N_THREADS), num_boost_round=XGB_TREES) te = np.flatnonzero(ry == Y) pred[te] = bst.inplace_predict(np.ascontiguousarray(X[te])) log(f"xgboost {Y}: {len(idx):,} training rows") return pred def nn_forecasts(keys, X_ind, y, mi, P, test_years): """MLP on [x, x*P1, x*P6, x*P12]; early stopping on the last 2 years, full-window refit; mean of 5 seeds.""" import torch torch.set_num_threads(N_THREADS) torch.use_deterministic_algorithms(True) Xt = torch.from_numpy(np.ascontiguousarray(X_ind)) Pt = torch.from_numpy(P) mit = torch.from_numpy(mi.astype(np.int64)) yt = torch.from_numpy(np.nan_to_num(y, nan=0.0).astype(np.float32)) valid = np.isfinite(y) month, ry = keys["month"].to_numpy(), keys["ret_year"].to_numpy() p_in = X_ind.shape[1] * 4 def inputs(rows): x = Xt[rows] return torch.cat([x, (x[:, None, :] * Pt[mit[rows]]).flatten(1)], 1) def net(): layers, w = [], p_in for h in (256, 128, 64): layers += [torch.nn.Linear(w, h), torch.nn.SiLU(), torch.nn.Dropout(NN_DROPOUT)] w = h return torch.nn.Sequential(*layers, torch.nn.Linear(w, 1)) def batches(rows, only_valid=True): if only_valid: rows = rows[valid[rows]] cuts = np.flatnonzero(np.diff(month[rows])) + 1 return [torch.from_numpy(r) for r in np.split(rows, cuts) if len(r) > 1] def rank_ic(pv, yv): a, b = pv.argsort().argsort().float(), yv.argsort().argsort().float() a, b = a - a.mean(), b - b.mean() return float((a @ b) / (a.norm() * b.norm() + 1e-12)) def train(rows_b, epochs, seed, val_b=None): torch.manual_seed(seed) np.random.seed(seed) model = net() opt = torch.optim.AdamW(model.parameters(), lr=NN_LR, weight_decay=NN_WD) g = np.random.default_rng(seed) best, best_ep, best_state, wait = -np.inf, epochs, None, 0 for ep in range(1, epochs + 1): model.train() for i in g.permutation(len(rows_b)): b = rows_b[i] loss = torch.mean((model(inputs(b)).squeeze(-1) - yt[b]) ** 2) opt.zero_grad() loss.backward() opt.step() if val_b is None: continue model.eval() with torch.no_grad(): ic = float(np.mean([rank_ic(model(inputs(b)).squeeze(-1), yt[b]) for b in val_b])) if ic > best: best, best_ep, wait = ic, ep, 0 best_state = {k: v.detach().clone() for k, v in model.state_dict().items()} else: wait += 1 if wait >= NN_PATIENCE: break if val_b is not None: model.load_state_dict(best_state) return model, best_ep preds = np.zeros(len(keys)) for Y in test_years: tr = np.flatnonzero(ry < Y - NN_VAL_YEARS) va = np.flatnonzero((ry >= Y - NN_VAL_YEARS) & (ry < Y)) te = np.flatnonzero(ry == Y) bte = batches(te, only_valid=False) out = [] for seed in NN_SEEDS: _, e_star = train(batches(tr), NN_MAX_EPOCHS, seed, val_b=batches(va)) model, _ = train(batches(np.concatenate([tr, va])), e_star, seed) # full-window refit model.eval() with torch.no_grad(): pr = torch.cat([model(inputs(b)).squeeze(-1) for b in bte]).numpy() rows = torch.cat(bte).numpy() z = pd.Series(pr).groupby(month[rows]).transform(lambda s: (s - s.mean()) / s.std()).to_numpy() o = np.zeros(len(keys)) o[rows] = z out.append(o) preds[te] = np.mean(out, axis=0)[te] log(f"nn {Y}: {len(tr):,} training rows, {len(NN_SEEDS)} seeds") preds[~np.isin(ry, test_years)] = np.nan return preds # ===================================================================================================================== # 5. Portfolio # ===================================================================================================================== def woodbury_solve(D, Sf, X, b): Sf = (Sf + Sf.T) / 2 Db, DX = b / D, X / D[:, None] M = np.linalg.inv(Sf) + X.T @ DX M = (M + M.T) / 2 return Db - DX @ np.linalg.solve(M, X.T @ Db) def zscore_by(v: np.ndarray, groups: np.ndarray) -> np.ndarray: return pd.Series(v).groupby(groups).transform(lambda s: (s - s.mean()) / s.std()).to_numpy() # ===================================================================================================================== # main # ===================================================================================================================== def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame: random.seed(1) np.random.seed(1) col = "features" if "features" in features.columns else features.columns[0] feats = sorted(set(features[col].astype(str)) & set(chars.columns)) if "eom_ret" not in chars.columns: chars = chars.assign(eom_ret=pd.to_datetime(chars[DATE_COL]) + pd.offsets.MonthEnd(1)) test = chars[TEST_COL].astype(str).str.strip().isin(["1", "True", "true"]) if TEST_COL in chars.columns \ else pd.Series(True, index=chars.index) test_keys = chars.loc[test, [ID_COL, DATE_COL]].assign(**{ID_COL: lambda d: d[ID_COL].astype(np.int64)}) log(f"{len(chars):,} rows, {len(feats)} features, {test.sum():,} test rows") keys, X_ind, X_mkt, X_fm, cols, n_ind = prepare(chars, feats) # CTF-FIX: eom can arrive as plain dates, the comparisons below need timestamps test_eoms = np.sort(pd.to_datetime(test_keys[DATE_COL]).unique()) test_years = sorted({(pd.Timestamp(e) + pd.offsets.MonthEnd(1)).year for e in test_eoms}) log(f"prepared; test return years {test_years[0]}-{test_years[-1]}") covs, res_vol = risk_model(daily_ret, keys, X_fm, n_ind, test_eoms) rank_res, spec = targets(keys, X_fm, cols, n_ind, res_vol) X_fe = fe_bundle(keys, X_ind, X_mkt, feats, daily_ret) del X_mkt mi, P = payoff_tables(keys, X_ind) log("targets, engineered features and payoffs built") g = keys["eom"].to_numpy() p_ridge = ridge_forecasts(keys, np.concatenate([X_ind, X_fe], axis=1), rank_res, test_years) p_xgb = xgb_forecasts(keys, X_ind, rank_res, test_years) p_nn = nn_forecasts(keys, X_ind, spec, mi, P, test_years) in_test = np.isin(keys["ret_year"].to_numpy(), test_years) mu = np.nanmean(np.column_stack([zscore_by(p, g) for p in (p_xgb, p_ridge, p_nn)]), axis=1) mu[~in_test] = np.nan out = [] for e in test_eoms: rows = np.flatnonzero((g == e) & np.isfinite(mu) & np.isfinite(res_vol)) if len(rows) < 2: continue X = X_fm[rows].astype(np.float64) D = res_vol[rows] ** 2 Sf = covs[pd.Timestamp(e)] w = woodbury_solve(D, Sf, X, mu[rows]) Xw = X.T @ w w *= np.sqrt(VOL_ANN ** 2 / 12 / ((np.sum(w ** 2 * D) + Xw @ Sf @ Xw) * 21)) out.append(pd.DataFrame({ID_COL: keys["id"].to_numpy()[rows], DATE_COL: e, "w": w})) w = pd.concat(out, ignore_index=True) w[DATE_COL] = pd.to_datetime(w[DATE_COL]) test_keys[DATE_COL] = pd.to_datetime(test_keys[DATE_COL]) w = w.merge(test_keys, on=[ID_COL, DATE_COL], how="inner") w[ID_COL] = w[ID_COL].astype("int64") w[DATE_COL] = w[DATE_COL].dt.strftime("%Y-%m-%d") w["w"] = pd.to_numeric(w["w"], errors="coerce").fillna(0.0).astype("float64") w = w[[ID_COL, DATE_COL, "w"]].sort_values([DATE_COL, ID_COL]).reset_index(drop=True) log(f"returning {len(w):,} rows across {w[DATE_COL].nunique()} months") return w if __name__ == "__main__": import time from pathlib import Path data = Path(__file__).resolve().parent / "Data" t0 = time.time() weights = main(pd.read_parquet(data / "ctff_chars.parquet"), pd.read_parquet(data / "ctff_features.parquet"), pd.read_parquet(data / "ctff_daily_ret.parquet")) out_path = Path(__file__).resolve().parent / "portfolio_weights_model.csv" weights.to_csv(out_path, index=False) print(f"Saved {len(weights):,} rows to {out_path} ({time.time() - t0:.0f}s)")