foxhunt/crates/ml-validation/src/statistical.rs

//! Statistical validation tests for backtesting robustness.
//!
//! Provides:
//! - Math helpers: normal CDF/PPF, Sharpe ratio, skewness, excess kurtosis
//! - Deflated Sharpe Ratio (Bailey & Lopez de Prado, 2014)
//! - Monte Carlo permutation test for strategy significance
//! - Probability of Backtest Overfitting (Bailey et al., 2017) via CSCV

use rand::seq::SliceRandom;
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use serde::{Deserialize, Serialize};

// ─── Section 1: Math helpers ────────────────────────────────────────────────

/// Standard normal CDF using Abramowitz & Stegun rational approximation.
///
/// Approximation error is less than 7.5e-8 across the entire domain.
/// Values outside [-8, 8] are clamped to 0.0 or 1.0 respectively.
pub fn normal_cdf(x: f64) -> f64 {
    if x <= -8.0 {
        return 0.0;
    }
    if x >= 8.0 {
        return 1.0;
    }

    // Abramowitz & Stegun approximation 26.2.17
    const B1: f64 = 0.319_381_53;
    const B2: f64 = -0.356_563_782;
    const B3: f64 = 1.781_477_937;
    const B4: f64 = -1.821_255_978;
    const B5: f64 = 1.330_274_429;

    let abs_x = x.abs();
    let t = 1.0 / (1.0 + 0.231_641_9 * abs_x);
    let t2 = t * t;
    let t3 = t2 * t;
    let t4 = t3 * t;
    let t5 = t4 * t;

    let poly = B1 * t + B2 * t2 + B3 * t3 + B4 * t4 + B5 * t5;

    // Standard normal PDF = (1/sqrt(2*pi)) * exp(-x^2/2)
    let pdf = (-0.5 * abs_x * abs_x).exp() / (2.0 * std::f64::consts::PI).sqrt();

    let cdf_right = 1.0 - pdf * poly;

    if x >= 0.0 {
        cdf_right
    } else {
        1.0 - cdf_right
    }
}

/// Inverse standard normal CDF (quantile function) using Peter Acklam's algorithm.
///
/// Handles edge cases: p <= 0 returns -8.0, p >= 1 returns 8.0, p = 0.5 returns 0.0.
pub fn normal_ppf(p: f64) -> f64 {
    if p <= 0.0 {
        return -8.0;
    }
    if p >= 1.0 {
        return 8.0;
    }
    if (p - 0.5).abs() < 1e-15 {
        return 0.0;
    }

    // Rational approximation coefficients (Peter Acklam)
    const A1: f64 = -3.969_683_028_665_376e1;
    const A2: f64 = 2.209_460_984_245_205e2;
    const A3: f64 = -2.759_285_104_469_687e2;
    const A4: f64 = 1.383_577_518_672_69e2;
    const A5: f64 = -3.066_479_806_614_716e1;
    const A6: f64 = 2.506_628_277_459_239;

    const B1: f64 = -5.447_609_879_822_406e1;
    const B2: f64 = 1.615_858_368_580_409e2;
    const B3: f64 = -1.556_989_798_598_866e2;
    const B4: f64 = 6.680_131_188_771_972e1;
    const B5: f64 = -1.328_068_155_288_572e1;

    const C1: f64 = -7.784_894_002_430_293e-3;
    const C2: f64 = -3.223_964_580_411_365e-1;
    const C3: f64 = -2.400_758_277_161_838;
    const C4: f64 = -2.549_732_539_343_734;
    const C5: f64 = 4.374_664_141_464_968;
    const C6: f64 = 2.938_163_982_698_783;

    const D1: f64 = 7.784_695_709_041_462e-3;
    const D2: f64 = 3.224_671_290_700_398e-1;
    const D3: f64 = 2.445_134_137_142_996;
    const D4: f64 = 3.754_408_661_907_416;

    const P_LOW: f64 = 0.024_25;
    const P_HIGH: f64 = 1.0 - P_LOW;

    if p < P_LOW {
        // Rational approximation for lower region
        let q = (-2.0 * p.ln()).sqrt();
        let num = ((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6;
        let den = (((D1 * q + D2) * q + D3) * q + D4) * q + 1.0;
        num / den
    } else if p <= P_HIGH {
        // Rational approximation for central region
        let q = p - 0.5;
        let r = q * q;
        let num = (((((A1 * r + A2) * r + A3) * r + A4) * r + A5) * r + A6) * q;
        let den = ((((B1 * r + B2) * r + B3) * r + B4) * r + B5) * r + 1.0;
        num / den
    } else {
        // Rational approximation for upper region
        let q = (-2.0 * (1.0 - p).ln()).sqrt();
        let num = ((((C1 * q + C2) * q + C3) * q + C4) * q + C5) * q + C6;
        let den = (((D1 * q + D2) * q + D3) * q + D4) * q + 1.0;
        -(num / den)
    }
}

/// Annualized Sharpe ratio (mean / sample-std, no annualization factor).
///
/// Returns 0.0 if fewer than 2 returns or if standard deviation is approximately zero.
pub fn sharpe_ratio(returns: &[f64]) -> f64 {
    if returns.len() < 2 {
        return 0.0;
    }

    let n = returns.len() as f64;
    let mean = returns.iter().sum::<f64>() / n;

    let variance = returns.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / (n - 1.0);
    let std = variance.sqrt();

    if std < 1e-15 {
        return 0.0;
    }

    mean / std
}

/// Adjusted Fisher-Pearson skewness coefficient.
///
/// Returns 0.0 if fewer than 3 observations.
pub fn skewness(returns: &[f64]) -> f64 {
    if returns.len() < 3 {
        return 0.0;
    }

    let n = returns.len() as f64;
    let mean = returns.iter().sum::<f64>() / n;

    let m2 = returns.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / n;
    let m3 = returns.iter().map(|r| (r - mean).powi(3)).sum::<f64>() / n;

    if m2 < 1e-30 {
        return 0.0;
    }

    let skew_raw = m3 / m2.powf(1.5);

    // Adjusted Fisher-Pearson: multiply by sqrt(n*(n-1)) / (n-2)
    let adjustment = (n * (n - 1.0)).sqrt() / (n - 2.0);
    skew_raw * adjustment
}

/// Excess kurtosis (subtracts 3 from the raw kurtosis).
///
/// Returns 0.0 if fewer than 4 observations.
pub fn excess_kurtosis(returns: &[f64]) -> f64 {
    if returns.len() < 4 {
        return 0.0;
    }

    let n = returns.len() as f64;
    let mean = returns.iter().sum::<f64>() / n;

    let m2 = returns.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / n;
    let m4 = returns.iter().map(|r| (r - mean).powi(4)).sum::<f64>() / n;

    if m2 < 1e-30 {
        return 0.0;
    }

    let kurt_raw = m4 / (m2 * m2);

    // Excess kurtosis: subtract 3
    // Apply sample correction: ((n-1)/((n-2)*(n-3))) * ((n+1)*kurt_raw - 3*(n-1))
    // This is the standard bias-corrected excess kurtosis formula
    let numerator = (n + 1.0) * (kurt_raw - 3.0) + 6.0;
    let denominator = (n - 2.0) * (n - 3.0);

    if denominator.abs() < 1e-30 {
        return 0.0;
    }

    (n - 1.0) * numerator / denominator
}

// ─── Section 2: Deflated Sharpe Ratio ───────────────────────────────────────

/// Result of the Deflated Sharpe Ratio test.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct DsrResult {
    /// The observed Sharpe ratio being tested.
    pub observed_sharpe: f64,
    /// Expected maximum Sharpe ratio under the null hypothesis.
    pub expected_max_sharpe: f64,
    /// Standard error of the Sharpe ratio estimate.
    pub sharpe_std_error: f64,
    /// The DSR test statistic.
    pub deflated_sharpe: f64,
    /// One-sided p-value (probability of observing this DSR under null).
    pub pvalue: f64,
}

/// Deflated Sharpe Ratio (Bailey & Lopez de Prado, 2014).
///
/// Tests whether an observed Sharpe ratio is statistically significant after
/// accounting for multiple testing (number of strategy trials).
///
/// # Arguments
/// * `observed_sharpe` - The Sharpe ratio to test
/// * `num_trials` - Number of independent strategy trials conducted
/// * `sharpe_variance` - Variance of Sharpe ratios across trials
/// * `skew` - Skewness of the return distribution
/// * `kurt` - Excess kurtosis of the return distribution
/// * `num_observations` - Number of return observations
pub fn deflated_sharpe_ratio(
    observed_sharpe: f64,
    num_trials: usize,
    sharpe_variance: f64,
    skew: f64,
    kurt: f64,
    num_observations: usize,
) -> DsrResult {
    // Euler-Mascheroni constant
    const GAMMA: f64 = 0.577_215_664_901_532_9;

    let n = num_trials.max(1) as f64;
    let obs = num_observations.max(2) as f64;

    // Expected maximum Sharpe under null:
    // SR* = sqrt(V[SR]) * ((1 - gamma) * ppf(1 - 1/N) + gamma * ppf(1 - 1/(N*e)))
    let sr_std = sharpe_variance.abs().sqrt();

    let term1 = if n > 1.0 {
        (1.0 - GAMMA) * normal_ppf(1.0 - 1.0 / n)
    } else {
        0.0
    };

    let n_e = n * std::f64::consts::E;
    let term2 = if n_e > 1.0 {
        GAMMA * normal_ppf(1.0 - 1.0 / n_e)
    } else {
        0.0
    };

    let expected_max_sharpe = sr_std * (term1 + term2);

    // Sharpe ratio standard error (Bailey & Lopez de Prado, 2014):
    // SE = sqrt((1 - γ₃·SR + (γ₄/4)·SR²) / (T-1))
    // where γ₃ = skewness, γ₄ = excess kurtosis
    let sr_sq = observed_sharpe * observed_sharpe;
    let se_inner = 1.0 - skew * observed_sharpe + kurt / 4.0 * sr_sq;
    let se_inner_clamped = se_inner.max(1e-10); // prevent negative sqrt
    let sharpe_std_error = (se_inner_clamped / (obs - 1.0)).sqrt();

    // DSR statistic
    let deflated_sharpe = if sharpe_std_error > 1e-15 {
        (observed_sharpe - expected_max_sharpe) / sharpe_std_error
    } else {
        0.0
    };

    // p-value: 1 - Phi(DSR)
    let pvalue = 1.0 - normal_cdf(deflated_sharpe);

    DsrResult {
        observed_sharpe,
        expected_max_sharpe,
        sharpe_std_error,
        deflated_sharpe,
        pvalue,
    }
}

// ─── Section 3: Monte Carlo Permutation Test ────────────────────────────────

/// Result of a Monte Carlo permutation test.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct PermutationResult {
    /// The observed Sharpe ratio of the original return series.
    pub observed_sharpe: f64,
    /// Proportion of permuted Sharpe ratios >= observed (empirical p-value).
    pub pvalue: f64,
    /// Number of permutations performed.
    pub num_permutations: usize,
    /// Mean of the null distribution of Sharpe ratios.
    pub null_mean: f64,
    /// Standard deviation of the null distribution.
    pub null_std: f64,
}

/// Monte Carlo permutation test for strategy significance.
///
/// Shuffles daily returns `num_permutations` times using a seeded ChaCha8Rng,
/// recomputes the Sharpe ratio each time, and returns the fraction of permuted
/// Sharpe ratios that are greater than or equal to the observed Sharpe.
///
/// # Arguments
/// * `daily_returns` - Slice of daily return values
/// * `num_permutations` - Number of random shuffles to perform
/// * `seed` - RNG seed for reproducibility
pub fn permutation_test(
    daily_returns: &[f64],
    num_permutations: usize,
    seed: u64,
) -> PermutationResult {
    let observed = sharpe_ratio(daily_returns);

    if daily_returns.len() < 2 || num_permutations == 0 {
        return PermutationResult {
            observed_sharpe: observed,
            pvalue: 1.0,
            num_permutations: 0,
            null_mean: 0.0,
            null_std: 0.0,
        };
    }

    let mut rng = ChaCha8Rng::seed_from_u64(seed);
    let mut shuffled: Vec<f64> = daily_returns.to_vec();
    let mut count_ge = 0_usize;
    let mut null_sharpes: Vec<f64> = Vec::with_capacity(num_permutations);

    for _ in 0..num_permutations {
        shuffled.copy_from_slice(daily_returns);
        shuffled.shuffle(&mut rng);
        let perm_sharpe = sharpe_ratio(&shuffled);
        null_sharpes.push(perm_sharpe);
        if perm_sharpe >= observed {
            count_ge = count_ge.saturating_add(1);
        }
    }

    let pvalue = count_ge as f64 / num_permutations as f64;

    let null_n = null_sharpes.len() as f64;
    let null_mean = if null_n > 0.0 {
        null_sharpes.iter().sum::<f64>() / null_n
    } else {
        0.0
    };

    let null_std = if null_n > 1.0 {
        let variance =
            null_sharpes.iter().map(|s| (s - null_mean).powi(2)).sum::<f64>() / (null_n - 1.0);
        variance.sqrt()
    } else {
        0.0
    };

    PermutationResult {
        observed_sharpe: observed,
        pvalue,
        num_permutations,
        null_mean,
        null_std,
    }
}

// ─── Section 4: Probability of Backtest Overfitting ─────────────────────────

/// Result of the Probability of Backtest Overfitting (PBO) test.
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct PboResult {
    /// Estimated probability of backtest overfitting [0, 1].
    pub pbo: f64,
    /// Number of IS/OOS combinations evaluated.
    pub num_combinations: usize,
    /// Distribution of logit values across combinations.
    pub logit_distribution: Vec<f64>,
}

/// Compute binomial coefficient C(n, k) using iterative multiplication.
///
/// Returns 0 if k > n.
pub fn binomial_coefficient(n: usize, k: usize) -> usize {
    if k > n {
        return 0;
    }
    if k == 0 || k == n {
        return 1;
    }

    // Use the smaller of k and n-k for efficiency
    let k = k.min(n.saturating_sub(k));
    let mut result: usize = 1;
    for i in 0..k {
        // result = result * (n - i) / (i + 1)
        // Do multiplication first, then division, to stay in integer domain
        result = result.saturating_mul(n.saturating_sub(i)) / (i.saturating_add(1));
    }
    result
}

/// Generate combinations of `k` items chosen from `n`.
///
/// If `C(n, k) <= 10000`, generates all combinations exhaustively.
/// Otherwise, randomly samples up to 10000 combinations using a seeded RNG.
fn generate_combinations(n: usize, k: usize) -> Vec<Vec<usize>> {
    let total = binomial_coefficient(n, k);

    if total == 0 || k == 0 || k > n {
        return Vec::new();
    }

    const MAX_COMBOS: usize = 10_000;

    if total <= MAX_COMBOS {
        // Exhaustive generation
        let mut combos = Vec::with_capacity(total);
        let mut current = Vec::with_capacity(k);
        generate_combos_recursive(n, k, 0, &mut current, &mut combos);
        combos
    } else {
        // Random sampling
        let mut rng = ChaCha8Rng::seed_from_u64(42);
        let mut combos = Vec::with_capacity(MAX_COMBOS);
        let mut seen = std::collections::HashSet::new();

        for _ in 0..MAX_COMBOS.saturating_mul(3) {
            if combos.len() >= MAX_COMBOS {
                break;
            }
            // Generate a random combination
            let mut indices: Vec<usize> = (0..n).collect();
            indices.shuffle(&mut rng);
            let mut combo: Vec<usize> = indices.into_iter().take(k).collect();
            combo.sort_unstable();

            if seen.insert(combo.clone()) {
                combos.push(combo);
            }
        }

        combos
    }
}

/// Recursive helper for exhaustive combination generation.
fn generate_combos_recursive(
    n: usize,
    k: usize,
    start: usize,
    current: &mut Vec<usize>,
    result: &mut Vec<Vec<usize>>,
) {
    if current.len() == k {
        result.push(current.clone());
        return;
    }

    let remaining = k.saturating_sub(current.len());
    for i in start..=n.saturating_sub(remaining) {
        current.push(i);
        generate_combos_recursive(n, k, i.saturating_add(1), current, result);
        current.pop();
    }
}

/// Probability of Backtest Overfitting (adapted from Bailey et al., 2017) via CSCV.
///
/// Single-strategy adaptation of Combinatorially Symmetric Cross-Validation:
/// 1. Takes N fold Sharpe ratios and generates C(N, N/2) combinations
/// 2. For each combination, one half is treated as in-sample (IS), the other as OOS
/// 3. Computes mean Sharpe for IS and OOS halves
/// 4. PBO = fraction of combinations where IS mean > OOS mean
///
/// **Interpretation for single strategies:**
/// - PBO ≈ 0.5: performance is consistent across time periods (no overfitting signal)
/// - PBO >> 0.5: performance degrades from IS to OOS (temporal overfitting)
/// - PBO << 0.5: performance improves from IS to OOS (unusual, possible regime shift)
///
/// **Note:** The original CSCV paper tests multiple strategy configurations against
/// each other. This adaptation tests a single strategy's consistency across time.
/// For single strategies, DSR and permutation tests are the primary significance measures;
/// PBO adds value when `num_trials > 1` in a hyperopt setting.
///
/// Requires even N >= 4. Returns `pbo = 0.5` for degenerate inputs.
pub fn probability_of_backtest_overfitting(per_fold_sharpes: &[f64]) -> PboResult {
    let n = per_fold_sharpes.len();

    // Need even N >= 4
    if n < 4 || n % 2 != 0 {
        return PboResult {
            pbo: 0.5,
            num_combinations: 0,
            logit_distribution: Vec::new(),
        };
    }

    let half = n / 2;
    let combinations = generate_combinations(n, half);

    if combinations.is_empty() {
        return PboResult {
            pbo: 0.5,
            num_combinations: 0,
            logit_distribution: Vec::new(),
        };
    }

    let num_combinations = combinations.len();
    let mut logit_distribution = Vec::with_capacity(num_combinations);
    let mut overfit_count = 0_usize;

    for combo in &combinations {
        let is_set: std::collections::HashSet<usize> = combo.iter().copied().collect();

        // Compute mean Sharpe for IS and OOS halves
        let is_sum: f64 = combo
            .iter()
            .filter_map(|&i| per_fold_sharpes.get(i).copied())
            .sum();
        let oos_sum: f64 = (0..n)
            .filter(|i| !is_set.contains(i))
            .filter_map(|i| per_fold_sharpes.get(i).copied())
            .sum();

        let is_count = combo.len() as f64;
        let oos_count = (n - combo.len()) as f64;

        if is_count < 1.0 || oos_count < 1.0 {
            continue;
        }

        let is_mean = is_sum / is_count;
        let oos_mean = oos_sum / oos_count;

        // Performance degradation: IS outperforming OOS signals overfitting
        let degradation = is_mean - oos_mean;

        // Logit of relative OOS performance.
        // ratio = OOS / (|IS| + |OOS|): 0.5 when equal, <0.5 when IS > OOS
        let denom = is_mean.abs() + oos_mean.abs();
        let ratio = if denom > 1e-15 {
            (oos_mean.abs() / denom).clamp(0.01, 0.99)
        } else {
            0.5 // both means ≈ 0
        };
        // Adjust sign: if OOS is negative but IS positive, ratio should reflect underperformance
        let signed_ratio = if degradation > 0.0 {
            // IS > OOS: overfitting signal, push ratio below 0.5
            ratio.min(0.49)
        } else {
            // OOS >= IS: genuine signal, push ratio above 0.5
            ratio.max(0.51)
        };
        let logit = (signed_ratio / (1.0 - signed_ratio)).ln();
        logit_distribution.push(logit);

        // Count as overfitting if IS mean > OOS mean
        if is_mean > oos_mean {
            overfit_count = overfit_count.saturating_add(1);
        }
    }

    let pbo = if !logit_distribution.is_empty() {
        overfit_count as f64 / logit_distribution.len() as f64
    } else {
        0.5
    };

    PboResult {
        pbo,
        num_combinations,
        logit_distribution,
    }
}

// ─── Tests ──────────────────────────────────────────────────────────────────

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

    // ── Math helper tests ───────────────────────────────────────────────

    #[test]
    fn test_normal_cdf_known_values() {
        // Phi(0) = 0.5
        assert!((normal_cdf(0.0) - 0.5).abs() < 1e-6, "Phi(0) should be 0.5");

        // Phi(1.96) ~ 0.975
        assert!(
            (normal_cdf(1.96) - 0.975).abs() < 0.002,
            "Phi(1.96) should be ~0.975, got {}",
            normal_cdf(1.96)
        );

        // Phi(-1.96) ~ 0.025
        assert!(
            (normal_cdf(-1.96) - 0.025).abs() < 0.002,
            "Phi(-1.96) should be ~0.025, got {}",
            normal_cdf(-1.96)
        );
    }

    #[test]
    fn test_normal_ppf_known_values() {
        // ppf(0.5) = 0
        assert!(
            normal_ppf(0.5).abs() < 1e-6,
            "ppf(0.5) should be 0, got {}",
            normal_ppf(0.5)
        );

        // ppf(0.975) ~ 1.96
        assert!(
            (normal_ppf(0.975) - 1.96).abs() < 0.01,
            "ppf(0.975) should be ~1.96, got {}",
            normal_ppf(0.975)
        );
    }

    #[test]
    fn test_normal_cdf_ppf_roundtrip() {
        for i in 1..=99 {
            let p = i as f64 / 100.0;
            let roundtrip = normal_cdf(normal_ppf(p));
            assert!(
                (roundtrip - p).abs() < 0.005,
                "Roundtrip failed for p={}: got {}",
                p,
                roundtrip
            );
        }
    }

    #[test]
    fn test_sharpe_ratio_positive_trend() {
        // All-positive returns should give SR > 0
        let returns: Vec<f64> = (1..=100).map(|i| 0.001 * i as f64).collect();
        let sr = sharpe_ratio(&returns);
        assert!(
            sr > 0.0,
            "Sharpe of all-positive returns should be > 0, got {}",
            sr
        );
    }

    #[test]
    fn test_sharpe_ratio_zero_returns() {
        let returns = vec![0.0; 50];
        let sr = sharpe_ratio(&returns);
        assert!(
            sr.abs() < 1e-10,
            "Sharpe of all-zero returns should be 0, got {}",
            sr
        );
    }

    #[test]
    fn test_skewness_symmetric() {
        // Symmetric returns: [-3, -2, -1, 0, 1, 2, 3] -> skew ~ 0
        let returns: Vec<f64> = (-30..=30).map(|i| i as f64 * 0.01).collect();
        let skew = skewness(&returns);
        assert!(
            skew.abs() < 0.5,
            "Skewness of symmetric returns should be near 0, got {}",
            skew
        );
    }

    #[test]
    fn test_excess_kurtosis_normal() {
        // For a uniform-ish distribution, kurtosis should be bounded
        let returns: Vec<f64> = (0..1000).map(|i| ((i * 7 + 13) % 100) as f64 / 100.0).collect();
        let kurt = excess_kurtosis(&returns);
        // Uniform distribution has excess kurtosis of -1.2
        // Our pseudo-uniform should have bounded kurtosis
        assert!(
            kurt.abs() < 5.0,
            "Excess kurtosis should be bounded, got {}",
            kurt
        );
    }

    // ── DSR tests ───────────────────────────────────────────────────────

    #[test]
    fn test_dsr_single_trial_not_penalized() {
        // SR=2.0, only 1 trial, normal returns (excess kurtosis=0)
        // -> no multiple testing penalty -> should be significant
        let result = deflated_sharpe_ratio(2.0, 1, 1.0, 0.0, 0.0, 252);
        assert!(
            result.pvalue < 0.05,
            "Single trial SR=2.0 should have p<0.05, got {}",
            result.pvalue
        );
    }

    #[test]
    fn test_dsr_many_trials_penalized() {
        // SR=1.0, 1000 trials, normal returns (excess kurtosis=0)
        // -> heavy penalty -> should NOT be significant
        let result = deflated_sharpe_ratio(1.0, 1000, 1.0, 0.0, 0.0, 252);
        assert!(
            result.pvalue > 0.05,
            "1000 trials SR=1.0 should have p>0.05, got {}",
            result.pvalue
        );
    }

    #[test]
    fn test_dsr_higher_sharpe_more_significant() {
        // SR=3 should be more significant than SR=1 with same parameters
        let r1 = deflated_sharpe_ratio(1.0, 100, 1.0, 0.0, 0.0, 252);
        let r3 = deflated_sharpe_ratio(3.0, 100, 1.0, 0.0, 0.0, 252);
        assert!(
            r3.pvalue < r1.pvalue,
            "SR=3 p-value ({}) should be < SR=1 p-value ({})",
            r3.pvalue,
            r1.pvalue
        );
    }

    #[test]
    fn test_dsr_se_formula_normal_returns() {
        // For normal returns (skew=0, excess_kurt=0), SE should simplify to:
        // SE = sqrt(1 / (T-1)) = 1/sqrt(T-1)
        let result = deflated_sharpe_ratio(1.0, 1, 1.0, 0.0, 0.0, 252);
        let expected_se = 1.0 / (251.0_f64).sqrt();
        assert!(
            (result.sharpe_std_error - expected_se).abs() < 1e-6,
            "SE for normal returns should be 1/sqrt(T-1)={}, got {}",
            expected_se,
            result.sharpe_std_error
        );
    }

    // ── Permutation test ────────────────────────────────────────────────

    #[test]
    fn test_permutation_test_random_returns_high_pvalue() {
        // Alternating +/- returns have no real signal
        let returns: Vec<f64> = (0..200)
            .map(|i| if i % 2 == 0 { 0.01 } else { -0.01 })
            .collect();
        let result = permutation_test(&returns, 500, 12345);
        assert!(
            result.pvalue > 0.05,
            "Alternating returns should have p>0.05, got {}",
            result.pvalue
        );
    }

    #[test]
    fn test_permutation_test_reproducible() {
        let returns: Vec<f64> = (0..100).map(|i| (i as f64).sin() * 0.01).collect();
        let r1 = permutation_test(&returns, 200, 42);
        let r2 = permutation_test(&returns, 200, 42);
        assert!(
            (r1.pvalue - r2.pvalue).abs() < 1e-10,
            "Same seed should produce identical p-values: {} vs {}",
            r1.pvalue,
            r2.pvalue
        );
        assert!(
            (r1.observed_sharpe - r2.observed_sharpe).abs() < 1e-10,
            "Same seed should produce identical observed Sharpe: {} vs {}",
            r1.observed_sharpe,
            r2.observed_sharpe
        );
    }

    // ── PBO tests ───────────────────────────────────────────────────────

    #[test]
    fn test_pbo_produces_valid_output() {
        // 8 folds with varied Sharpes -- verify structural correctness
        let sharpes = vec![1.0, 1.1, 1.05, 0.95, 1.02, 0.98, 1.08, 1.03];
        let result = probability_of_backtest_overfitting(&sharpes);

        assert!(
            (0.0..=1.0).contains(&result.pbo),
            "PBO must be in [0,1], got {}",
            result.pbo
        );
        // C(8,4) = 70 combinations
        assert_eq!(result.num_combinations, 70);
        assert!(
            !result.logit_distribution.is_empty(),
            "Should produce logit distribution"
        );
        // All logits should be finite
        assert!(
            result.logit_distribution.iter().all(|l| l.is_finite()),
            "All logits should be finite"
        );
    }

    #[test]
    fn test_pbo_too_few_folds() {
        // N=2 is too few -> degenerate result
        let sharpes = vec![1.0, 0.5];
        let result = probability_of_backtest_overfitting(&sharpes);
        assert!(
            (result.pbo - 0.5).abs() < 1e-10,
            "N=2 should return pbo=0.5, got {}",
            result.pbo
        );
        assert_eq!(result.num_combinations, 0);
    }

    #[test]
    fn test_pbo_odd_folds_degenerate() {
        // N=5 (odd) -> degenerate result
        let sharpes = vec![1.0, 2.0, 3.0, 4.0, 5.0];
        let result = probability_of_backtest_overfitting(&sharpes);
        assert!(
            (result.pbo - 0.5).abs() < 1e-10,
            "Odd N should return pbo=0.5, got {}",
            result.pbo
        );
        assert_eq!(result.num_combinations, 0);
    }

    #[test]
    fn test_pbo_identical_sharpes_symmetric() {
        // All identical Sharpes → IS_mean = OOS_mean for every combination
        // Neither IS > OOS nor OOS > IS → PBO = 0.0 (no overfit count)
        let sharpes = vec![1.0, 1.0, 1.0, 1.0];
        let result = probability_of_backtest_overfitting(&sharpes);
        assert!(
            result.pbo < 0.01,
            "Identical Sharpes should have PBO ≈ 0.0 (no IS > OOS), got {}",
            result.pbo
        );
    }

    #[test]
    fn test_binomial_coefficient() {
        assert_eq!(binomial_coefficient(8, 4), 70);
        assert_eq!(binomial_coefficient(6, 3), 20);
        assert_eq!(binomial_coefficient(0, 0), 1);
        assert_eq!(binomial_coefficient(5, 0), 1);
        assert_eq!(binomial_coefficient(5, 5), 1);
        assert_eq!(binomial_coefficient(3, 5), 0);
    }
}