diff --git a/changes.md b/changes.md new file mode 100644 index 0000000..cf6bd90 --- /dev/null +++ b/changes.md @@ -0,0 +1,140 @@ +# Changes + +Branch: `fix/dual-noise-units-and-compose-cholesky` +Commit: `2265b84` + +Two bugs fixed. Both were about **units and matrix algebra**, not about logic +or features. Nothing was added or removed — two calculations were wrong and are +now right. + +--- + +## Bug 1 — the jump counter wasn't a rate + +**File:** `state/noise.py` + +**What it should do:** measure how often a market jumps, in *jumps per day*. + +**What it did:** multiplied the jump count by the bar width instead of dividing +by elapsed time. That's the upside-down version of a rate. + +```python +lambda_eta = jump_mask.sum() * dt # before — a count, scaled wrongly +lambda_eta = jump_mask.sum() / (n*dt) # after — a rate +``` + +The volatility estimate `sigma_tau` had the same problem: it measured +volatility *per bar* instead of *per day*. + +**Why it mattered:** the repo's headline result (Theorem III.1, the +Cramér-Rao bound) **adds** these two numbers together. If one is per-bar and +the other is per-day, adding them is meaningless — like adding 5 kilometres to +3 hours. + +**Measured effect.** A test market with a known 3 jumps/day and 2% daily +volatility: + +| | true | before | after | +|---|---|---|---| +| jumps per day | 3.0 | 0.38 | **3.00** | +| daily volatility | 0.020 | 0.0022 | **0.0198** | + +Worse than being wrong, it was *inconsistently* wrong: feeding the same market +in as 5-minute bars vs 1-minute bars gave answers 5x apart. The Cramér-Rao +bound is supposed to describe a market, not your data feed. It now agrees +within 1% across both. + +**One caller fixed:** `demo/run_egamec.py` passed `dt=1/252` for daily bars, +which told the code each bar was 1/252 of a day. Corrected to `dt=1.0`. + +--- + +## Bug 2 — event composition crashed on market crises + +**File:** `events/operators.py` + +**What it should do:** combine two market events into one (e.g. "rate hike +after a systemic crisis") and work out the combined uncertainty. + +**What it did:** `Sigma_w` holds noise in a form where the real covariance is +`S x S-transpose`. Every other part of the codebase reads it that way. +`compose()` used `S x S` — no transpose. + +```python +op1.Sigma_w @ op1.Sigma_w # before +op1.Sigma_w @ op1.Sigma_w.T # after +``` + +**Why it mattered:** the two forms only agree when the matrix is diagonal. +Exactly one event type has a non-diagonal one — `systemic_crisis_operator`, +where the whole point is that asset prices fall *together*. + +For that operator the wrong form produces a matrix that isn't a valid +covariance at all, so the next line crashed: + +``` +LinAlgError: Matrix is not positive definite +``` + +**In plain terms: you could not combine a market crisis with any other event.** +Not "you got a wrong number" — the code stopped. + +Nothing in the codebase composed a crisis operator, and every existing +composition test used diagonal-noise events, so this was never hit. + +--- + +## Files changed + +| File | What | +|---|---| +| `state/noise.py` | 2 lines of math + docstrings stating the time units | +| `events/operators.py` | 1 line of math + a comment on why | +| `demo/run_egamec.py` | 1 line — corrected `dt` for daily bars | +| `tests/test_noise.py` | +6 tests (new, additive) | +| `tests/test_events.py` | +5 tests (new, additive) | + +Total: 268 added, 9 removed. Only **4 lines** are real logic; the rest is +tests and comments. + +--- + +## Testing + +Suite: **133 -> 144 tests, all passing.** + +Each new test was checked against the old code to confirm it actually catches +the bug: **5 of 6** new noise tests and **4 of 5** new event tests fail on the +previous commit. The ones that pass on old code are intentional — each group +keeps one diagonal-noise case to make sure previously-correct behaviour didn't +move. + +The tests check *meaning*, not memorised numbers. The strongest one holds the +jump **count** fixed while stretching the time window — 6 jumps in 2 days is a +different rate than 6 jumps in 20 days. A count-shaped estimator returns the +same number for both and gets caught. + +--- + +## What still works + +- All 144 tests pass +- `demo/run_egamec.py`, `demo/denoised_price_2026.py`, `demo/hindcast_2008.py` + all run clean +- CI is unaffected (it runs `pytest tests/` plus 3 import checks) +- `compose()` has zero non-test callers, so that fix can't regress anything — + it only makes a previously-impossible operation possible + +## What this affects downstream + +`notebooks/day03_dual_noise.ipynb` prints the estimate next to a hardcoded +"true" value. That line now visibly disagrees. + +The estimate is the part that's now **correct**: it returns 0.190 against that +notebook's true annual volatility of 0.20 (its `dt` is in years, so the answer +comes out annualised). The stale part is the hardcoded string, which is a +per-bar number labelled "/day". + +Two-line notebook fix. Not included in this commit. + +--- diff --git a/demo/run_egamec.py b/demo/run_egamec.py index dda57ef..8e7b192 100644 --- a/demo/run_egamec.py +++ b/demo/run_egamec.py @@ -36,7 +36,9 @@ results = {} for ticker in tickers[:10]: # demo: first 10 tickers rets = df.xs(ticker, level="ticker")["returns"].values - params = calibrator.calibrate(rets, dt=1/252) + # synthetic_market generates DAILY bars, so dt = 1 day. (dt=1/252 would + # declare each bar to be 1/252 of a day and inflate every rate by 252×.) + params = calibrator.calibrate(rets, dt=1.0) results[ticker] = params avg_sigma_tau = np.mean([p.sigma_tau for p in results.values()]) diff --git a/events/operators.py b/events/operators.py index 5d06f5a..1c63785 100644 --- a/events/operators.py +++ b/events/operators.py @@ -852,10 +852,19 @@ def compose(op1: EventOperator, op2: EventOperator) -> EventOperator: A_comp = op1.A_w @ op2.A_w b_comp = op1.A_w @ op2.b_w + op1.b_w - # Uncertainty propagation: Var(A1(A2 s + ε2) + ε1) = A1 Σ2 A1' + Σ1 - Sigma_comp_sq = op1.Sigma_w @ op1.Sigma_w + op1.A_w @ (op2.Sigma_w @ op2.Sigma_w) @ op1.A_w.T - # Take element-wise sqrt to get back to Cholesky-scale - Sigma_comp = np.linalg.cholesky(Sigma_comp_sq + 1e-9 * np.eye(Sigma_comp_sq.shape[0])) + # Uncertainty propagation. Sigma_w is a Cholesky FACTOR, so the covariance + # it represents is Σ Σᵀ — not Σ Σ. The two coincide only for diagonal Σ, + # which is why this went unnoticed: systemic_crisis_operator is the one + # constructor with off-diagonal noise, and nothing composed it. For that + # operator `Σ Σ` is not even PSD, so the cholesky() below raised + # LinAlgError — composing a systemic crisis with anything was impossible. + # Cov(A1(A2 s + L2 ε2) + L1 ε1) = L1 L1ᵀ + A1 (L2 L2ᵀ) A1ᵀ + Cov_comp = ( + op1.Sigma_w @ op1.Sigma_w.T + + op1.A_w @ (op2.Sigma_w @ op2.Sigma_w.T) @ op1.A_w.T + ) + # Back to a Cholesky factor, so the result composes again under this same law. + Sigma_comp = np.linalg.cholesky(Cov_comp + 1e-9 * np.eye(Cov_comp.shape[0])) return EventOperator( name=f"({op1.name}) ∘ ({op2.name})", diff --git a/state/noise.py b/state/noise.py index 39db030..b02e83b 100644 --- a/state/noise.py +++ b/state/noise.py @@ -19,8 +19,16 @@ @dataclass class DualNoiseParams: - sigma_tau: float # physical noise volatility (σ_τ) - lambda_eta: float # behavioral jump intensity (λ_η) + """ + Calibrated dual-noise parameters. + + σ_τ and λ_η are both *per unit time* on the same clock (per day when + calibrated from bars whose dt is expressed in days). tau_t and + cramer_rao_bound add them, so a mismatch there silently rescales + Theorem III.1. + """ + sigma_tau: float # physical noise volatility (σ_τ), per √(unit time) + lambda_eta: float # behavioral jump intensity (λ_η), jumps per unit time m2_eta: float # second moment of jump size ∫z² ν^η(dz) alpha_levy: float = 1.5 # Lévy tail index ∈ (1, 2) @@ -67,6 +75,11 @@ def detect_jumps( Test statistic: |r_t| / sigma_hat vs c_alpha, where sigma_hat = sqrt(BPV/n) is the per-bar volatility estimate. + + `dt` is accepted for signature symmetry with calibrate() but is not + used: both sides of the comparison are per-bar quantities, so the bar + width cancels. Rescaling to a per-unit-time clock happens in + calibrate(), not here. """ if len(returns) == 0: return np.zeros(0, dtype=bool) @@ -79,12 +92,31 @@ def calibrate( intraday_returns: NDArray[np.float64], dt: float = 1 / 78, ) -> DualNoiseParams: + """ + Estimate (σ_τ, λ_η, m₂^η) from a sample of `n` bars of width `dt`. + + Units. Every rate below is per unit of the clock `dt` is measured in + (dt = 1/78 → bars are 1/78 of a day → rates are per day). Theorem III.1 + adds σ_τ² and λ_η·m₂^η, so the two must share that clock: + + T = n · dt total sample time + σ_τ² = BPV / T integrated variance ÷ elapsed time + λ_η = N_jumps / T jump count ÷ elapsed time + m₂^η = mean(z²) jump-size second moment (no time unit) + + Passing daily bars therefore means dt=1.0, not dt=1/252. + """ + n = len(intraday_returns) + if n == 0 or dt <= 0: + return DualNoiseParams(sigma_tau=0.0, lambda_eta=0.0, m2_eta=1e-6) + bpv = self.estimate_bpv(intraday_returns) jump_mask = self.detect_jumps(intraday_returns, bpv, dt) jump_sizes = intraday_returns[jump_mask] - sigma_tau = np.sqrt(bpv / len(intraday_returns)) - lambda_eta = jump_mask.sum() * dt if dt > 0 else 0.0 + elapsed = n * dt + sigma_tau = float(np.sqrt(bpv / elapsed)) + lambda_eta = float(jump_mask.sum() / elapsed) m2_eta = float(np.mean(jump_sizes**2)) if len(jump_sizes) > 0 else 1e-6 return DualNoiseParams( diff --git a/tests/test_events.py b/tests/test_events.py index 2d0eab3..cd0be3e 100644 --- a/tests/test_events.py +++ b/tests/test_events.py @@ -313,5 +313,100 @@ def test_systemic_crisis_single_asset(self): assert op.Sigma_w.shape == (5, 5) +class TestCompositionNoisePropagation: + """ + Regression tests for the Sigma_w propagation law in compose(). + + Sigma_w is documented as a Cholesky FACTOR, and every other test in this + file reads its covariance as `S @ S.T`. compose() used `S @ S` — which + agrees only when S is diagonal. Every composition test above uses + diagonal-noise operators, and systemic_crisis_operator (the one + constructor with off-diagonal noise) was never composed, so the gap + stayed invisible. + + The failure mode was not merely a wrong covariance: for the crisis + operator `S @ S` is not PSD, so np.linalg.cholesky raised LinAlgError and + a systemic crisis could not be composed with any other event at all. + """ + + @staticmethod + def _cov(op): + return op.Sigma_w @ op.Sigma_w.T + + def test_composition_covariance_matches_propagation_law(self): + """ + Cov(A1(A2 s + L2 ε2) + L1 ε1) = L1 L1ᵀ + A1 (L2 L2ᵀ) A1ᵀ. + + Uses systemic_crisis as the INNER operator so its off-diagonal block + is the thing being propagated. + """ + n = 3 + op_crisis = systemic_crisis_operator(severity=0.8, n=n) + op_rate = rate_change_operator(change_bps=50, n=n) + + comp = compose(op_rate, op_crisis) # crisis first, then rate change + + expected = ( + self._cov(op_rate) + + op_rate.A_w @ self._cov(op_crisis) @ op_rate.A_w.T + ) + np.testing.assert_allclose(self._cov(comp), expected, atol=1e-8) + + def test_composition_preserves_cross_asset_correlation(self): + """ + The economic content: composing a macro event onto a systemic crisis + must not destroy the crisis's correlated-price structure. `S @ S` + silently reshapes those off-diagonals. + """ + n = 4 + comp = compose( + rate_change_operator(change_bps=25, n=n), + systemic_crisis_operator(severity=0.9, n=n), + ) + Cov = self._cov(comp) + pidx = [i * 5 + P for i in range(n)] + for a in range(n): + for b in range(a + 1, n): + assert Cov[pidx[a], pidx[b]] > 1e-6, ( + f"price correlation between assets {a},{b} lost in composition" + ) + + def test_composed_sigma_is_a_valid_cholesky_factor(self): + """Result must be lower-triangular so it composes again under the same law.""" + comp = compose( + systemic_crisis_operator(severity=0.6, n=3), + rate_change_operator(change_bps=25, n=3), + ) + S = comp.Sigma_w + np.testing.assert_allclose(S, np.tril(S), atol=1e-12) + assert np.linalg.eigvalsh(self._cov(comp)).min() >= -1e-10 + + def test_composition_is_associative_in_covariance(self): + """ + Three-way composition must give the same covariance either way it is + bracketed — the property that makes the monoid claim in the module + docstring meaningful, and the one an S@S law breaks. + """ + n = 3 + a = rate_change_operator(change_bps=25, n=n) + b = systemic_crisis_operator(severity=0.7, n=n) + c = earnings_shock_operator(surprise_pct=8.0, asset_idx=1, n=n) + + left = compose(compose(a, b), c) + right = compose(a, compose(b, c)) + np.testing.assert_allclose(self._cov(left), self._cov(right), atol=1e-7) + + def test_diagonal_noise_case_is_unchanged(self): + """ + Guards against over-correction: for diagonal Sigma the old and new laws + agree, so previously-correct behaviour must not move. + """ + op1 = stock_split_operator(ratio=2.0, n=1) + op2 = dividend_operator(div_yield=0.02, asset_idx=0, n=1) + comp = compose(op1, op2) + expected = self._cov(op1) + op1.A_w @ self._cov(op2) @ op1.A_w.T + np.testing.assert_allclose(self._cov(comp), expected, atol=1e-9) + + if __name__ == "__main__": pytest.main([__file__, "-v"]) diff --git a/tests/test_noise.py b/tests/test_noise.py index aa50887..0060338 100644 --- a/tests/test_noise.py +++ b/tests/test_noise.py @@ -138,6 +138,127 @@ def test_higher_vol_gives_higher_sigma_tau(self): "Higher vol should give higher sigma_tau" +class TestCalibrationUnits: + """ + Regression tests for the per-unit-time scaling in calibrate(). + + The prior code computed `lambda_eta = N_jumps * dt`, which is the + reciprocal of a rate: it under-reported intensity by a factor of n·dt² + and — the giveaway — did not change when the sample got longer. Nothing + caught it because the only assertion on lambda_eta was `>= 0`, which a + dimensionally inverted quantity satisfies just fine. + + These tests pin the SEMANTICS (λ is a rate), not a magic number. + """ + + @staticmethod + def _planted(n_days, bars_per_day, jumps_per_day, seed=7): + """Quiet Brownian bars with `jumps_per_day` unmistakable jumps planted.""" + rng = np.random.default_rng(seed) + n = n_days * bars_per_day + r = rng.normal(0, 0.0005, n) + idx = np.linspace(0, n - 1, n_days * jumps_per_day).astype(int) + r[idx] += 0.05 + return r, len(set(idx)) + + def test_lambda_eta_recovers_planted_intensity(self): + """λ_η must come back in jumps per day, not jumps × dt.""" + bars = 78 + r, n_jumps = self._planted(n_days=10, bars_per_day=bars, jumps_per_day=3) + + params = DualNoiseCalibrator().calibrate(r, dt=1 / bars) + + # 10 days, 3 planted jumps/day → ~3/day. Detector may miss a few, so + # allow a band — but the old code returned 0.038 here, 80x low. + assert 2.0 < params.lambda_eta < 4.0, ( + f"λ_η = {params.lambda_eta:.4f}, expected ≈ {n_jumps / 10:.1f} jumps/day" + ) + + def test_lambda_eta_is_invariant_to_sample_length(self): + """A rate is a property of the process, not of how long you watched it.""" + bars = 78 + short, _ = self._planted(n_days=2, bars_per_day=bars, jumps_per_day=3) + long, _ = self._planted(n_days=20, bars_per_day=bars, jumps_per_day=3) + + cal = DualNoiseCalibrator() + lam_short = cal.calibrate(short, dt=1 / bars).lambda_eta + lam_long = cal.calibrate(long, dt=1 / bars).lambda_eta + + assert abs(lam_short - lam_long) < 1.0, ( + f"intensity drifted with sample length: {lam_short:.3f} vs {lam_long:.3f}" + ) + + def test_same_jump_count_over_longer_window_is_a_lower_rate(self): + """ + Holds jump COUNT fixed and varies elapsed time — the one comparison a + count-shaped estimator cannot fake. 6 jumps in 2 days is 3/day; + 6 jumps in 20 days is 0.3/day. `N_jumps * dt` returns the same number + for both, because neither N nor dt changed. + """ + bars = 78 + rng = np.random.default_rng(3) + + def six_jumps_over(n_days): + n = n_days * bars + r = rng.normal(0, 0.0005, n) + r[np.linspace(0, n - 1, 6).astype(int)] += 0.05 + return r + + cal = DualNoiseCalibrator() + lam_dense = cal.calibrate(six_jumps_over(2), dt=1 / bars).lambda_eta + lam_sparse = cal.calibrate(six_jumps_over(20), dt=1 / bars).lambda_eta + + assert lam_dense > 5 * lam_sparse, ( + f"same 6 jumps over 2d vs 20d gave {lam_dense:.3f} vs {lam_sparse:.3f} " + f"— λ_η is tracking count, not rate" + ) + + def test_sigma_tau_is_per_unit_time_not_per_bar(self): + """ + σ_τ shares the clock with λ_η, so it scales as BPV/(n·dt). Sampling the + SAME process at 5-min vs 1-min bars must give the same daily σ_τ; the + old per-bar estimate `sqrt(BPV/n)` differed by √5 between them. + """ + rng = np.random.default_rng(11) + daily_sigma = 0.02 + + cal = DualNoiseCalibrator() + est = {} + for bars in (78, 390): # 5-min and 1-min bars over one day + r = rng.normal(0, daily_sigma / np.sqrt(bars), bars * 5) + est[bars] = cal.calibrate(r, dt=1 / bars).sigma_tau + + np.testing.assert_allclose(est[78], est[390], rtol=0.25) + np.testing.assert_allclose(est[78], daily_sigma, rtol=0.25) + + def test_cramer_rao_bound_is_invariant_to_sampling_frequency(self): + """ + Theorem III.1 is a statement about a PROCESS, so observing that process + on 5-min vs 1-min bars must yield the same 1-day bound. This is the + end-to-end check that σ_τ² and λ_η·m₂^η share a clock: under the old + scaling σ_τ² moved by 5x and λ_η by 25x between these two samplings, + in opposite directions. + """ + cal = DualNoiseCalibrator() + bound = {} + for bars in (78, 390): + rng = np.random.default_rng(5) + n = bars * 10 + r = rng.normal(0, 0.02 / np.sqrt(bars), n) + r[np.linspace(0, n - 1, 30).astype(int)] += 0.05 # 3 jumps/day + bound[bars] = cal.calibrate(r, dt=1 / bars).cramer_rao_bound(h=1.0) + + np.testing.assert_allclose(bound[78], bound[390], rtol=0.35) + + def test_empty_and_degenerate_input(self): + """Guard the early return added alongside the rescaling.""" + cal = DualNoiseCalibrator() + for returns, dt in [(np.array([]), 1 / 78), (np.array([0.01, 0.02]), 0.0)]: + p = cal.calibrate(returns, dt=dt) + assert np.isfinite(p.sigma_tau) and np.isfinite(p.lambda_eta) + assert p.m2_eta > 0 + + class TestCramerRaoBound: """Regression tests for DualNoiseParams.cramer_rao_bound (Bug #3)."""