"""PortfolioLens — Phase 3 data-preparation & feature-engineering helpers. Turns the **committed snapshot** (``data/snapshot/``) into a leakage-aware, stationary feature matrix for modeling. Like ``poc_quality`` it is a reusable module the Quarto chapter imports, so the chapter stays thin. Everything here operates on cached data — **no network, no API keys**. Design (per the indicator-council steer and strategy.md §8): * **Align before you engineer.** A business-day master calendar; monthly/weekly macro are *lagged to their publication date* (``lag_release``) so a reference-period value is never used before it was actually known, then forward-filled as a step function. * **Two layers.** A *regime/exposure layer* (credit, curve, financial conditions, Sahm, VIX) gates gross exposure via a small, fixed, **equal-weighted diffusion vote** — no fitted weights. A thin *return-ranking layer* (volatility-scaled 12-1 momentum) ranks the tradable proxies. * **Stationary, comparable features.** Levels become rolling z-scores / percentile ranks / diffusion shares / spreads-and-ratios. * **Labels are forward-looking by construction** and kept strictly separate from features (``forward_risk_adjusted_return``) — the project's risk-adjusted target. Run ``python scripts/poc_prepare.py`` to (re)build and write ``data/prepared/feature_matrix.parquet`` for the Phase-4 hand-off. """ from __future__ import annotations import argparse import os import sys from pathlib import Path import numpy as np import pandas as pd # make sibling modules importable whether run standalone or imported by the chapter sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import poc_quality as q # noqa: E402 REPO_ROOT = Path(__file__).resolve().parents[1] SNAPSHOT_DIR = REPO_ROOT / "data" / "snapshot" PREPARED_DIR = REPO_ROOT / "data" / "prepared" # The buyable proxies we actually rank/score (^VIX is an index, not investable, so it is # a regime input only — excluded here). DEFAULT_TRADABLE = ["SPY", "XLE", "ITA", "DBA", "BDRY", "GLD", "UUP", "BTC-USD", "ETH-USD", "LMT", "COPPER"] # ~business days per calendar period, for rolling windows BDAYS_MONTH, BDAYS_QTR, BDAYS_YEAR = 21, 63, 252 Z_WINDOW = 756 # ~3y rolling window for z-scores (long enough to span a regime) # ───────────────────────────────────────────────────────────────────────────── # 3a/3d — calendar, release lag, alignment (CRISP-DM: clean / integrate / format) # ───────────────────────────────────────────────────────────────────────────── def master_calendar(panel: pd.DataFrame) -> pd.DatetimeIndex: """Business-day index spanning the panel (the modeling clock).""" return pd.bdate_range(panel.index.min(), panel.index.max()) # How long after its reference period each cadence is actually published. FRED stamps monthly # series at the FIRST of the reference month, but many monthly macro prints (durable-goods orders, # consumer sentiment, industrial production) are released mid-to-late the FOLLOWING month — so a # 2-month forward shift conservatively guarantees no value is used before it was published (erring # toward staleness, never look-ahead). True per-release ALFRED-vintage timestamps are the documented # Phase-4 upgrade; these uniform buffers are the PoC approximation. _RELEASE_LAG = {"monthly": dict(months=2), "quarterly": dict(months=4), "weekly": dict(days=7)} def lag_release(s: pd.Series, **offset) -> pd.Series: """Shift a series' timestamps forward to approximate its publication date. A value stamped by *reference period* is only knowable on the later *release date*; shifting the index forward prevents look-ahead. ``offset`` is passed to ``pd.DateOffset`` (e.g. ``months=1``). """ out = s.copy() out.index = s.index + pd.DateOffset(**offset) return out def align_features(panel: pd.DataFrame, manifest: pd.DataFrame, cal: pd.DatetimeIndex | None = None) -> pd.DataFrame: """Reindex every series onto the business-day master, lag-aware. Macro (monthly/weekly/quarterly) is lagged to publication then forward-filled as a step function (never interpolated, never filled across the release frontier). Price series are reindexed and lightly forward-filled across holidays. """ cal = master_calendar(panel) if cal is None else cal meta = manifest.set_index("filekey") out: dict[str, pd.Series] = {} for col in panel.columns: s = panel[col].dropna() if len(s) == 0: out[col] = pd.Series(np.nan, index=cal) continue freq = meta.loc[col, "frequency"] if col in meta.index else "business-day" if freq in _RELEASE_LAG: s = lag_release(s, **_RELEASE_LAG[freq]) # forward-fill from actual observation dates onto the business-day calendar s = s.reindex(cal.union(s.index)).sort_index().ffill().reindex(cal) out[col] = s return pd.DataFrame(out, index=cal) # ───────────────────────────────────────────────────────────────────────────── # stationary transforms # ───────────────────────────────────────────────────────────────────────────── def rolling_z(df: pd.DataFrame, window: int = Z_WINDOW, min_periods: int | None = None) -> pd.DataFrame: """Rolling z-score (value minus rolling mean, over rolling std).""" mp = min_periods or max(BDAYS_MONTH, window // 4) mu = df.rolling(window, min_periods=mp).mean() sd = df.rolling(window, min_periods=mp).std() return (df - mu) / sd.replace(0, np.nan) def rolling_rank(df: pd.DataFrame, window: int = Z_WINDOW, min_periods: int | None = None) -> pd.DataFrame: """Rolling percentile rank of the current value within its trailing window [0,1].""" mp = min_periods or max(BDAYS_MONTH, window // 4) def _pct(x: np.ndarray) -> float: return (np.sum(x <= x[-1]) - 1) / (len(x) - 1) if len(x) > 1 else np.nan return df.rolling(window, min_periods=mp).apply(_pct, raw=True) # ───────────────────────────────────────────────────────────────────────────── # 3c — construct: return-ranking layer (momentum) # ───────────────────────────────────────────────────────────────────────────── def momentum(prices: pd.DataFrame, lookback: int = BDAYS_YEAR, skip: int = BDAYS_MONTH, vol_scale: bool = True, vol_window: int = BDAYS_QTR) -> pd.DataFrame: """12-1 month momentum: return from t-lookback to t-skip (skip the last month to avoid short-term reversal), optionally volatility-scaled to tame momentum crashes.""" p = prices.sort_index() mom = p.shift(skip) / p.shift(lookback) - 1.0 if vol_scale: rets = np.log(p / p.shift(1)) vol = rets.rolling(vol_window).std() * np.sqrt(BDAYS_YEAR) mom = mom / vol.replace(0, np.nan) return mom def cross_sectional_rank(df: pd.DataFrame) -> pd.DataFrame: """Rank each row across columns into [0,1] — who leads *today* (1 = strongest).""" return df.rank(axis=1, pct=True) # ───────────────────────────────────────────────────────────────────────────── # 3c — construct: regime / exposure layer # ───────────────────────────────────────────────────────────────────────────── def curve_slope(panel: pd.DataFrame) -> pd.DataFrame: """Yield-curve slope features (negative = inverted = recession lead).""" out = pd.DataFrame(index=panel.index) if {"DGS10", "DGS2"} <= set(panel.columns): out["curve_10y2y"] = panel["DGS10"] - panel["DGS2"] if "T10Y3M" in panel.columns: out["curve_10y3m"] = panel["T10Y3M"] return out def copper_gold(panel: pd.DataFrame, window: int = Z_WINDOW) -> pd.DataFrame: """Dr. Copper / gold ratio — a growth & risk-appetite barometer (+ its z-score).""" out = pd.DataFrame(index=panel.index) if {"COPPER", "GLD"} <= set(panel.columns): ratio = (panel["COPPER"] / panel["GLD"]).rename("copper_gold") out["copper_gold"] = ratio out["copper_gold_z"] = rolling_z(ratio.to_frame(), window)["copper_gold"] return out def quadrant(panel: pd.DataFrame, window: int = BDAYS_QTR) -> pd.DataFrame: """Reflation/stagflation quadrant from 3-month changes in breakevens (T10YIE) and real yields (DFII10) — routes the conflict-sensitive sleeve (strategy.md §7).""" out = pd.DataFrame(index=panel.index) if {"T10YIE", "DFII10"} <= set(panel.columns): d_be = panel["T10YIE"].diff(window) d_rr = panel["DFII10"].diff(window) out["breakeven_chg_3m"] = d_be out["realrate_chg_3m"] = d_rr lab = pd.Series(np.nan, index=panel.index, dtype="object") lab[(d_be > 0) & (d_rr < 0)] = "reflation" # commodities/gold risk-on lab[(d_be > 0) & (d_rr >= 0)] = "overheat" lab[(d_be <= 0) & (d_rr < 0)] = "goldilocks" lab[(d_be <= 0) & (d_rr >= 0)] = "deflation/tightening" out["quadrant"] = lab return out def lei_proxy(panel: pd.DataFrame, window: int = 6 * BDAYS_MONTH) -> pd.DataFrame: """Free homemade diffusion-LEI: share of available leading components rising over ~6 months. Components: AWHAEMAN, ICSA (inverted), PERMIT, NEWORDER, UMCSENT, SPY, and the 10y-fed-funds term spread — 7 of the 10 Conference-Board LEI inputs.""" comps: dict[str, pd.Series] = {} if "AWHAEMAN" in panel.columns: comps["hours"] = panel["AWHAEMAN"] if "ICSA" in panel.columns: comps["claims_inv"] = -panel["ICSA"] # fewer claims = improving if "PERMIT" in panel.columns: comps["permits"] = panel["PERMIT"] if "NEWORDER" in panel.columns: comps["capex_orders"] = panel["NEWORDER"] if "UMCSENT" in panel.columns: comps["consumer_exp"] = panel["UMCSENT"] if "SPY" in panel.columns: comps["sp500"] = panel["SPY"] if {"DGS10", "DFF"} <= set(panel.columns): comps["term_spread"] = panel["DGS10"] - panel["DFF"] out = pd.DataFrame(index=panel.index) if comps: C = pd.DataFrame(comps) d = C.diff(window) # Keep NaN where a component has no value yet (e.g. AWHAEMAN before 2006) so the diffusion # averages only the AVAILABLE components, instead of counting "no data" as a falling vote # (which would bias the index low early in the sample and freeze lei_n_components at 7). rising = d.gt(0).astype(float).where(d.notna()) out["lei_diffusion"] = rising.mean(axis=1) # share rising among available, [0,1] out["lei_n_components"] = rising.notna().sum(axis=1) # how many components are live return out def regime_score(panel: pd.DataFrame, window: int = Z_WINDOW) -> pd.DataFrame: """Equal-weighted diffusion vote over robust stress signals → a gross-exposure dial. Each component casts a 0/1 risk-off vote; the score is the share firing (1 = max stress). No fitted weights — the council's anti-overfitting steer. The dial scales gross exposure from 1.0 (calm) down toward 0.5 (full stress); it is an *exposure* knob, not a buy/sell trigger. """ votes = pd.DataFrame(index=panel.index) if "BAA10Y" in panel.columns: # credit stress votes["credit"] = (rolling_z(panel[["BAA10Y"]], window)["BAA10Y"] > 0).astype(float) if "T10Y3M" in panel.columns: # curve inversion votes["curve"] = (panel["T10Y3M"] < 0).astype(float) if "NFCI" in panel.columns: # tight conditions votes["nfci"] = (panel["NFCI"] > 0).astype(float) if "SAHM" in panel.columns: # recession onset votes["sahm"] = (panel["SAHM"] >= 0.5).astype(float) if "VIX" in panel.columns: # elevated volatility votes["vix"] = (rolling_z(panel[["VIX"]], window)["VIX"] > 0).astype(float) out = votes.copy() score = votes.mean(axis=1) out["regime_score"] = score out["gross_exposure_dial"] = 1.0 - 0.5 * score return out # ───────────────────────────────────────────────────────────────────────────── # label (forward-looking — NEVER a feature) # ───────────────────────────────────────────────────────────────────────────── def forward_risk_adjusted_return(prices: pd.DataFrame, horizon: int = BDAYS_MONTH, vol_window: int = BDAYS_QTR) -> pd.DataFrame: """Forward h-day log return divided by trailing-vol*sqrt(h) — the risk-adjusted target. Uses ``shift(-horizon)`` so it is look-ahead *by design*; for training only.""" p = prices.sort_index() logp = np.log(p) fwd = logp.shift(-horizon) - logp daily_vol = np.log(p / p.shift(1)).rolling(vol_window).std() return fwd / (daily_vol * np.sqrt(horizon)).replace(0, np.nan) # ───────────────────────────────────────────────────────────────────────────── # assembly # ───────────────────────────────────────────────────────────────────────────── def build_feature_matrix(snapshot_dir: str | Path = SNAPSHOT_DIR, tradables: list[str] | None = None, horizons: tuple[int, ...] = (BDAYS_MONTH, BDAYS_QTR)): """Return ``(X, Y)`` — the aligned feature matrix and the forward-return labels.""" manifest = q.load_manifest(snapshot_dir) panel = q.load_panel(snapshot_dir) cal = master_calendar(panel) aligned = align_features(panel, manifest, cal) trad = [c for c in (tradables or DEFAULT_TRADABLE) if c in aligned.columns] feats: dict[str, pd.Series] = {} # return-ranking layer mom = momentum(aligned[trad]) rank = cross_sectional_rank(mom) for c in trad: feats[f"mom_{c}"] = mom[c] feats[f"momrank_{c}"] = rank[c] # regime / exposure layer reg = regime_score(aligned) feats["regime_score"] = reg["regime_score"] feats["gross_exposure_dial"] = reg["gross_exposure_dial"] cs = curve_slope(aligned) for col in cs.columns: feats[col] = cs[col] for c in ["BAA10Y", "NFCI", "ANFCI", "VIX"]: if c in aligned.columns: feats[f"{c}_z"] = rolling_z(aligned[[c]])[c] if "SAHM" in aligned.columns: feats["sahm_flag"] = (aligned["SAHM"] >= 0.5).astype(float) # quadrant + channel routing (compute each frame once) qd, cg, lp = quadrant(aligned), copper_gold(aligned), lei_proxy(aligned) for frame in (qd, cg, lp): for col in frame.columns: feats[col] = frame[col] X = pd.DataFrame(feats, index=cal) labels: dict[str, pd.Series] = {} for h in horizons: fra = forward_risk_adjusted_return(aligned[trad], horizon=h) for c in trad: labels[f"fwd_rar_{h}d_{c}"] = fra[c] Y = pd.DataFrame(labels, index=cal) return X, Y def run(snapshot_dir: str | Path = SNAPSHOT_DIR, prepared_dir: str | Path = PREPARED_DIR) -> pd.DataFrame: """Build the feature matrix and write it to ``data/prepared/`` for Phase 4.""" prepared_dir = Path(prepared_dir) prepared_dir.mkdir(parents=True, exist_ok=True) X, Y = build_feature_matrix(snapshot_dir) full = X.join(Y) full.index.name = "date" full.to_parquet(prepared_dir / "feature_matrix.parquet") print(f"Feature matrix: {full.shape[0]} business days x {X.shape[1]} features " f"+ {Y.shape[1]} labels -> {prepared_dir / 'feature_matrix.parquet'}") return full if __name__ == "__main__": ap = argparse.ArgumentParser(description="PortfolioLens Phase-3 feature build") ap.add_argument("--snapshot", default=str(SNAPSHOT_DIR)) ap.add_argument("--prepared", default=str(PREPARED_DIR)) args = ap.parse_args() run(args.snapshot, args.prepared)