f/fundlab/decompose.py

"""Strategy decomposition: which pool components best explain a fund's returns.

For each actively-managed fund we:

  1. take its total-return (Adj Close) daily returns,
  2. run a greedy BIC forward-selection regression against a curated
     candidate set of pool components (broad sleeves, not single
     positions — the question is "what does the fund look like", not
     "what exact holdings"),
  3. report the chosen weights (betas), fit (R^2, alpha/t), and a
     rolling-window check for time-varying weights,

and hand a verdict to the report: *static sleeve mix* (weights stable),
*mostly stable with some timing*, or *time-varying / not a sleeve mix*.

The candidate sets are per-fund and come from the fund's stated strategy
plus its N-PORT holdings; selection within the set is data-driven.

NumPy only (no statsmodels in the venv). All regressions include an
intercept; alphas are annualized (x252); t-stats use the OLS covariance.
"""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np
import pandas as pd

from fundlab import pool as _pool

DATA = Path.home() / "prog/fin/stocks"
RESULTS = Path(__file__).parent / "decompose_results.json"

# ------------------------------------------------------------------ data
def adj_close(sym: str, root: Path = DATA) -> pd.Series | None:
    f = root / f"{sym}-history.csv"
    if not f.exists():
        return None
    try:
        s = (pd.read_csv(f, parse_dates=["Date"], index_col="Date")
             ["Adj Close"].dropna())
        s = s[~s.index.duplicated(keep="last")].sort_index()
    except Exception:
        return None
    if len(s) < 252:
        return None
    s.name = sym
    return s


def returns_panel(symbols: list[str], start: str | None = None) -> pd.DataFrame:
    """Daily simple returns, outer-joined on dates (NaN where a series
    has no data). Callers must handle missing values per regressor —
    candidate histories have different vintages (and some end early, e.g.
    finux stops in 2017), so a global inner join can be empty."""
    cols = []
    seen = set()
    for s in symbols:                      # de-dup: a symbol may appear
        if s in seen:                      # in both the benchmark mix and
            continue                       # the sleeve set
        seen.add(s)
        a = adj_close(s)
        if a is not None:
            cols.append(a)
    p = pd.concat(cols, axis=1, sort=False)
    if start is not None:
        p = p[p.index >= pd.Timestamp(start)]
    return p.pct_change()


# ------------------------------------------------------------------ stats
def ols(y: np.ndarray, X: np.ndarray) -> dict:
    """OLS of y on X (X must already include the intercept column)."""
    n, k = X.shape
    beta, *_ = np.linalg.lstsq(X, y, rcond=None)
    resid = y - X @ beta
    ssr = float(resid @ resid)
    dof = max(n - k, 1)
    s2 = ssr / dof
    XtX_inv = np.linalg.pinv(X.T @ X)
    se = np.sqrt(np.maximum(np.diag(XtX_inv) * s2, 0.0))
    t = np.where(se > 0, beta / np.where(se > 0, se, 1.0), 0.0)
    yss = float(((y - y.mean()) ** 2).sum())
    r2 = 1.0 - ssr / yss if yss > 0 else 0.0
    adj_r2 = 1.0 - (1.0 - r2) * (n - 1) / dof
    bic = n * np.log(ssr / n) + k * np.log(n)
    dw = float(((resid[1:] - resid[:-1]) ** 2).sum() / ssr) if ssr > 0 else 0.0
    return {"beta": beta, "se": se, "t": t, "r2": r2, "adj_r2": adj_r2,
            "bic": bic, "dw": dw, "ssr": ssr, "n": n, "resid": resid}


def forward_select(y: np.ndarray, cand: dict[str, np.ndarray],
                   max_k: int = 6, min_d_bic: float = 2.0,
                   y_ok: np.ndarray | None = None) -> tuple[list[str], dict, np.ndarray]:
    """Greedy forward selection on BIC, with per-model complete cases.

    ``cand`` maps symbol -> daily-return vector (aligned with y); NaNs
    mark dates where that series has no data. Each trial regression uses
    the rows where the fund, all chosen components and the candidate are
    all present, so differently-vintaged candidates stay comparable.
    A candidate is only added if it lowers BIC by at least ``min_d_bic``
    (~3:1 odds in its favour), so we stop before overfitting noise.

    Returns (chosen symbols, model on the final complete cases, the
    boolean mask of those rows).
    """
    if y_ok is None:
        y_ok = ~np.isnan(y)
    ones0 = np.ones(len(y))
    chosen: list[str] = []

    def fit(ch: list[str]) -> tuple[dict, np.ndarray]:
        ok = y_ok.copy()
        for c in ch:
            ok &= ~np.isnan(cand[c])
        X = np.column_stack([ones0[ok], *[cand[c][ok] for c in ch]])
        return ols(y[ok], X), ok

    cur, cur_ok = fit(chosen)
    for _ in range(max_k):
        best_j, best_b, best_m = None, cur["bic"], None
        for j in cand:
            if j in chosen:
                continue
            ok = cur_ok & ~np.isnan(cand[j])
            X = np.column_stack([ones0[ok], *[cand[c][ok] for c in chosen],
                                 cand[j][ok]])
            if X.shape[0] < 252 or np.linalg.matrix_rank(X) < X.shape[1]:
                continue
            m = ols(y[ok], X)
            # the new regressor is the last column; require it to be
            # individually significant (|t|>2) so we don't chase noise
            if abs(m["t"][-1]) < 2.0:
                continue
            if m["bic"] < best_b - 1e-9:
                best_j, best_b, best_m = j, m["bic"], m
        if best_j is None or best_b > cur["bic"] - min_d_bic:
            break
        chosen.append(best_j)
        cur, cur_ok = fit(chosen)
    return chosen, cur, cur_ok


def rolling_beta(y: np.ndarray, X: np.ndarray, w: int = 252,
                 step: int = 21, dates: pd.Index | None = None) -> pd.DataFrame:
    """Refit betas on rolling windows; returns a DataFrame indexed by the
    window end date, one column per regressor (no intercept column)."""
    n = len(y)
    rows = []
    idx = []
    i = w
    while i <= n:
        Xm = X[i - w:i]
        ym = y[i - w:i]
        if np.linalg.matrix_rank(Xm) < Xm.shape[1]:
            i += step
            continue
        b, *_ = np.linalg.lstsq(Xm, ym, rcond=None)
        rows.append(b[1:])
        idx.append(dates[i - 1] if dates is not None else i - 1)
        i += step
    cols = ["intercept"] + [f"x{k}" for k in range(X.shape[1] - 1)]
    return pd.DataFrame(rows, index=idx, columns=cols[1:])


# ------------------------------------------------------------------ funds
# per-fund candidate sets: broad pool sleeves consistent with the fund's
# stated strategy and N-PORT holdings. Pool labels for the report come
# from fundlab.pool.
# Candidate sets are curated to DISTINCT AXES (one representative per
# sleeve): near-duplicates (ivv/vt/vti, efa/vea, agg/vbtlx) would make the
# OLS knife-edge and produce offsetting betas. Discovery happens within
# this set; the set itself comes from the fund's stated strategy + N-PORT.
CANDIDATES: dict[str, list[str]] = {
    # QQQ 65% + SPY 29% + cash, options overlay
    "atesx": ["qqq", "ivv", "iwm", "shv", "bil"],
    # systematic short-duration IG credit + cash
    "atrfx": ["vstbx", "vicbx", "vweax", "ief", "ivv", "shv", "bil"],
    # market neutral: long equity, short credit
    "cvsix": ["ivv", "ive", "vweax", "vicbx", "tlt", "shv", "bil"],
    # US large-cap core plus
    "jlpsx": ["ivv", "ijt", "iwm", "efa", "vwo", "shv", "bil"],
    # global multi-asset
    "pmaix": ["ivv", "efa", "vwo", "vnq", "agg", "vweax", "djp", "gld",
              "fxe", "tlt", "shv", "bil"],
    # long/short mortgage & ABS
    "pmorx": ["vmbix", "vweax", "vicbx", "finux", "ief", "tlt", "shv",
              "bil"],
    # market-neutral style premia
    "qspnx": ["ive", "ivw", "iwm", "ijt", "ijs", "efa", "shv", "bil"],
    # low-volatility equity fund of funds
    "svarx": ["ivv", "ive", "ijk", "efa", "vwo", "shv", "bil", "agg"],
    # full credit-spectrum strategic income
    "cosix": ["agg", "vweax", "vmbix", "finux", "fghnx", "ief", "tlt",
              "shv", "bil"],
    # hedge fund (multi-strategy)
    "mbxix": ["ivv", "efa", "tlt", "ief", "vweax", "djp", "gld", "fxe",
              "shv", "bil"],
    # global macro, sovereign-centric
    "eagmx": ["tlt", "ief", "vweax", "finux", "fxe", "fxb", "fxy", "djp",
              "gld", "shv", "bil"],
    # Leuthold Core: fund of ETFs (no return history - holdings only)
    "lcorx": ["ivv", "efa", "vwo", "agg", "djp", "shv", "bil"],
    # US dividend-growth equity
    "lamhx": ["ivv", "ive", "ijk", "ijt", "iwm", "efa", "shv", "bil"],
}
# share classes: same underlying fund
ALIAS = {"pmfkx": "pmaix", "lcrix": "lcorx", "egrsx": "eagmx"}
FULL_WINDOW = "1990-01-01"   # effectively all history
RECENT_WINDOW = "2021-01-01"  # last ~5 years


def _drift(rb: pd.DataFrame, full_beta: np.ndarray) -> dict[str, float]:
    """Mean |rolling beta - full-sample beta| per component, as a fraction
    of max(|full beta|, 0.25) — 0 means perfectly stable weights."""
    out = {}
    for k, col in enumerate(rb.columns):
        b = full_beta[k + 1]
        d = float((rb[col] - b).abs().mean())
        out[col] = d / max(abs(b), 0.25)
    return out


def decompose(fund: str, start: str = FULL_WINDOW,
              candidates: dict[str, list[str]] | None = None) -> dict:
    cand_syms = (candidates or CANDIDATES).get(fund, [])
    r = returns_panel([fund] + cand_syms, start=start)
    if fund not in r.columns:
        return {"fund": fund, "error": "no return history in the data set",
                "verdict": "no return history", "start": start}
    y = r[fund].to_numpy()
    dates = r.index
    cand = {s: r[s].to_numpy() for s in r.columns if s != fund}
    y_ok = ~np.isnan(y)
    if y_ok.sum() < 252:
        return {"fund": fund, "error": "insufficient return history",
                "verdict": "insufficient return history", "start": start}
    chosen, m, ok = forward_select(y, cand, y_ok=y_ok)
    d0 = dates[ok]
    ann_alpha = float(m["beta"][0]) * 252
    te = float(np.std(m["resid"]) * np.sqrt(252))
    port_ann = float(np.mean(y[ok]) * 252)
    out = {
        "fund": fund,
        "start": str(d0[0].date()),
        "end": str(d0[-1].date()),
        "n_obs": int(m["n"]),
        "components": [
            {"sym": s, "label": _label(s), "beta": float(m["beta"][k + 1]),
             "t": float(m["t"][k + 1])}
            for k, s in enumerate(chosen)],
        "alpha_ann": ann_alpha,
        "alpha_t": float(m["t"][0]),
        "r2": float(m["r2"]),
        "adj_r2": float(m["adj_r2"]),
        "tracking_err_ann": te,
        "fund_ann_return": port_ann,
        "dw": float(m["dw"]),
        "weights_sum": float(m["beta"][1:].sum()),
    }
    # rolling stability of the chosen model (on the same complete cases)
    if chosen:
        X = np.column_stack([np.ones(int(ok.sum())),
                             *[r[s].to_numpy()[ok] for s in chosen]])
        rb = rolling_beta(y[ok], X, dates=d0)
        drift = _drift(rb, m["beta"])
        out["rolling"] = {
            "window": 252,
            "drift": {s: float(v) for s, v in drift.items()},
            "max_drift": float(max(drift.values())) if drift else 0.0,
            "n_windows": int(len(rb)),
        }
    out["verdict"] = _verdict(out)
    return out


# per-fund cross-check notes: what the N-PORT holdings + strategy say, and
# how the returns decomposition should be read against that
NOTES: dict[str, str] = {
    "atesx": ("Holdings (May 2026): QQQ 65% + SPY 29% + MMF 0.6%, with 4.9% "
              "'other assets in excess of liabilities' — an options overlay. "
              "But the rolling beta to those SAME holdings stays 0.13–0.89 "
              "(median 0.30, never above 1): the 'risk managed' in the name is "
              "real — a systematic equity de-risking overlay. Decomposition: one "
              "TACTICAL US-equity sleeve, not a static mix."),
    "atrfx": ("Systematic alpha over short-duration IG credit + cash. Returns are "
              "dominated by idiosyncratic credit/derivatives P&L (R² ≤ 0.22 vs bond "
              "sleeves) and the best-fit weights are knife-edge. Read as: "
              "cash-like carry + systematic alpha, no meaningful static sleeve."),
    "cvsix": ("Market neutral (long US equity, short credit). Full sample (since "
              "1990) is unexplainable — the strategy has changed over 36 years; the "
              "last 5 years show a small net equity/credit tilt (ivv +0.14, "
              "vweax +0.06) explaining 74%. The rest is spread/option alpha "
              "(full-sample annualized alpha +5.5%, t=6.7)."),
    "jlpsx": ("US large-cap core plus: essentially 1.04x the S&P 500 "
              "(R² 0.96 over 5y, stable). The 'plus' is small optionality "
              "(tiny ijt/vwo tilts in the 5y fit). The cleanest fund on the list."),
    "pmaix": ("Global multi-asset fund of funds (N-PORT: 99.5% in unaffiliated "
              "underlying funds/loans). Returns decompose into high-yield credit "
              "(vweax +0.62), intl equity (efa +0.23), commodities (+0.05), "
              "bonds (−0.15): R² 0.68, stable weights, alpha +3.5%/yr (t=3.3). "
              "The sleeves show through the underlying funds."),
    "pmorx": ("Long/short mortgage & ABS — returns mostly idiosyncratic (R² 0.10). "
              "5y direction is long MBS (vmbix +0.31) / short intermediate rates "
              "(ief −0.39), consistent with a carry/relative-value mortgage "
              "strategy. Not a static sleeve."),
    "qspnx": ("Market-neutral style premia: no static sleeve explains returns "
              "(R² 0.18 5y). Alpha vs a cash-like benchmark: +12.8%/yr full sample "
              "(t=4.0). Decomposition = pure factor harvesting (value/size/style "
              "tilts in both books); the 'exposures' in the table are residuals, "
              "not sleeves."),
    "svarx": ("Low-volatility equity fund of funds. Small but positive net market "
              "(efa +0.07, agg +0.10; R² 0.24, low drift). The edge is in "
              "volatility selection, not the mix: +5.2%/yr alpha over that small "
              "sleeve (t=5.1)."),
    "cosix": ("Strategic income across the credit spectrum. Full sample (since "
              "1990) unexplainable — vintage; the last 5 years are the honest "
              "current mix: high-yield +0.30, MBS +0.29, IG core +0.18 (R² 0.86)."),
    "mbxix": ("Multi-strategy hedge fund: 53% explained over a decade "
              "(ivv +0.39, ief −0.67, fxe −0.28, tlt +0.17, djp +0.08) — equity "
              "long, duration short, FX/commodity tilts, large active residual. "
              "Caveat: newest N-PORT on file is Sep 2024 — the fund may have "
              "changed strategy or stopped filing."),
    "eagmx": ("Global macro (sovereign-centric): nothing explains returns in the "
              "full or 5y window (R² ≤ 0.05) — textbook macro, positions are "
              "tactical and asset-agnostic. The whole story is +5.1%/yr (t=8.0) "
              "over a flat benchmark."),
    "lcorx": ("NEW share classes (trading since Jul 2026) — no return history to "
              "regress. Holdings (Dec 2025 N-PORT): 91.7% ETFs + 8.4% money "
              "market; the strategy is Leuthold core multi-asset via ETFs. "
              "Re-run the decomposition after a year of NAV accumulates."),
    "lamhx": ("Dividend growth: R² 0.95; S&P 500 + value/mid tilt "
              "(ivv +0.62, ive +0.26, ijk +0.20, iwm −0.14 over 5y), stable "
              "weights. Closest to a passive fund with an overlay on this list."),
}


def _label(sym: str) -> str:
    try:
        for e in _pool.load_pool():
            if e.symbol == sym:
                return e.label
    except Exception:
        pass
    return sym


def _verdict(d: dict) -> str:
    if d.get("error"):
        return d["error"]
    md = d.get("rolling", {}).get("max_drift", 0.0)
    r2 = d["r2"]
    if r2 >= 0.95 and md < 0.5:
        v = "static sleeve mix (weights stable, ~fully explained)"
    elif r2 >= 0.85:
        v = "mostly stable sleeve mix" + ("" if md < 1.0 else " with drifting weights")
    elif r2 >= 0.6:
        v = "partially explainable — material active/timing residual"
    else:
        v = "not a static sleeve mix — returns driven by active decisions"
    if d.get("components") is None or not d.get("components"):
        v = "no component explains returns (market neutral or cash-like)"
    return v


def run_all() -> dict[str, dict]:
    out = {}
    done = set()
    for fund in CANDIDATES:
        if fund in ALIAS or fund in done:   # share class: parent covers it
            continue
        full = decompose(fund)
        full["recent"] = decompose(fund, start=RECENT_WINDOW)
        full["note"] = NOTES.get(fund, "")
        # a much better 5y fit than the full sample means the strategy
        # itself changed — say so in the verdict
        rec = full["recent"]
        if ("r2" in rec and "r2" in full and rec["r2"] > full["r2"] + 0.2):
            full["verdict"] += (f"  [strategy evolved — 5y R² = "
                                f"{rec['r2']:.2f}]")
        out[fund] = full
        done.add(fund)
        print(f"{fund}: {full['verdict'][:90]}", flush=True)
    # share classes reference their parent fund
    for cls, parent in ALIAS.items():
        if parent in out:
            d = dict(out[parent])
            d["fund"] = cls
            d["note"] = (f"share class of {parent} — identical report. "
                         + d.get("note", ""))
            out[cls] = d
    RESULTS.write_text(json.dumps(out, indent=1))
    return out


if __name__ == "__main__":
    run_all()