"""
CTF submission -- MGT 924 Problem Set 3 (Abhay Kumar & Thomas Robichaud, Yale SOM)

Model: "Three-learner ensemble with a clustered Barra-style risk model"

  1. Characteristics are rank-transformed within each month to [-0.5, 0.5]
     (missing -> 0, the cross-sectional median).
  2. Three learners predict the cross-sectionally standardised next-month
     excess return, re-estimated every December on an expanding window:
        (a) ridge regression                      (linear)
        (b) ridge on random Fourier features      (smooth non-linear)
        (c) LightGBM                              (trees / interactions)
     The ridge penalties are chosen by a time-ordered hold-out (the most
     recent months of the training window). Predictions are z-scored within
     the month and combined with non-negative stacking weights estimated from
     the learners' own past out-of-sample forecasts (equal weights until 36
     months of such forecasts exist).
  3. Risk model: characteristics are grouped by hierarchical clustering of
     their correlation matrix; daily cross-sectional regressions of stock
     returns on 12 industry dummies and the cluster exposures give factor
     returns and residuals. The model's own forecast is added as a factor so
     that the optimiser cannot hide risk in the direction of the signal.
     Sigma = B Omega B' + D  (EWMA factor covariance, EWMA idiosyncratic var).
  4. Portfolio: dollar-neutral mean-variance weights (Woodbury solve), scaled
     to a constant ex-ante volatility. The amount of idiosyncratic-variance
     inflation (a ridge penalty on the optimiser) is chosen each December
     from the realised Sharpe ratio of past out-of-sample portfolios.

Every quantity used for the portfolio formed at month t depends only on data
dated t or earlier: targets enter training only once realised, the refit
calendar is anchored to the calendar (not to the end of the sample), and all
tuning uses past out-of-sample months.

Disclosure: the use of LightGBM on this dataset and the "tangency minus
minimum-variance" way of imposing dollar neutrality follow public leaderboard
submissions (P. Geertsema). The risk model follows Jensen, Kelly, Malamud and
Pedersen, "Machine Learning and the Implementable Efficient Frontier"; the
clustering follows Jensen, Kelly and Pedersen, "Is There a Replication Crisis
in Finance?". AI assistance was used to write the code.

"""
import os
from collections import deque

import numpy as np
import pandas as pd
from scipy.cluster.hierarchy import fcluster, linkage
from scipy.spatial.distance import squareform
from scipy.stats import rankdata

SEED = 42

CFG = dict(
    min_train=60,          # months of realised targets before the first fit
    val_months=60,         # hold-out length for the ridge penalties (capped at 1/3 of sample)
    lam_grid=np.logspace(-3, 3, 13),   # ridge penalty, in units of mean Gram eigenvalue
    n_rff=512,             # random frequencies (-> 1024 sin/cos features)
    rff_scales=(0.5, 1.0, 2.0),
    y_clip=3.0,            # winsorise the standardised target
    use_lgbm=True,
    lgbm=dict(n_estimators=400, learning_rate=0.03, num_leaves=63,
              colsample_bytree=0.5, subsample=0.7, subsample_freq=1,
              reg_lambda=10.0, max_bin=63),
    lgbm_min_child_frac=0.001,   # min_child_samples as a share of training rows
    stack=True,            # learn the ensemble weights from past out-of-sample forecasts
    stack_window=120,      # months of out-of-sample forecasts used for the weights
    stack_min_hist=36,
    n_clusters=80,         # style factors from clustered characteristics
    use_industry=True,     # 12 industry factors from SIC codes (else a single market factor)
    signal_factor=True,    # the model's own forecast is a risk factor (aligns alpha and risk)
    vol_calibration=False, # rescale by trailing realised/forecast vol (tested; not needed once the
                           # forecast is a risk factor, and it lowered the pre-test Sharpe ratio)
    bias_window=36,
    bias_min_hist=12,
    risk_window=60,        # months of daily factor returns kept
    fac_halflife=252,      # days
    idio_window=24,        # months
    idio_halflife=6,       # months
    idio_min_days=40,
    kappa_grid=(0.02, 0.1, 0.25, 1.0, 4.0),
    kappa_default=1.0,
    kappa_min_hist=36,     # months of out-of-sample history before kappa is tuned
    vol_target=0.10 / np.sqrt(12.0),
    sig_clip=3.0,
)


# --------------------------------------------------------------------------- helpers
def _as_bool(col):
    """ctff_test arrives as the strings '0'/'1' from WRDS; accept bool, numeric or string flags."""
    a = np.asarray(col)
    if a.dtype == bool:
        return a
    if np.issubdtype(a.dtype, np.number):
        return np.nan_to_num(a.astype(np.float64)) != 0
    return np.isin(np.char.lower(a.astype(str)), ["1", "true", "t", "1.0"])


def _zscore(x):
    x = np.asarray(x, dtype=np.float64)
    s = x.std()
    if not np.isfinite(s) or s < 1e-12:
        return np.zeros_like(x)
    return (x - x.mean()) / s


def _rank_block(a):
    """Within-month rank transform of an (n, K) block to [-0.5, 0.5]; NaN -> 0."""
    a = np.asarray(a, dtype=np.float64)
    out = np.zeros(a.shape, dtype=np.float32)
    cnt = np.sum(~np.isnan(a), axis=0)
    ok = cnt > 1
    if ok.any():
        r = rankdata(a[:, ok], axis=0, nan_policy="omit", method="average")
        r = (r - 0.5) / cnt[ok] - 0.5
        out[:, ok] = np.nan_to_num(r, nan=0.0)
    return out


def _ridge_path(G, c, G_val, c_val, lam_grid):
    """Pick the ridge penalty that minimises squared error on the hold-out Gram."""
    p = G.shape[0]
    scale = max(np.trace(G) / p, 1e-12)
    G_tr, c_tr = G - G_val, c - c_val
    ev, Q = np.linalg.eigh(G_tr)
    ev = np.clip(ev, 0.0, None)
    qc = Q.T @ c_tr
    tr_scale = max(np.trace(G_tr) / p, 1e-12)
    best, best_lam = np.inf, lam_grid[len(lam_grid) // 2]
    for lam in lam_grid:
        b = Q @ (qc / (ev + lam * tr_scale))
        loss = b @ G_val @ b - 2.0 * b @ c_val
        if loss < best - 1e-15:
            best, best_lam = loss, lam
    ev, Q = np.linalg.eigh(G)
    ev = np.clip(ev, 0.0, None)
    beta = Q @ ((Q.T @ c) / (ev + best_lam * scale))
    return beta, float(best_lam)


def _nnls_small(S, c):
    """Non-negative least squares for a handful of variables (enumerate the active sets)."""
    m = len(c)
    best, best_w = 0.0, np.zeros(m)
    for mask in range(1, 2 ** m):
        idx = [i for i in range(m) if mask >> i & 1]
        try:
            w = np.linalg.solve(S[np.ix_(idx, idx)] + 1e-10 * np.eye(len(idx)), c[idx])
        except np.linalg.LinAlgError:
            continue
        if (w < 0).any():
            continue
        obj = w @ S[np.ix_(idx, idx)] @ w - 2.0 * w @ c[idx]
        if obj < best - 1e-15:
            best, best_w = obj, np.zeros(m)
            best_w[idx] = w
    return best_w


def _clusters(G, n_clusters):
    """Hierarchical clustering of characteristics on 1 - correlation. Returns K x C averaging matrix."""
    K = G.shape[0]
    d = np.sqrt(np.clip(np.diag(G), 0.0, None))
    ok = d > 1e-8
    C = np.zeros((K, K))
    C[np.ix_(ok, ok)] = G[np.ix_(ok, ok)] / np.outer(d[ok], d[ok])
    np.fill_diagonal(C, 1.0)
    D = np.clip(1.0 - C, 0.0, 2.0)
    D = (D + D.T) / 2.0
    np.fill_diagonal(D, 0.0)
    nc = int(min(n_clusters, max(1, K // 2)))
    lab = fcluster(linkage(squareform(D, checks=False), method="average"), nc, criterion="maxclust")
    labs = np.unique(lab)
    A = np.zeros((K, len(labs)), dtype=np.float32)
    for j, l in enumerate(labs):
        m = (lab == l) & ok
        if m.any():
            A[m, j] = 1.0 / m.sum()
    return A[:, A.sum(0) > 0]


# Fama-French-12-style industry groups from SIC ranges (group 11 = everything else / missing)
_SIC_GROUPS = [
    (0, [(100, 999), (2000, 2399), (2700, 2749), (2770, 2799), (3100, 3199), (3940, 3989)]),      # non-durables
    (1, [(2500, 2519), (2590, 2599), (3630, 3659), (3710, 3711), (3714, 3714), (3716, 3716),
         (3750, 3751), (3792, 3792), (3900, 3939), (3990, 3999)]),                                # durables
    (2, [(2520, 2589), (2600, 2699), (2750, 2769), (3000, 3099), (3200, 3569), (3580, 3629),
         (3700, 3709), (3712, 3713), (3715, 3715), (3717, 3749), (3752, 3791), (3793, 3799),
         (3830, 3839), (3860, 3899)]),                                                            # manufacturing
    (3, [(1200, 1399), (2900, 2999)]),                                                            # energy
    (4, [(2800, 2829), (2840, 2899)]),                                                            # chemicals
    (5, [(3570, 3579), (3660, 3692), (3694, 3699), (3810, 3829), (7370, 7379)]),                  # business equipment
    (6, [(4800, 4899)]),                                                                          # telecoms
    (7, [(4900, 4949)]),                                                                          # utilities
    (8, [(5000, 5999), (7200, 7299), (7600, 7699)]),                                              # shops
    (9, [(2830, 2839), (3693, 3693), (3840, 3859), (8000, 8099)]),                                # health
    (10, [(6000, 6999)]),                                                                         # finance
]
_N_IND = 12


def _industry(sic):
    code = pd.to_numeric(pd.Series(np.asarray(sic)), errors="coerce").fillna(-1).to_numpy()
    out = np.full(len(code), _N_IND - 1, dtype=np.int8)
    for g, ranges in _SIC_GROUPS:
        for lo, hi in ranges:
            out[(code >= lo) & (code <= hi)] = g
    return out


def _exposures(Xb, A, ind=None):
    """[industry dummies (or a market intercept), standardised cluster scores] for one month."""
    S = Xb @ A
    S = S - S.mean(0)
    sd = S.std(0)
    sd[sd < 1e-8] = 1.0
    if ind is None:
        lead = np.ones((len(S), 1))
    else:
        lead = np.zeros((len(S), _N_IND))
        lead[np.arange(len(S)), ind] = 1.0
    return np.column_stack([lead, S / sd]).astype(np.float64)


def _mvo(mu, B, Om, d):
    """Dollar-neutral mean-variance direction under Sigma = B Om B' + diag(d) (Woodbury)."""
    di = 1.0 / d
    BtDi = B.T * di
    M = np.linalg.inv(Om) + BtDi @ B
    rhs = np.column_stack([mu, np.ones(len(mu))])
    y = di[:, None] * rhs
    sol = y - BtDi.T @ np.linalg.solve(M, BtDi @ rhs)
    a, b = sol[:, 0], sol[:, 1]
    sb = b.sum()
    return a - (a.sum() / sb) * b if abs(sb) > 1e-12 else a - a.mean()


def _exante_vol(w, B, Om, d):
    f = B.T @ w
    return float(np.sqrt(max(f @ Om @ f + np.sum(d * w * w), 0.0)))


def _scale(w, B, Om, d, target):
    v = _exante_vol(w, B, Om, d)
    return w * (target / v) if v > 1e-12 else np.zeros_like(w)


def _rank_ls(sig):
    r = rankdata(sig) - (len(sig) + 1) / 2.0
    s = np.abs(r).sum()
    return r / s * 2.0 if s > 0 else r


# --------------------------------------------------------------------------- pipeline
def run(chars, features, daily_ret, diagnostics=False, cfg=None, verbose=True):
    cfg = dict(CFG, **(cfg or {}))
    np.random.seed(SEED)
    rng = np.random.default_rng(SEED)

    def log(*a):
        if verbose:
            print(*a, flush=True)

    feats = [f for f in features["features"].tolist() if f in chars.columns]
    K = len(feats)

    # ---- panel layout (sorted by month, then id) ---------------------------------
    eom_dt = pd.to_datetime(chars["eom"])
    ym_row = (eom_dt.dt.year * 12 + eom_dt.dt.month - 1).to_numpy()
    order = np.lexsort((chars["id"].to_numpy(), ym_row))
    ym_row = ym_row[order]
    ids_row = chars["id"].to_numpy()[order]
    eom_out = chars["eom"].to_numpy()[order]
    ret = chars["ret_exc_lead1m"].to_numpy(dtype=np.float64)[order]
    N = len(order)
    ind_row = None
    if cfg["use_industry"] and "sic" in chars.columns:
        ind_row = _industry(chars["sic"].to_numpy()[order])
    n_lead = _N_IND if ind_row is not None else 1

    n_lead += int(bool(cfg["signal_factor"]))
    sig_cache = {}
    ens_w = None                      # stacking weights (None -> equal weights)
    live_prev = None                  # (month, forecasts actually used that month) for stacking
    stack_hist = deque(maxlen=cfg["stack_window"])

    def signals(t_):
        """z-scored forecasts of the current models for month t_ (cached per model vintage)."""
        if t_ not in sig_cache:
            Xb_ = X[starts[t_]:ends[t_]]
            b1_, b2_, booster_ = models
            p = {"ridge": _zscore(Xb_ @ b1_), "rff": _zscore(rff(Xb_) @ b2_)}
            if booster_ is not None:
                p["lgbm"] = _zscore(booster_.predict(Xb_))
            names_ = list(p)
            p["equal"] = _zscore(np.mean([p[k] for k in names_], axis=0))
            wts_ = ens_w if ens_w is not None else np.ones(len(names_)) / len(names_)
            p["ensemble"] = _zscore(sum(w_ * p[k] for w_, k in zip(wts_, names_)))
            sig_cache[t_] = p
        return sig_cache[t_]

    def expo(t_):
        sl_ = slice(starts[t_], ends[t_])
        B_ = _exposures(X[sl_], A, None if ind_row is None else ind_row[sl_])
        if cfg["signal_factor"]:
            B_ = np.column_stack([B_, signals(t_)["ensemble"]])
        return B_

    yms, starts = np.unique(ym_row, return_index=True)
    ends = np.append(starts[1:], N)
    T = len(yms)
    uid, gid = np.unique(ids_row, return_inverse=True)
    n_ids = len(uid)
    log(f"[ctf] rows={N:,} months={T} ids={n_ids:,} features={K}")

    # ---- features: within-month ranks, float32 ------------------------------------
    X = np.empty((N, K), dtype=np.float32)
    step = 64
    for c0 in range(0, K, step):
        cols = feats[c0:c0 + step]
        raw = chars[cols].to_numpy(dtype=np.float32)[order]
        for t in range(T):
            X[starts[t]:ends[t], c0:c0 + len(cols)] = _rank_block(raw[starts[t]:ends[t]])
        del raw
    log("[ctf] features ranked")

    # ---- target: cross-sectionally standardised, winsorised -----------------------
    y = np.zeros(N, dtype=np.float32)
    yok = np.zeros(N, dtype=bool)
    for t in range(T):
        r = ret[starts[t]:ends[t]]
        m = np.isfinite(r)
        if m.sum() > 10:
            z = np.zeros(len(r))
            z[m] = np.clip(_zscore(r[m]), -cfg["y_clip"], cfg["y_clip"])
            y[starts[t]:ends[t]] = z
            yok[starts[t]:ends[t]] = m

    # ---- daily returns -> per-month (stocks of previous month-end) x days ----------
    d_dt = pd.to_datetime(daily_ret["date"])
    d_ym = (d_dt.dt.year.to_numpy(dtype=np.int32) * 12 + d_dt.dt.month.to_numpy(dtype=np.int32) - 1)
    d_day = d_dt.to_numpy().astype("datetime64[D]").astype(np.int32)
    del d_dt
    d_id = daily_ret["id"].to_numpy()
    d_gid = np.clip(np.searchsorted(uid, d_id), 0, n_ids - 1).astype(np.int32)
    d_ret = daily_ret["ret_exc"].to_numpy(dtype=np.float32)
    keep = (uid[d_gid] == d_id) & np.isfinite(d_ret)
    del d_id
    o = np.flatnonzero(keep)
    o = o[np.argsort(d_ym[o], kind="stable")]
    d_ym, d_day, d_gid, d_ret = d_ym[o], d_day[o], d_gid[o], d_ret[o]
    del o, keep

    daily_R = {}   # month index j -> (n_{j-1}, days_j) returns during month j for stocks held at j-1
    pos = np.full(n_ids, -1, dtype=np.int64)
    for j in range(1, T):
        if yms[j] != yms[j - 1] + 1:
            continue
        lo, hi = np.searchsorted(d_ym, yms[j]), np.searchsorted(d_ym, yms[j], side="right")
        if hi - lo == 0:
            continue
        g_prev = gid[starts[j - 1]:ends[j - 1]]
        pos[g_prev] = np.arange(len(g_prev))
        p = pos[d_gid[lo:hi]]
        m = p >= 0
        days, di = np.unique(d_day[lo:hi][m], return_inverse=True)
        R = np.full((len(g_prev), len(days)), np.nan, dtype=np.float32)
        R[p[m], di] = d_ret[lo:hi][m]
        daily_R[j] = R
        pos[g_prev] = -1
    del d_ym, d_day, d_gid, d_ret
    log("[ctf] daily returns aligned")

    # ---- random Fourier features ---------------------------------------------------
    P2 = int(cfg["n_rff"])
    W = rng.standard_normal((K, P2)).astype(np.float32)
    sc = np.asarray(cfg["rff_scales"], dtype=np.float32)[np.arange(P2) % len(cfg["rff_scales"])]
    W *= sc / (0.2887 * np.sqrt(max(K, 1)))   # ranks have s.d. ~0.289 -> projection s.d. ~ scale

    def rff(Xb):
        proj = Xb @ W
        Z = np.concatenate([np.sin(proj), np.cos(proj)], axis=1)
        return Z - Z.mean(0)

    # ---- state ---------------------------------------------------------------------
    G1 = np.zeros((K, K)); c1 = np.zeros(K)
    G2 = np.zeros((2 * P2, 2 * P2)); c2 = np.zeros(2 * P2)
    GX = np.zeros((K, K))                      # Gram of all features seen (for clustering)
    hist1, hist2 = deque(maxlen=cfg["val_months"]), deque(maxlen=cfg["val_months"])
    n_tgt_months = 0
    models = None
    A = None
    fac, ssr, nob = {}, {}, {}
    kappas = list(cfg["kappa_grid"])
    kappa = cfg["kappa_default"]
    cand_hist = []                              # (month index, [return per kappa])
    w_out = np.zeros(N, dtype=np.float64)
    has_w = np.zeros(N, dtype=bool)
    diag_rows = []

    def risk_month(j):
        """Daily factor regressions for month j, exposures from month j-1."""
        R = daily_R.get(j)
        if R is None or A is None:
            fac.pop(j, None); ssr.pop(j, None); nob.pop(j, None)
            return
        B = expo(j - 1)
        full = ~np.isnan(R).any(1)
        if full.sum() < 3 * B.shape[1]:
            fac.pop(j, None); ssr.pop(j, None); nob.pop(j, None)
            return
        F = np.linalg.lstsq(B[full], R[full].astype(np.float64), rcond=None)[0]   # (k, days)
        E = R - B @ F
        g_prev = gid[starts[j - 1]:ends[j - 1]]
        s = np.zeros(n_ids, dtype=np.float32); n = np.zeros(n_ids, dtype=np.float32)
        s[g_prev] = np.nansum(E * E, axis=1)
        n[g_prev] = np.sum(~np.isnan(E), axis=1)
        fac[j], ssr[j], nob[j] = F.T, s, n

    def risk_at(t, g):
        """Monthly factor covariance and idiosyncratic variances known at the end of month t."""
        k = A.shape[1] + n_lead
        Fs = [fac[j] for j in range(max(1, t - cfg["risk_window"] + 1), t + 1) if j in fac]
        if Fs:
            Fd = np.vstack(Fs)
            wts = 0.5 ** (np.arange(len(Fd))[::-1] / cfg["fac_halflife"])
            Om = (Fd * wts[:, None]).T @ Fd / wts.sum() * 21.0
            if len(Fd) < 3 * k:
                Om = np.diag(np.diag(Om))
        else:
            Om = np.eye(k) * 1e-4
        Om = (Om + Om.T) / 2.0 + np.eye(k) * max(np.trace(Om) / k, 1e-10) * 1e-4
        num = np.zeros(len(g)); den = np.zeros(len(g)); cnt = np.zeros(len(g))
        for j in range(max(1, t - cfg["idio_window"] + 1), t + 1):
            if j in ssr:
                wj = 0.5 ** ((t - j) / cfg["idio_halflife"])
                num += wj * ssr[j][g]; den += wj * nob[j][g]; cnt += nob[j][g]
        okv = cnt >= cfg["idio_min_days"]
        d = np.full(len(g), np.nan)
        d[okv] = num[okv] / den[okv] * 21.0
        if okv.sum() >= 10:
            lo_, med, hi_ = np.percentile(d[okv], [10, 50, 99])
            d = np.clip(np.where(okv, d, med), max(lo_, 1e-8), hi_)
        else:
            d = np.full(len(g), 0.01)
        return Om, d

    # ---- walk forward ----------------------------------------------------------------
    for t in range(T):
        s_t, e_t = starts[t], ends[t]
        cal_month = int(yms[t] % 12) + 1

        # (1) month t-1's targets are realised by the end of month t
        if t >= 1:
            sl = slice(starts[t - 1], ends[t - 1])
            Xp, yp, mk = X[sl], y[sl].astype(np.float64), yok[sl]
            GX += (Xp.T @ Xp).astype(np.float64)
            if mk.sum() > 10:
                Xm = Xp[mk].astype(np.float64)
                g1, k1 = Xm.T @ Xm, Xm.T @ yp[mk]
                Zm = rff(Xp)[mk]
                g2 = (Zm.T @ Zm).astype(np.float32)
                k2 = (Zm.T @ yp[mk].astype(np.float32)).astype(np.float64)
                G1 += g1; c1 += k1; G2 += g2; c2 += k2
                hist1.append((g1.astype(np.float32), k1)); hist2.append((g2, k2))
                n_tgt_months += 1
                if live_prev is not None and live_prev[0] == t - 1:
                    S_ = live_prev[1][mk]
                    stack_hist.append((S_.T @ S_, S_.T @ yp[mk]))

        # (2) refit every December (and at the first month with enough history)
        refit = n_tgt_months >= cfg["min_train"] and (models is None or cal_month == 12)
        if refit:
            V = int(min(cfg["val_months"], n_tgt_months // 3))
            gv1 = sum(h[0].astype(np.float64) for h in list(hist1)[-V:]); cv1 = sum(h[1] for h in list(hist1)[-V:])
            gv2 = sum(h[0].astype(np.float64) for h in list(hist2)[-V:]); cv2 = sum(h[1] for h in list(hist2)[-V:])
            b1, lam1 = _ridge_path(G1, c1, gv1, cv1, cfg["lam_grid"])
            b2, lam2 = _ridge_path(G2, c2, gv2, cv2, cfg["lam_grid"])
            booster = None
            if cfg["use_lgbm"]:
                import lightgbm as lgb
                n_tr = starts[t]                       # rows of months <= t-1
                mk = yok[:n_tr]
                params = dict(cfg["lgbm"])
                params["min_child_samples"] = int(max(50, cfg["lgbm_min_child_frac"] * mk.sum()))
                reg = lgb.LGBMRegressor(objective="regression", random_state=SEED, deterministic=True,
                                        force_row_wise=True, n_jobs=int(min(32, os.cpu_count() or 1)),
                                        verbosity=-1, **params)
                if mk.all():
                    reg.fit(X[:n_tr], y[:n_tr])
                else:
                    reg.fit(X[:n_tr], y[:n_tr], sample_weight=mk.astype(np.float32))
                booster = reg.booster_
            models = (b1.astype(np.float32), b2.astype(np.float32), booster)
            if cfg["stack"] and len(stack_hist) >= cfg["stack_min_hist"]:
                wst = _nnls_small(sum(h[0] for h in stack_hist), sum(h[1] for h in stack_hist))
                ens_w = wst / wst.sum() if wst.sum() > 1e-12 else None
            sig_cache.clear()
            A = _clusters(GX, cfg["n_clusters"])
            for j in range(max(1, t - cfg["risk_window"] - 12), t + 1):
                risk_month(j)
            for j in [j for j in fac if j < t - cfg["risk_window"] - 12]:
                fac.pop(j, None); ssr.pop(j, None); nob.pop(j, None)


            past = np.array([r for (m, r) in cand_hist if m < t])
            if len(past) >= cfg["kappa_min_hist"]:
                sd = past.std(0)
                sr = np.where(sd > 1e-12, past.mean(0) / np.where(sd > 1e-12, sd, 1.0), -np.inf)
                kappa = kappas[int(np.argmax(sr))]
            log(f"[ctf] refit {yms[t] // 12}-{cal_month:02d}: train_months={n_tgt_months} "
                f"lam_ridge={lam1:.3g} lam_rff={lam2:.3g} kappa={kappa} "
                f"stack={'equal' if ens_w is None else np.round(ens_w, 2).tolist()} factors={A.shape[1] + n_lead}")
        elif A is not None:
            risk_month(t)

        if models is None or (e_t - s_t) < 20:
            continue

        # (3) predict, combine, build the portfolio
        for old in [m for m in sig_cache if m < t - cfg["risk_window"] - 14]:
            del sig_cache[old]
        preds = signals(t)
        ens = preds["ensemble"]
        comp = [k for k in preds if k not in ("ensemble", "equal")]
        live_prev = (t, np.column_stack([preds[k] for k in comp]))

        g = gid[s_t:e_t]
        B = expo(t)
        Om, d = risk_at(t, g)
        dbar = d.mean()
        clip = cfg["sig_clip"]
        mu = np.clip(ens, -clip, clip)
        rr = np.nan_to_num(ret[s_t:e_t], nan=0.0)

        # forecast-risk calibration: realised / forecast volatility of this strategy over the
        # trailing window (returns of months <= t-1 are known at t)
        past = np.array([r for (m, r) in cand_hist if m < t][-cfg["bias_window"]:])
        bias = np.ones(len(kappas))
        if cfg["vol_calibration"] and len(past) >= cfg["bias_min_hist"]:
            bias = np.clip(past.std(axis=0, ddof=1) / cfg["vol_target"], 0.5, 4.0)
        cand_w = [_scale(_mvo(mu, B, Om, d + kp * dbar), B, Om, d, cfg["vol_target"]) for kp in kappas]
        if np.isfinite(ret[s_t:e_t]).any():
            cand_hist.append((t, [float(w @ rr) for w in cand_w]))
        ki = kappas.index(kappa)
        w = cand_w[ki] / bias[ki]
        w_out[s_t:e_t] = w
        has_w[s_t:e_t] = True

        if diagnostics:
            row = dict(ym=int(yms[t]), n=int(e_t - s_t), kappa=float(kappa),
                       exante_vol=_exante_vol(w, B, Om, d) * float(bias[ki]), bias=float(bias[ki]),
                       gross=float(np.abs(w).sum()),
                       has_ret=bool(np.isfinite(ret[s_t:e_t]).any()), ret_final=float(w @ rr))
            yy = y[s_t:e_t]
            for name, p_ in preds.items():
                row[f"ic_{name}"] = float(np.corrcoef(p_, yy)[0, 1]) if yy.std() > 0 else np.nan
                row[f"ret_rankls_{name}"] = float(_rank_ls(p_) @ rr)
                if name != "ensemble":
                    wm = _scale(_mvo(np.clip(p_, -clip, clip), B, Om, d + kappa * dbar), B, Om, d, cfg["vol_target"])
                    row[f"ret_mvo_{name}"] = float(wm @ rr)
            for kp, wc in zip(kappas, cand_w):
                row[f"ret_kappa_{kp:g}"] = float(wc @ rr)
            for k_, w_ in zip(comp, ens_w if ens_w is not None else np.ones(len(comp)) / len(comp)):
                row[f"stackw_{k_}"] = float(w_)
            k = B.shape[1]
            wd = _scale(_mvo(mu, B[:, :1] * 0.0, np.eye(1), d + kappa * dbar), B, Om, d, cfg["vol_target"])
            row["ret_diag_only"] = float(wd @ rr)
            wu = _mvo(mu, B, Om, d + kappa * dbar)
            wu = wu / max(np.abs(wu).sum(), 1e-12)
            row["ret_mvo_unit_gross"] = float(wu @ rr)
            row["n_factors"] = int(k)
            zr = np.random.default_rng(SEED + t).standard_normal(len(mu))
            wr = _scale(zr - zr.mean(), B, Om, d, cfg["vol_target"])
            row["chk_random_ret"] = float(wr @ rr)
            row["chk_mvo_raw_ret"] = float(cand_w[ki] @ rr)
            diag_rows.append(row)

    out = pd.DataFrame({"id": ids_row.astype(np.int64), "eom": eom_out, "w": w_out})
    if "ctff_test" in chars.columns:
        test = _as_bool(chars["ctff_test"])[order]
        out = out[has_w | test]
    out = out.reset_index(drop=True)
    if diagnostics:
        dg = pd.DataFrame(diag_rows)
        if len(dg):
            dg["eom"] = pd.to_datetime(dict(year=dg["ym"] // 12, month=dg["ym"] % 12 + 1, day=1)) + pd.offsets.MonthEnd(0)
            if "ctff_test" in chars.columns:
                tm = pd.Series(_as_bool(chars["ctff_test"])[order]).groupby(ym_row).any()
                dg["test"] = dg["ym"].map(tm).to_numpy()
        return out, dg
    return out


def main(chars: pd.DataFrame, features: pd.DataFrame, daily_ret: pd.DataFrame) -> pd.DataFrame:
    """Competition entry point: returns a DataFrame with columns id, eom, w."""
    np.random.seed(SEED)
    return run(chars, features, daily_ret, diagnostics=False)


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")
    pf = main(chars, features, daily_ret)
    pf.to_csv("output.csv", index=False)