/// The profile of the objective in one parameter.
///
/// For each value, the best fit achievable with that parameter pinned there and
/// the other two re-optimised. This is the honest picture of what a quote set
/// determines: a profile with a sharp minimum means the parameter is pinned, and
/// a flat one means it is not, however confidently the optimiser reported it.
pub fn profile(surface: &Surface, fitted: &Sabr, p: Parameter, values: &[f64]) -> Vec<f64> {
    let free: Vec<Parameter> = [Parameter::Alpha, Parameter::Rho, Parameter::Nu]
        .into_iter()
        .filter(|&q| q != p)
        .collect();

    values
        .iter()
        .map(|&v| {
            let start = p.set(fitted, v);
            let refitted = fit(surface, &start, &free);
            surface.rms_error(&refitted)
        })
        .collect()
}

/// How far a parameter can move, in each direction, while the refitted surface
/// stays within `tolerance` volatility points of the quotes.
///
/// The number a calibration report should carry beside every fitted parameter
/// and never does.
pub fn identifiable_range(
    surface: &Surface,
    fitted: &Sabr,
    p: Parameter,
    tolerance: f64,
) -> (f64, f64) {
    let centre = p.get(fitted);
    let (lo_bound, hi_bound) = bounds_of(p);
    let base = surface.rms_error(fitted);

    let acceptable = |v: f64| {
        let free: Vec<Parameter> = [Parameter::Alpha, Parameter::Rho, Parameter::Nu]
            .into_iter()
            .filter(|&q| q != p)
            .collect();
        let refitted = fit(surface, &p.set(fitted, v), &free);
        surface.rms_error(&refitted) - base <= tolerance
    };

    // Bisect on each side for the last value that still fits.
    let walk = |limit: f64| {
        if acceptable(limit) {
            return limit;
        }
        let (mut good, mut bad) = (centre, limit);
        for _ in 0..40 {
            let mid = 0.5 * (good + bad);
            if acceptable(mid) {
                good = mid;
            } else {
                bad = mid;
            }
        }
        good
    };

    (walk(lo_bound), walk(hi_bound))
}

/// Fit an Ornstein-Uhlenbeck process to an observed path, by least squares.
///
/// The discrete transition is an AR(1),
///
/// ```text
///     X[i+1] = c + phi X[i] + noise,   phi = exp(-kappa dt),
/// ```
///
/// so regressing each observation on the previous one recovers everything. This
/// is also the maximum likelihood estimator, because the transition is Gaussian
/// with constant variance and least squares and MLE coincide there.
///
/// It is the estimator every desk uses, and the fitting and testing chapter is
/// largely about the trouble it causes.
pub fn fit_ou(path: &[f64], dt: f64) -> OuFit {
    let n = path.len().saturating_sub(1);
    if n < 2 || dt <= 0.0 {
        return OuFit { kappa: f64::NAN, theta: f64::NAN, sigma: f64::NAN };
    }

    let (x, y) = (&path[..n], &path[1..]);
    let count = n as f64;
    let mean_x = x.iter().sum::<f64>() / count;
    let mean_y = y.iter().sum::<f64>() / count;

    let mut sxx = 0.0;
    let mut sxy = 0.0;
    for i in 0..n {
        let dx = x[i] - mean_x;
        sxx += dx * dx;
        sxy += dx * (y[i] - mean_y);
    }
    if sxx <= 0.0 {
        return OuFit { kappa: f64::NAN, theta: f64::NAN, sigma: f64::NAN };
    }

    let phi = sxy / sxx;
    let intercept = mean_y - phi * mean_x;

    // Residual variance, which inverts to the diffusion coefficient.
    let mut residual = 0.0;
    for i in 0..n {
        let e = y[i] - intercept - phi * x[i];
        residual += e * e;
    }
    let residual_variance = residual / (count - 2.0);

    // A path that failed to mean revert at all comes back as phi >= 1, which has
    // no kappa behind it. Reported as infinite rather than silently clamped.
    if phi <= 0.0 {
        return OuFit { kappa: f64::INFINITY, theta: mean_x, sigma: f64::NAN };
    }
    let kappa = -phi.ln() / dt;
    let theta = intercept / (1.0 - phi);
    let sigma = if kappa > 0.0 && phi < 1.0 {
        (residual_variance * 2.0 * kappa / (1.0 - phi * phi)).sqrt()
    } else {
        (residual_variance / dt).sqrt()
    };

    OuFit { kappa, theta, sigma }
}

/// The average estimated mean reversion over many independent histories.
///
/// `window` is the length of each history in years and `per_year` how often it
/// is sampled, so the number of observations is their product. The fitting and
/// testing chapter's claim is that the answer depends on `window` and
/// essentially not on `per_year`, which is the opposite of the usual instinct.
pub fn mean_reversion_estimate(
    true_kappa: f64,
    window: f64,
    per_year: f64,
    trials: usize,
    seed: u64,
) -> f64 {
    let dt = 1.0 / per_year;
    let steps = (window * per_year).round() as usize;
    if steps < 3 {
        return f64::NAN;
    }

    let mut rng = Rng::new(seed);
    let mut total = 0.0;
    let mut counted = 0.0;
    for _ in 0..trials {
        // Started in the stationary distribution, so nothing here is a burn-in
        // artefact.
        let stationary_sd = 1.0 / (2.0 * true_kappa).sqrt();
        let start = stationary_sd * rng.next_normal();
        let path = ou_path(true_kappa, 0.0, 1.0, start, dt, steps, &mut rng);
        let fit = fit_ou(&path, dt);
        if fit.kappa.is_finite() {
            total += fit.kappa;
            counted += 1.0;
        }
    }
    total / counted
}

/// Repeat the experiment: simulate, test, estimate, and summarise.
pub fn inference_quality(
    kappa: f64,
    years: f64,
    dt: f64,
    trials: usize,
    seed: u64,
) -> InferenceQuality {
    let mut rejections = 0usize;
    let mut half_lives = Vec::with_capacity(trials);

    for i in 0..trials {
        let path = stationary_ou_sample(kappa, 0.01, years, dt, seed.wrapping_add(i as u64 * 7919));
        let test = unit_root_test(&path, dt);
        if test.rejects_random_walk {
            rejections += 1;
        }
        if test.implied_kappa > 0.0 {
            half_lives.push(std::f64::consts::LN_2 / test.implied_kappa);
        } else {
            // A non-positive estimate means no reversion was found at all.
            half_lives.push(f64::INFINITY);
        }
    }

    half_lives.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let at = |q: f64| half_lives[((trials as f64 - 1.0) * q).round() as usize];

    InferenceQuality {
        power: rejections as f64 / trials as f64,
        half_life_interval: (at(0.05), at(0.95)),
        median_half_life: at(0.5),
    }
}

pub fn cir_pit_shape(
    kappa: f64,
    theta: f64,
    xi: f64,
    mu: f64,
    dt: f64,
    n: usize,
    bins: usize,
    seed: u64,
) -> PitShape {
    let path = cir_path(kappa, theta, xi, mu, dt, n, seed);
    let returns: Vec<f64> = path.windows(2).map(|w| (w[1] / w[0]).ln()).collect();
    let realised = (returns.iter().map(|r| r * r).sum::<f64>() / returns.len() as f64 / dt).sqrt();
    let pit = lognormal_pit(&path, dt, mu, realised);

    let mut counts = vec![0usize; bins];
    for &u in &pit {
        let b = ((u * bins as f64) as usize).min(bins - 1);
        counts[b] += 1;
    }
    let scale = bins as f64 / pit.len() as f64;
    let density = counts.iter().map(|&c| c as f64 * scale).collect();

    let tail_below = pit.iter().filter(|&&u| u < 0.01).count();
    let tail_above = pit.iter().filter(|&&u| u > 0.99).count();

    PitShape {
        density,
        ks: ks_statistic(&pit),
        tail_below,
        tail_above,
        tail_expected: pit.len() as f64 / 100.0,
    }
}

/// The Kolmogorov-Smirnov statistic scaled by `sqrt(n)`.
///
/// Comparable across sample sizes. The asymptotic five percent critical value is
/// `1.358`, so anything much beyond that is a model the data rejects.
pub fn ks_statistic(values: &[f64]) -> f64 {
    ks_uniform(values) * (values.len() as f64).sqrt()
}