"""Factor-space EPO. CTF submission, freeze v1.0. Each provided characteristic defines a monthly long-short rank portfolio, precision-weighted by inverse idiosyncratic volatility (ivol_capm_252d, clipped cross-sectionally at the 5th and 95th percentiles, power 1). The combination across factor portfolios is shrunk GLS on the factor correlation matrix, (1-w)C + wI: the Simple EPO of Pedersen, Babu and Levine (2021), equivalent to the ridge of Kozak, Nagel and Santosh (2020) at g = w/(1-w). Factor moments are exponentially weighted, halflife 120 months. w is chosen algorithmically by an out-of-sample sweep over formations preceding the first ctff_test month. Stock weights are the factor combination mapped back through the rank portfolios, scaled to a fixed ex-ante volatility target. Causality: formation t uses characteristics through t and realised returns through month t (the lead return labelled at formation t-1). Factor admission requires complete return history from panel start through t-1. daily_ret is accepted per the signature but unused; covariance is factor-level. Determinism: no randomness; sorted groupings; single code path. Dependencies: pre-installed packages only (numpy, pandas, polars, scipy). """ import time import numpy as np import pandas as pd import polars as pl from scipy.linalg import eigh OMEGA_GRID = (0.05, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 0.95) HALFLIFE = 120.0 IV_COL, IV_POWER = "ivol_capm_252d", 1.0 CLIP_LO, CLIP_HI = 0.05, 0.95 VOL_TARGET = 0.10 / np.sqrt(12.0) MIN_SWEEP, DEFAULT_OMEGA = 24, 0.5 MIN_BURN, MAX_BURN = 24, 180 def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame: t_start = time.time() feats = [f for f in features[features.columns[0]].tolist() if f in chars.columns] keep = ["id", "eom", "ret_exc_lead1m"] + ([IV_COL] if IV_COL in chars.columns else []) + \ (["ctff_test"] if "ctff_test" in chars.columns else []) + feats keep = list(dict.fromkeys(keep)) work = chars[keep].copy() work["eom"] = pd.to_datetime(work["eom"]) + pd.offsets.MonthEnd(0) lf = pl.from_pandas(work) del work print("rows", lf.height, "features", len(feats), "prep_s", round(time.time() - t_start, 1), flush=True) # precision weight: within-month clipped inverse ivol, median fill, else 1 if IV_COL in lf.columns: iv = pl.col(IV_COL).cast(pl.Float64) lo = iv.quantile(CLIP_LO).over("eom") hi = iv.quantile(CLIP_HI).over("eom") lf = lf.with_columns( (1.0 / pl.max_horizontal(pl.min_horizontal(iv, hi), lo) ** IV_POWER).alias("p_r") ).with_columns( pl.col("p_r").fill_null(pl.col("p_r").median().over("eom")).fill_null(1.0).alias("p") ) else: lf = lf.with_columns(pl.lit(1.0).alias("p")) # rank * precision, unit gross per factor per month lf = lf.with_columns([ ((pl.col(f).rank("average") - (pl.col(f).count() + 1) / 2).over("eom") .fill_null(0.0) * pl.col("p")).cast(pl.Float32).alias(f) for f in feats ]).with_columns([ pl.when(pl.col(f).abs().sum().over("eom") > 0) .then(pl.col(f) / pl.col(f).abs().sum().over("eom")) .otherwise(0.0).alias(f) for f in feats ]) lead = pl.col("ret_exc_lead1m").cast(pl.Float64) agg = [((pl.col(f) * lead).sum() / pl.when(pl.col(f).abs().sum() > 0).then(pl.col(f).abs().sum())).alias(f) for f in feats] flag = ((pl.col("ctff_test").cast(pl.Utf8).str.to_lowercase() .is_in(["true", "1", "1.0"]).any()) if "ctff_test" in lf.columns else pl.lit(False)).alias("is_test") panel = lf.group_by("eom").agg(agg + [flag]).sort("eom") X = panel.select(feats).to_numpy().astype(np.float64) alive = ~np.isnan(X) Xz = np.nan_to_num(X) is_test = panel["is_test"].to_numpy() eoms = panel["eom"] T = X.shape[0] ft = int(np.argmax(is_test)) if is_test.any() else T if ft >= 2 * MIN_BURN: t0 = min(MAX_BURN, ft // 2) else: t0 = max(2, min(MIN_BURN, max(ft, 4) // 2)) test_ts = [t for t in range(max(t0, 2), T) if is_test[t]] or [T - 1] print("months", T, "first_test", ft, "burn", t0, "test_formations", len(test_ts), flush=True) emap = pl.DataFrame({"eom": eoms, "t": pl.Series(np.arange(T, dtype=np.int64))}) frames = { (k[0] if isinstance(k, tuple) else k): g for k, g in lf.join(emap, on="eom", how="inner") .filter(pl.col("t").is_in([int(t) for t in test_ts])) .sort("id") .partition_by("t", as_dict=True).items() } rho = 0.5 ** (1.0 / HALFLIFE) Wm, W2 = 0.0, 0.0 s, S = np.zeros(X.shape[1]), np.zeros((X.shape[1], X.shape[1])) for t in range(t0): Wm, W2 = rho * Wm + 1.0, rho * rho * W2 + 1.0 s, S = rho * s + Xz[t], rho * S + np.outer(Xz[t], Xz[t]) grid = np.array(OMEGA_GRID) R_pre, omega_star, results, n_done = [], None, [], 0 for t in range(t0, T): take_pre = t < ft take_test = t in frames if take_pre or take_test: idx = np.where(alive[:t].all(0))[0] mu = s[idx] / Wm C = (S[np.ix_(idx, idx)] / Wm - np.outer(mu, mu)) * (Wm * Wm / max(Wm * Wm - W2, 1e-12)) D = np.maximum(np.sqrt(np.diag(C)), 1e-12) m = mu / D lam, V = eigh(C / np.outer(D, D)) b = V.T @ m if take_pre: y = V.T @ (Xz[t, idx] / D) row = np.empty(len(grid)) for i, w_ in enumerate(grid): ls = (1 - w_) * lam + w_ al = b / ls row[i] = (al @ y) / np.sqrt((al * al * ls).sum()) R_pre.append(row) if take_test: if omega_star is None: if len(R_pre) >= MIN_SWEEP: RP = np.asarray(R_pre) sr = RP.mean(0) / RP.std(0, ddof=1) * np.sqrt(12) omega_star = float(grid[int(np.argmax(sr))]) print("sweep", len(R_pre), "formations | SR", [f"{g:.2f}:{v:.2f}" for g, v in zip(grid, sr)], flush=True) else: omega_star = DEFAULT_OMEGA print("sweep too short, default omega", omega_star, flush=True) print("omega*", omega_star, flush=True) ls = (1 - omega_star) * lam + omega_star al = b / ls theta = (V @ al) / D / np.sqrt((al * al * ls).sum()) * VOL_TARGET g = frames[t] Wn = g.select([feats[i] for i in idx]).to_numpy() wv = Wn @ theta if np.isfinite(wv).all() and np.abs(wv).sum() > 0: results.append(pd.DataFrame( {"id": g["id"].to_pandas(), "eom": eoms[t], "w": wv})) n_done += 1 if n_done % 100 == 0: print("solved", n_done, "of", len(test_ts), "live", len(idx), "min", round((time.time() - t_start) / 60.0, 1), flush=True) Wm, W2 = rho * Wm + 1.0, rho * rho * W2 + 1.0 s, S = rho * s + Xz[t], rho * S + np.outer(Xz[t], Xz[t]) if not results: raise ValueError("no weights produced") out = pd.concat(results, ignore_index=True) out["id"] = pd.to_numeric(out["id"]).astype("int64") out["eom"] = pd.to_datetime(out["eom"]).dt.strftime("%Y-%m-%d") out["w"] = out["w"].astype(float) print("rows", len(out), "months", out["eom"].nunique(), "total_min", round((time.time() - t_start) / 60.0, 1), flush=True) return out[["id", "eom", "w"]] def report_test_sharpe(out, chars): """Evaluation convention used throughout development; not called by the pipeline.""" r = chars[["id", "eom", "ret_exc_lead1m"]].copy() r["eom"] = pd.to_datetime(r["eom"]).dt.strftime("%Y-%m-%d") j = out.merge(r, on=["id", "eom"], how="inner") ser = (j["w"] * j["ret_exc_lead1m"]).groupby(j["eom"]).sum().sort_index() sr = ser.mean() / ser.std(ddof=1) * np.sqrt(12) print("TEST months", len(ser), "SR", round(float(sr), 3), "vol", round(float(ser.std(ddof=1)) * np.sqrt(12), 4), flush=True) if __name__ == "__main__": chars = pd.read_parquet("ctff_chars.parquet") features = pd.read_parquet("ctff_features.parquet") daily_ret = pd.read_parquet("ctff_daily_ret.parquet") out = main(chars, features, daily_ret) out.to_csv("output.csv", index=False) print("wrote output.csv rows", len(out), "months", out["eom"].nunique(), flush=True)