diff --git a/docs/experiments/2026-08-09-cts-factor-collinearity.md b/docs/experiments/2026-08-09-cts-factor-collinearity.md new file mode 100644 index 00000000..d371074c --- /dev/null +++ b/docs/experiments/2026-08-09-cts-factor-collinearity.md @@ -0,0 +1,354 @@ +# CTS factor collinearity — the suspected momentum cluster is not there, and the real defect is a different one + +**Date:** 2026-08-09 +**Issue:** #208 (from the QuantCrawler teardown #193, §1.5 — "confluence, not consensus") +**Harness:** `keel/research/cts_factors.py` — the reusable, importable, unit-tested measurement +library, alongside `research/independence.py` and `research/cscv.py`. It is library code, so it +lives in the package and ships in the wheel; `scripts/` is release/operator tooling and is +excluded from the wheel, which is the wrong home for something `tests/research/` imports. +**Script:** `docs/experiments/2026-08-09-cts-factor-collinearity.py` — the pre-declared +configuration and the report that produced every number below, and *only* that. Same split as +`2026-08-08-between-family-independence.py` (driver) against `research/independence.py` +(library), and `2026-08-05-coinbase-asset-class-probe.py`. +**Status:** measurement only. **No weight, factor, gate, threshold or rule changed.** +**Ledger:** three `diagnostic_only` rows, session `cts-factor-collinearity-2026-08-09` — +`...-unconditional-daily-`, `...-unconditional-hourly-` and `...-conditional-fired-2026-08-09`. +Each carries its own `hypotheses_tested` (45 / 36 / 45) so the multiple-testing budget is +charged per arm rather than pooled. This run influenced no shipped decision. + +**Verdict: the suspicion in #208 is refuted for momentum and confirmed, weakly, for trend.** + +The three factors the issue names as one momentum axis — `rsi_extreme`, `rsi_divergence`, +`deceleration` — have a **mean within-cluster φ of −0.018** on the primary sample, against +**+0.025** for the rest of the matrix. They are not merely uncorrelated; the sign is negative, +and the two RSI factors are **mutually exclusive**: over 6,822 observations they co-occur +**zero times** (Jaccard 0.000, lift 0.00). Sharing an input array is not the same as sharing an +answer, and this is the case where it was not. + +The trend pair *does* cluster: `condition_aligned` × `ema_fan_aligned` **φ = +0.190**, against ++0.018 for the rest of the matrix — an order of magnitude above background, stable on all five +assets (0.110 to 0.218) and reproduced at **+0.201** on a 27× larger hourly sample. It is the +one pre-declared cluster that is real. It is also small: φ = 0.190 is 3.6% shared variance. + +⭐ **The finding worth acting on is not collinearity at all.** `round_number_proximity` +(weight 1) is **present on 100.0% of BTC-USD, ETH-USD and PAXG-USD bars** and ~20% of ADA/XLM +bars. `levels.is_round_number(price, step=Decimal("0.005"))` treats `step` as an **absolute** +0.005, and Coinbase quotes those three products to two decimals — so every price is an exact +multiple of 0.005 and the check can never fail. On three of five live assets that factor is a +**constant +1 added to every CTS score**, which is worse than redundancy: a redundant factor at +least varies. + +| the question #208 asked | answer | number | +|---|---|---:| +| momentum cluster is collinear? | **no — refuted** | mean within φ **−0.018** | +| trend cluster is collinear? | yes, mildly | mean within φ **+0.190** | +| does collinearity inflate the score? | **barely** | Var ratio **1.161** (σ 7.8% wide) | +| largest collinearity anywhere in the matrix | not a pre-declared cluster | `deceleration` × `candlestick_pattern` **φ +0.254** | +| is any factor broken? | **yes, and it is not a correlation problem** | `round_number_proximity` P(present) = **1.000** on 3 of 5 assets | + +## What this is, and what it explicitly is not + +The teardown's stated rule is one indicator per category with **3-of-4 agreement**. ⚠️ **That +ratio is unvalidated SEO-adjacent content and is used here for nothing.** No threshold below is +theirs; no result is compared against it; the number 3, the number 4 and the four categories +appear in no computation. What was taken seriously is the *mechanism*, which is not theirs and +not in dispute: an **additive** score presumes its terms are close to independent, because +`total = Σ wᵢxᵢ` treats every point as new evidence. keel's CTS puts 4 of 14 raw points on three +momentum reads and another 4 on two trend reads. Whether those reads actually move together was +never measured. Now it is. + +This is a measurement, not a redesign. `DEFAULT_WEIGHTS` is read and never written; nothing in +`keel/research/cts_factors.py` is imported by the live path. + +## Method + +### There was no sample, so one had to be reconstructed + +The issue proposes dumping per-factor contributions from the `signals` table. That table has +persisted a full `cts_factors` breakdown since P3 Task 1 (`engine._persist_signal`) — and it +holds **one row** in `keel.db` and **zero** in `keel-live.db`. The live book is five daily +turtles that fire a handful of times a year, and only gate-cleared signals are ever written. A +correlation matrix from one observation is not a weak measurement, it is not a measurement. + +It is recoverable because `engine.assemble_cts_context` is a **pure function of +`(setup, candles)`** — no repo, no clock, no network — and `keel.db.candles` holds 611,176 bars +back to 2021-07-18. So the harness replays the real scoring path over history: shipped +`assemble_cts_context`, shipped `indicators_cts.score`, nothing reimplemented. + +### Both samples, because the obvious one is a collider + +| arm | population | granularity | window | N | +|---|---|---|---|---:| +| **1 (headline)** | every bar | ONE_DAY | expanding | **6,822** | +| 2 | every bar | ONE_HOUR | rolling 500 | **186,725** | +| 3 | every 20th bar, BTC only | ONE_HOUR | 250 / 500 / 1000 | 2,207 each | +| 4 | fired, gate-cleared signals | ONE_DAY | expanding | **77** | + +**N for arm 1 = 6,822** = 1,647 (BTC) + 1,647 (ETH) + 258 (PAXG) + 1,635 (ADA) + 1,635 (XLM) +daily bars, being every cached bar after a 200-bar warm-up (the EMA fan's longest period), on +the five live-allowlist assets. Each observation is one bar with a synthetic `Setup` priced at +that bar's close. Only `setup.entry` is read by the context assembly — `stop`, `target` and +`context` are not consulted — so the synthetic setup invents no risk model, it only supplies a +price for the three factors that need one. + +**Restricting to fired signals would have been selecting on the outcome.** `Rule.detect()` and +the engine's choppy / higher-TF-bias / kill-zone gates are functions of the same regime and +momentum state several factors read — the choppy gate admits a bar only when +`regime.detect_condition` is tradeable, which is most of what `condition_aligned` measures. +Conditioning on a common descendant of the variables being correlated is a collider, and it +moves the correlations by an unknown amount and sign. Arm 4 measures it anyway and shows exactly +that (below). The unconditional arm carries the conclusion. + +### The window is not a free parameter, and it is the live path's own + +Several `analysis.*` calls behind the context read the **whole list handed to them**: +`levels.find_levels` scans every pivot in the series, `regime.detect_phase` compares the last +close against `candles[0]`. Factor presence therefore depends on how much history the caller +passes. The live path — `agent.run_once` → `repo.get_candles(product, gran)` with no bounds → +`engine.evaluate` — passes the **entire cached series**, which arm 1's expanding window +reproduces exactly. Arm 2 needs a fixed window only because expanding is O(n²) and infeasible +over 44k hourly bars; arm 3 measures what that substitution costs (nothing, for these clusters). + +### Measures + +`independence.compare()` does **not** fit and is not used: its five §80.16 measurements assume +daily calendar-aligned in-market vectors plus a P&L series and derive entries from rising edges +of a held position. Factor presence is a per-bar flag with no holding period and no P&L, so +`pnl_correlation`, `entry_distances` and `median_entry_distance` have no analog. Two primitives +transfer and are imported by name: + +- **φ** — Pearson correlation of two {0,1} vectors, computed exactly from the 2×2 contingency + table (integer arithmetic and one square root, which matters at 186,725 × 36). It is + *identical* to `independence.pearson` on the same input, and + `tests/research/test_cts_factors.py::test_phi_agrees_with_shipped_pearson` asserts that on + five seeds against the shipped function as oracle. That oracle is why `_pearson` was + **deliberately promoted to `pearson`** — an oracle behind a private name is one refactor from + vanishing. +- **Jaccard** — `independence.jaccard`, unchanged. Robust where φ is not: a factor present on 2% + of bars makes the joint-absence cell dominate φ's denominator. +- **lift** — P(both)/(P(a)P(b)). The one measure whose scale does not shrink with rarity, so it + catches "rare, but always together". + +Clusters were **pre-declared** in `cts_factors.SUSPECTED_CLUSTERS`, taken verbatim from #208 +before the run, so no cluster here was found by staring at a matrix. + +### Multiple testing — controlled, and then discounted + +10 of the 11 factors vary in arm 1 (`seasonality` is weighted 0 and hardcoded `False`), giving +**45 pairwise tests**. Family-wise error is controlled by **Holm–Bonferroni at α = 0.05** over +exactly the pairs tested — not over 11·10/2, which would charge the budget for pairs no test was +run on. Arm 2 has 9 varying factors (`sr_touches` is constant, see below) and so **36 tests**. + +⚠️ **Then ignore the p-values, and this is said plainly rather than dressed up.** At N = 6,822, +φ = 0.041 clears Holm; at N = 186,725, φ = 0.054 lands at p = 1.9×10⁻¹¹⁷. **23 of 45 pairs are +"significant" in arm 1 and 12 of the 12 printed in arm 2.** Significance here answers "is this +correlation exactly zero", which nobody asked. The question asked is "is it large enough to +distort an additive score", which is an effect-size question. The correction is reported to show +the budget was accounted for, not because it carries any argument. **No conclusion in this +document rests on a p-value.** + +## Result 1 — the momentum cluster does not exist + +Arm 1, pooled, N = 6,822: + +| pair | φ | Jaccard | lift | +|---|---:|---:|---:| +| `rsi_extreme` × `rsi_divergence` | **−0.049** | **0.000** | **0.00** | +| `deceleration` × `rsi_extreme` | −0.019 | 0.015 | 0.73 | +| `deceleration` × `rsi_divergence` | +0.014 | 0.070 | 1.09 | +| **mean within cluster** | **−0.018** | 0.028 | | +| *mean over every other pair* | *+0.025* | | | + +**The within-cluster correlation is below background, and negative.** The cluster shares 4 of 14 +raw points and behaves like three unrelated factors that occasionally exclude one another. + +⭐ **`rsi_extreme` and `rsi_divergence` never co-occur — not "rarely", zero times in 6,822 +bars.** This is structural once looked at: `rsi_divergence` scores only when +`indicators.rsi_divergence(...) == "bullish"`, which requires price making a lower low while RSI +makes a higher low, and that higher low is by construction *off* the oversold extreme that +`rsi_extreme` requires. Two factors reading the same array, wired to fire in disjoint states. +Same array, opposite question. + +Consistent on every asset separately, so the pooling is not hiding a cancellation: + +| within-cluster mean φ | BTC | ETH | PAXG | ADA | XLM | +|---|---:|---:|---:|---:|---:| +| momentum | +0.007 | −0.024 | −0.032 | −0.018 | −0.036 | +| trend | +0.187 | +0.155 | +0.110 | +0.218 | +0.199 | + +And reproduced on the 27× larger hourly sample (arm 2, N = 186,725): momentum mean within φ +**−0.004**, largest-magnitude within-cluster pair **−0.011**; trend **+0.201** against +0.009 +background. + +⚠️ **This is a negative result and it is reported as one.** #208's structural argument — two +factors off the same `indicators.rsi(closes)` array, a third momentum read on the same closes — +is a correct description of the wiring and a wrong prediction about the output. The measurement +was worth making precisely because the wiring looked damning. + +## Result 2 — the trend pair is the only real cluster, and it is small + +`condition_aligned` × `ema_fan_aligned`: **φ = +0.190**, Jaccard 0.268, lift 1.47 (arm 1); +**φ = +0.201**, Jaccard 0.270, lift 1.47 (arm 2). Background for the rest of the matrix is ++0.018 and +0.009 respectively — so the pair sits **10–22× above background**, on both samples, +on every asset. + +φ = 0.190 is **3.6% shared variance**. Those two factors carry 4 of 14 raw points. That is a +real effect, it is the only pre-declared cluster that survived, and it is not large enough to +justify a code change on its own. + +## Result 3 — the strongest pair in the matrix was not predicted by anyone + +| pair | φ | Jaccard | lift | +|---|---:|---:|---:| +| **`deceleration` × `candlestick_pattern`** | **+0.254** | 0.244 | **2.10** | +| `condition_aligned` × `rsi_divergence` | −0.195 | **0.000** | **0.00** | +| `ema_fan_aligned` × `rsi_extreme` | +0.191 | 0.066 | 2.90 | +| `condition_aligned` × `ema_fan_aligned` | +0.190 | 0.268 | 1.47 | +| `in_pullback` × `fib_confluence` | +0.168 | 0.377 | 1.10 | + +`deceleration` × `candlestick_pattern` is the largest |φ| anywhere in the matrix on **all three** +unconditional configurations (0.254 daily, 0.257 hourly, 0.252 at every window length tested) — +and it crosses the very category boundary the teardown's taxonomy would draw, "momentum" against +"candlestick". Both are in fact reads of *recent bar geometry*: a decelerating leg and a +long-wicked reversal candle are two descriptions of the same exhaustion bar. **If any pair in +CTS is double-counting, it is this one, and it is not one #208 named.** It is worth 2 of 14 raw +points, and at φ = 0.254 (6.5% shared variance) it is still not large. + +⭐ **A second structural pattern shows up that additivity handles worse than redundancy: +mutual exclusion.** `condition_aligned` × `rsi_divergence` co-occur **zero** times (lift 0.00), +as do `rsi_extreme` × `rsi_divergence`. A pair that can never both score means the ceiling of +14 raw points is unreachable — the achievable maximum is lower than the nominal one, so the +`entry_technique` thresholds (`low=5`, `high=8`) sit at a different place on the real +distribution than on the nominal one. Observed mean total is **5.14** with σ 2.03, against a +nominal ceiling of 14. This is unmeasured elsewhere and is flagged, not fixed. + +## Result 4 — how much does any of this actually inflate the score? + +The single number that answers #208's worry. Additivity assumes the cross terms vanish; this is +how much they do not. `independent` is the variance the same factors would produce at the same +base rates with zero correlation (`Σ wᵢ² pᵢ(1−pᵢ)`, each factor Bernoulli): + +| arm | N | mean total | Var observed | Var independent | **ratio** | +|---|---:|---:|---:|---:|---:| +| 1 — unconditional, daily | 6,822 | 5.14 | 4.11 | 3.54 | **1.161** | +| 2 — unconditional, hourly | 186,725 | 4.05 | 3.14 | 2.88 | **1.089** | +| 4 — fired signals only | 77 | 6.97 | 1.68 | 2.54 | **0.663** | + +**Correlation widens the CTS total's standard deviation by 7.8%** (√1.161) on the primary +sample, 4.4% on the hourly one. Reading the ratio as the equicorrelation +`n_eff = n/(1+(n−1)ρ̄)` gives **≈8.6 effective independent factors out of the 10 that vary**. +⚠️ That is the standard equicorrelated approximation reported for interpretation only — it is +**not** §78.2's `N̂ = ρ̂ + (1−ρ̂)·M`, which corrects a trials count and is a different quantity, +and it is fed into no gate and no MinBTL computation. + +**8.6 of 10 is not the "confluence, not consensus" failure the issue anticipated.** + +## Result 5 — the conditional sample, and why it is not the headline + +Arm 4 drives the real `engine.evaluate` with the shipped turtle (byte-identical to +`keel-live.db` rules 1–5), gates and all: **77 gate-cleared signals** over 7,822 daily bars +(BTC 21, ETH 27, PAXG 2, ADA 14, XLM 13). Conditioning does exactly what the collider argument +predicts: + +| factor | P(present), every bar | P(present), fired | +|---|---:|---:| +| `condition_aligned` | 0.275 | **0.818** | +| `ema_fan_aligned` | 0.301 | **0.909** | +| `rsi_extreme` | 0.023 | **0.416** | +| `deceleration` | 0.174 | **0.013** | +| `candlestick_pattern` | 0.202 | **0.013** | + +The gates admit bars where the trend factors are already true and the exhaustion factors are +already false, so the variance ratio falls **below 1** (0.663) — conditioning has truncated the +distribution, not revealed independence. Every correlation in this arm is a conditional one and +none of them estimates the quantity #208 asked about. Reporting only this table would have been +selecting on the outcome. Both clusters point the same way here anyway (momentum −0.047, trend ++0.202), and only **3 of 45** pairs survive Holm at N = 77. + +⚠️ Two live factors are all but dead on the population that actually gets scored: +`deceleration` and `candlestick_pattern` each fire on **1 of 77** gate-cleared signals. They +contribute 2 of 14 nominal raw points and ~0.03 in practice. + +## Result 6 — window sensitivity, and one factor that vanishes with it + +BTC-USD hourly, every 20th bar, N = 2,207 each: + +| window | momentum mean φ | trend mean φ | Var ratio | P(`in_pullback`) | strongest pair | +|---:|---:|---:|---:|---:|---| +| 250 | −0.014 | +0.201 | 1.069 | 0.842 | `deceleration`×`candlestick_pattern` 0.252 | +| 500 | −0.014 | +0.201 | 1.065 | 0.851 | `deceleration`×`candlestick_pattern` 0.252 | +| 1000 | −0.014 | +0.201 | 1.033 | 0.857 | `deceleration`×`candlestick_pattern` 0.252 | + +Both cluster columns are **identical by construction, not by luck**: every factor in either +cluster reads a bounded lookback (`detect_condition` 20 bars, RSI 14, the fan's longest period +200), so past the warm-up the window cannot reach them. `P(in_pullback)` is the control proving +the window is applied at all — `detect_phase` compares against `candles[0]`, so it moves. + +⚠️ **`sr_touches` is present on 79.9% of daily bars and 0.0% of hourly bars under a 500-bar +window** — constant, hence dropped from arm 2's testing family entirely. `find_levels` needs +enough history to accumulate 3 distinct touches at a level, and 500 hourly bars never provide +it. A factor worth 2 of 14 points is therefore load-bearing or completely inert depending on how +much history the caller happens to pass, and the live caller passes an **ever-growing** series. +Not a correlation finding, but it fell out of this one. + +## What this changes + +**0. Nothing, in code, in this PR.** The issue is explicit that this is research and not a +drive-by tweak to the scoring weights, and no weight, factor, threshold or gate was touched. +Everything below is a recommendation or a follow-up issue. + +**1. #208's headline hypothesis is refuted and should be closed as such.** The momentum cluster +is not a cluster (mean φ −0.018, and the two RSI factors are mutually exclusive). Collapsing +`{rsi_extreme, rsi_divergence, deceleration}` to one representative — the remedy the issue +proposes — would **destroy information**, not remove duplication. Do not do it. + +**2. Recommend a follow-up issue on `levels.is_round_number`, which is the real defect this run +surfaced.** `step=Decimal("0.005")` is an absolute half-cent, not a fraction of price, so on any +product Coinbase quotes to 2 decimals the check is `price % 0.005 == 0` → always true. Verified +directly: P(present) = **1.0000** on BTC-USD, ETH-USD and PAXG-USD (2dp), 0.2195 on ADA-USD +(5dp), 0.1901 on XLM-USD (6dp). Three of five live assets receive an unconditional +1 on every +CTS score, and the same factor means something entirely different on the other two. This is a +correctness bug in a shipped analysis primitive, not a weighting question, and it is out of +scope here by the issue's own terms. + +**3. The collapse candidate, if there ever is one, is `deceleration` × `candlestick_pattern` +(φ 0.254) — not the momentum trio.** Recommend it be *recorded* and not acted on: 6.5% shared +variance over 2 of 14 raw points does not justify spending trials budget, and arm 4 shows both +factors fire on 1 of 77 real signals anyway. The cheaper question is why they almost never fire +on gate-cleared bars at all. + +**4. Additive CTS is defensible on this evidence.** ≈8.6 effective independent factors out of +10 active, 7.8% variance inflation. The teardown's critique is a correct general principle that +keel's implementation happens to survive. ⚠️ **That is a statement about independence, not about +edge.** Under §73.5 independence is necessary and never sufficient: whether any of these 11 +factors *predicts* anything is a separate, unasked and unanswered question, and nothing here +should be read as validating CTS. + +**5. `seasonality` remains structurally untestable and should stay out of every future family +count.** Weight 0 and hardcoded `False` in `assemble_cts_context` — it is not a factor scoring 0, +it is a factor that is never computed. + +## Caveats + +- In-sample, one window, no out-of-sample split, no promotion gate. A diagnostic. +- **Correlation, not causation, and not edge.** Nothing here says a factor is useful. +- Bars are **serially dependent**, and the 45 tests are not independent of one another. The + Holm correction handles the multiplicity across pairs; it does **not** correct for + autocorrelation within a series, which inflates effective N and makes every p-value optimistic. + Since no conclusion rests on a p-value this does not move the finding, but a future run + quoting significance would need a block bootstrap. Arm 3's `step=20` thinning was a cost + measure, not a remedy. +- Arm 1 pools five assets on the assumption the factor relationship is the same on each; the + per-asset table is printed for exactly that reason and the clusters agree across all five. + `round_number_proximity` is the one factor that emphatically does **not** pool, and it is + reported per-asset for that reason. +- PAXG-USD contributes 258 of 6,822 daily observations (listed 2025-05-08) and 2 of 77 fired + signals. Its per-asset column is the weakest in every table. +- Arm 4's N = 77 supports almost nothing on its own; it is here to demonstrate the selection + effect, which it does, and not to estimate a correlation. +- The synthetic setup prices entry at the bar's **close**. Three factors + (`round_number_proximity`, `fib_confluence`, `sr_touches` via `nearest_level`) are measured + against that price, so their base rates would shift under a different entry convention. The + close is the neutral choice — it is what a market order fills at — but it is a choice. +- All five assets are crypto over one broadly-correlated window. §79.16's warning about rising + cross-market correlation applies to factor structure as much as to returns. diff --git a/docs/experiments/2026-08-09-cts-factor-collinearity.py b/docs/experiments/2026-08-09-cts-factor-collinearity.py new file mode 100644 index 00000000..979112ca --- /dev/null +++ b/docs/experiments/2026-08-09-cts-factor-collinearity.py @@ -0,0 +1,358 @@ +#!/usr/bin/env python +"""Are keel's 11 CTS confluence factors collinear? -- issue #208, QuantCrawler teardown §1.5. + +Every empirical claim in `docs/experiments/2026-08-09-cts-factor-collinearity.md` comes from this +script. Strictly read-only: it opens the candle cache with `mode=ro`, drives no broker, writes +nothing but stdout, and changes no weight, gate or rule. The measurement machinery lives in +`keel/research/cts_factors.py` (importable, unit-tested); this file is only the pre-declared +configuration and the report. + +THE CRITIQUE UNDER TEST. The teardown's stated design rule for multi-indicator scoring is one +indicator per category with 3-of-4 agreement, plus a warning that stacking same-category +indicators manufactures false confirmation. ⚠️ **Their 3-of-4 ratio is unvalidated SEO-adjacent +content and is used here for exactly nothing** -- no threshold below is theirs, and no result is +compared against it. What transfers is the mechanism: an ADDITIVE score is only meaningful if its +terms are close to independent, and keel's CTS awards 4 of 14 raw points to three momentum reads +and another 4 to two trend reads. Whether that is a real problem is an empirical question about +keel's own code, and it had never been asked. + +PRE-DECLARED BEFORE THE RUN (§78.5), so nothing below is a cluster found by staring at a matrix: + + clusters `cts_factors.SUSPECTED_CLUSTERS`, taken verbatim from issue #208 + assets the live allowlist (`config.live-sandbox.yaml`) + primary arm UNCONDITIONAL, ONE_DAY, expanding window + alpha 0.05 family-wise, Holm-Bonferroni over the pairs actually tested + +FOUR ARMS: + + 1 UNCONDITIONAL / ONE_DAY / expanding window -- the headline. Expanding from the first + cached bar is not a convenience: it is exactly what the live path does + (`agent.run_once` -> `repo.get_candles(product, gran)` with no bounds -> `engine.evaluate`), + and `levels.find_levels` / `regime.detect_phase` both read the whole list they are given. + 2 UNCONDITIONAL / ONE_HOUR / rolling window -- large-N robustness, ~25x arm 1's sample. + Expanding is O(n^2) and infeasible over 44k hourly bars, so the window is fixed and + arm 3 checks what that costs. + 3 WINDOW SENSITIVITY -- arm 2's structure at three window lengths. + 4 CONDITIONAL / fired signals -- the shipped turtle through the real + `engine.evaluate`, gates and all. Reported for comparison and explicitly NOT the headline: + detection and the gates are functions of the same regime state several factors read, so + conditioning on them is a collider and moves the correlations by an unknown amount. + + .venv/bin/python docs/experiments/2026-08-09-cts-factor-collinearity.py + .venv/bin/python docs/experiments/2026-08-09-cts-factor-collinearity.py --db path/to.db +""" + +from __future__ import annotations + +import argparse +import logging +import sqlite3 +from decimal import Decimal + +from keel.research.cts_factors import ( + FACTOR_NAMES, + SUSPECTED_CLUSTERS, + FactorSample, + PairStat, + cluster_report, + holm_adjust, + pair_stats, + pool, + replay_every_bar, + replay_fired, + variance_report, +) +from keel.strategy.indicators_cts import DEFAULT_WEIGHTS +from keel.strategy.rules.turtle_breakout import TurtleBreakout +from keel.types import Candle, Granularity + +# -- PRE-DECLARED CONFIGURATION ------------------------------------------------------------------- + +DEFAULT_DB = "/Users/elmehdiaitbrahim/keel/keel.db" + +#: The live allowlist (`config.live-sandbox.yaml`), matching every prior experiment in this +#: directory so the sample is comparable across write-ups. +ASSETS = ("BTC-USD", "ETH-USD", "PAXG-USD", "ADA-USD", "XLM-USD") + +PRIMARY_GRANULARITY = Granularity.ONE_DAY +HOURLY_GRANULARITY = Granularity.ONE_HOUR + +#: Arm 2's fixed window, and arm 3's ladder around it. +HOURLY_WINDOW = 500 +SENSITIVITY_WINDOWS = (250, 500, 1000) +SENSITIVITY_STEP = 20 + +ALPHA = 0.05 + +#: Arm 4's rule: the SHIPPED turtle, byte-identical to `keel-live.db` rules 1-5 and to +#: `2026-08-08-between-family-independence.py`'s arm A. +TURTLE_PARAMS: dict[str, object] = { + "entry_lookback": 40, + "exit_lookback": 20, + "atr_period": 20, + "atr_stop_mult": Decimal("2"), + "target_rr": Decimal("6"), + "adx_period": 14, + "adx_threshold": 25.0, + "s1_filter": False, + "use_macd_confirm": False, + "min_volume_filter": False, + "volume_ma_period": 20, + "volume_mult": 1.2, +} + + +def load_candles(db_path: str, product_id: str, granularity: Granularity) -> list[Candle]: + """Ascending candles for one product/granularity, read-only. + + Same shape as `2026-08-08-between-family-independence.py`'s loader, and deliberately stdlib + `sqlite3` rather than `Repository`: an experiment must not be able to write to the cache it + is reading, and `mode=ro` is the cheapest way to make that structural. + """ + connection = sqlite3.connect(f"file:{db_path}?mode=ro", uri=True) + try: + rows = connection.execute( + "SELECT ts, o, h, l, c, v FROM candles " + "WHERE product_id = ? AND granularity = ? ORDER BY ts", + (product_id, granularity.value), + ).fetchall() + finally: + connection.close() + return [ + Candle( + ts=ts, + open=Decimal(o), + high=Decimal(h), + low=Decimal(low), + close=Decimal(c), + volume=Decimal(v), + ) + for ts, o, h, low, c, v in rows + ] + + +# -- report ----------------------------------------------------------------------------------- + + +def _fmt(value: Decimal | float, places: str = "0.001") -> str: + return str(Decimal(str(value)).quantize(Decimal(places))) + + +def print_base_rates(title: str, samples: dict[str, FactorSample], pooled: FactorSample) -> None: + print(f"\n{title} -- base rate P(factor present)") + header = f"{'factor':24} {'wt':>3} {'pooled':>8} " + " ".join( + f"{a.removesuffix('-USD'):>7}" for a in samples + ) + print(header) + print("-" * len(header)) + for name in FACTOR_NAMES: + row = f"{name:24} {DEFAULT_WEIGHTS[name]:3d} {_fmt(pooled.base_rate(name)):>8} " + row += " ".join(f"{_fmt(s.base_rate(name)):>7}" for s in samples.values()) + print(row) + print(f"{'N (observations)':24} {'':3} {pooled.n:8d} " + " ".join( + f"{s.n:7d}" for s in samples.values() + )) + + +def print_pairs(title: str, stats: list[PairStat], limit: int | None = None) -> None: + print(f"\n{title} -- pairwise co-occurrence, |phi| descending") + header = ( + f"{'factor A':24} {'factor B':24} {'phi':>7} {'jaccard':>8} {'lift':>6} " + f"{'p (Holm)':>10} {'FWER':>5}" + ) + print(header) + print("-" * len(header)) + shown = stats if limit is None else stats[:limit] + for stat in shown: + print( + f"{stat.a:24} {stat.b:24} {_fmt(stat.phi):>7} {_fmt(stat.jaccard):>8} " + f"{_fmt(stat.lift, '0.01'):>6} {stat.p_holm:10.2e} " + f"{('sig' if stat.significant else '-'):>5}" + ) + if limit is not None and len(stats) > limit: + print(f" ... {len(stats) - limit} further pairs, all |phi| <= {_fmt(shown[-1].phi)}") + + +def print_clusters(title: str, sample: FactorSample, stats: list[PairStat]) -> None: + print(f"\n{title} -- do the pre-declared clusters cluster?") + header = ( + f"{'cluster':10} {'members':>8} {'pairs':>6} {'mean phi':>9} {'max phi':>9} " + f"{'mean phi':>9} {'mean J':>7} {'raw pts':>8}" + ) + print(header) + print(f"{'':10} {'':8} {'':6} {'(within)':>9} {'(within)':>9} {'(other)':>9} {'(within)':>7} " + f"{'(share)':>8}") + print("-" * len(header)) + for report in cluster_report(sample, stats): + print( + f"{report.name:10} {len(report.members):8d} {report.within_pairs:6d} " + f"{_fmt(report.mean_within_phi):>9} {_fmt(report.max_within_phi):>9} " + f"{_fmt(report.mean_other_phi):>9} {_fmt(report.mean_within_jaccard):>7} " + f"{_fmt(report.weight_share):>8}" + ) + + +def print_variance(title: str, sample: FactorSample) -> None: + report = variance_report(sample) + print(f"\n{title} -- variance of the CTS total, observed vs. under independence") + print(f" mean CTS total {_fmt(report.mean_total, '0.01')}") + print(f" Var(total), observed {_fmt(report.observed, '0.01')}") + print(f" Var(total), independent same-marginals {_fmt(report.independent, '0.01')}") + print(f" inflation ratio {_fmt(report.ratio, '0.001')}") + + +def analyse(title: str, samples: dict[str, FactorSample], pair_limit: int | None = None) -> None: + """Print every table for one arm: base rates, pairs, clusters, variance.""" + pooled = pool(samples.values()) + stats = holm_adjust(pair_stats(pooled), alpha=ALPHA) + varying = pooled.varying() + print(f"\n{'=' * 96}\n{title}\n{'=' * 96}") + print( + f"factors with variance: {len(varying)} of {len(FACTOR_NAMES)} " + f"-> {len(stats)} pairwise tests in the family; Holm-Bonferroni at alpha={ALPHA}" + ) + constant = [n for n in FACTOR_NAMES if n not in varying] + if constant: + print( + "constant in this sample (untestable, excluded from the family): " + + ", ".join(constant) + ) + print_base_rates(title, samples, pooled) + print_pairs(title, stats, limit=pair_limit) + print_clusters(title, pooled, stats) + print_variance(title, pooled) + + +def arm1(db_path: str) -> None: + print("\nARM 1 -- UNCONDITIONAL, ONE_DAY, expanding window (the live path's own window)") + samples: dict[str, FactorSample] = {} + for product_id in ASSETS: + candles = load_candles(db_path, product_id, PRIMARY_GRANULARITY) + if not candles: + print(f" {product_id}: no {PRIMARY_GRANULARITY.value} candles cached -- skipped") + continue + samples[product_id] = replay_every_bar(product_id, candles, window=None) + if samples: + analyse("ARM 1: unconditional, ONE_DAY, expanding window", samples) + _per_asset_clusters(samples) + + +def _per_asset_clusters(samples: dict[str, FactorSample]) -> None: + """The pooled matrix assumes the relationship is the same on every asset. Check it.""" + print("\nPER-ASSET within-cluster mean phi (pooling is only honest if these agree)") + header = f"{'cluster':10} " + " ".join(f"{a.removesuffix('-USD'):>9}" for a in samples) + print(header) + print("-" * len(header)) + per_asset = { + name: holm_adjust(pair_stats(sample), alpha=ALPHA) for name, sample in samples.items() + } + for cluster in SUSPECTED_CLUSTERS: + cells = [] + for product_id, stats in per_asset.items(): + reports = {r.name: r for r in cluster_report(samples[product_id], stats)} + report = reports.get(cluster) + cells.append(_fmt(report.mean_within_phi) if report else "n/a") + print(f"{cluster:10} " + " ".join(f"{c:>9}" for c in cells)) + + +def arm2(db_path: str) -> None: + print(f"\nARM 2 -- UNCONDITIONAL, ONE_HOUR, rolling window={HOURLY_WINDOW}") + samples: dict[str, FactorSample] = {} + for product_id in ASSETS: + candles = load_candles(db_path, product_id, HOURLY_GRANULARITY) + if not candles: + continue + samples[product_id] = replay_every_bar(product_id, candles, window=HOURLY_WINDOW) + if samples: + analyse( + f"ARM 2: unconditional, ONE_HOUR, rolling window={HOURLY_WINDOW}", + samples, + pair_limit=12, + ) + + +def arm3(db_path: str) -> None: + print(f"\n{'=' * 96}") + print("ARM 3 -- window sensitivity (BTC-USD ONE_HOUR, every " + f"{SENSITIVITY_STEP}th bar)") + print(f"{'=' * 96}") + candles = load_candles(db_path, "BTC-USD", HOURLY_GRANULARITY) + if not candles: + return + header = ( + f"{'window':>7} {'N':>7} {'momentum':>10} {'trend':>10} {'var ratio':>10} " + f"{'P(pullback)':>12} {'strongest pair':>44}" + ) + print(header) + print(f"{'(bars)':>7} {'':7} {'(mean phi)':>10} {'(mean phi)':>10} {'':10} {'':12} " + f"{'(|phi|, any of the 45)':>44}") + print("-" * len(header)) + for window in SENSITIVITY_WINDOWS: + sample = replay_every_bar("BTC-USD", candles, window=window, step=SENSITIVITY_STEP) + stats = holm_adjust(pair_stats(sample), alpha=ALPHA) + reports = {r.name: r for r in cluster_report(sample, stats)} + top = stats[0] + print( + f"{window:7d} {sample.n:7d} " + f"{_fmt(reports['momentum'].mean_within_phi) if 'momentum' in reports else 'n/a':>10} " + f"{_fmt(reports['trend'].mean_within_phi) if 'trend' in reports else 'n/a':>10} " + f"{_fmt(variance_report(sample).ratio):>10} " + f"{_fmt(sample.base_rate('in_pullback')):>12} " + f"{(top.a + ' x ' + top.b + ' ' + _fmt(top.phi)):>44}" + ) + print( + "\n The three cluster columns are IDENTICAL by construction, not luck: every factor\n" + " either cluster reads a bounded lookback (`regime.detect_condition` 20 bars, RSI 14,\n" + " the EMA fan's longest period 200), so past the 200-bar warm-up the window cannot reach\n" + " them. `P(pullback)` is the control that proves the window is being applied at all --\n" + " `regime.detect_phase` compares against `candles[0]`, so it moves with window length,\n" + " and `levels.find_levels` scans every bar it is handed." + ) + + +def arm4(db_path: str) -> None: + print(f"\n{'=' * 96}") + print("ARM 4 -- CONDITIONAL: fired signals only (shipped turtle through engine.evaluate)") + print("⚠️ Conditioned on rule detection AND the choppy / higher-TF / kill-zone gates, which") + print(" read the same regime state several factors read. A collider. NOT the headline.") + print(f"{'=' * 96}") + samples: dict[str, FactorSample] = {} + for product_id in ASSETS: + candles = load_candles(db_path, product_id, PRIMARY_GRANULARITY) + if not candles: + continue + rule = TurtleBreakout(product_id=product_id, **TURTLE_PARAMS) # type: ignore[arg-type] + sample = replay_fired(rule, PRIMARY_GRANULARITY, candles, window=None) + print(f" {product_id}: {sample.n} gate-cleared signals over {len(candles)} bars") + if sample.n: + samples[product_id] = sample + if not samples: + print(" no fired signals on any asset -- nothing to correlate") + return + analyse("ARM 4: conditional on a fired, gate-cleared signal", samples, pair_limit=12) + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--db", default=DEFAULT_DB, help=f"Candle cache (default: {DEFAULT_DB})") + parser.add_argument("--arms", default="1,2,3,4", help="Comma-separated arm numbers to run") + args = parser.parse_args() + + # `engine.evaluate` emits an INFO telemetry event per non-firing rule per bar; arm 4 drives it + # ~9,000 times. Silencing INFO keeps the report readable and costs nothing measured. + logging.disable(logging.INFO) + + print("CTS factor collinearity -- issue #208 (QuantCrawler teardown §1.5)") + print(f"db={args.db} assets={','.join(ASSETS)}") + print("⚠️ RESEARCH ONLY. No weight, factor, gate or rule is changed by this run.") + print("⚠️ The teardown's '3-of-4 agreement' ratio is unvalidated and is used for nothing") + print(" here; only its mechanism (additivity presumes independence) is under test.") + + selected = {arm.strip() for arm in args.arms.split(",")} + for name, fn in (("1", arm1), ("2", arm2), ("3", arm3), ("4", arm4)): + if name in selected: + fn(args.db) + + +if __name__ == "__main__": + main() diff --git a/docs/experiments/trials-ledger.jsonl b/docs/experiments/trials-ledger.jsonl index acb5d0a6..c83bfd6f 100644 --- a/docs/experiments/trials-ledger.jsonl +++ b/docs/experiments/trials-ledger.jsonl @@ -70,3 +70,6 @@ {"decision":"diagnostic_only","kind":"ablation","params":{"arm_a":"turtle_breakout 40/20 ATR(20) 2N rr6 ADX(14)>25","arm_b":"rsi_meanrev constructor defaults, timeframe=ONE_DAY","assets":["BTC-USD","ETH-USD","PAXG-USD","ADA-USD","XLM-USD"],"fee_pct":"0.006","granularity":"ONE_DAY","measurement":"KB-80.16 five measurements","not_a_candidacy_test":"rsi_meanrev settled-refuted per KB-74.3","outcome":"degenerate - arm B took 0 trades on all 5 assets; no relationship measurable","slippage_pct":"0.0005"},"per_bar_pnl":[],"per_trade_pnl":[],"prev_hash":"ef7708484ef5a3dbfac307f09f153f9d0de6250596a3c90edbc593d901c32ea6","provenance":"a_priori","row_hash":"9e68efdd976787e92a572f9bdfee7010642334cbd8d0cbb719ffcd1123df1313","rule":"turtle_breakout|rsi_meanrev","series_missing":true,"session":"between-family-independence-2026-08-08","summary":{},"timestamp":1786185738,"trial_id":"between-family-turtle-vs-rsi-2026-08-08"} {"decision":"diagnostic_only","kind":"ablation","params":{"arm_a":"turtle_breakout 40/20","arm_b":"turtle_breakout 80/40","assets":["BTC-USD","ETH-USD","PAXG-USD"],"finding":"PR #103 used closed-trade P&L attribution, not mark-to-market","granularity":"ONE_DAY","measurement":"harness calibration vs 2026-07-20-horizon-independence.md","reproduced_closed_convention":{"BTC-USD":"0.802 vs published 0.802","ETH-USD":"0.934 vs published 0.934"}},"per_bar_pnl":[],"per_trade_pnl":[],"prev_hash":"9e68efdd976787e92a572f9bdfee7010642334cbd8d0cbb719ffcd1123df1313","provenance":"a_priori","row_hash":"eca7915eab07b63acd3503a8dfb515f6df99ab3114bd1214e03e6aed0ccbc512","rule":"turtle_breakout","series_missing":true,"session":"between-family-independence-2026-08-08","summary":{},"timestamp":1786185747,"trial_id":"harness-calibration-crosshorizon-2026-08-08"} {"decision":"diagnostic_only","kind":"ablation","params":{"alignment":"common timestamps","arms":"shipped turtle 40/20 vs ITSELF on another asset, identical params","assets":["BTC-USD","ETH-USD","PAXG-USD","ADA-USD","XLM-USD"],"cross_horizon_comparator":{"jaccard":"0.510","pnl_closed":"0.508","position":"0.585"},"finding":"asset breadth is ~an order of magnitude more independent than horizon breadth; supports the minbtl allowlist-expansion case","granularity":"ONE_DAY","mean_jaccard":"0.144","mean_pnl_corr_closed":"0.011","mean_pnl_corr_mtm":"0.056","mean_position_corr":"0.175","measurement":"KB-80.16 five measurements, CROSS-ASSET","pairs":10},"per_bar_pnl":[],"per_trade_pnl":[],"prev_hash":"eca7915eab07b63acd3503a8dfb515f6df99ab3114bd1214e03e6aed0ccbc512","provenance":"a_priori","row_hash":"5e1984ea47926c00a94b4a7da578f32deea4900b7a085199b4153a6fa2fe4dce","rule":"turtle_breakout","series_missing":true,"session":"between-family-independence-2026-08-08","summary":{},"timestamp":1786188041,"trial_id":"cross-asset-shipped-turtle-2026-08-08"} +{"decision":"diagnostic_only","kind":"ablation","params":{"arm":"1 - UNCONDITIONAL, every bar","assets":["BTC-USD","ETH-USD","PAXG-USD","ADA-USD","XLM-USD"],"factors_varying":10,"finding":"momentum cluster REFUTED (mean within phi -0.018 vs +0.025 background; the two RSI factors co-occur 0 times in 6822 bars). trend cluster CONFIRMED but small (+0.190). strongest pair in the matrix is deceleration x candlestick_pattern +0.254, not pre-declared. CTS total variance inflation 1.161 -> ~8.6 effective independent factors of 10.","granularity":"ONE_DAY","hypotheses_tested":45,"measurement":"pairwise phi/Jaccard/lift over CTS factor-presence vectors","multiple_testing":"Holm-Bonferroni, alpha=0.05, family = the 45 pairs actually tested","n_observations":6822,"no_code_change":"DEFAULT_WEIGHTS, factors, gates and thresholds all unchanged","p_values_load_bearing":false,"pre_declared_clusters":{"momentum":["rsi_extreme","rsi_divergence","deceleration"],"trend":["condition_aligned","ema_fan_aligned"]},"significant_pairs":23,"source_ratio_not_used":"QuantCrawler 3-of-4 agreement is unvalidated; used in no computation","window":"expanding from first cached bar (reproduces the live path exactly)"},"per_bar_pnl":[],"per_trade_pnl":[],"prev_hash":"5e1984ea47926c00a94b4a7da578f32deea4900b7a085199b4153a6fa2fe4dce","provenance":"a_priori","row_hash":"6d9f3f9ee286b45270f4b3527b6a54f7981f8e3915cb81d5f080619412fe5b8a","rule":"n/a (factor-level, not rule-level)","series_missing":true,"session":"cts-factor-collinearity-2026-08-09","summary":{"max_abs_phi_any_pair":"0.254","mean_cts_total":"5.14","mean_other_phi":"0.025","mean_within_phi_momentum":"-0.018","mean_within_phi_trend":"0.190","variance_ratio":"1.161"},"timestamp":1786446356,"trial_id":"cts-factor-collinearity-unconditional-daily-2026-08-09"} +{"decision":"diagnostic_only","kind":"ablation","params":{"arm":"2 - UNCONDITIONAL, every bar","assets":["BTC-USD","ETH-USD","PAXG-USD","ADA-USD","XLM-USD"],"factors_varying":9,"finding":"reproduces arm 1: momentum -0.004, trend +0.201, strongest pair still deceleration x candlestick_pattern +0.257. sr_touches is CONSTANT (0.000) under a 500-bar hourly window and was dropped from the family.","granularity":"ONE_HOUR","hypotheses_tested":36,"measurement":"replication of arm 1 at 27x the sample","multiple_testing":"Holm-Bonferroni, alpha=0.05; at this N significance is uninformative","n_observations":186725,"window":"rolling 500 bars (expanding is O(n^2) over 44k bars)"},"per_bar_pnl":[],"per_trade_pnl":[],"prev_hash":"6d9f3f9ee286b45270f4b3527b6a54f7981f8e3915cb81d5f080619412fe5b8a","provenance":"a_priori","row_hash":"40756b17a5dbc00630dfb5fa815ee03e8d5532acaa2ab09f97c9595e6a1783e4","rule":"n/a (factor-level, not rule-level)","series_missing":true,"session":"cts-factor-collinearity-2026-08-09","summary":{"mean_other_phi":"0.009","mean_within_phi_momentum":"-0.004","mean_within_phi_trend":"0.201","variance_ratio":"1.089"},"timestamp":1786446356,"trial_id":"cts-factor-collinearity-unconditional-hourly-2026-08-09"} +{"decision":"diagnostic_only","kind":"ablation","params":{"arm":"4 - CONDITIONAL (collider; reported for contrast, never as the headline)","assets":["BTC-USD","ETH-USD","PAXG-USD","ADA-USD","XLM-USD"],"finding":"demonstrates the selection effect it was run to demonstrate: P(condition_aligned) 0.275 -> 0.818, P(ema_fan_aligned) 0.301 -> 0.909, and the variance ratio falls BELOW 1 (0.663) because conditioning truncates rather than decorrelates.","granularity":"ONE_DAY","hypotheses_tested":45,"measurement":"the SAME matrix conditioned on a fired, gate-cleared signal","n_observations":77,"rule":"turtle_breakout 40/20 ATR(20) 2N rr6 ADX(14)>25 (shipped, keel-live rules 1-5)","significant_pairs":3,"underpowered":"N=77 estimates nothing on its own","window":"expanding"},"per_bar_pnl":[],"per_trade_pnl":[],"prev_hash":"40756b17a5dbc00630dfb5fa815ee03e8d5532acaa2ab09f97c9595e6a1783e4","provenance":"a_priori","row_hash":"4b196d5369087d0e8e689d9138e6a9a9d186daf1369e29158ac1a215d1221e2b","rule":"n/a (factor-level, not rule-level)","series_missing":true,"session":"cts-factor-collinearity-2026-08-09","summary":{"mean_within_phi_momentum":"-0.047","mean_within_phi_trend":"0.202","variance_ratio":"0.663"},"timestamp":1786446356,"trial_id":"cts-factor-collinearity-conditional-fired-2026-08-09"} diff --git a/keel/research/cts_factors.py b/keel/research/cts_factors.py new file mode 100644 index 00000000..67d5a889 --- /dev/null +++ b/keel/research/cts_factors.py @@ -0,0 +1,530 @@ +"""Do keel's 11 CTS confluence factors carry independent evidence? -- issue #208. + +CTS (`strategy/indicators_cts.py`) is an ADDITIVE score: 14 raw points spread over 11 named +factors, summed, and mapped to one of three entry techniques. Additivity is only defensible if +the factors are approximately independent. If three of them are three reads of the same momentum +series, then when they agree they agree BY CONSTRUCTION, and the score reports "strong +confluence" on what is really one piece of evidence counted three times. That is the +"confluence, not consensus" critique from the QuantCrawler teardown (#193 §1.5), turned into a +measurement. + +⚠️ **This module measures. It does not score, gate, or change anything.** Nothing here is +imported by the live path; `DEFAULT_WEIGHTS` is read, never written. + +WHY THIS NEEDED A REPLAY HARNESS AT ALL. The obvious data source is the `signals` table, which +has persisted a per-factor `cts_factors` breakdown since P3 Task 1 (`engine._persist_signal`). +It holds **one row** in `keel.db` and **zero** in `keel-live.db` -- the live book's rules are +daily turtles that fire a handful of times a year, and the audit trail only records signals that +cleared every gate. A correlation matrix needs thousands of observations, so the sample has to +be reconstructed. It can be: `engine.assemble_cts_context` is a pure function of +`(setup, candles)`, and `keel.db.candles` holds ~611k bars back to 2021. + +TWO SAMPLES, BOTH REPORTED, BECAUSE THE OBVIOUS ONE IS BIASED: + +- `replay_every_bar` -- UNCONDITIONAL. Every bar in the series gets a synthetic `Setup` priced + at that bar's close, and the full context is assembled. This is the population the factors + actually live in. +- `replay_fired` -- CONDITIONAL. Only bars where a real `Rule.detect()` produced a `Setup` that + then cleared `engine.evaluate`'s choppy / higher-TF-bias / kill-zone gates. This is the exact + population the live scorer sees. + +The conditional sample is the tempting one ("correlate the factors on real signals"), and it is +the wrong one to draw the headline from. Rule detection and the engine's gates are functions of +the same regime and momentum state the factors read -- the choppy gate in particular admits a +bar only when `regime.detect_condition` is tradeable, which is most of what `condition_aligned` +measures. Conditioning on it truncates the joint distribution and moves every correlation +involving those factors by an unknown amount and sign (Berkson's paradox / collider +stratification). Reporting only the conditional matrix would be selecting on the outcome. Both +run; the unconditional one carries the conclusion. + +THE WINDOW MATTERS, AND IT IS NOT A FREE PARAMETER. Several `analysis.*` calls behind the +context read the WHOLE list handed to them rather than a fixed lookback -- `levels.find_levels` +scans every pivot in the series, and `regime.detect_phase` compares the last close against +`candles[0]`. So factor presence depends on how much history the caller passes. The live path +(`agent.run_once` -> `repo.get_candles(product, gran)` with no bounds -> `engine.evaluate`) +passes the entire cached series, which `window=None` reproduces exactly (expanding window from +the first cached bar). A fixed `window` is offered because expanding is O(n^2) and infeasible on +44k hourly bars; when one is used, say so and check the answer against a second value. + +MEASURES. For two boolean factor vectors the pairwise statistics are: + + phi Pearson correlation of the two {0,1} vectors, computed exactly from the 2x2 + contingency counts. Identical to `research.independence.pearson` on the same input + -- which is what `tests/research/test_cts_factors.py` asserts, using the shipped + function as the oracle rather than trusting this module's arithmetic. + Jaccard `research.independence.jaccard`, unchanged: #{both present} / #{either present}. + Robust where phi is not, because a factor present on 2% of bars makes the + joint-absence cell dominate phi's denominator. + lift P(both) / (P(a)P(b)). 1.0 is independence; it is the one measure whose scale does + not shrink with rarity, so it catches "these two are rare, but always co-occur". + +`independence.compare()` is deliberately NOT used. Its five §80.16 measurements assume daily +calendar-aligned in-market vectors plus a P&L series, and derive entries from the rising edges +of a held position. Factor presence is a per-bar flag with no holding period and no P&L, so +`pnl_correlation`, `entry_distances` and `median_entry_distance` have no analog. Only the two +primitives above transfer, and they are imported by name. +""" + +from __future__ import annotations + +import math +from collections.abc import Iterable, Sequence +from dataclasses import dataclass, field +from decimal import Decimal, localcontext + +from keel.research.independence import jaccard +from keel.strategy import engine, indicators_cts +from keel.strategy.rules.base import Rule, Setup +from keel.types import Candle, Granularity + +__all__ = [ + "FACTOR_NAMES", + "SUSPECTED_CLUSTERS", + "ClusterReport", + "FactorSample", + "PairStat", + "VarianceReport", + "cluster_report", + "contingency", + "factor_presence", + "holm_adjust", + "pair_stats", + "phi", + "pool", + "replay_every_bar", + "replay_fired", + "variance_report", +] + +#: The 11 factor names, in `DEFAULT_WEIGHTS` order. Read from the shipped table so a factor +#: added or removed there shows up here without an edit -- this module must never carry its own +#: copy of the scorer's shape. +FACTOR_NAMES: tuple[str, ...] = tuple(indicators_cts.DEFAULT_WEIGHTS) + +#: The clusters issue #208 suspects, PRE-DECLARED before the run (§78.5) so the result cannot be +#: a cluster discovered by staring at the matrix. `momentum` is three reads downstream of one +#: price/velocity series -- `rsi_extreme` and `rsi_divergence` share the very same +#: `indicators.rsi(closes)` array (`engine.assemble_cts_context`), and `deceleration` is a third +#: momentum read on the same closes. `trend` is two structure reads. 4 of 14 raw points each. +SUSPECTED_CLUSTERS: dict[str, tuple[str, ...]] = { + "momentum": ("rsi_extreme", "rsi_divergence", "deceleration"), + "trend": ("condition_aligned", "ema_fan_aligned"), +} + +#: Bars discarded at the head of every replay. The context's slowest input is the EMA fan +#: (`indicators.ema_fan`, longest period 200), so a shorter warm-up would score the opening +#: stretch of every series against a fan that has not converged. +DEFAULT_WARMUP = 200 + + +# --------------------------------------------------------------------------- +# The sample +# --------------------------------------------------------------------------- + + +@dataclass +class FactorSample: + """Per-bar presence vectors for every CTS factor, plus the CTS total on each bar. + + `vectors[name][i]` is 1 when factor `name` was present on observation `i`. Every vector has + length `n`; `totals[i]` is the score `indicators_cts.score` awarded that observation under + `DEFAULT_WEIGHTS`, kept so variance can be measured on the real total rather than on a + reconstruction of it. + """ + + vectors: dict[str, list[int]] = field(default_factory=dict) + totals: list[int] = field(default_factory=list) + labels: list[str] = field(default_factory=list) + + @property + def n(self) -> int: + return len(self.totals) + + def base_rate(self, name: str) -> Decimal: + """Fraction of observations on which `name` was present; 0 for an empty sample.""" + vector = self.vectors.get(name, []) + return Decimal(sum(vector)) / Decimal(len(vector)) if vector else Decimal(0) + + def varying(self) -> tuple[str, ...]: + """Factor names that are neither always-present nor always-absent in this sample. + + A constant vector has no variance, so every correlation involving it is undefined and + `phi` returns 0 by convention. Those pairs must be excluded from the multiple-testing + budget rather than counted as 45 tests of which several were structurally guaranteed to + return zero. `seasonality` is always in here: it is weighted 0 and hardcoded `False`. + """ + out = [] + for name in FACTOR_NAMES: + vector = self.vectors.get(name, []) + if vector and 0 < sum(vector) < len(vector): + out.append(name) + return tuple(out) + + +def pool(samples: Iterable[FactorSample]) -> FactorSample: + """Concatenate per-asset samples into one. + + Pooling across assets is safe in a way pooling P&L is not: presence is a dimensionless + boolean, so there is no units problem (cf. `backtest()`'s per-unit quote-currency P&L, which + `2026-08-08-between-family-independence.md` refuses to pool). It does assume the factor + relationship is the same across assets, which is why the per-asset matrices are printed too. + """ + merged = FactorSample(vectors={name: [] for name in FACTOR_NAMES}) + for sample in samples: + for name in FACTOR_NAMES: + merged.vectors[name].extend(sample.vectors.get(name, [])) + merged.totals.extend(sample.totals) + merged.labels.extend(sample.labels) + return merged + + +# --------------------------------------------------------------------------- +# Replay +# --------------------------------------------------------------------------- + + +def factor_presence(setup: Setup, candles: list[Candle]) -> tuple[dict[str, bool], int]: + """`({factor name: present}, cts total)` for one `(setup, candles)` pair. + + Delegates to the SHIPPED `engine.assemble_cts_context` + `indicators_cts.score` -- nothing + about the scoring path is reimplemented here, so a change to either is a change to this + measurement, which is the point. + """ + context = engine.assemble_cts_context(setup, candles) + result = indicators_cts.score(context) + return {f.name: f.present for f in result.factors}, result.total + + +def _synthetic_setup(product_id: str, candle: Candle) -> Setup: + """A `Setup` standing in for "some rule proposed an entry at this bar's close". + + `assemble_cts_context` reads ONLY `setup.entry` (documented in its docstring), so the stop + and target below are structurally inert -- they exist because `Setup` requires them. They + are set to a plain 2:1 so the object is not a nonsense risk profile if anything downstream + ever looks. Using the close as the entry is the neutral choice: it is the price a + market-order rule would get on that bar, and it is what `round_number_proximity`, + `fib_confluence` and the nearest-S/R-level lookup are measured against. + """ + close = candle.close + return Setup( + product_id=product_id, + direction="long", + entry=close, + stop=close * Decimal("0.98"), + target=close * Decimal("1.04"), + context={}, + ts=candle.ts, + ) + + +def _window(candles: list[Candle], index: int, window: int | None) -> list[Candle]: + """History visible at bar `index`: expanding from bar 0 when `window` is None.""" + start = 0 if window is None else max(0, index - window + 1) + return candles[start : index + 1] + + +def replay_every_bar( + product_id: str, + candles: list[Candle], + *, + warmup: int = DEFAULT_WARMUP, + window: int | None = None, + step: int = 1, +) -> FactorSample: + """UNCONDITIONAL sample: assemble the CTS context on every bar from `warmup` onward. + + No rule is consulted and no gate is applied, so the resulting matrix describes the factors + themselves rather than the factors as filtered by whatever fired. `step > 1` subsamples the + bar index (thinning autocorrelation and cost at once); `window` is documented on the module. + """ + sample = FactorSample(vectors={name: [] for name in FACTOR_NAMES}) + for index in range(warmup, len(candles), step): + visible = _window(candles, index, window) + presence, total = factor_presence(_synthetic_setup(product_id, candles[index]), visible) + for name in FACTOR_NAMES: + sample.vectors[name].append(1 if presence.get(name) else 0) + sample.totals.append(total) + sample.labels.append(f"{product_id}@{candles[index].ts}") + return sample + + +def replay_fired( + rule: Rule, + granularity: Granularity, + candles: list[Candle], + *, + warmup: int = DEFAULT_WARMUP, + window: int | None = None, + step: int = 1, +) -> FactorSample: + """CONDITIONAL sample: only bars on which `engine.evaluate` actually emitted a `Signal`. + + Drives the real `engine.evaluate` rather than reimplementing its gate ladder, so the choppy, + higher-TF-bias and kill-zone gates apply exactly as they do in production, and a DCA-class + setup is exempted from them exactly as it is in production. The context is then re-assembled + from the emitted signal's setup -- the same move `sim.portfolio_sim._record_cts_telemetry` + makes, and necessary because `Signal` carries the CTS total but not the per-factor breakdown. + + ⚠️ See the module docstring: this sample is conditioned on quantities correlated with the + factors being measured. It is reported for completeness and comparison, never alone. + """ + sample = FactorSample(vectors={name: [] for name in FACTOR_NAMES}) + for index in range(warmup, len(candles), step): + visible = _window(candles, index, window) + for signal in engine.evaluate([rule], {granularity: visible}): + if signal.setup is None: + continue + presence, total = factor_presence(signal.setup, visible) + for name in FACTOR_NAMES: + sample.vectors[name].append(1 if presence.get(name) else 0) + sample.totals.append(total) + sample.labels.append(f"{rule.name}@{candles[index].ts}") + return sample + + +# --------------------------------------------------------------------------- +# Pairwise statistics +# --------------------------------------------------------------------------- + + +def contingency(xs: Sequence[int], ys: Sequence[int]) -> tuple[int, int, int, int]: + """The 2x2 table `(n11, n10, n01, n00)` over two aligned boolean vectors.""" + n11 = n10 = n01 = n00 = 0 + for x, y in zip(xs, ys): + if x: + if y: + n11 += 1 + else: + n10 += 1 + elif y: + n01 += 1 + else: + n00 += 1 + return n11, n10, n01, n00 + + +def phi(xs: Sequence[int], ys: Sequence[int]) -> Decimal: + """Phi coefficient -- Pearson correlation of two {0,1} vectors, from the 2x2 table. + + Exact integer arithmetic up to one square root, and O(n) with no Decimal work in the loop, + which matters at ~200k observations x 45 pairs. Returns 0 when either vector is constant, + matching `independence.pearson`'s documented convention for a variance-free series (and for + the same reason: one degenerate factor must not abort a whole matrix). + """ + n11, n10, n01, n00 = contingency(xs, ys) + with localcontext() as ctx: + ctx.prec = 50 + numerator = Decimal(n11 * n00 - n10 * n01) + product = Decimal((n11 + n10) * (n01 + n00) * (n11 + n01) * (n10 + n00)) + if product <= 0: + return Decimal(0) + return +(numerator / product.sqrt()) + + +def _two_sided_p(correlation: Decimal, n: int) -> float: + """p-value for `phi` under H0 of independence, from the chi-square-of-1-df identity. + + `chi2 = n * phi^2` for a 2x2 table, and `sqrt(chi2_1)` is a standard normal deviate, so the + two-sided p is `erfc(|phi| * sqrt(n/2))`. Asymptotic; at the sample sizes here the normal + approximation is not what limits the conclusion (see `holm_adjust`). + """ + if n < 2: + return 1.0 + return math.erfc(abs(float(correlation)) * math.sqrt(n / 2)) + + +@dataclass(frozen=True) +class PairStat: + """One factor pair's co-occurrence, with the multiple-testing bookkeeping attached.""" + + a: str + b: str + n: int + n11: int + n10: int + n01: int + n00: int + phi: Decimal + jaccard: Decimal + lift: Decimal + p_value: float + p_holm: float = 1.0 + significant: bool = False + + +def pair_stats(sample: FactorSample, names: Sequence[str] | None = None) -> list[PairStat]: + """Every unordered pair of `names` (default: the sample's non-constant factors). + + Ordered by descending |phi|, which is a presentation choice only -- the multiple-testing + correction in `holm_adjust` is applied over the whole family regardless of print order. + """ + chosen = tuple(names) if names is not None else sample.varying() + out: list[PairStat] = [] + for index, first in enumerate(chosen): + for second in chosen[index + 1 :]: + xs = sample.vectors.get(first, []) + ys = sample.vectors.get(second, []) + n11, n10, n01, n00 = contingency(xs, ys) + n = n11 + n10 + n01 + n00 + correlation = phi(xs, ys) + expected = Decimal((n11 + n10) * (n11 + n01)) + lift = (Decimal(n11) * Decimal(n)) / expected if expected > 0 else Decimal(0) + out.append( + PairStat( + a=first, + b=second, + n=n, + n11=n11, + n10=n10, + n01=n01, + n00=n00, + phi=correlation, + jaccard=jaccard(xs, ys), + lift=lift, + p_value=_two_sided_p(correlation, n), + ) + ) + return sorted(out, key=lambda stat: abs(stat.phi), reverse=True) + + +def holm_adjust(stats: Sequence[PairStat], alpha: float = 0.05) -> list[PairStat]: + """Holm-Bonferroni step-down over the whole family of pairwise tests. + + Holm rather than plain Bonferroni: it controls the same family-wise error rate and is + uniformly more powerful, and there is no reason to accept the weaker test. The family is + every pair TESTED, `k(k-1)/2` over the non-constant factors -- not `11*10/2`, which would + charge the budget for pairs no test was run on. + + ⚠️ Read the honest caveat in the write-up before quoting any of these. At n in the tens or + hundreds of thousands, `phi = 0.01` clears any threshold; significance here answers "is this + correlation exactly zero", which nobody asked, while the question actually asked is "is it + big enough to distort an additive score". That is an effect-size question, and the p-values + are reported to show the correction was applied, not to carry the argument. + """ + ordered = sorted(stats, key=lambda stat: stat.p_value) + total = len(ordered) + adjusted: list[PairStat] = [] + running = 0.0 + still_rejecting = True + for index, stat in enumerate(ordered): + running = max(running, min(1.0, stat.p_value * (total - index))) + still_rejecting = still_rejecting and running <= alpha + adjusted.append( + PairStat( + **{ + **stat.__dict__, + "p_holm": running, + "significant": still_rejecting, + } + ) + ) + return sorted(adjusted, key=lambda stat: abs(stat.phi), reverse=True) + + +# --------------------------------------------------------------------------- +# Cluster + variance summaries +# --------------------------------------------------------------------------- + + +@dataclass(frozen=True) +class ClusterReport: + """Within-cluster vs. rest-of-matrix correlation for one pre-declared cluster.""" + + name: str + members: tuple[str, ...] + within_pairs: int + mean_within_phi: Decimal + max_within_phi: Decimal + mean_other_phi: Decimal + mean_within_jaccard: Decimal + weight_share: Decimal + + +def cluster_report( + sample: FactorSample, + stats: Sequence[PairStat], + clusters: dict[str, tuple[str, ...]] | None = None, + weights: dict[str, int] | None = None, +) -> list[ClusterReport]: + """Does each pre-declared cluster actually cluster? + + The comparison that matters is `mean_within_phi` against `mean_other_phi` (every pair with + at most one member in the cluster). A cluster whose internal correlation is indistinguishable + from the rest of the matrix is not a cluster, however plausible its wiring diagram looks. + Means are taken over SIGNED phi: taking |phi| would report a cluster whose members + systematically EXCLUDE one another as tightly coupled, which is a different finding and must + not be laundered into this one. `max_within_phi` is signed for the same reason. + """ + clusters = SUSPECTED_CLUSTERS if clusters is None else clusters + weights = indicators_cts.DEFAULT_WEIGHTS if weights is None else weights + total_weight = sum(weights.values()) + out: list[ClusterReport] = [] + + for name, members in clusters.items(): + inside = [s for s in stats if s.a in members and s.b in members] + outside = [s for s in stats if not (s.a in members and s.b in members)] + if not inside: + continue + out.append( + ClusterReport( + name=name, + members=members, + within_pairs=len(inside), + mean_within_phi=_mean(s.phi for s in inside), + max_within_phi=max(inside, key=lambda s: abs(s.phi)).phi, + mean_other_phi=_mean(s.phi for s in outside), + mean_within_jaccard=_mean(s.jaccard for s in inside), + weight_share=( + Decimal(sum(weights.get(m, 0) for m in members)) / Decimal(total_weight) + if total_weight + else Decimal(0) + ), + ) + ) + # `sample` is accepted for symmetry with the other reporters and to keep call sites uniform; + # every figure above is derivable from `stats` alone. + del sample + return out + + +@dataclass(frozen=True) +class VarianceReport: + """How much correlation inflates the spread of the CTS total. + + `observed` is the sample variance of the real `indicators_cts.score` totals. `independent` + is what that variance would be if the factors were mutually independent with the SAME base + rates -- `sum(w_i^2 * p_i * (1 - p_i))`, since each factor is Bernoulli. Their ratio is the + single number that answers issue #208's actual worry: additivity assumes the cross terms + vanish, and this is how much they do not. + """ + + observed: Decimal + independent: Decimal + ratio: Decimal + mean_total: Decimal + + +def variance_report(sample: FactorSample, weights: dict[str, int] | None = None) -> VarianceReport: + """Observed vs. independence-implied variance of the CTS total.""" + weights = indicators_cts.DEFAULT_WEIGHTS if weights is None else weights + n = sample.n + if n < 2: + return VarianceReport(Decimal(0), Decimal(0), Decimal(0), Decimal(0)) + + mean = Decimal(sum(sample.totals)) / Decimal(n) + observed = sum(((Decimal(t) - mean) ** 2 for t in sample.totals), Decimal(0)) / Decimal(n - 1) + + independent = Decimal(0) + for name, weight in weights.items(): + rate = sample.base_rate(name) + independent += Decimal(weight) ** 2 * rate * (Decimal(1) - rate) + + ratio = observed / independent if independent > 0 else Decimal(0) + return VarianceReport( + observed=observed, independent=independent, ratio=ratio, mean_total=mean + ) + + +def _mean(values: Iterable[Decimal]) -> Decimal: + collected = list(values) + return sum(collected, Decimal(0)) / Decimal(len(collected)) if collected else Decimal(0) diff --git a/keel/research/independence.py b/keel/research/independence.py index 082bc1b9..55918fd9 100644 --- a/keel/research/independence.py +++ b/keel/research/independence.py @@ -43,12 +43,21 @@ class IndependenceReport: median_entry_distance: int | None -def _pearson(xs: Sequence[Decimal], ys: Sequence[Decimal]) -> Decimal: +def pearson(xs: Sequence[Decimal], ys: Sequence[Decimal]) -> Decimal: """Pearson correlation. Returns 0 when either series is constant. A constant series has no variance to correlate against, so the coefficient is undefined rather than zero -- but every caller here wants "no measurable relationship", and raising would abort a whole matrix over one degenerate pair. Documented, not silent. + + **Public deliberately (issue #208).** It was `_pearson` while `compare()` was its only + caller. `research/cts_factors.py` measures correlation between BOOLEAN factor-presence + vectors, where the fast exact form is the phi coefficient computed from 2x2 contingency + counts rather than a running covariance -- but phi and Pearson-on-{0,1} are the same + number, so this function is that module's test ORACLE. An oracle imported through a + private name is an oracle one refactor away from silently disappearing, and importing + `_pearson` across modules is exactly the wart `sim/portfolio_sim.py` already carries + against `engine._assemble_cts_context`. Promoted rather than duplicated. """ n = len(xs) if n < 2 or n != len(ys): @@ -119,10 +128,10 @@ def compare( b_active=sum(b_positions), both_active=sum(1 for x, y in zip(a_positions, b_positions) if x and y), jaccard=jaccard(a_positions, b_positions), - position_correlation=_pearson( + position_correlation=pearson( [Decimal(v) for v in a_positions], [Decimal(v) for v in b_positions] ), - pnl_correlation=_pearson(a_pnl, b_pnl), + pnl_correlation=pearson(a_pnl, b_pnl), entry_distances=distances, median_entry_distance=_median(distances), ) diff --git a/keel/sim/portfolio_sim.py b/keel/sim/portfolio_sim.py index b0e36581..95af50b4 100644 --- a/keel/sim/portfolio_sim.py +++ b/keel/sim/portfolio_sim.py @@ -29,7 +29,7 @@ - **`cts_factor_populated` / `rejected_for_missing_input`** (both `dict[str, int]`, keyed by CTS context-key name -- see `strategy.indicators_cts.DEFAULT_WEIGHTS`): for *every* ENTER signal `evaluate()` emits (whether or not it ends up opened), the engine's own CTS context assembly is - reused verbatim (`engine._assemble_cts_context` + `indicators_cts.score`, not reimplemented) to + reused verbatim (`engine.assemble_cts_context` + `indicators_cts.score`, not reimplemented) to determine, per factor, whether it was present or absent on that bar. Present factors increment `cts_factor_populated[name]`; absent ones increment `rejected_for_missing_input[name]`. This is symmetric by construction (`populated[k] + missing[k]` == the number of signals a given `k` was @@ -483,7 +483,7 @@ def _record_cts_telemetry( setup = signal.setup if setup is None: continue - cts_result = indicators_cts.score(engine._assemble_cts_context(setup, window_1h)) + cts_result = indicators_cts.score(engine.assemble_cts_context(setup, window_1h)) for factor in cts_result.factors: bucket = ( telemetry.cts_factor_populated diff --git a/keel/strategy/engine.py b/keel/strategy/engine.py index 3b3f04a1..68fe1a52 100644 --- a/keel/strategy/engine.py +++ b/keel/strategy/engine.py @@ -139,7 +139,7 @@ def evaluate( ) continue - context = _assemble_cts_context(setup, trading_candles) + context = assemble_cts_context(setup, trading_candles) cts_result = indicators_cts.score(context, weights) technique = indicators_cts.entry_technique(cts_result.total) @@ -242,7 +242,7 @@ def _trading_granularity(rule: Rule, candles_by_tf: dict[Granularity, list[Candl # --------------------------------------------------------------------------- -def _assemble_cts_context(setup: Setup, candles: list[Candle]) -> dict[str, Any]: +def assemble_cts_context(setup: Setup, candles: list[Candle]) -> dict[str, Any]: """Build the 11-key `indicators_cts.score()` context from `analysis.*` + `setup`. Every flag is computed fresh from `analysis.regime`/`analysis.indicators`/ @@ -251,6 +251,15 @@ def _assemble_cts_context(setup: Setup, candles: list[Candle]) -> dict[str, Any] checks -- not lifted verbatim from `setup.context`, whose keys are each rule's own explainability shape (e.g. `pattern` vs. `candlestick_pattern`) and aren't guaranteed to line up with the CTS context's documented key names (see `indicators_cts.py`'s docstring). + + **Public deliberately (issue #208).** This is a PURE function of `(setup, candles)` -- + no I/O, no repo, no clock -- which makes it the one piece of the live scoring path that + can be replayed offline over historical candles. Two modules outside `strategy/` do + exactly that: `sim/portfolio_sim._record_cts_telemetry` (which was already reaching + through the `_` prefix) and `research/cts_factors`, which replays it bar-by-bar to test + whether the 11 factors carry independent evidence (issue #208). Only `setup.entry` is + read from the setup; `stop`/`target`/`context` are not consulted, so an offline replay + can present a synthetic entry price without inventing a risk model. """ if len(candles) < 2: return {} diff --git a/tests/research/test_cts_factors.py b/tests/research/test_cts_factors.py new file mode 100644 index 00000000..8f640fb9 --- /dev/null +++ b/tests/research/test_cts_factors.py @@ -0,0 +1,333 @@ +"""CTS factor co-occurrence harness (issue #208). + +The load-bearing test is `test_phi_agrees_with_shipped_pearson`: `phi` is a fast contingency-count +reimplementation of Pearson on {0,1} vectors, and the only reason to trust it is that it agrees +with `independence.pearson`, which `compare()` has shipped since PR #103. That is the oracle, and +it is why `pearson` was promoted out of `_pearson`. +""" + +from __future__ import annotations + +import random +from decimal import Decimal + +import pytest + +from keel.research.cts_factors import ( + FACTOR_NAMES, + SUSPECTED_CLUSTERS, + FactorSample, + cluster_report, + contingency, + factor_presence, + holm_adjust, + pair_stats, + phi, + pool, + replay_every_bar, + replay_fired, + variance_report, +) +from keel.research.independence import pearson +from keel.strategy.indicators_cts import DEFAULT_WEIGHTS +from keel.strategy.rules.base import Rule, Setup +from keel.types import Candle, Granularity + + +def _candles(count: int, seed: int = 7) -> list[Candle]: + """A pseudo-random but deterministic OHLCV walk, 2dp like real Coinbase USD candles.""" + rng = random.Random(seed) + price = 100.0 + out: list[Candle] = [] + for index in range(count): + price = max(1.0, price * (1 + rng.uniform(-0.03, 0.03))) + high = price * (1 + abs(rng.uniform(0, 0.02))) + low = price * (1 - abs(rng.uniform(0, 0.02))) + out.append( + Candle( + ts=1_600_000_000 + index * 86_400, + open=Decimal(f"{rng.uniform(low, high):.2f}"), + high=Decimal(f"{high:.2f}"), + low=Decimal(f"{low:.2f}"), + close=Decimal(f"{price:.2f}"), + volume=Decimal("1000"), + ) + ) + return out + + +def _sample(vectors: dict[str, list[int]]) -> FactorSample: + """A `FactorSample` from explicit vectors, with totals implied by `DEFAULT_WEIGHTS`.""" + length = len(next(iter(vectors.values()))) + filled = {name: vectors.get(name, [0] * length) for name in FACTOR_NAMES} + totals = [ + sum(DEFAULT_WEIGHTS[name] * filled[name][i] for name in FACTOR_NAMES) + for i in range(length) + ] + n = length + return FactorSample(vectors=filled, totals=totals, labels=[str(i) for i in range(n)]) + + +# -- phi, against the shipped oracle --------------------------------------------------------- + + +@pytest.mark.parametrize("seed", [1, 2, 3, 11, 101]) +def test_phi_agrees_with_shipped_pearson(seed: int) -> None: + """`phi` from 2x2 counts == `independence.pearson` on the same {0,1} vectors.""" + rng = random.Random(seed) + # Correlate b with a so the pairs are not all near-zero: 70% copy, 30% independent coin. + xs = [rng.randint(0, 1) for _ in range(400)] + ys = [x if rng.random() < 0.7 else rng.randint(0, 1) for x in xs] + + oracle = pearson([Decimal(v) for v in xs], [Decimal(v) for v in ys]) + assert phi(xs, ys) == pytest.approx(oracle, abs=Decimal("1e-20")) + + +def test_phi_is_one_for_identical_vectors_and_minus_one_for_inverted() -> None: + xs = [1, 0, 1, 1, 0, 0, 1] + assert phi(xs, xs) == pytest.approx(Decimal(1), abs=Decimal("1e-20")) + assert phi(xs, [1 - x for x in xs]) == pytest.approx(Decimal(-1), abs=Decimal("1e-20")) + + +def test_phi_returns_zero_for_a_constant_vector() -> None: + """Documented convention, shared with `independence.pearson`: no variance -> 0, not a raise.""" + assert phi([1, 1, 1, 1], [0, 1, 0, 1]) == Decimal(0) + assert phi([0, 0, 0, 0], [0, 0, 0, 0]) == Decimal(0) + + +def test_contingency_counts_all_four_cells() -> None: + n11, n10, n01, n00 = contingency([1, 1, 0, 0, 1], [1, 0, 1, 0, 0]) + assert (n11, n10, n01, n00) == (1, 2, 1, 1) + + +# -- sample bookkeeping ---------------------------------------------------------------------- + + +def test_varying_excludes_constant_factors_from_the_testing_family() -> None: + sample = _sample( + { + "condition_aligned": [1, 0, 1, 0], + "ema_fan_aligned": [1, 1, 0, 0], + "rsi_extreme": [0, 0, 0, 0], # constant absent + "in_pullback": [1, 1, 1, 1], # constant present + } + ) + varying = sample.varying() + assert set(varying) == {"condition_aligned", "ema_fan_aligned"} + # 2 varying factors -> exactly 1 pair charged to the multiple-testing budget, not 55. + assert len(pair_stats(sample)) == 1 + + +def test_base_rate_and_pooling() -> None: + a = _sample({"condition_aligned": [1, 1, 0, 0]}) + b = _sample({"condition_aligned": [1, 1, 1, 1]}) + merged = pool([a, b]) + assert merged.n == 8 + assert merged.base_rate("condition_aligned") == Decimal("0.75") + assert len(merged.labels) == 8 + + +def test_pair_stats_reports_mutual_exclusion_as_zero_jaccard_and_zero_lift() -> None: + sample = _sample( + {"condition_aligned": [1, 1, 0, 0, 0, 0], "rsi_divergence": [0, 0, 1, 1, 0, 0]} + ) + (stat,) = pair_stats(sample) + assert stat.n11 == 0 + assert stat.jaccard == Decimal(0) + assert stat.lift == Decimal(0) + assert stat.phi < Decimal(0) + + +# -- multiple testing ------------------------------------------------------------------------ + + +def test_holm_is_step_down_and_monotone() -> None: + sample = _sample( + { + "condition_aligned": [1, 0] * 50, + "ema_fan_aligned": [1, 0] * 50, # identical -> tiny p + "rsi_extreme": [1, 1, 0, 0] * 25, # independent-ish -> large p + "in_pullback": [0, 1, 1, 0] * 25, + } + ) + stats = holm_adjust(pair_stats(sample), alpha=0.05) + by_raw = sorted(stats, key=lambda s: s.p_value) + assert [s.p_holm for s in by_raw] == sorted(s.p_holm for s in by_raw) # non-decreasing + assert all(s.p_holm >= s.p_value for s in stats) # never anti-conservative + assert all(s.p_holm <= 1.0 for s in stats) + # Step-down stops at the first failure: no significant flag may follow a non-significant one. + flags = [s.significant for s in by_raw] + assert flags == sorted(flags, reverse=True) + + +def test_holm_charges_the_family_size_actually_tested() -> None: + """A pair count of k(k-1)/2 over VARYING factors, not over all 11.""" + sample = _sample({"condition_aligned": [1, 0] * 20, "ema_fan_aligned": [1, 0] * 20}) + stats = holm_adjust(pair_stats(sample)) + assert len(stats) == 1 + assert stats[0].p_holm == pytest.approx(stats[0].p_value) + + +# -- variance ------------------------------------------------------------------------------- + + +def test_variance_ratio_is_one_when_factors_are_independent_by_construction() -> None: + """Two factors crossed on a full 2x2 design have exactly zero covariance.""" + sample = _sample( + {"condition_aligned": [1, 1, 0, 0] * 40, "ema_fan_aligned": [1, 0, 1, 0] * 40} + ) + report = variance_report(sample) + assert report.ratio == pytest.approx(Decimal(1), abs=Decimal("0.02")) + + +def test_variance_ratio_exceeds_one_when_two_factors_are_duplicates() -> None: + """The failure mode #208 is about: one piece of evidence scored twice.""" + duplicated = [1, 0, 1, 1, 0, 0, 1, 0] * 20 + sample = _sample({"condition_aligned": duplicated, "ema_fan_aligned": list(duplicated)}) + assert variance_report(sample).ratio > Decimal("1.9") + + +def test_variance_report_is_degenerate_safe() -> None: + assert variance_report(_sample({"condition_aligned": [1]})).ratio == Decimal(0) + + +# -- clusters ------------------------------------------------------------------------------- + + +def test_cluster_report_separates_within_from_the_rest_of_the_matrix() -> None: + tight = [1, 0, 1, 1, 0, 0, 1, 0] * 20 + sample = _sample( + { + "condition_aligned": tight, + "ema_fan_aligned": list(tight), # trend cluster, perfectly coupled + "rsi_extreme": [1, 1, 0, 0] * 40, # momentum, unrelated to the above + "rsi_divergence": [1, 0, 0, 1] * 40, + "deceleration": [0, 1, 0, 1] * 40, + } + ) + reports = {r.name: r for r in cluster_report(sample, pair_stats(sample))} + assert reports["trend"].mean_within_phi > Decimal("0.99") + assert reports["trend"].mean_within_phi > reports["trend"].mean_other_phi + assert abs(reports["momentum"].mean_within_phi) < Decimal("0.2") + # Weight share is read off the shipped table, never hardcoded here. + assert reports["momentum"].weight_share == Decimal( + sum(DEFAULT_WEIGHTS[m] for m in SUSPECTED_CLUSTERS["momentum"]) + ) / Decimal(sum(DEFAULT_WEIGHTS.values())) + + +def test_cluster_means_are_signed_so_mutual_exclusion_is_not_read_as_coupling() -> None: + """Averaging |phi| would report an anti-correlated cluster as a tight one. It must not.""" + alternating = [1, 0] * 60 + sample = _sample( + { + "condition_aligned": alternating, + "ema_fan_aligned": [1 - v for v in alternating], # perfectly ANTI-correlated + } + ) + (report,) = [r for r in cluster_report(sample, pair_stats(sample)) if r.name == "trend"] + assert report.mean_within_phi < Decimal("-0.99") + + +# -- replay --------------------------------------------------------------------------------- + + +def test_replay_every_bar_produces_one_observation_per_bar_past_warmup() -> None: + candles = _candles(260) + sample = replay_every_bar("BTC-USD", candles, warmup=250) + assert sample.n == 10 + assert all(len(vector) == 10 for vector in sample.vectors.values()) + assert set(sample.vectors) == set(FACTOR_NAMES) + assert all(value in (0, 1) for value in sample.vectors["condition_aligned"]) + assert sample.labels[0].startswith("BTC-USD@") + + +def test_replay_step_subsamples_the_bar_index() -> None: + candles = _candles(260) + assert replay_every_bar("BTC-USD", candles, warmup=250, step=5).n == 2 + + +def test_replay_totals_match_the_shipped_scorer() -> None: + """`totals` must be `indicators_cts.score`'s number, not a re-derivation from the vectors.""" + candles = _candles(230) + sample = replay_every_bar("BTC-USD", candles, warmup=220) + for index in range(sample.n): + rebuilt = sum( + DEFAULT_WEIGHTS[name] * sample.vectors[name][index] for name in FACTOR_NAMES + ) + assert sample.totals[index] == rebuilt + + +def test_expanding_and_rolling_windows_differ_where_the_window_reaches() -> None: + """`regime.detect_phase` reads `candles[0]`, so window length is not cosmetic.""" + candles = _candles(600, seed=3) + expanding = replay_every_bar("BTC-USD", candles, warmup=210, window=None) + rolling = replay_every_bar("BTC-USD", candles, warmup=210, window=210) + assert expanding.n == rolling.n + assert expanding.vectors["in_pullback"] != rolling.vectors["in_pullback"] + + +class _AlwaysFires(Rule): + """A rule that proposes a 2:1 long on every bar -- isolates `replay_fired`'s gate ladder.""" + + name = "always_fires" + granularity = Granularity.ONE_DAY + + def __init__(self) -> None: + self.params = {} + self.product_id = "BTC-USD" + + def detect(self, candles_by_tf): # type: ignore[no-untyped-def] + candles = candles_by_tf[Granularity.ONE_DAY] + close = candles[-1].close + return Setup( + product_id=self.product_id, + direction="long", + entry=close, + stop=close * Decimal("0.98"), + target=close * Decimal("1.04"), + context={}, + ts=candles[-1].ts, + ) + + def exit_signal(self, held, candles_by_tf) -> bool: # type: ignore[no-untyped-def] + return False + + def describe(self) -> dict: + return {"name": self.name} + + +def test_replay_fired_is_a_strict_subset_of_every_bar() -> None: + """The gates can only remove observations, never add or invent them.""" + candles = _candles(400, seed=5) + unconditional = replay_every_bar("BTC-USD", candles, warmup=210) + fired = replay_fired(_AlwaysFires(), Granularity.ONE_DAY, candles, warmup=210) + assert 0 < fired.n <= unconditional.n + assert set(fired.vectors) == set(FACTOR_NAMES) + + +def test_replay_fired_yields_nothing_when_the_rule_never_detects() -> None: + class _NeverFires(_AlwaysFires): + def detect(self, candles_by_tf): # type: ignore[no-untyped-def] + return None + + fired = replay_fired(_NeverFires(), Granularity.ONE_DAY, _candles(260), warmup=250) + assert fired.n == 0 + assert fired.varying() == () + + +def test_factor_presence_delegates_to_the_shipped_scoring_path() -> None: + """Every key `indicators_cts` knows about comes back, and the total is its own.""" + candles = _candles(240) + setup = Setup( + product_id="BTC-USD", + direction="long", + entry=candles[-1].close, + stop=candles[-1].close * Decimal("0.98"), + target=candles[-1].close * Decimal("1.04"), + context={}, + ts=candles[-1].ts, + ) + presence, total = factor_presence(setup, candles) + assert set(presence) == set(FACTOR_NAMES) + assert presence["seasonality"] is False # weighted 0 and hardcoded absent (spec §9) + assert total == sum(DEFAULT_WEIGHTS[name] for name, present in presence.items() if present)