fix: correct merger premium indexing, crisis correlation loss, and noise bound API

events/operators.py

@@ -488,13 +488,15 @@ def systemic_crisis_operator(severity: float, n: int = 1, d: int = 5) -> EventOp
         b[i * d + I] = -severity * 0.20        # opacity spike
 
     # Near-singular covariance: correlations → 1 in crisis
+    # Build the target covariance matrix for price components, then Cholesky-factor it.
     rho = 0.5 + 0.4 * severity
-    Sigma = np.zeros((n * d, n * d))
+    sigma_p = severity * 0.08  # per-asset price std dev
+    Cov = np.eye(n * d) * sigma_p ** 2
     for i in range(n):
         for j in range(n):
-            Sigma[i * d + P, j * d + P] = (severity * 0.08) ** 2 * (rho if i != j else 1.0)
-    Sigma = np.sqrt(np.abs(Sigma)) * np.sign(Sigma)  # back to std dev scale
-    Sigma = np.diag(np.maximum(np.diag(Sigma), severity * 0.08))  # ensure positive diagonal
+            if i != j:
+                Cov[i * d + P, j * d + P] = sigma_p ** 2 * rho
+    Sigma = np.linalg.cholesky(Cov)
 
     return EventOperator(
         name=f"systemic_crisis_severity{severity:.1f}",
@@ -624,13 +626,19 @@ def merger_operator(acquirer_idx: int, target_idx: int, n: int = 2, d: int = 5,
             A[out_idx * d:(out_idx + 1) * d, i * d:(i + 1) * d] = np.eye(d)
             out_idx += 1
 
+    # Compute acquirer's output row index
+    acq_out_idx = 0
+    for i in range(acquirer_idx):
+        if i != target_idx:
+            acq_out_idx += 1
+
     b = np.zeros(m * d)
-    b[0 * d + P] = np.log(1 + premium_pct / 100)  # deal premium on acquirer row
-    b[0 * d + V] = np.log(1.5)  # volume addition (approximate log-sum-exp)
-    b[0 * d + L] = 0.05  # slight leverage increase from deal financing
+    b[acq_out_idx * d + P] = np.log(1 + premium_pct / 100)  # deal premium on acquirer row
+    b[acq_out_idx * d + V] = np.log(1.5)  # volume addition (approximate log-sum-exp)
+    b[acq_out_idx * d + L] = 0.05  # slight leverage increase from deal financing
 
     Sigma = np.eye(m * d) * 0.05
-    Sigma[0 * d + P, 0 * d + P] = 0.08  # acquirer price most uncertain
+    Sigma[acq_out_idx * d + P, acq_out_idx * d + P] = 0.08  # acquirer price most uncertain
 
     return EventOperator(
         name=f"merger_acq{acquirer_idx}_tgt{target_idx}_prem{premium_pct:.0f}pct",

state/noise.py

@@ -14,7 +14,7 @@
 """
 import numpy as np
 from dataclasses import dataclass
-from typing import Optional
+from numpy.typing import NDArray
 
 
 @dataclass
@@ -27,9 +27,8 @@ class DualNoiseParams:
     @property
     def tau_t(self) -> float:
         """Composite temperature parameter from §III.5: τ_t = √(σ_τ² + λ_η · m₂^η)"""
-        return np.sqrt(self.sigma_tau**2 + self.lambda_eta * self.m2_eta)
+        return float(np.sqrt(self.sigma_tau**2 + self.lambda_eta * self.m2_eta))
 
-    @property
     def cramer_rao_bound(self, h: float = 1.0) -> float:
         """Prediction variance lower bound for horizon h (Theorem III.1)."""
         return (self.sigma_tau**2 + self.lambda_eta * self.m2_eta) * h
@@ -49,14 +48,19 @@ class DualNoiseCalibrator:
     def __init__(self, alpha_lm: float = 0.001):
         self.alpha_lm = alpha_lm  # Lee-Mykland significance level
 
-    def estimate_bpv(self, returns: np.ndarray) -> float:
+    def estimate_bpv(self, returns: NDArray[np.float64]) -> float:
         """Bipower variation estimate of integrated physical variance."""
         if len(returns) < 2:
             return 0.0
         bpv = (np.pi / 2) * np.sum(np.abs(returns[1:]) * np.abs(returns[:-1]))
-        return bpv
-
-    def detect_jumps(self, returns: np.ndarray, bpv: float, dt: float = 1/78) -> np.ndarray:
+        return float(bpv)
+
+    def detect_jumps(
+        self,
+        returns: NDArray[np.float64],
+        bpv: float,
+        dt: float = 1 / 78,
+    ) -> NDArray[np.bool_]:
         """
         Lee-Mykland (2008) jump detection. Returns boolean mask of jump indices.
         dt = bar interval as fraction of day (e.g. 1/78 for 5-min bars).
@@ -68,9 +72,13 @@ def detect_jumps(self, returns: np.ndarray, bpv: float, dt: float = 1/78) -> np.
             return np.zeros(0, dtype=bool)
         sigma_hat = np.sqrt(bpv / len(returns))  # per-bar vol estimate
         critical = sigma_hat * self._lee_mykland_critical(len(returns))
-        return np.abs(returns) > critical
+        return np.asarray(np.abs(returns) > critical, dtype=np.bool_)
 
-    def calibrate(self, intraday_returns: np.ndarray, dt: float = 1/78) -> DualNoiseParams:
+    def calibrate(
+        self,
+        intraday_returns: NDArray[np.float64],
+        dt: float = 1 / 78,
+    ) -> DualNoiseParams:
         bpv = self.estimate_bpv(intraday_returns)
         jump_mask = self.detect_jumps(intraday_returns, bpv, dt)
         jump_sizes = intraday_returns[jump_mask]
@@ -89,4 +97,9 @@ def calibrate(self, intraday_returns: np.ndarray, dt: float = 1/78) -> DualNoise
     def _lee_mykland_critical(n: int, alpha: float = 0.001) -> float:
         """Approximate critical value for Lee-Mykland test."""
         c = np.sqrt(2 * np.log(n))
-        return c + (np.log(np.log(n)) + np.log(4 * np.pi) - 2 * np.log(-np.log(1 - alpha))) / (2 * c)
+        critical = c + (
+            np.log(np.log(n))
+            + np.log(4 * np.pi)
+            - 2 * np.log(-np.log(1 - alpha))
+        ) / (2 * c)
+        return float(critical)

tests/test_events.py

@@ -214,5 +214,104 @@ def test_event_sequence_dimension_tracking(self):
         assert len(log) == len(operators)
 
 
+class TestMergerPremiumIndexing:
+    """Regression tests for merger_operator b_w indexing (Bug #1)."""
+
+    def test_merger_premium_follows_acquirer_index(self):
+        """Deal premium must be applied at the acquirer's output row, not row 0."""
+        n, d = 5, 5
+        acquirer_idx, target_idx = 2, 4
+        premium_pct = 30.0
+        op = merger_operator(
+            acquirer_idx=acquirer_idx, target_idx=target_idx,
+            n=n, premium_pct=premium_pct,
+        )
+        # Acquirer (idx 2) should map to output row 2:
+        # output rows: 0->0, 1->1, 2(acq)->2, 3->3, 4(tgt)->skipped
+        acq_out_idx = 2
+        expected_premium = np.log(1 + premium_pct / 100)
+
+        # Premium should be at acquirer output row
+        np.testing.assert_allclose(
+            op.b_w[acq_out_idx * d + P], expected_premium, rtol=1e-6,
+            err_msg="Premium should be at acquirer output row",
+        )
+        # Row 0 should have zero price shift (it's an uninvolved pass-through asset)
+        np.testing.assert_allclose(
+            op.b_w[0 * d + P], 0.0, atol=1e-10,
+            err_msg="Row 0 should have no premium when acquirer_idx != 0",
+        )
+
+    def test_merger_premium_at_zero_unchanged(self):
+        """When acquirer_idx=0, behavior is unchanged (backward compatibility)."""
+        op = merger_operator(acquirer_idx=0, target_idx=1, n=3, premium_pct=25.0)
+        expected = np.log(1.25)
+        np.testing.assert_allclose(op.b_w[0 * 5 + P], expected, rtol=1e-6)
+
+    def test_merger_volume_and_leverage_follow_acquirer(self):
+        """Volume spike and leverage shift must also follow the acquirer row."""
+        n, d = 4, 5
+        op = merger_operator(acquirer_idx=2, target_idx=0, n=n)
+        # Output mapping: 0(tgt)->skip, 1->0, 2(acq)->1, 3->2
+        acq_out_idx = 1
+        assert op.b_w[acq_out_idx * d + V] > 0, "Volume spike should be at acquirer row"
+        assert op.b_w[acq_out_idx * d + L] > 0, "Leverage shift should be at acquirer row"
+        # Other rows should have no volume/leverage shifts from merger
+        for row in [0, 2]:
+            assert op.b_w[row * d + V] == 0.0, f"Row {row} should have no volume shift"
+            assert op.b_w[row * d + L] == 0.0, f"Row {row} should have no leverage shift"
+
+
+class TestSystemicCrisisCorrelation:
+    """Regression tests for systemic_crisis_operator Sigma_w (Bug #2)."""
+
+    def test_systemic_crisis_has_cross_asset_correlation(self):
+        """Crisis Sigma_w must have non-zero off-diagonal entries in the price block."""
+        n = 4
+        op = systemic_crisis_operator(severity=0.8, n=n)
+        S = op.Sigma_w
+        # The implied covariance Cov = S @ S.T should have non-zero off-diags
+        Cov = S @ S.T
+        price_indices = [i * 5 + P for i in range(n)]
+        for i in price_indices:
+            for j in price_indices:
+                if i != j:
+                    assert abs(Cov[i, j]) > 1e-6, \
+                        f"Cov[{i},{j}] should be non-zero (crisis = correlated moves)"
+
+    def test_systemic_crisis_covariance_is_psd(self):
+        """Implied covariance must be PSD and preserve component variances."""
+        severity = 1.0
+        op = systemic_crisis_operator(severity=severity, n=5)
+        Cov = op.Sigma_w @ op.Sigma_w.T
+        eigs = np.linalg.eigvalsh(Cov)
+        assert eigs.min() >= -1e-10, f"Min eigenvalue {eigs.min()} — covariance not PSD"
+        np.testing.assert_allclose(np.diag(Cov), (severity * 0.08) ** 2)
+
+    def test_systemic_crisis_correlation_increases_with_severity(self):
+        """Higher severity should produce higher cross-asset correlation."""
+        n = 3
+        def avg_price_corr(severity):
+            op = systemic_crisis_operator(severity=severity, n=n)
+            Cov = op.Sigma_w @ op.Sigma_w.T
+            pidx = [i * 5 + P for i in range(n)]
+            stds = np.sqrt(np.array([Cov[i, i] for i in pidx]))
+            corrs = []
+            for a in range(n):
+                for b in range(a + 1, n):
+                    corrs.append(Cov[pidx[a], pidx[b]] / (stds[a] * stds[b] + 1e-12))
+            return np.mean(corrs)
+
+        rho_low = avg_price_corr(0.3)
+        rho_high = avg_price_corr(0.9)
+        assert rho_high > rho_low, \
+            f"Higher severity should give higher correlation: {rho_high:.3f} vs {rho_low:.3f}"
+
+    def test_systemic_crisis_single_asset(self):
+        """n=1 should not crash (no off-diagonal terms needed)."""
+        op = systemic_crisis_operator(severity=0.5, n=1)
+        assert op.Sigma_w.shape == (5, 5)
+
+
 if __name__ == "__main__":
     pytest.main([__file__, "-v"])

tests/test_noise.py

@@ -138,6 +138,40 @@ def test_higher_vol_gives_higher_sigma_tau(self):
             "Higher vol should give higher sigma_tau"
 
 
+class TestCramerRaoBound:
+    """Regression tests for DualNoiseParams.cramer_rao_bound (Bug #3)."""
+
+    def test_cramer_rao_bound_callable_with_h(self):
+        """cramer_rao_bound must accept explicit horizon parameter h."""
+        params = DualNoiseParams(sigma_tau=0.01, lambda_eta=2.0, m2_eta=0.001)
+        # Should not raise TypeError
+        result = params.cramer_rao_bound(h=5.0)
+        assert isinstance(result, float)
+        assert result > 0
+
+    def test_cramer_rao_bound_scales_with_horizon(self):
+        """Bound at horizon h should be h times the bound at horizon 1."""
+        params = DualNoiseParams(sigma_tau=0.01, lambda_eta=2.0, m2_eta=0.001)
+        bound_1 = params.cramer_rao_bound(h=1.0)
+        bound_5 = params.cramer_rao_bound(h=5.0)
+        np.testing.assert_allclose(bound_5, 5.0 * bound_1, rtol=1e-10)
+
+    def test_cramer_rao_bound_default_h(self):
+        """Default h=1.0 should produce the same result as explicit h=1.0."""
+        params = DualNoiseParams(sigma_tau=0.02, lambda_eta=1.0, m2_eta=0.005)
+        np.testing.assert_allclose(
+            params.cramer_rao_bound(),
+            params.cramer_rao_bound(h=1.0),
+        )
+
+    def test_cramer_rao_bound_matches_formula(self):
+        """Verify against Theorem III.1: (σ_τ² + λ_η · m₂^η) · h."""
+        params = DualNoiseParams(sigma_tau=0.03, lambda_eta=5.0, m2_eta=0.01)
+        h = 3.0
+        expected = (0.03**2 + 5.0 * 0.01) * h
+        np.testing.assert_allclose(params.cramer_rao_bound(h=h), expected, rtol=1e-10)
+
+
 if __name__ == "__main__":
     import pytest
     pytest.main([__file__, "-v"])