""" Markowitz-ML for the JKP Common Task Framework (MGT 924, PS3). Every month hold w = Sigma^-1 mu (optionally blended with the forecast-only portfolio w = mu) mu : ridge-regression forecast of each stock's next-month return, relative to the market Sigma : Barra-style factor covariance = industry factors + K characteristic "themes" + idiosyncratic What is in here and why (each point validated on the full data, see docs/step2_v2_validation.md): + one risk factor per 2-digit SIC industry instead of a single market factor (+0.4 test Sharpe) + K = 40 characteristic themes in the risk model (+0.4) + "blend" rule: half Markowitz, half forecast portfolio (same Sharpe, less leverage) - industry-demeaned forecast target: redundant once Sigma has industry factors (dropped) - lasso / elastic net / ensemble forecasts: no out-of-sample gain over ridge (dropped) Sections, in the order things happen: 0. KNOBS what you can change 1. prepare data -> one polars DataFrame (ranked characteristics, industry, idio vol, target) 2. split TRAIN | VALIDATION | TEST (three consecutive blocks of months) 3. fit_mu ridge regression -> mu_t = X_t @ beta 4. fit_risk factor model -> Sigma_t = B_t F B_t' + D_t 5. portfolios weight rules + backtest 6. run grid over the knobs on VALIDATION (prints Sharpe), refit, run TEST 7. main competition entry point: main(chars, features, daily_ret) -> DataFrame[id, eom, w] Nothing used at month t comes from after month t: characteristics are what was known at t, and all fitting and knob selection happen on months before the block being evaluated. Run locally: python code/markowitz_ml.py CTF Admin Modifications (2026-10-08): -------------------------------------- 1. In prepare(), the conversion of sic to a number no longer stops on a value that is not a number. Reason: the full dataset has blank sic values, which stopped the run. A blank is now treated as missing and gets the same "unknown industry" code (-1) that missing values already got. """ from datetime import date from pathlib import Path import numpy as np import polars as pl from scipy.cluster.hierarchy import fcluster, linkage # ============================================================================= # 0. KNOBS # ============================================================================= VAL_START = date(1980, 1, 1) # TRAIN = months before this, VALIDATION = from here to TEST # TEST = months flagged ctff_test == 1 (Dec-1989 .. Nov-2023) LAMBDAS = [1e-3, 1e-2] # ridge penalties to try (characteristics are ranks, variance ~ 1/12) N_FACTORS = [20, 40] # number of characteristic themes in the risk model to try INDUSTRY_RISK = [False, True] # risk model: single market factor, or one factor per industry BLEND = 0.5 # weight on the forecast-only portfolio inside the "blend" rule SELECT_AMONG = ["Sigma^-1 mu", "blend"] # rules the validation step may choose between TARGET_VOL = 0.10 # annual ex-ante volatility every portfolio is scaled to SEED = 0 ROOT = Path(__file__).resolve().parents[1] FEATS: list[str] = [] # characteristic names, filled by prepare() # ============================================================================= # 1. PREPARE -> one polars DataFrame, one row per stock-month # ============================================================================= def prepare(chars, features) -> pl.DataFrame: """chars = ctff_chars table, features = table with a 'features' column (polars or pandas).""" chars = chars if isinstance(chars, pl.DataFrame) else pl.from_pandas(chars) features = features if isinstance(features, pl.DataFrame) else pl.from_pandas(features) FEATS[:] = [f for f in features["features"].to_list() if f in chars.columns] r = pl.col("ret_exc_lead1m") panel = ( chars.lazy() .select(["id", "eom", "ret_exc_lead1m", "ctff_test", "sic"] + FEATS) .with_columns(pl.col(pl.Float32, pl.Float64).fill_nan(None)) # NaN and null both = missing .with_columns( id=pl.col("id").cast(pl.Int64), eom=pl.col("eom").cast(pl.Date), is_test=pl.col("ctff_test").cast(pl.String).is_in(["1", "true", "True"]), # CTF-FIX: a blank sic becomes missing (then -1) instead of stopping the run ind=(pl.col("sic").cast(pl.Int64, strict=False) // 100).fill_null(-1), # 2-digit SIC industry ivol=pl.col("ivol_capm_252d").cast(pl.Float64), # raw trailing idio vol (daily) # target = next-month return relative to the month's average: we only need to rank stocks y=r - r.mean().over("eom"), ) # rank every characteristic within the month to [-0.5, 0.5]; missing -> 0 (= the median) .with_columns( [((pl.col(f).rank().over("eom") - 0.5) / pl.col(f).count().over("eom") - 0.5).fill_null(0.0) for f in FEATS] ) # missing idio vol -> month's median (overall median if the whole month is missing), never zero .with_columns( ivol=pl.col("ivol").fill_null(pl.col("ivol").median().over("eom")).fill_null(pl.col("ivol").median()) .fill_null(0.02).clip(1e-4) ) .collect() .sort(["eom", "id"]) ) return panel # ============================================================================= # 2. SPLIT -> three consecutive blocks of months # ============================================================================= def split(panel: pl.DataFrame): """TRAIN and VALIDATION need the next-month return (to fit / score); TEST keeps every row.""" pre = panel.filter(~pl.col("is_test") & pl.col("ret_exc_lead1m").is_not_null()) train = pre.filter(pl.col("eom") < VAL_START) val = pre.filter(pl.col("eom") >= VAL_START) test = panel.filter(pl.col("is_test")) return train, val, test # ============================================================================= # 3. EXPECTED RETURNS -> ridge regression (slides eq. 5) # ============================================================================= def fit_mu(panel: pl.DataFrame, lam: float) -> np.ndarray: """beta = (X'X/n + lam*I)^-1 X'y/n. No intercept: X and y are demeaned within each month.""" X = panel.select(FEATS).to_numpy() y = panel["y"].to_numpy() n = len(y) return np.linalg.solve(X.T @ X / n + lam * np.eye(X.shape[1]), X.T @ y / n) # ============================================================================= # 4. RISK MODEL -> Sigma_t = B_t F B_t' + D_t # B_t : exposures = industry dummies (or a single intercept) + K theme exposures # F : covariance of the monthly factor returns # D_t : idiosyncratic variance of each stock = c * ivol_t^2 # A fitted risk model is a dict {"L", "codes", "F", "c"}. # ============================================================================= def make_themes(panel: pl.DataFrame, K: int) -> np.ndarray: """ Group the characteristics into K themes by hierarchical clustering on their correlation. Returns a p x K matrix L so that a stock's theme exposures are X @ L = its average rank within each theme. """ corr = np.corrcoef(panel.select(FEATS).to_numpy(), rowvar=False) dist = np.clip(1 - corr[np.triu_indices(len(FEATS), k=1)], 0, None) labels = fcluster(linkage(dist, method="average"), K, criterion="maxclust") L = np.zeros((len(FEATS), K)) for k in range(K): members = labels == k + 1 L[members, k] = 1 / members.sum() return L def exposures(X: np.ndarray, ind: np.ndarray, risk: dict) -> np.ndarray: """B = [industry dummies (or 1), X L].""" themes = X @ risk["L"] if risk["codes"] is None: # single market factor return np.column_stack([np.ones(len(X)), themes]) ind = np.where(np.isin(ind, risk["codes"]), ind, -1) # industries unseen when fitting -> "unknown" return np.column_stack([(ind[:, None] == risk["codes"][None, :]).astype(float), themes]) def fit_risk(panel: pl.DataFrame, K: int, industry: bool) -> dict: """ Each month regress realized returns on exposures: r_t = B_t f_t + e_t. F = cov(f_t) factor covariance c = typical e^2 / ivol^2 maps trailing daily idio vol^2 into monthly residual variance """ codes = np.union1d(panel["ind"].unique().to_numpy(), [-1]) if industry else None risk = {"L": make_themes(panel, K), "codes": codes} f_hist, ratio = [], [] for (eom,), g in panel.group_by("eom", maintain_order=True): B = exposures(g.select(FEATS).to_numpy(), g["ind"].to_numpy(), risk) r = g["ret_exc_lead1m"].to_numpy() f = np.linalg.lstsq(B, r, rcond=None)[0] # this month's factor returns e = r - B @ f # this month's idiosyncratic returns f_hist.append(f) ratio.append(np.median(e**2) / np.median(g["ivol"].to_numpy() ** 2)) F = np.cov(np.array(f_hist).T) risk["F"] = F + 1e-6 * np.mean(np.diag(F)) * np.eye(len(F)) # tiny ridge keeps F invertible risk["c"] = float(np.median(ratio)) return risk def month_risk(g: pl.DataFrame, risk: dict): """(B, F, D) for one month.""" B = exposures(g.select(FEATS).to_numpy(), g["ind"].to_numpy(), risk) D = risk["c"] * g["ivol"].to_numpy() ** 2 return B, risk["F"], D def sigma_solve(B, F, D, v): """ Sigma^-1 v for Sigma = B F B' + diag(D). Same result as np.linalg.solve(B @ F @ B.T + np.diag(D), v) but via the Woodbury identity, so the N x N matrix is never built (N ~ 2,000 stocks, ~100 factors -> ~100x faster). """ DinvB = B / D[:, None] small = np.linalg.inv(F) + B.T @ DinvB # (#factors x #factors) return v / D - DinvB @ np.linalg.solve(small, B.T @ (v / D)) def scale_to_vol(w, B, F, D): """Rescale w so that its ex-ante volatility sqrt(w' Sigma w) equals TARGET_VOL (annual).""" Bw = B.T @ w var = Bw @ F @ Bw + np.sum(D * w**2) return w * (TARGET_VOL / np.sqrt(12)) / np.sqrt(var) if var > 0 else w # ============================================================================= # 5. PORTFOLIOS AND BACKTEST # A rule maps (mu, B, F, D) -> weights scaled to TARGET_VOL. Add a rule to the dict to try it. # ============================================================================= def rule_mu(mu, B, F, D): return scale_to_vol(mu, B, F, D) # forecasts only def rule_minvar(mu, B, F, D): return scale_to_vol(sigma_solve(B, F, D, np.ones_like(mu)), B, F, D) # risk only def rule_markowitz(mu, B, F, D): return scale_to_vol(sigma_solve(B, F, D, mu), B, F, D) # both def rule_blend(mu, B, F, D): # mix of the two above: hedges against estimation error in Sigma return scale_to_vol((1 - BLEND) * rule_markowitz(mu, B, F, D) + BLEND * rule_mu(mu, B, F, D), B, F, D) PORTFOLIOS = {"I mu": rule_mu, "Sigma^-1 1": rule_minvar, "Sigma^-1 mu": rule_markowitz, "blend": rule_blend} def backtest(panel: pl.DataFrame, beta: np.ndarray, risk: dict, keep: str | None = None): """Run every rule through every month. Returns monthly returns (+ the weights of rule `keep`).""" rows, weights = [], [] for (eom,), g in panel.group_by("eom", maintain_order=True): mu = g.select(FEATS).to_numpy() @ beta # forecasts for this month B, F, D = month_risk(g, risk) # risk model for this month r = np.nan_to_num(g["ret_exc_lead1m"].to_numpy()) # realized next-month returns (0 if unknown) row = {"eom": eom} for name, rule in PORTFOLIOS.items(): w = rule(mu, B, F, D) row[name] = w @ r if name == keep: weights.append(pl.DataFrame({"id": g["id"], "w": w}).with_columns(eom=pl.lit(eom))) rows.append(row) rets = pl.DataFrame(rows) return (rets, pl.concat(weights).select(["id", "eom", "w"])) if keep else rets def sharpes(rets: pl.DataFrame) -> dict: """Annualized Sharpe ratio of every rule's monthly returns.""" return {p: float(rets[p].mean() / rets[p].std() * np.sqrt(12)) for p in PORTFOLIOS} def show(s: dict) -> str: return " ".join(f"{p} = {v:5.2f}" for p, v in s.items()) # ============================================================================= # 6. RUN -> select knobs on VALIDATION, refit on TRAIN + VALIDATION, produce TEST results # ============================================================================= def run(panel: pl.DataFrame, verbose: bool = True): np.random.seed(SEED) train, val, test = split(panel) if verbose: for name, block in [("train", train), ("val", val), ("test", test)]: print(f"{name:5s}: {block['eom'].n_unique():4d} months {block.height:9,} stock-months " f"{block['eom'].min()} .. {block['eom'].max()}") print("\n== VALIDATION Sharpe (fit on train) ==") results = [] for industry in INDUSTRY_RISK: for K in N_FACTORS: risk = fit_risk(train, K, industry) for lam in LAMBDAS: beta = fit_mu(train, lam) s = sharpes(backtest(val, beta, risk)) results.append({"industry": industry, "K": K, "lambda": lam, **s}) if verbose: print(f"industry={str(industry):5s} K={K:2d} lambda={lam:<6g} {show(s)}") # best (knobs, rule) by validation Sharpe among the rules we allow best = max(((row, rule) for row in results for rule in SELECT_AMONG), key=lambda rr: rr[0][rr[1]]) knobs, rule = best if verbose: print(f"\nselected: industry={knobs['industry']} K={knobs['K']} lambda={knobs['lambda']} rule={rule}") # refit once on TRAIN + VALIDATION with the chosen knobs (static model), then score TEST trainval = pl.concat([train, val]) risk = fit_risk(trainval, knobs["K"], knobs["industry"]) beta = fit_mu(trainval, knobs["lambda"]) rets, weights = backtest(test, beta, risk, keep=rule) if verbose: print(f"\n== TEST Sharpe ({rets.height} months) ==\n{show(sharpes(rets))}") return rets, weights, pl.DataFrame(results), {**knobs, "rule": rule} # ============================================================================= # 7. COMPETITION ENTRY POINT # ============================================================================= def main(chars, features, daily_ret=None): """Competition contract: pandas DataFrames in, pandas DataFrame [id, eom, w] out.""" import pandas as pd panel = prepare(chars, features) _, weights, _, _ = run(panel, verbose=False) out = weights.to_pandas() out["id"] = out["id"].astype("int64") out["eom"] = pd.to_datetime(out["eom"]).dt.date out["w"] = out["w"].astype("float64").fillna(0.0) return out[["id", "eom", "w"]] if __name__ == "__main__": OUT = ROOT / "output" OUT.mkdir(exist_ok=True) panel = prepare(pl.read_parquet(ROOT / "data/raw/ctff_chars.parquet"), pl.read_parquet(ROOT / "data/raw/ctff_features.parquet")) rets, weights, results, chosen = run(panel) rets.write_csv(OUT / "step3_monthly_returns.csv") weights.write_parquet(OUT / "step3_weights.parquet") results.write_csv(OUT / "step3_validation_grid.csv") print(f"\nsaved to {OUT}")