"""MGT 924 PS3 / JKP Common Task Framework submission. Youyang (Yolanda) Deng and Jiayu Meng. main(chars, features, daily_ret) -> DataFrame[id, eom, w] Three candidate models, chosen with CANDIDATE below: "A" Ridge return forecasts + 1-factor dollar-neutral mean-variance (Jiayu) "B" XGBoost rank forecasts + 20-factor PCA risk model Markowitz (Yolanda) "AB" 50/50 blend of A and B, each scaled by its trailing realised volatility No look-ahead: every quantity used for the portfolio at eom t is computable at t. - Features are ranked / standardised within each month only. - A row with eom s carries the return over month s+1, so models used at t are trained on rows with eom <= t - 1 month only (re-estimated once a year). - Risk models use daily returns dated <= t. - Blend weights use realised strategy returns of eoms < t. Deterministic: fixed seeds and a fixed XGBoost thread count. CTF Admin Modifications (2026-10-07): -------------------------------------- 1. Added [CTF-DEBUG] progress print statements (with flush=True) in main(), _xgb_predictions() and _pca_markowitz(). Reason: the script printed nothing, so a long run left an empty job log. The statements only report progress and do not change any result. """ from datetime import date, timedelta import numpy as np import pandas as pd import polars as pl import xgboost as xgb CANDIDATE = "B" # chosen by the pre-registered holdout rule (output/holdout_selection_rule.md) SEED = 20251004 XGB_THREADS = 8 # fixed so results do not depend on the machine VAL_YEARS = 5 # XGBoost validation window (years) XGB_GRID = [2, 4, 6] # max_depth candidates, re-chosen every year XGB_PARAMS = { "objective": "reg:squarederror", "tree_method": "hist", "max_bin": 64, "learning_rate": 0.05, "subsample": 0.5, "colsample_bytree": 0.5, "min_child_weight": 100, "reg_lambda": 1.0, "seed": SEED, "nthread": XGB_THREADS, } XGB_MAX_ROUNDS, XGB_EARLY_STOP, XGB_FALLBACK = 1000, 50, (4, 100) PCA_FACTORS, PCA_SHRINK = 20, 0.5 RISK_WINDOW, RISK_HALF_LIFE, RISK_MIN_OBS = 252, 126, 200 VOL_TARGET = 0.10 # annualised ex-ante volatility of Model B each month BLEND_BURN_IN = 36 # months of history used to scale the AB blend BLEND_MIN_OBS = 12 # ---------------------------------------------------------------------------- # Data preparation # ---------------------------------------------------------------------------- def _month_before(d: date) -> date: return d.replace(day=1) - timedelta(days=1) def _prepare(chars: pd.DataFrame, feats: list[str]) -> pl.DataFrame: df = pl.from_pandas(chars[["id", "eom", "ret_exc_lead1m", "ctff_test"] + feats]) return df.with_columns( pl.col("id").cast(pl.Int64), pl.col("eom").cast(pl.Date), pl.col("ret_exc_lead1m").cast(pl.Float64), pl.col("ctff_test").cast(pl.Utf8).str.to_lowercase().is_in(["1", "true", "t", "1.0"]), *[pl.col(f).cast(pl.Float64) for f in feats], ).sort(["eom", "id"]) def _rank_features(df: pl.DataFrame, feats: list[str]) -> pl.DataFrame: """Cross-sectional rank of each feature within the month, in [-0.5, 0.5].""" def ranked(f): r, n = pl.col(f).rank("average").over("eom"), pl.col(f).count().over("eom") return pl.when(n > 1).then((r - 1) / (n - 1) - 0.5).otherwise(0.0).cast(pl.Float32).alias(f) return df.with_columns([ranked(f) for f in feats]) def _zscore_features(df: pl.DataFrame, feats: list[str]) -> pl.DataFrame: """Model A: impute the monthly median, then z-score within the month.""" df = df.with_columns([pl.col(f).fill_null(pl.col(f).median().over("eom")).fill_null(0.0) for f in feats]) def z(f): # A column that is constant within the month (e.g. one observed value then # median-imputed) must map to 0; rounding can leave a ~1e-15 sd, so test # relative to the column's scale rather than for exact zero. sd, mean = pl.col(f).std().over("eom"), pl.col(f).mean().over("eom") tiny = sd.is_null() | (sd <= 1e-10 * (mean.abs() + 1.0)) return pl.when(tiny).then(0.0).otherwise((pl.col(f) - mean) / sd).alias(f) return df.with_columns([z(f) for f in feats]) def _prepare_daily(daily_ret: pd.DataFrame, last: date) -> pl.DataFrame: return (pl.from_pandas(daily_ret[["id", "date", "ret_exc"]]) .with_columns(pl.col("id").cast(pl.Int64), pl.col("date").cast(pl.Date), pl.col("ret_exc").cast(pl.Float64)) .filter((pl.col("date") <= last) & pl.col("ret_exc").is_not_null()) .sort("date")) # ---------------------------------------------------------------------------- # Model B: XGBoost forecasts # ---------------------------------------------------------------------------- def _qdm(df: pl.DataFrame, feats: list[str]) -> xgb.QuantileDMatrix: """Quantised matrix: same splits as a DMatrix under tree_method='hist' (checked: identical predictions) but ~1 byte per cell, so the full sample fits in memory.""" X = df.select(feats).to_numpy() return xgb.QuantileDMatrix(X, label=df["y"].to_numpy(), max_bin=XGB_PARAMS["max_bin"], nthread=XGB_THREADS) def _xgb_fit(labelled: pl.DataFrame, i_val: int, i_end: int, feats: list[str]): """Rows [0, i_val) train, [i_val, i_end) validate (labelled is sorted by eom). Pick max_depth / rounds on val (lowest MSE), then refit on rows [0, i_end).""" if i_end == 0: return None depth, rounds = XGB_FALLBACK if 0 < i_val < i_end: dtr = _qdm(labelled.slice(0, i_val), feats) # plain DMatrix for the (small) validation set: XGBoost caches its predictions # across boosting rounds, which it does not do for a QuantileDMatrix eval set val = labelled.slice(i_val, i_end - i_val) dva = xgb.DMatrix(val.select(feats).to_numpy(), label=val["y"].to_numpy(), nthread=XGB_THREADS) best = None for d in XGB_GRID: bst = xgb.train({**XGB_PARAMS, "max_depth": d}, dtr, XGB_MAX_ROUNDS, evals=[(dva, "val")], early_stopping_rounds=XGB_EARLY_STOP, verbose_eval=False) if best is None or bst.best_score < best[0]: best = (bst.best_score, d, bst.best_iteration + 1) _, depth, rounds = best del dtr, dva, bst dfull = _qdm(labelled.slice(0, i_end), feats) return xgb.train({**XGB_PARAMS, "max_depth": depth}, dfull, rounds) def _xgb_predictions(ranked: pl.DataFrame, feats: list[str], eoms: list[date]) -> pl.DataFrame: """Re-estimated every calendar year Y with labels known by Jan Y (eom <= Dec Y-1).""" r = pl.col("ret_exc_lead1m") labelled = (ranked.filter(r.is_not_null()) .with_columns(((r.rank("average").over("eom") - 1) / (pl.len().over("eom") - 1) - 0.5).cast(pl.Float32).alias("y")) .sort("eom", "id")) eom_col = labelled["eom"] out = [] for Y in sorted({e.year for e in eoms}): cutoff, val_start = date(Y - 1, 12, 31), date(Y - VAL_YEARS, 1, 1) i_val = int(eom_col.search_sorted(val_start, side="left")) # first eom >= val_start i_end = int(eom_col.search_sorted(cutoff, side="right")) # past last eom <= cutoff # CTF-FIX: progress output for the job log print(f"[CTF-DEBUG] XGBoost refit for {Y}: {i_end:,} training rows", flush=True) model = _xgb_fit(labelled, i_val, i_end, feats) test = ranked.filter(pl.col("eom").is_in([e for e in eoms if e.year == Y])) if model is None or test.is_empty(): continue pred = model.predict(xgb.DMatrix(test.select(feats).to_numpy(), nthread=XGB_THREADS)) out.append(test.select("id", "eom").with_columns(pred=pl.Series(pred, dtype=pl.Float64))) return pl.concat(out) if out else pl.DataFrame(schema={"id": pl.Int64, "eom": pl.Date, "pred": pl.Float64}) # ---------------------------------------------------------------------------- # Model B: PCA factor risk model + Markowitz # ---------------------------------------------------------------------------- def _pca_markowitz(preds: pl.DataFrame, daily: pl.DataFrame) -> pl.DataFrame: dates = daily["date"].to_numpy() trading_days = np.unique(dates) out = [] for (eom,), grp in preds.sort("eom").group_by("eom", maintain_order=True): last_day = np.searchsorted(trading_days, np.datetime64(eom), side="right") days = trading_days[max(0, last_day - RISK_WINDOW):last_day] if len(days) < 2: continue start = np.searchsorted(dates, days[0], side="left") end = np.searchsorted(dates, days[-1], side="right") ids = grp["id"].to_numpy() hist = (daily.slice(start, end - start) .filter(pl.col("id").is_in(ids)) .with_columns(pl.col("ret_exc").clip(-0.5, 0.5))) day_idx = np.searchsorted(days, hist["date"].to_numpy()) pos = {int(i): j for j, i in enumerate(ids)} col_idx = np.array([pos[int(i)] for i in hist["id"].to_numpy()], dtype=np.int64) R = np.zeros((len(days), len(ids))) obs = np.zeros_like(R, dtype=bool) R[day_idx, col_idx] = hist["ret_exc"].to_numpy() obs[day_idx, col_idx] = True R = np.where(obs, R - R.sum(0) / np.maximum(obs.sum(0), 1), 0.0) # exponentially weighted PCA factor model: Sigma = B diag(f) B' + diag(d) T = len(days) wts = 0.5 ** (np.arange(T)[::-1] / RISK_HALF_LIFE) wts /= wts.sum() Rw = R * np.sqrt(wts)[:, None] U, s, Vt = np.linalg.svd(Rw, full_matrices=False) k = min(PCA_FACTORS, len(s)) B, f = Vt[:k].T, s[:k] ** 2 * 21 resid = Rw - (U[:, :k] * s[:k]) @ Vt[:k] d = (resid ** 2).sum(0) / np.maximum((obs * wts[:, None]).sum(0), 1e-12) enough = obs.sum(0) >= RISK_MIN_OBS med = np.median(d[enough]) if enough.any() else np.median(d) d = np.maximum(np.where(enough, d, med), 0.1 * med) * 21 if not np.all(np.isfinite(d)) or med <= 0: continue def vol(w): bw = B.T @ w return float(np.sqrt((f * bw ** 2).sum() + (d * w ** 2).sum())) mu = grp["pred"].to_numpy() mu = mu - mu.mean() if not np.any(mu): continue DinvB = B / d[:, None] w = mu / d - DinvB @ np.linalg.solve(np.diag(1 / f) + B.T @ DinvB, B.T @ (mu / d)) w = (1 - PCA_SHRINK) * w / vol(w) + PCA_SHRINK * mu / vol(mu) w = w * (VOL_TARGET / np.sqrt(12)) / vol(w) out.append(pl.DataFrame({"id": ids, "eom": [eom] * len(ids), "w": w})) # CTF-FIX: progress output for the job log if len(out) % 60 == 0: print(f"[CTF-DEBUG] Portfolio weights built through {eom} ({len(out)} months)", flush=True) return pl.concat(out) if out else pl.DataFrame(schema={"id": pl.Int64, "eom": pl.Date, "w": pl.Float64}) # ---------------------------------------------------------------------------- # Model A: Ridge + 1-factor mean-variance (port of Jiayu's ps3_model.R) # ---------------------------------------------------------------------------- def _ridge_fit(X, y, months, n_folds=5, n_lambda=100): """cv.glmnet(alpha = 0, standardize = FALSE) with 5 consecutive month folds; returns the lambda.min fit as a prediction function.""" n, p = X.shape xm, ym = X.mean(0), y.mean() lam_max = np.abs((X - xm).T @ (y - ym)).max() / n / 1e-3 lams = lam_max * np.logspace(0, -4, n_lambda) mid = np.unique(months, return_inverse=True)[1] fold = np.minimum((mid * n_folds) // (mid.max() + 1), n_folds - 1) err = np.zeros(n_lambda) for k in range(n_folds): tr, te = fold != k, fold == k if tr.sum() == 0 or te.sum() == 0: continue xmt, ymt = X[tr].mean(0), y[tr].mean() Xc = X[tr] - xmt ev, V = np.linalg.eigh(Xc.T @ Xc / tr.sum()) proj = V.T @ (Xc.T @ (y[tr] - ymt) / tr.sum()) coefs = proj[:, None] / (ev[:, None] + lams[None, :]) # p x n_lambda, eigenbasis preds = ymt + ((X[te] - xmt) @ V) @ coefs # all lambdas at once err += ((y[te][:, None] - preds) ** 2).sum(0) lam = lams[np.argmin(err)] Xc = X - xm b = np.linalg.solve(Xc.T @ Xc / n + lam * np.eye(p), Xc.T @ (y - ym) / n) return lambda Z: ym + (Z - xm) @ b def _model_year_a(eom: date) -> int: """Jiayu's prediction year yr covers eom in [Dec yr-1, Nov yr].""" return eom.year + (1 if eom.month == 12 else 0) def _ridge_predictions(z: pl.DataFrame, feats: list[str], eoms: list[date]) -> pl.DataFrame: labelled = z.filter(pl.col("ret_exc_lead1m").is_not_null()) out = [] for yr in sorted({_model_year_a(e) for e in eoms}): train = labelled.filter(pl.col("eom") < date(yr - 1, 12, 31)) test = z.filter(pl.col("eom").is_in([e for e in eoms if _model_year_a(e) == yr])) if train["eom"].n_unique() < 2 or test.is_empty(): continue predict = _ridge_fit(train.select(feats).to_numpy(), train["ret_exc_lead1m"].to_numpy(), train["eom"].to_numpy()) out.append(test.select("id", "eom").with_columns( pred=pl.Series(predict(test.select(feats).to_numpy())))) return pl.concat(out) if out else pl.DataFrame(schema={"id": pl.Int64, "eom": pl.Date, "pred": pl.Float64}) def _one_factor_mv(preds: pl.DataFrame, daily: pl.DataFrame) -> pl.DataFrame: """Previous calendar year's daily returns, equal-weighted market factor, dollar-neutral Sigma^-1 mu scaled to gross exposure 1.""" preds = preds.with_columns((pl.col("eom").dt.year() + (pl.col("eom").dt.month() == 12).cast(pl.Int32)).alias("yr")) out = [] for (yr,), grp in preds.group_by("yr", maintain_order=True): d = daily.filter(pl.col("date").is_between(date(yr - 1, 1, 1), date(yr - 1, 12, 31))) if d.is_empty(): continue mkt = d.group_by("date").agg(pl.col("ret_exc").mean().alias("mkt")) mkt_var = float(mkt["mkt"].var()) risk = (d.join(mkt, on="date").group_by("id") .agg(pl.len().alias("n"), (pl.cov("ret_exc", "mkt") / mkt_var).alias("beta"), pl.col("ret_exc").var().alias("var_r")) .with_columns(pl.max_horizontal(pl.col("var_r") - pl.col("beta") ** 2 * mkt_var, pl.lit(1e-6)).alias("resid_var")) .filter((pl.col("n") >= 60) & pl.col("beta").is_finite() & pl.col("resid_var").is_finite())) g = grp.join(risk.select("id", "beta", "resid_var"), on="id", how="inner") for (eom,), m in g.group_by("eom", maintain_order=True): mu, b, dinv = m["pred"].to_numpy(), m["beta"].to_numpy(), 1 / m["resid_var"].to_numpy() den = 1 + mkt_var * (b ** 2 * dinv).sum() inv_mu = dinv * mu - dinv * b * (mkt_var * (b * dinv * mu).sum() / den) inv_one = dinv - dinv * b * (mkt_var * (b * dinv).sum() / den) w = inv_mu - (inv_mu.sum() / inv_one.sum()) * inv_one if np.abs(w).sum() > 0: out.append(pl.DataFrame({"id": m["id"], "eom": m["eom"], "w": w / np.abs(w).sum()})) return pl.concat(out) if out else pl.DataFrame(schema={"id": pl.Int64, "eom": pl.Date, "w": pl.Float64}) # ---------------------------------------------------------------------------- # AB blend # ---------------------------------------------------------------------------- def _blend(wa: pl.DataFrame, wb: pl.DataFrame, labels: pl.DataFrame) -> pl.DataFrame: """w = 0.5 w_A / vol_A + 0.5 w_B / vol_B, where vol_X at eom t is the std of strategy X's realised monthly returns for portfolios formed before t.""" def scale(w): r = (w.join(labels, on=["id", "eom"], how="inner") .group_by("eom").agg((pl.col("w") * pl.col("ret_exc_lead1m")).sum().alias("r")) .sort("eom")) return r.with_columns(pl.col("r").shift(1).rolling_std(window_size=10_000, min_samples=BLEND_MIN_OBS) .alias("vol")).select("eom", "vol") sa, sb = scale(wa), scale(wb) a = wa.join(sa, on="eom").with_columns((0.5 * pl.col("w") / pl.col("vol")).alias("wa")) b = wb.join(sb, on="eom").with_columns((0.5 * pl.col("w") / pl.col("vol")).alias("wb")) both = (a.select("id", "eom", "wa").join(b.select("id", "eom", "wb"), on=["id", "eom"], how="full", coalesce=True) .with_columns(pl.col("wa").fill_null(0.0), pl.col("wb").fill_null(0.0))) valid = sa.filter(pl.col("vol").is_not_null()).join(sb.filter(pl.col("vol").is_not_null()), on="eom")["eom"] return (both.filter(pl.col("eom").is_in(valid)) .with_columns((pl.col("wa") + pl.col("wb")).alias("w")).select("id", "eom", "w")) # ---------------------------------------------------------------------------- # Entry point # ---------------------------------------------------------------------------- def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame: feats = [f for f in features["features"].tolist() if f in chars.columns] data = _prepare(chars, feats) test_eoms = data.filter(pl.col("ctff_test"))["eom"].unique().sort().to_list() if not test_eoms: raise ValueError("No ctff_test rows in input") daily = _prepare_daily(daily_ret, test_eoms[-1]) # CTF-FIX: progress output for the job log print(f"[CTF-DEBUG] Data prepared: {data.height:,} rows, {len(feats)} features, " f"{len(test_eoms)} test months, candidate {CANDIDATE}", flush=True) # AB needs some months before the test period to estimate strategy volatilities all_eoms = data["eom"].unique().sort().to_list() first = all_eoms.index(test_eoms[0]) need = all_eoms[max(0, first - BLEND_BURN_IN):] if CANDIDATE == "AB" else test_eoms weights = {} if CANDIDATE in ("B", "AB"): preds = _xgb_predictions(_rank_features(data, feats), feats, need) weights["B"] = _pca_markowitz(preds, daily) if CANDIDATE in ("A", "AB"): preds = _ridge_predictions(_zscore_features(data, feats), feats, need) weights["A"] = _one_factor_mv(preds, daily) if CANDIDATE == "AB": labels = data.filter(pl.col("ret_exc_lead1m").is_not_null()).select("id", "eom", "ret_exc_lead1m") w = _blend(weights["A"], weights["B"], labels) else: w = weights[CANDIDATE] w = (w.filter(pl.col("eom").is_in(test_eoms) & pl.col("w").is_finite()) .sort(["eom", "id"])) out = w.to_pandas() out["id"] = out["id"].astype("int64") out["eom"] = pd.to_datetime(out["eom"]) out["w"] = out["w"].astype("float64") # CTF-FIX: progress output for the job log print(f"[CTF-DEBUG] Returning {len(out):,} weights for {out['eom'].nunique()} months", flush=True) return out[["id", "eom", "w"]]