From ff7ff0d5ba113083bdc5d7503313df64a9a9f763 Mon Sep 17 00:00:00 2001 From: Shorya Mishra Date: Wed, 29 Jul 2026 23:01:54 +0530 Subject: [PATCH] fix: correct merger premium indexing, crisis correlation loss, and noise bound API --- events/operators.py | 24 +++++++---- state/noise.py | 33 ++++++++++----- tests/test_events.py | 99 ++++++++++++++++++++++++++++++++++++++++++++ tests/test_noise.py | 34 +++++++++++++++ 4 files changed, 172 insertions(+), 18 deletions(-) diff --git a/events/operators.py b/events/operators.py index a1c008e..5d06f5a 100644 --- a/events/operators.py +++ b/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", diff --git a/state/noise.py b/state/noise.py index fa4adac..39db030 100644 --- a/state/noise.py +++ b/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) diff --git a/tests/test_events.py b/tests/test_events.py index bfe27de..2d0eab3 100644 --- a/tests/test_events.py +++ b/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"]) diff --git a/tests/test_noise.py b/tests/test_noise.py index 25064f8..aa50887 100644 --- a/tests/test_noise.py +++ b/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"])