Fix/dual noise units and compose cholesky

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.
+
+---

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()])

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})",

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(

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"])