sharpebench-core 0.28.0

//! Regime-conditional distributional comparison: where does the edge actually live?
//!
//! The board deflates and tests significance on *pooled* returns. Pooling is a
//! real blind spot: two strategies can share a pooled distribution and still
//! behave nothing alike once you condition on the market state. One earns its
//! whole edge in a crisis and bleeds the rest of the time; the other is the
//! mirror image. Pooled, they look like twins. [`crate::decay`] gestures at the
//! same worry ("is it a one-regime fluke?") but answers it only through an IC
//! half-life, which is a time axis, not a state axis.
//!
//! So: given per-period returns for two strategies plus a regime label per
//! period, this module compares them *within* each regime and compares
//! **distributions**, not just means. A mean gap that reverses sign between
//! regimes is the headline finding, and
//! [`RegimeDistributionReport::pooled_hides_reversal`] says so outright.
//!
//! Adapted from "Regime-Conditional Distributional Comparison of Trading
//! Strategies: A GAMLSS ZAGA Framework Applied to the S&P 500".
//!
//! # What the ZAGA idea buys, and what is actually implemented here
//!
//! The load-bearing modelling claim in that paper is that a strategy return
//! series is a **mixture**: a point mass at (or near) zero for the periods with
//! no position and no trade, plus a continuous part for the periods that
//! actually traded. Smoothing that mixture into one continuous density is the
//! modelling error the ZAGA (zero-adjusted gamma) form exists to avoid: the
//! near-zero mass drags the fitted mean and variance toward zero and makes an
//! aggressive, selective strategy look like a timid one.
//!
//! **Implemented here:** the mixture split itself. Per regime, per strategy, the
//! near-zero-return mass is separated from the continuous part; both are compared.
//!
//! **It is a return mass, not a no-trade mass.** The input is returns and
//! nothing else: no position, order or trade flag reaches this module. A period
//! whose return is within `zero_tol` may be a period the strategy sat out, or a
//! period it held a position that went nowhere, or one whose gain was eaten by
//! fees. Those are different behaviours with the same observable here, so the
//! split is reported as a near-zero-return mass and reading it as inactivity is
//! the caller's inference, not this module's measurement. Supply trade or
//! position flags elsewhere if participation itself is the quantity of
//! interest.
//! The continuous part is summarised by its mean, standard deviation and median,
//! by the sign split of its non-zero returns, and by a method-of-moments gamma
//! (shape, rate) fitted to the *magnitudes* of those returns. The two continuous
//! parts are also compared non-parametrically with a two-sample
//! Kolmogorov-Smirnov statistic, so a difference in shape shows up even when the
//! means agree.
//!
//! **Deliberately NOT implemented:** a GAMLSS fit. There are no link functions,
//! no penalised splines, no smooth covariate terms for mu / sigma / nu, no
//! iterative backfitting, no likelihood maximisation, and no standard errors or
//! p-values on fitted parameters. The gamma parameters here are moment matches
//! on magnitudes, not maximum-likelihood estimates, and they are descriptive
//! summaries rather than a fitted model. Nobody should read this module as
//! evidence that a GAMLSS ZAGA fit lives in the crate. If you need the real
//! thing, fit it elsewhere and feed the result in.
//!
//! Regime labels are an **input**. This module does not infer regimes and does
//! not contain a classifier: whoever labels the periods owns that judgement, and
//! burying a classifier here would let a regime definition quietly become part
//! of the scoring kernel.

use serde::{Deserialize, Serialize};

use crate::stats::{mean, std_dev, variance};

/// Knobs for [`compare_by_regime`].
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct RegimeCompareOpts {
    /// A return whose absolute value is at or below this lands in the near-zero
    /// mass rather than the continuous part. Exact zeros are the common case,
    /// but a fee-only or rounding-dust period belongs with them rather than
    /// being modelled as a tiny continuous return. This is a threshold on the
    /// **return**, not a detection of inactivity: see the module header.
    pub zero_tol: f64,
    /// Regimes with fewer than this many periods are still reported but are
    /// excluded from the reversal verdict: a two-period "regime" flipping sign
    /// is noise, not a finding.
    pub min_periods: usize,
    /// Mean gaps at or below this magnitude are called a tie rather than an edge.
    pub tie_tol: f64,
}

impl Default for RegimeCompareOpts {
    fn default() -> Self {
        Self {
            zero_tol: 1e-9,
            min_periods: 8,
            tie_tol: 1e-12,
        }
    }
}

/// The zero-mass / continuous-part split for one strategy inside one regime.
///
/// This is the ZAGA decomposition and nothing more: see the module header for
/// what is and is not fitted.
#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
pub struct ZagaSplit {
    /// Periods observed in this regime.
    pub n: usize,
    /// Fraction of periods whose return is within `zero_tol`, in [0, 1]. This is
    /// the `nu` role in a ZAGA. It is a near-zero-**return** mass: with no trade
    /// or position flag in the input it cannot separate sitting out from holding
    /// a position that went nowhere, so it is not a participation rate.
    ///
    /// Serialized as `near_zero_return_mass`. The pre-rename key `zero_mass` is
    /// still accepted on input and is deprecated for removal in a later major.
    #[serde(alias = "zero_mass")]
    pub near_zero_return_mass: f64,
    /// Periods in the continuous part.
    pub n_nonzero: usize,
    /// Fraction of the continuous part that is positive, in [0, 1]. A gamma is
    /// non-negative, so the sign has to be carried separately; this is where it
    /// lives.
    pub positive_share: f64,
    /// Mean of the continuous part (signed).
    pub cont_mean: f64,
    /// Standard deviation of the continuous part (signed returns, sample `n-1`).
    pub cont_sd: f64,
    /// Median of the continuous part (signed).
    pub cont_median: f64,
    /// Method-of-moments gamma shape fitted to the *magnitudes* of the
    /// continuous part. 0.0 when undefined (no dispersion, or nothing to fit).
    pub gamma_shape: f64,
    /// Method-of-moments gamma rate for the same magnitudes. 0.0 when undefined.
    pub gamma_rate: f64,
    /// Mean over all periods in the regime, zeros included. Kept alongside
    /// `cont_mean` because the gap between the two is exactly the distortion the
    /// mixture split exists to expose.
    pub pooled_mean: f64,
}

/// Head-to-head comparison of two strategies inside one regime.
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
pub struct RegimeComparison {
    /// The caller's regime label.
    pub regime: String,
    /// Periods carrying this label.
    pub n_periods: usize,
    /// Strategy A's mixture split in this regime.
    pub a: ZagaSplit,
    /// Strategy B's mixture split in this regime.
    pub b: ZagaSplit,
    /// `a.near_zero_return_mass - b.near_zero_return_mass`. Positive means A had
    /// more near-zero-return periods here, which is a genuine behavioural
    /// difference even when the means agree. It is not by itself a difference in
    /// how often each traded.
    ///
    /// Serialized as `near_zero_return_mass_gap`. The pre-rename key
    /// `zero_mass_gap` is still accepted on input and is deprecated for removal in
    /// a later major.
    #[serde(alias = "zero_mass_gap")]
    pub near_zero_return_mass_gap: f64,
    /// `a.pooled_mean - b.pooled_mean` inside this regime.
    pub mean_gap: f64,
    /// `a.cont_mean - b.cont_mean`: the gap once the near-zero-return periods are
    /// taken out. Diverging from `mean_gap` means the pooled comparison was mostly
    /// measuring how often each strategy moved at all, not how well it moved.
    pub cont_mean_gap: f64,
    /// Two-sample Kolmogorov-Smirnov statistic between the two continuous parts,
    /// in [0, 1]. Non-parametric, so a pure shape difference (same mean, fatter
    /// left tail) still registers.
    pub ks_statistic: f64,
    /// `+1` if A leads on `mean_gap` here, `-1` if B leads, `0` within `tie_tol`.
    pub edge_sign: i8,
    /// Whether this regime cleared `min_periods` and therefore counts toward the
    /// reversal verdict.
    pub counted: bool,
}

/// Full regime-conditional comparison of two strategies.
#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
pub struct RegimeDistributionReport {
    /// One entry per distinct regime label, ordered lexicographically by label so
    /// the report is byte-identical on every recompute.
    pub regimes: Vec<RegimeComparison>,
    /// `mean(a) - mean(b)` over every aligned period, ignoring regime. This is
    /// the number the rest of the board would have seen.
    pub pooled_mean_gap: f64,
    /// `+1` / `-1` / `0` for `pooled_mean_gap` under `tie_tol`.
    pub pooled_edge_sign: i8,
    /// The counted regimes the pooled verdict does not represent.
    ///
    /// When `pooled_edge_sign` is non-zero these are the counted regimes whose
    /// own `edge_sign` is non-zero and opposite to it. When the pooled gap is a
    /// tie, a pooled sign cannot be contradicted, yet the sharpest reversal of
    /// all lands there: a regime where A leads and one where B leads, cancelling
    /// exactly. So a tied pooled verdict lists every counted regime with a sign
    /// whenever both signs are present, and lists none otherwise (one regime
    /// carrying the whole edge is concentration, reported by `edge_dispersion`,
    /// not a reversal).
    pub reversal_regimes: Vec<String>,
    /// True when `reversal_regimes` is non-empty: counted regimes disagree with
    /// each other or with the pooled sign. The pooled number is then not so much
    /// wrong as unusable, since it is averaging over a sign change. A tie in
    /// `pooled_mean_gap` does not suppress this.
    pub pooled_hides_reversal: bool,
    /// Spread between the best and worst counted regime `mean_gap`. A large
    /// spread with no reversal still says the edge is concentrated, not general.
    pub edge_dispersion: f64,
}

fn sign_of(x: f64, tie_tol: f64) -> i8 {
    if x > tie_tol {
        1
    } else if x < -tie_tol {
        -1
    } else {
        0
    }
}

fn sorted_copy(xs: &[f64]) -> Vec<f64> {
    let mut v = xs.to_vec();
    v.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
    v
}

fn median_of(xs: &[f64]) -> f64 {
    let n = xs.len();
    if n == 0 {
        return 0.0;
    }
    let s = sorted_copy(xs);
    if n % 2 == 1 {
        s[n / 2]
    } else {
        0.5 * (s[n / 2 - 1] + s[n / 2])
    }
}

/// Two-sample Kolmogorov-Smirnov statistic: the largest gap between the two
/// empirical CDFs. Ties are advanced together so repeated values do not inflate
/// the statistic. 0.0 when either sample is empty.
fn ks_two_sample(a: &[f64], b: &[f64]) -> f64 {
    if a.is_empty() || b.is_empty() {
        return 0.0;
    }
    let xa = sorted_copy(a);
    let xb = sorted_copy(b);
    let na = xa.len() as f64;
    let nb = xb.len() as f64;
    let (mut i, mut j) = (0usize, 0usize);
    let mut d = 0.0_f64;
    loop {
        let v = match (xa.get(i), xb.get(j)) {
            (Some(&va), Some(&vb)) => va.min(vb),
            (Some(&va), None) => va,
            (None, Some(&vb)) => vb,
            (None, None) => break,
        };
        while i < xa.len() && xa[i] <= v {
            i += 1;
        }
        while j < xb.len() && xb[j] <= v {
            j += 1;
        }
        let diff = (i as f64 / na - j as f64 / nb).abs();
        if diff > d {
            d = diff;
        }
    }
    d
}

/// Split one strategy's returns inside one regime into the near-zero-return mass
/// plus the continuous part, and summarise both.
fn zaga_split(xs: &[f64], zero_tol: f64) -> ZagaSplit {
    let n = xs.len();
    let cont: Vec<f64> = xs.iter().copied().filter(|r| r.abs() > zero_tol).collect();
    let n_nonzero = cont.len();
    let near_zero_return_mass = if n == 0 {
        0.0
    } else {
        (n - n_nonzero) as f64 / n as f64
    };
    let positive_share = if n_nonzero == 0 {
        0.0
    } else {
        cont.iter().filter(|&&r| r > 0.0).count() as f64 / n_nonzero as f64
    };

    // Gamma by moment matching on magnitudes. The gamma is the continuous half of
    // a ZAGA and is defined on the positive line, so the sign is carried by
    // `positive_share` rather than folded in here.
    let mags: Vec<f64> = cont.iter().map(|r| r.abs()).collect();
    let mm = mean(&mags);
    let mv = variance(&mags);
    let (gamma_shape, gamma_rate) = if mm > 0.0 && mv > 0.0 {
        (mm * mm / mv, mm / mv)
    } else {
        (0.0, 0.0)
    };

    ZagaSplit {
        n,
        near_zero_return_mass,
        n_nonzero,
        positive_share,
        cont_mean: mean(&cont),
        cont_sd: std_dev(&cont),
        cont_median: median_of(&cont),
        gamma_shape,
        gamma_rate,
        pooled_mean: mean(xs),
    }
}

/// Compare two strategies regime by regime.
///
/// `a` and `b` are per-period returns for the two strategies and `regimes` is
/// the regime label for the same periods; all three are aligned and the shortest
/// of them sets the comparison length, so a caller who trims one series does not
/// silently compare misaligned periods. A regime with no periods simply does not
/// appear.
///
/// Pure and deterministic: no randomness, no clock, and regimes come out in
/// lexicographic label order.
pub fn compare_by_regime(
    a: &[f64],
    b: &[f64],
    regimes: &[&str],
    opts: RegimeCompareOpts,
) -> RegimeDistributionReport {
    let n = a.len().min(b.len()).min(regimes.len());
    if n == 0 {
        return RegimeDistributionReport {
            regimes: Vec::new(),
            pooled_mean_gap: 0.0,
            pooled_edge_sign: 0,
            reversal_regimes: Vec::new(),
            pooled_hides_reversal: false,
            edge_dispersion: 0.0,
        };
    }

    let mut labels: Vec<String> = regimes[..n].iter().map(|s| (*s).to_string()).collect();
    labels.sort();
    labels.dedup();

    let pooled_mean_gap = mean(&a[..n]) - mean(&b[..n]);
    let pooled_edge_sign = sign_of(pooled_mean_gap, opts.tie_tol);

    let mut out: Vec<RegimeComparison> = Vec::with_capacity(labels.len());
    for label in labels {
        let idx: Vec<usize> = (0..n).filter(|&i| regimes[i] == label).collect();
        let ra: Vec<f64> = idx.iter().map(|&i| a[i]).collect();
        let rb: Vec<f64> = idx.iter().map(|&i| b[i]).collect();
        let sa = zaga_split(&ra, opts.zero_tol);
        let sb = zaga_split(&rb, opts.zero_tol);
        let mean_gap = sa.pooled_mean - sb.pooled_mean;
        let cont_a: Vec<f64> = ra
            .iter()
            .copied()
            .filter(|r| r.abs() > opts.zero_tol)
            .collect();
        let cont_b: Vec<f64> = rb
            .iter()
            .copied()
            .filter(|r| r.abs() > opts.zero_tol)
            .collect();
        out.push(RegimeComparison {
            regime: label,
            n_periods: idx.len(),
            near_zero_return_mass_gap: sa.near_zero_return_mass - sb.near_zero_return_mass,
            mean_gap,
            cont_mean_gap: sa.cont_mean - sb.cont_mean,
            ks_statistic: ks_two_sample(&cont_a, &cont_b),
            edge_sign: sign_of(mean_gap, opts.tie_tol),
            counted: idx.len() >= opts.min_periods,
            a: sa,
            b: sb,
        });
    }

    // A tied pooled gap is the case where a reversal hides best: equal and
    // opposite regime edges average to nothing, so there is no pooled sign left
    // to contradict. Report the disagreement on its own terms there.
    let counted_signs: Vec<i8> = out
        .iter()
        .filter(|r| r.counted)
        .map(|r| r.edge_sign)
        .collect();
    let cross_regime_reversal = counted_signs.contains(&1) && counted_signs.contains(&-1);
    let reversal_regimes: Vec<String> = out
        .iter()
        .filter(|r| {
            r.counted
                && r.edge_sign != 0
                && if pooled_edge_sign == 0 {
                    cross_regime_reversal
                } else {
                    r.edge_sign != pooled_edge_sign
                }
        })
        .map(|r| r.regime.clone())
        .collect();

    let counted_gaps: Vec<f64> = out
        .iter()
        .filter(|r| r.counted)
        .map(|r| r.mean_gap)
        .collect();
    let edge_dispersion = if counted_gaps.len() < 2 {
        0.0
    } else {
        let hi = counted_gaps
            .iter()
            .copied()
            .fold(f64::NEG_INFINITY, f64::max);
        let lo = counted_gaps.iter().copied().fold(f64::INFINITY, f64::min);
        hi - lo
    };

    RegimeDistributionReport {
        regimes: out,
        pooled_mean_gap,
        pooled_edge_sign,
        pooled_hides_reversal: !reversal_regimes.is_empty(),
        reversal_regimes,
        edge_dispersion,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    /// Alternating labels so both regimes get a deterministic, equal-sized slice.
    fn labels(n: usize) -> Vec<&'static str> {
        (0..n)
            .map(|i| if i % 2 == 0 { "bull" } else { "bear" })
            .collect()
    }

    #[test]
    fn pooling_hides_a_sign_reversal() {
        // A wins in bull by exactly as much as it loses in bear, so the pooled
        // gap is ~0 and the pooled comparison says "identical". It is not.
        let n = 120;
        let lab = labels(n);
        let a: Vec<f64> = (0..n)
            .map(|i| if i % 2 == 0 { 0.010 } else { -0.010 })
            .collect();
        let b: Vec<f64> = (0..n)
            .map(|i| if i % 2 == 0 { -0.010 } else { 0.010 })
            .collect();
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert!(
            rep.pooled_mean_gap.abs() < 1e-12,
            "pooled gap should vanish: {}",
            rep.pooled_mean_gap
        );
        assert_eq!(rep.regimes.len(), 2);
        let bull = rep.regimes.iter().find(|r| r.regime == "bull").unwrap();
        let bear = rep.regimes.iter().find(|r| r.regime == "bear").unwrap();
        assert_eq!(bull.edge_sign, 1);
        assert_eq!(bear.edge_sign, -1);
        assert!(rep.edge_dispersion > 0.03, "{}", rep.edge_dispersion);
    }

    #[test]
    fn tied_pooled_gap_still_reports_the_regime_reversal() {
        // The exact-cancellation case: A wins in bull by what it loses in bear,
        // so the pooled gap is 0 and there is no pooled sign to contradict. That
        // is the sharpest reversal there is, and it has to be reported as one.
        let n = 120;
        let lab = labels(n);
        let a: Vec<f64> = (0..n)
            .map(|i| if i % 2 == 0 { 0.010 } else { -0.010 })
            .collect();
        let b: Vec<f64> = (0..n)
            .map(|i| if i % 2 == 0 { -0.010 } else { 0.010 })
            .collect();
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert_eq!(rep.pooled_edge_sign, 0, "the pooled gap is a tie");
        assert!(
            rep.pooled_hides_reversal,
            "a tie must not suppress the reversal: {rep:?}"
        );
        assert_eq!(
            rep.reversal_regimes,
            vec!["bear".to_string(), "bull".to_string()]
        );
    }

    #[test]
    fn tied_pooled_gap_without_opposing_regimes_reports_no_reversal() {
        // The pooled gap is a tie because a four-period "shock" regime cancels a
        // forty-period "calm" one, but only "calm" clears `min_periods`. A single
        // signed counted regime is concentration, not reversal.
        let mut a = Vec::new();
        let mut b = Vec::new();
        let mut lab: Vec<&'static str> = Vec::new();
        for _ in 0..40 {
            a.push(0.010);
            b.push(0.000);
            lab.push("calm");
        }
        for _ in 0..4 {
            a.push(-0.100);
            b.push(0.000);
            lab.push("shock");
        }
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert_eq!(rep.pooled_edge_sign, 0, "the pooled gap cancels: {rep:?}");
        let shock = rep.regimes.iter().find(|r| r.regime == "shock").unwrap();
        assert!(!shock.counted, "the four-period regime must not count");
        assert!(
            !rep.pooled_hides_reversal,
            "only one counted regime has a sign: {rep:?}"
        );
        assert!(rep.reversal_regimes.is_empty());
    }

    #[test]
    fn reversal_is_flagged_against_a_signed_pooled_verdict() {
        // A is better overall, but strictly worse in "crisis".
        let mut a = Vec::new();
        let mut b = Vec::new();
        let mut lab: Vec<&'static str> = Vec::new();
        for _ in 0..60 {
            a.push(0.01);
            b.push(0.001);
            lab.push("calm");
        }
        for _ in 0..20 {
            a.push(-0.02);
            b.push(0.005);
            lab.push("crisis");
        }
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert_eq!(rep.pooled_edge_sign, 1, "A leads pooled");
        assert!(rep.pooled_hides_reversal);
        assert_eq!(rep.reversal_regimes, vec!["crisis".to_string()]);
    }

    #[test]
    fn consistent_edge_reports_no_reversal() {
        let n = 100;
        let lab = labels(n);
        let a: Vec<f64> = (0..n).map(|i| 0.004 + 0.001 * (i as f64).sin()).collect();
        let b: Vec<f64> = (0..n).map(|i| 0.001 + 0.001 * (i as f64).sin()).collect();
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert!(!rep.pooled_hides_reversal);
        assert!(rep.reversal_regimes.is_empty());
        assert!(rep.regimes.iter().all(|r| r.edge_sign == 1));
    }

    #[test]
    fn zero_mass_separates_a_selective_strategy_from_a_timid_one() {
        // A trades one period in four at 4%; B trades every period at 1%. Same
        // pooled mean, completely different behaviour: only the mixture split
        // sees it.
        let n = 80;
        let lab = vec!["all"; n];
        let a: Vec<f64> = (0..n)
            .map(|i| if i % 4 == 0 { 0.04 } else { 0.0 })
            .collect();
        let b: Vec<f64> = vec![0.01; n];
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        let r = &rep.regimes[0];
        assert!(
            r.mean_gap.abs() < 1e-12,
            "pooled means match: {}",
            r.mean_gap
        );
        assert!(
            (r.a.near_zero_return_mass - 0.75).abs() < 1e-12,
            "A sits out 3 of 4: {}",
            r.a.near_zero_return_mass
        );
        assert!(r.b.near_zero_return_mass.abs() < 1e-12, "B always trades");
        assert!(
            r.near_zero_return_mass_gap > 0.7,
            "the split should shout: {}",
            r.near_zero_return_mass_gap
        );
        assert!(
            (r.cont_mean_gap - 0.03).abs() < 1e-12,
            "per-trade gap is 4% vs 1%: {}",
            r.cont_mean_gap
        );
        assert!(
            r.ks_statistic > 0.9,
            "disjoint supports: {}",
            r.ks_statistic
        );
    }

    #[test]
    fn ks_catches_a_shape_difference_at_equal_means() {
        // Same mean, same zero mass, different dispersion.
        let n = 100;
        let lab = vec!["all"; n];
        let a: Vec<f64> = (0..n)
            .map(|i| if i % 2 == 0 { 0.001 } else { -0.001 })
            .collect();
        let b: Vec<f64> = (0..n)
            .map(|i| if i % 2 == 0 { 0.05 } else { -0.05 })
            .collect();
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        let r = &rep.regimes[0];
        assert!(r.mean_gap.abs() < 1e-12);
        assert!(
            r.ks_statistic > 0.4,
            "shape difference should register: {}",
            r.ks_statistic
        );
        assert!(r.b.cont_sd > r.a.cont_sd);
    }

    #[test]
    fn gamma_moments_match_a_known_magnitude_sample() {
        // Magnitudes {1, 2, 3, 4}: mean 2.5, sample variance 5/3.
        // shape = 2.5^2 / (5/3) = 3.75, rate = 2.5 / (5/3) = 1.5.
        let a = vec![1.0, -2.0, 3.0, -4.0];
        let b = vec![1.0, -2.0, 3.0, -4.0];
        let lab = vec!["all"; 4];
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        let s = rep.regimes[0].a;
        assert!((s.gamma_shape - 3.75).abs() < 1e-9, "{}", s.gamma_shape);
        assert!((s.gamma_rate - 1.5).abs() < 1e-9, "{}", s.gamma_rate);
        assert!((s.positive_share - 0.5).abs() < 1e-12);
    }

    #[test]
    fn short_regimes_do_not_drive_the_reversal_verdict() {
        // A three-period "crisis" reversal is noise and must not be counted.
        let mut a = vec![0.01; 60];
        let mut b = vec![0.001; 60];
        let mut lab: Vec<&'static str> = vec!["calm"; 60];
        a.extend([-0.05, -0.05, -0.05]);
        b.extend([0.02, 0.02, 0.02]);
        lab.extend(["crisis", "crisis", "crisis"]);
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        let crisis = rep.regimes.iter().find(|r| r.regime == "crisis").unwrap();
        assert_eq!(crisis.edge_sign, -1, "the reversal is still reported");
        assert!(!crisis.counted, "but it is not counted");
        assert!(!rep.pooled_hides_reversal);
    }

    #[test]
    fn misaligned_inputs_truncate_to_the_shortest() {
        let a = vec![0.01; 10];
        let b = vec![0.0; 4];
        let lab = vec!["all"; 7];
        let rep = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert_eq!(rep.regimes[0].n_periods, 4);
    }

    #[test]
    fn empty_input_is_inert() {
        let rep = compare_by_regime(&[], &[], &[], RegimeCompareOpts::default());
        assert!(rep.regimes.is_empty());
        assert_eq!(rep.pooled_edge_sign, 0);
        assert!(!rep.pooled_hides_reversal);
    }

    #[test]
    fn output_order_is_lexicographic_and_reproducible() {
        let lab = vec!["zulu", "alpha", "mike", "alpha", "zulu", "mike"];
        let a = vec![0.01, 0.02, 0.03, 0.04, 0.05, 0.06];
        let b = vec![0.0; 6];
        let first = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        let second = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());
        assert_eq!(first, second, "recompute must be identical");
        let order: Vec<&str> = first.regimes.iter().map(|r| r.regime.as_str()).collect();
        assert_eq!(order, vec!["alpha", "mike", "zulu"]);
    }

    #[test]
    fn near_zero_return_mass_serializes_under_the_new_key_and_reads_the_old_one() {
        let lab = vec!["all"; 8];
        let a = vec![0.04, 0.0, 0.0, 0.0, 0.04, 0.0, 0.0, 0.0];
        let b = vec![0.01; 8];
        let report = compare_by_regime(&a, &b, &lab, RegimeCompareOpts::default());

        let json = serde_json::to_string(&report).unwrap();
        assert!(json.contains("\"near_zero_return_mass\":"), "{json}");
        assert!(json.contains("\"near_zero_return_mass_gap\":"), "{json}");
        assert!(
            !json.contains("\"zero_mass\":"),
            "old key must not serialize"
        );
        assert!(
            !json.contains("\"zero_mass_gap\":"),
            "old key must not serialize"
        );

        let round_trip: RegimeDistributionReport = serde_json::from_str(&json).unwrap();
        assert_eq!(round_trip, report);

        // A report saved before the rename still parses: `zero_mass` and
        // `zero_mass_gap` are accepted as aliases of the new keys.
        let legacy = json
            .replace("\"near_zero_return_mass_gap\":", "\"zero_mass_gap\":")
            .replace("\"near_zero_return_mass\":", "\"zero_mass\":");
        assert!(legacy.contains("\"zero_mass\":") && !legacy.contains("near_zero"));
        let from_legacy: RegimeDistributionReport = serde_json::from_str(&legacy).unwrap();
        assert_eq!(from_legacy, report);
    }
}