"""Frozen, paired residual-volatility normalization screen. Offline only. See docs/research/residual-normalization.md. A failed screen terminates the experiment; a passing screen allows the separately specified daily book replay. """ from __future__ import annotations import argparse import contextlib import hashlib import importlib.metadata import json import math import sqlite3 import statistics import sys from collections import Counter from datetime import date, datetime, timezone from pathlib import Path import numpy as np ROOT = Path(__file__).resolve().parents[1] if str(ROOT) not in sys.path: sys.path.insert(0, str(ROOT)) from app.services.backtest_service import _quintile_spread, _spearman # noqa: E402 from app.services.momentum_service import compute_residual_12_1_momentum # noqa: E402 START = date(2016, 7, 1) ERA_SPLIT = date(2021, 1, 1) LOOKBACK = 252 SKIP = 21 PAIRS = LOOKBACK - SKIP HORIZON = 30 MIN_CROSS = 20 MIN_WINDOWS = 12 IC_FLOOR = 0.03 SEED = 20260908 REPLICATES = 5000 IC_BLOCK = 4 ARMS = ("unrestricted_control", "paired_control", "normalized") PROTOCOL = ROOT / "docs/research/residual-normalization.md" def sha256(path: Path) -> str: with path.open("rb") as source: return hashlib.file_digest(source, "sha256").hexdigest() def git_blob_id(path: Path) -> str: """The id ``git hash-object`` gives this file. Line endings are normalized first, as git does on commit, so a CRLF checkout and an LF checkout of the same source record the same id — one that can be looked up in history.""" data = path.read_bytes().replace(b"\r\n", b"\n") return hashlib.sha1(b"blob %d\0" % len(data) + data).hexdigest() def read_snapshot(path: Path) -> tuple[dict, dict, dict]: """Fail closed on incomplete data; never open a live application database.""" path = path.resolve(strict=True) manifest = json.loads(Path(str(path) + ".manifest.json").read_text(encoding="utf-8")) if manifest.get("complete") is not True: raise ValueError("Snapshot completion manifest is not complete") if any(Path(str(path) + suffix).exists() for suffix in ("-wal", "-journal")): raise ValueError("Snapshot has a journal; use a completed, quiescent copy") # sqlite3's own context manager only ends the transaction; closing releases # the file handle, which on Windows otherwise pins the snapshot for the run. with contextlib.closing(sqlite3.connect(path.as_uri() + "?mode=ro", uri=True)) as con: con.execute("PRAGMA query_only=ON") counts = { key: con.execute(f"SELECT COUNT(*) FROM {table}").fetchone()[0] for key, table in ( ("ticker_count", "tickers"), ("ohlcv_row_count", "ohlcv_records"), ("rank_only_count", "research_rank_only"), ) } if any(manifest.get(k) != v for k, v in counts.items()): raise ValueError(f"Snapshot completion counts do not match: {counts}") tickers = con.execute( "SELECT id, symbol FROM tickers WHERE symbol NOT IN " "(SELECT symbol FROM research_rank_only) ORDER BY symbol" ).fetchall() benchmark = { date.fromisoformat(d[:10]): float(c) for d, c in con.execute( "SELECT date, close FROM benchmark_prices WHERE symbol='SPY' ORDER BY date" ) } if not benchmark or any(not math.isfinite(c) or c <= 0 for c in benchmark.values()): raise ValueError("Missing or invalid SPY history") prices = {} for ticker_id, symbol in tickers: rows = con.execute( "SELECT date, open, high, low, close, volume FROM ohlcv_records " "WHERE ticker_id=? ORDER BY date", (ticker_id,), ).fetchall() dates = [date.fromisoformat(r[0][:10]) for r in rows] if len(dates) != len(set(dates)): raise ValueError(f"Duplicate stock dates: {symbol}") # A missing price must stop the run: read as 0 it would put a -100% # return into the unrestricted control while the paired arms drop it. if any(r[k] is None for r in rows for k in range(1, 5)): raise ValueError(f"Missing OHLC price: {symbol}") columns = ( ([d.toordinal() for d in dates],) + tuple([float(r[k]) for r in rows] for k in range(1, 5)) + ([float(r[5] or 0) for r in rows],) ) prices[symbol] = columns return prices, benchmark, {**manifest, "verified_counts": counts} def signal_at( dates: list[date], closes: list[float], i: int, benchmark: dict[date, float], benchmark_index: dict[date, int], ) -> dict: """Production numerator, paired with a strictly aligned sample denominator.""" empty = {"unrestricted_control": None, "paired_control": None, "normalized": None, "residual_sd": None, "reason": "short_history"} if i < LOOKBACK: return empty window_dates = dates[i - LOOKBACK:i + 1] window_closes = closes[i - LOOKBACK:i + 1] raw = compute_residual_12_1_momentum(window_dates, window_closes, benchmark) if raw is not None and not math.isfinite(raw): raise ValueError("Non-finite production signal") empty["unrestricted_control"] = raw formation_dates = window_dates[:PAIRS + 1] formation_closes = window_closes[:PAIRS + 1] indices = [benchmark_index.get(d) for d in formation_dates] if any(x is None for x in indices) or any( b != a + 1 for a, b in zip(indices, indices[1:]) ): return {**empty, "reason": "unaligned_formation"} if any(not math.isfinite(c) or c <= 0 for c in formation_closes): return {**empty, "reason": "invalid_formation_price"} market = [benchmark[d] for d in formation_dates] stock_returns = [b / a - 1 for a, b in zip(formation_closes, formation_closes[1:])] market_returns = [b / a - 1 for a, b in zip(market, market[1:])] mean_s, mean_m = statistics.fmean(stock_returns), statistics.fmean(market_returns) var_m = sum((r - mean_m) ** 2 for r in market_returns) if var_m <= 0 or raw is None: return {**empty, "reason": "unavailable_beta"} beta = sum((s - mean_s) * (m - mean_m) for s, m in zip(stock_returns, market_returns)) / var_m residuals = [s - beta * m for s, m in zip(stock_returns, market_returns)] # Independent arithmetic verifies that the denominator uses the numerator's # exact return window. It must never silently redefine production momentum. if not math.isclose(sum(residuals), raw, rel_tol=1e-10, abs_tol=1e-12): raise AssertionError("Production numerator parity failed") sd = statistics.stdev(residuals) # Treat machine-precision residuals of a perfect market tracker as zero, # not as infinite signal strength. This is a numerical degeneracy check, # not a tunable volatility floor. numerical_zero = 64 * np.finfo(float).eps * max(1.0, max(map(abs, stock_returns))) if sd <= numerical_zero: return {**empty, "reason": "zero_residual_variance"} normalized = raw / sd if not math.isfinite(normalized): raise ValueError("Non-finite normalized signal") return {"unrestricted_control": raw, "paired_control": raw, "normalized": normalized, "residual_sd": sd, "reason": "eligible"} def weekly_indices(calendar: list[date]) -> list[int]: last = {} for i, day in enumerate(calendar): last[day.isocalendar()[:2]] = i return sorted(i for i in last.values() if calendar[i] >= START and i + HORIZON < len(calendar)) def nonoverlapping_indices(indices: list[int], horizon: int = HORIZON) -> list[int]: kept = [] for i in sorted(indices): if not kept or i >= kept[-1] + horizon: kept.append(i) return kept def interval(values) -> list[float] | None: if len(values) == 0: return None return [float(v) for v in np.quantile(values, [0.025, 0.975])] def bootstrap_means(matrix: np.ndarray) -> dict: """Shared moving-block indices preserve both the dates and arm pairing.""" n = len(matrix) if n < IC_BLOCK: return {"reason": "fewer_windows_than_block", "arms": {}, "paired_delta_ci95": None} rng = np.random.default_rng(SEED) starts = rng.integers(0, n - IC_BLOCK + 1, (REPLICATES, math.ceil(n / IC_BLOCK))) indices = (starts[:, :, None] + np.arange(IC_BLOCK)).reshape(REPLICATES, -1)[:, :n] means = matrix[indices].mean(axis=1) return { "arms": {arm: interval(means[:, i]) for i, arm in enumerate(ARMS)}, "paired_delta_ci95": interval(means[:, 2] - means[:, 1]), } def summarize(windows: list[dict]) -> dict: if not windows: return {"windows": 0, "arms": {}, "paired_mean_ic_delta": None} matrix = np.array([[w["arms"][arm]["ic"] for arm in ARMS] for w in windows]) uncertainty = bootstrap_means(matrix) arms = {} for i, arm in enumerate(ARMS): values = matrix[:, i] se = float(values.std(ddof=1) / math.sqrt(len(values))) if len(values) > 1 else None arms[arm] = { "mean_ic": float(values.mean()), "ic_t_stat_descriptive": float(values.mean() / se) if se and se > 0 else None, "ic_ci95": uncertainty["arms"].get(arm), "positive_window_fraction": float((values > 0).mean()), "avg_cross_section": statistics.fmean(w["arms"][arm]["n"] for w in windows), "mean_quintile_spread": statistics.fmean( w["arms"][arm]["quintile_spread"] for w in windows ), "reliable_window_count": len(windows) >= MIN_WINDOWS, } return {"windows": len(windows), "first_asof": windows[0]["asof"], "last_asof": windows[-1]["asof"], "last_forward_end": windows[-1]["forward_end"], "arms": arms, "paired_mean_ic_delta": float((matrix[:, 2] - matrix[:, 1]).mean()), "paired_ic_delta_ci95": uncertainty["paired_delta_ci95"]} def screen_pass(summaries: dict) -> dict: full = summaries["full"] means = {era: summary.get("arms", {}).get("normalized", {}).get("mean_ic") for era, summary in summaries.items()} checks = { "mean_ic_ge_0_03": means["full"] is not None and means["full"] >= IC_FLOOR, "at_least_12_windows": full["windows"] >= MIN_WINDOWS, "pre_2021_positive": means["pre_2021"] is not None and means["pre_2021"] > 0, "post_2021_positive": means["post_2021"] is not None and means["post_2021"] > 0, } return {"passed": all(checks.values()), "checks": checks, "failed_checks": [k for k, v in checks.items() if not v]} def run_screen(prices: dict, benchmark: dict) -> dict: calendar = sorted(benchmark) bench_index = {d: i for i, d in enumerate(calendar)} weekly = weekly_indices(calendar) observations = {i: [] for i in weekly} reasons = Counter() for count, (symbol, columns) in enumerate(sorted(prices.items()), 1): dates = [date.fromordinal(int(o)) for o in columns[0]] closes = columns[4] index = {d: i for i, d in enumerate(dates)} for asof in weekly: i, forward = index.get(calendar[asof]), index.get(calendar[asof + HORIZON]) if i is None or forward is None: reasons["missing_forward_or_asof"] += 1 continue if any(not math.isfinite(closes[j]) or closes[j] <= 0 for j in (i, forward)): reasons["invalid_forward_or_asof"] += 1 continue signal = signal_at(dates, closes, i, benchmark, bench_index) reasons[signal["reason"]] += 1 observations[asof].append({"symbol": symbol, **signal, "forward": closes[forward] / closes[i] - 1}) if count % 50 == 0 or count == len(prices): print(f"signal screen: {count}/{len(prices)} symbols", flush=True) usable = [i for i in weekly if sum(o["normalized"] is not None for o in observations[i]) >= MIN_CROSS] kept = nonoverlapping_indices(usable) windows = [] for i in kept: rows = observations[i] arms = {} for arm in ARMS: pairs = [(o[arm], o["forward"]) for o in rows if o[arm] is not None] ic = _spearman([p[0] for p in pairs], [p[1] for p in pairs]) if ic is None: raise ValueError(f"Degenerate {arm} cross-section on {calendar[i]}") arms[arm] = {"n": len(pairs), "ic": ic, "quintile_spread": _quintile_spread(pairs)} # Raised, not asserted: python -O must not strip the pairing guarantee. if arms["paired_control"]["n"] != arms["normalized"]["n"]: raise AssertionError(f"Paired arms scored different stocks on {calendar[i]}") windows.append({"asof": calendar[i].isoformat(), "forward_end": calendar[i + HORIZON].isoformat(), "arms": arms}) summaries = { "full": summarize(windows), "pre_2021": summarize([w for w in windows if w["forward_end"] < ERA_SPLIT.isoformat()]), "post_2021": summarize([w for w in windows if w["asof"] >= ERA_SPLIT.isoformat()]), } return {"summaries": summaries, "verdict": screen_pass(summaries), "windows": windows, "coverage": {"registry_symbols": len(prices), "weekly_asofs_considered": len(weekly), "usable_weekly_asofs": len(usable), "all_weekly_cell_reasons": dict(reasons), "retained_unrestricted_cells": sum(w["arms"][ARMS[0]]["n"] for w in windows), "retained_paired_cells": sum(w["arms"][ARMS[1]]["n"] for w in windows)}, "bootstrap": {"seed": SEED, "replicates": REPLICATES, "window_block": IC_BLOCK, "generator": "NumPy PCG64", "interval": "percentile 95%"}} def render_report(report: dict) -> str: screen = report["screen"] lines = ["# Residual normalization — signal screen", "", f"Generated: {report['generated_at']}", "", f"**Verdict: {'PASS' if screen['verdict']['passed'] else 'FAIL'}**", "", "Previously inspected historical data; no production change or new out-of-sample claim.", "", "| Era | Arm | Windows | Avg stocks | Mean IC | IC 95% interval | Quintile spread |", "|---|---|---:|---:|---:|---|---:|"] for era, summary in screen["summaries"].items(): for arm, row in summary["arms"].items(): ci = row["ic_ci95"] ci_text = f"[{ci[0]:.4f}, {ci[1]:.4f}]" if ci else "unavailable" lines.append(f"| {era} | {arm} | {summary['windows']} | {row['avg_cross_section']:.1f} " f"| {row['mean_ic']:.4f} | {ci_text} | {100 * row['mean_quintile_spread']:.2f}% |") lines.extend(["", "Paired normalized-minus-control IC differences:", ""]) for era, summary in screen["summaries"].items(): if summary["windows"]: lines.append(f"- {era}: {summary['paired_mean_ic_delta']:.4f}; " f"95% interval {summary['paired_ic_delta_ci95']}.") lines.extend(["", f"Screen checks: `{json.dumps(screen['verdict']['checks'])}`", "", f"Portfolio stage: **{report['portfolio_stage']}**.", "", "Signal IC and quintile spreads are predictive diagnostics, not simulated strategy returns.", "Intervals resample paired blocks of four non-overlapping forward windows; they remain exploratory.", "The unrestricted control shows the effect of requiring complete formation history.", "", "[Frozen protocol](../docs/research/residual-normalization.md)", ""]) return "\n".join(lines) def main() -> int: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--snapshot", type=Path, default=ROOT / "backtest_snapshots/research.sqlite") parser.add_argument("--out", type=Path, default=ROOT / "reports/residual-normalization-20260908.json") args = parser.parse_args() snapshot = args.snapshot.resolve(strict=True) before = sha256(snapshot) code_paths = [Path(__file__).resolve(), PROTOCOL, ROOT / "app/services/momentum_service.py", ROOT / "app/services/backtest_service.py"] hashes = {p.relative_to(ROOT).as_posix(): git_blob_id(p) for p in code_paths} prices, benchmark, manifest = read_snapshot(snapshot) print(f"manifest verified: {len(prices)} production-universe names; snapshot {before[:16]}", flush=True) screen = run_screen(prices, benchmark) if sha256(snapshot) != before or any(git_blob_id(ROOT / p) != h for p, h in hashes.items()): raise RuntimeError("Snapshot, protocol, or implementation changed during the run") report = { "schema_version": 1, "generated_at": datetime.now(timezone.utc).isoformat(), "status": "complete_signal_screen", "snapshot": str(snapshot), "snapshot_sha256": before, "source_git_blobs": hashes, "manifest": manifest, "python": sys.version, "packages": {p: importlib.metadata.version(p) for p in ("numpy", "sqlalchemy")}, "screen": screen, "portfolio_stage": "eligible_pending_replay" if screen["verdict"]["passed"] else "not_run_failed_signal_screen", "production_changed": False, } args.out.parent.mkdir(parents=True, exist_ok=True) args.out.write_text(json.dumps(report, indent=2, allow_nan=False) + "\n", encoding="utf-8") args.out.with_suffix(".md").write_text(render_report(report), encoding="utf-8") print(json.dumps({"verdict": screen["verdict"], "summaries": screen["summaries"], "report": str(args.out)}, indent=2), flush=True) return 0 if __name__ == "__main__": raise SystemExit(main())