feat(research): commit the regime-monitor replay harness, and reproduce v2/v3

v3 was calibrated by replaying the series offline, but that harness was never
committed -- so its published numbers could not be re-derived, and a v4 cut would
have had to choose anchors by argument rather than measurement. This is that
harness, and it reproduces the published figures.

scripts/run_regime_monitor_calibration.py replays State/Warning session by
session from the same inputs the live job uses (Alpaca for all 33 symbols, FRED
for VIX and HY OAS), with no database: breadth and divergence come from
breadth_service's pure helpers. It never reimplements an unchanged live sensor --
_compute_index, _score_pillars, P2, P4 and the Warning sensors are imported and
called. Only candidate formulas (proposed v4) and retired ones (v2, gone from the
codebase) are defined here and patched onto the module for a variant's duration.

Reproduction of the 408 sessions ending 2026-07-24, against the figures in
docs/research/regime-monitor-v3.md:

  v2 State avg      22.6   ->  22.68
  v2 State p80      35.1   ->  35.1     exact
  v2 State max      91.2   ->  91.2     exact
  v2 P3 pegged        39   ->    39     exact
  v2 W1 live         108   ->   108     exact
  v3 State max      87.4   ->  87.4     exact
  v3 band shares  73.3/15.0/8.3/3.4 -> 73.0/15.4/8.1/3.4

Three things the harness had to get right to reach that, each of which was
initially wrong and caught by a gate rather than by inspection:

  - "W1 live 108" counts NONZERO sessions, not non-null ones. v2's divergence
    gate returned 0.0 during any decline (v3 tapers instead), so the retired
    divergence formula had to be reconstructed too.
  - v2 sliced HY_OAS_REFERENCE_YEARS = 10.0 per session, not v3's 700 days. The
    percentile leg ranks against that window, so replaying it short shifted the
    middle of the distribution while leaving the max exact.
  - The published v2 numbers correspond to FULL OAS coverage. Replaying v2 with
    the 400-calendar-day fetch it shipped with yields max 100.0 and 133
    credit-less sessions -- so that truncation was not in force when the figures
    were taken. Recorded rather than assumed.

The script refuses to emit a band recommendation unless every hard gate passes
(33 symbols fetched, per-symbol warm-up and final bar, full basket on every
session, calendar anchors, 100% coverage, and a row-wise state_v4 <= state_v3
invariant), and exits non-zero. It is meant to be structurally impossible to read
a calibration result out of a run whose pipeline did not validate. No v4 code
ships in this commit.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

tests/unit/test_regime_calibration_script.py

@@ -0,0 +1,140 @@
+"""The calibration harness's pure helpers. No network, no data.
+
+The harness's whole value is that its numbers can be trusted, so the parts that
+decide whether a run is trustworthy — the variant patching and the gate
+evaluation — are worth pinning even though the script is research-only.
+"""
+
+import importlib.util
+from pathlib import Path
+
+import pytest
+
+from app.services import regime_monitor_service as rms
+
+_SPEC = importlib.util.spec_from_file_location(
+    "regime_calibration",
+    Path(__file__).resolve().parents[2] / "scripts" / "run_regime_monitor_calibration.py",
+)
+calib = importlib.util.module_from_spec(_SPEC)
+_SPEC.loader.exec_module(calib)
+
+
+class TestVariantPatching:
+    def test_patching_restores_module_state_exactly(self):
+        """A leaked patch would silently contaminate every later variant."""
+        before = {
+            name: getattr(rms, name)
+            for name in ("p1_trend_break", "p5_volatility", "p3_drawdown",
+                         "f2_credit_spreads", "HY_OAS_WINDOW_DAYS")
+        }
+        with calib.patched(calib.VARIANTS["v2_reconstruction"]):
+            assert rms.p3_drawdown is not before["p3_drawdown"]
+            assert rms.HY_OAS_WINDOW_DAYS == 3653
+        for name, original in before.items():
+            assert getattr(rms, name) is original or getattr(rms, name) == original
+
+    def test_patching_restores_even_when_the_body_raises(self):
+        original = rms.p5_volatility
+        with pytest.raises(RuntimeError):
+            with calib.patched(calib.VARIANTS["v4"]):
+                raise RuntimeError("boom")
+        assert rms.p5_volatility is original
+
+    def test_v3_variant_patches_nothing(self):
+        """The reproduction gate must run against live code, not a copy."""
+        assert calib.VARIANTS["v3"] == {}
+
+
+class TestCandidateFormulas:
+    def test_graduated_trend_break_never_exceeds_the_binary_one(self):
+        """Half of the provable state_v4 <= state_v3 invariant."""
+        closes = [100.0] * 200
+        for last in (105.0, 100.0, 99.0, 92.0, 80.0, 50.0):
+            series = closes[:-1] + [last]
+            v4 = calib._candidate_under_200(series)
+            v3 = rms._under_200(series)
+            assert v4 <= v3, f"close={last}: v4 {v4} > v3 {v3}"
+
+    def test_anchored_vix_never_exceeds_the_live_formula(self):
+        p5 = calib._candidate_p5(calib.P5_VIX_ANCHORS_A)
+        for vix in (10, 15, 17, 20, 25, 30, 40, 55, 82):
+            assert p5(vix) <= rms.p5_volatility(vix), f"vix={vix}"
+
+    def test_the_vix_table_keeps_resolving_past_thirty(self):
+        p5 = calib._candidate_p5(calib.P5_VIX_ANCHORS_A)
+        assert p5(30) < p5(40) < p5(50) < p5(55) == 100.0
+        assert rms.p5_volatility(30) == rms.p5_volatility(82) == 100.0  # the defect
+
+    def test_a_shallow_break_no_longer_pegs(self):
+        closes = [100.0] * 199 + [98.0]  # ~2% below a flat 200-DMA
+        assert rms._under_200(closes) == 100.0
+        assert calib._candidate_under_200(closes) < 40.0
+
+    def test_candidate_tables_are_well_formed(self):
+        for table in (calib.P5_VIX_ANCHORS_A, calib.P5_VIX_ANCHORS_B,
+                      calib.P1_TREND_BREAK_ANCHORS):
+            xs = [x for x, _ in table]
+            ys = [y for _, y in table]
+            assert xs == sorted(xs) and len(set(xs)) == len(xs)
+            assert ys == sorted(ys)
+            assert 0.0 <= min(ys) and max(ys) <= 100.0
+
+
+class TestGates:
+    def _row(self, date_str, state, coverage=100.0, w1=5.0):
+        return {
+            "date": date_str,
+            "state": {"score": state, "coverage": coverage, "pillars": []},
+            "warning": {"score": 10.0, "coverage": 100.0, "pillars": [
+                {"id": "breadth_divergence", "sensors": [{"id": "W1", "score": w1}]}
+            ]},
+            "data_quality": {"stale_inputs": []},
+        }
+
+    def test_the_invariant_gate_catches_a_v4_row_above_its_v3_row(self):
+        v3 = [self._row("2026-01-02", 40.0), self._row("2026-01-05", 50.0)]
+        v4 = [self._row("2026-01-02", 38.0), self._row("2026-01-05", 55.0)]
+        gate = calib._invariant_gate(v3, v4)
+        assert gate["passed"] is False
+        assert gate["detail"][0]["date"] == "2026-01-05"
+
+    def test_the_invariant_gate_passes_when_v4_is_never_higher(self):
+        v3 = [self._row("2026-01-02", 40.0)]
+        v4 = [self._row("2026-01-02", 40.0)]
+        assert calib._invariant_gate(v3, v4)["passed"] is True
+
+    def test_sensor_lookup_searches_both_axes(self):
+        """W1 lives under warning; a state-only search silently returns None and
+        made the W1 census read 0."""
+        row = self._row("2026-01-02", 40.0, w1=7.5)
+        assert calib._sensor_score(row, "breadth_divergence", "W1") == 7.5
+
+
+class TestBandShares:
+    def test_shares_use_the_candidate_bands_not_the_imported_default(self):
+        """band_for binds bands=STATE_BANDS as an import-time default, so the
+        grid must pass candidates explicitly or every row scores identically."""
+        scores = [10.0, 30.0, 55.0, 75.0]
+        loose = calib._band_shares(scores, (20.0, 50.0, 80.0))
+        tight = calib._band_shares(scores, (20.0, 50.0, 60.0))
+        assert loose["breaking"] == 0.0
+        assert tight["breaking"] == 25.0
+
+    def test_shares_sum_to_one_hundred(self):
+        shares = calib._band_shares([5.0, 25.0, 55.0, 85.0], (20.0, 50.0, 80.0))
+        assert sum(shares.values()) == pytest.approx(100.0)
+
+
+class TestScenarios:
+    def test_the_decisive_scenario_is_computed_not_asserted(self):
+        rows = {s["scenario"].split()[0]: s for s in
+                calib._scenarios(calib.P5_VIX_ANCHORS_A, calib.P1_TREND_BREAK_ANCHORS)}
+        # A 2022-style AI/tech drawdown with genuinely calm credit. This is the
+        # meaning anchor for the breaking threshold, so it is machine-checked.
+        assert rows["S3a"]["C1"] == 0.0
+        assert rows["S3a"]["state"] == pytest.approx(70.33, abs=0.01)
+        assert rows["S3b"]["state"] == pytest.approx(74.00, abs=0.01)
+        # ...and it must clear the chosen threshold under either P2 assumption.
+        assert min(rows["S3a"]["state"], rows["S3b"]["state"]) > 65.0
+        assert rows["S1"]["state"] < 20.0