/// Estimate [`Backbone`] from a daily rate series.
///
/// `window` is the number of business days in each block, and blocks do not
/// overlap. Rates are absolute, so `0.0472` for `4.72%`.
pub fn backbone(rates: &[f64], window: usize) -> Backbone {
    assert!(window >= 5, "a window of {window} days is too short to measure a volatility");

    let (mut xs, mut ys) = (Vec::new(), Vec::new());
    let (mut low, mut high) = (f64::INFINITY, f64::NEG_INFINITY);
    let mut start = 0;
    while start + window + 1 <= rates.len() {
        let block = &rates[start..start + window + 1];
        let level = block.iter().sum::<f64>() / block.len() as f64;
        let changes: Vec<f64> = block.windows(2).map(|w| w[1] - w[0]).collect();
        let mean = changes.iter().sum::<f64>() / changes.len() as f64;
        let variance =
            changes.iter().map(|c| (c - mean).powi(2)).sum::<f64>() / changes.len() as f64;
        let vol = (variance * 252.0).sqrt();
        if level > 0.0 && vol > 0.0 {
            xs.push(level.ln());
            ys.push(vol.ln());
            low = low.min(level);
            high = high.max(level);
        }
        start += window;
    }

    let n = xs.len();
    assert!(n > 2, "only {n} usable windows");
    let mx = xs.iter().sum::<f64>() / n as f64;
    let my = ys.iter().sum::<f64>() / n as f64;
    let sxx: f64 = xs.iter().map(|x| (x - mx).powi(2)).sum();
    let sxy: f64 = xs.iter().zip(&ys).map(|(x, y)| (x - mx) * (y - my)).sum();
    let beta = sxy / sxx;
    let residual: f64 = xs
        .iter()
        .zip(&ys)
        .map(|(x, y)| (y - (my + beta * (x - mx))).powi(2))
        .sum();
    let standard_error = (residual / (n as f64 - 2.0) / sxx).sqrt();

    Backbone { beta, standard_error, windows: n, level_low: low, level_high: high }
}

/// Hagan's lognormal (Black) implied volatility.
///
/// The standard 2002 expansion. Written in the usual two pieces: a
/// leading-order term that carries the smile's shape, and a bracket of
/// first-order corrections in the expiry.
///
/// The at-the-money case is handled separately because the general formula
/// has a removable singularity there — `z / x(z)` is `0/0` when the strike
/// equals the forward. Taking the limit rather than nudging the strike is
/// what keeps the at-the-money volatility exact, and the at-the-money point
/// is the one the whole surface is anchored on.
pub fn implied_vol(&self, f: f64, k: f64, t: f64) -> Option<f64> {
    if !(f > 0.0) || !(k > 0.0) || !(t > 0.0) {
        return None;
    }
    let Sabr { alpha, beta, rho, nu } = *self;
    if !(alpha > 0.0) || !(0.0..=1.0).contains(&beta) || rho.abs() >= 1.0 || nu < 0.0 {
        return None;
    }

    let one_m_b = 1.0 - beta;
    let fk = (f * k).powf(one_m_b);
    let log_fk = (f / k).ln();

    // The correction bracket, common to both branches.
    let term1 = one_m_b * one_m_b / 24.0 * alpha * alpha / fk;
    let term2 = 0.25 * rho * beta * nu * alpha / fk.sqrt();
    let term3 = (2.0 - 3.0 * rho * rho) / 24.0 * nu * nu;
    let correction = 1.0 + (term1 + term2 + term3) * t;

    // At the money, in the limit k -> f.
    if log_fk.abs() < 1e-9 {
        return Some(alpha / f.powf(one_m_b) * correction);
    }

    let z = nu / alpha * fk.sqrt() * log_fk;
    // x(z), the integral that converts the log-moneyness into the
    // volatility's own coordinate.
    let x = (((1.0 - 2.0 * rho * z + z * z).sqrt() + z - rho) / (1.0 - rho)).ln();
    if !x.is_finite() || x.abs() < 1e-12 {
        return Some(alpha / f.powf(one_m_b) * correction);
    }

    let denom = fk.sqrt()
        * (1.0
            + one_m_b * one_m_b / 24.0 * log_fk * log_fk
            + one_m_b.powi(4) / 1920.0 * log_fk.powi(4));

    let vol = alpha / denom * (z / x) * correction;
    if vol.is_finite() && vol > 0.0 {
        Some(vol)
    } else {
        None
    }
}

/// The at-the-money volatility's sensitivity to the forward: the backbone.
///
/// This is the quantity the smile dynamics chapter is about. It is not the
/// smile's slope in the strike — that is the skew, measured at one instant
/// with the forward held still. The backbone is what happens to the at-the-
/// money volatility when the market actually moves, and it is what
/// determines whether the model's delta is right.
///
/// Computed by moving the forward and asking the model again, rather than by
/// differentiating the formula, so that it measures the model's behaviour
/// and not our algebra.
pub fn backbone(&self, f: f64, t: f64, h: f64) -> Option<f64> {
    let up = self.atm_vol(f * (1.0 + h), t)?;
    let dn = self.atm_vol(f * (1.0 - h), t)?;
    Some((up - dn) / (2.0 * h * f))
}

/// Fit `alpha` and `rho` so that the model reproduces a given at-the-money
/// volatility and at-the-money skew, with `beta` and `nu` held where they were
/// put.
///
/// This is the operation the smile dynamics chapter turns into its main point.
/// `beta` is not determined by the smile: for any `beta` one likes, there is a
/// `rho` that reproduces today's skew, and the resulting models are near enough
/// indistinguishable from today's quotes. They then disagree completely about
/// what happens when the forward moves, because `beta` is the backbone and
/// `rho` is not. So calibrating perfectly to the market leaves the dynamics
/// undetermined, and the choice has to be made on other grounds.
///
/// Alternates two one-dimensional solves rather than searching in two
/// dimensions at once. It converges quickly because the coupling is weak in the
/// right direction: `alpha` mostly sets the level and `rho` mostly sets the
/// slope.
pub fn calibrate_alpha_rho(
    target_atm: f64,
    target_skew: f64,
    beta: f64,
    nu: f64,
    f: f64,
    t: f64,
) -> Option<Sabr> {
    let h = 0.01;
    let mut s = Sabr { alpha: target_atm * f.powf(1.0 - beta), beta, rho: 0.0, nu };

    for _ in 0..60 {
        // alpha for the level: the at-the-money volatility is increasing in it.
        let (mut lo, mut hi) = (1e-6, 10.0 * target_atm.max(1.0) * f.powf(1.0 - beta));
        for _ in 0..100 {
            let mid = 0.5 * (lo + hi);
            let trial = Sabr { alpha: mid, ..s };
            match trial.atm_vol(f, t) {
                Some(v) if v < target_atm => lo = mid,
                Some(_) => hi = mid,
                None => lo = mid,
            }
        }
        s.alpha = 0.5 * (lo + hi);

        // rho for the slope: the skew is increasing in it.
        let (mut lo, mut hi) = (-0.999, 0.999);
        for _ in 0..100 {
            let mid = 0.5 * (lo + hi);
            let trial = Sabr { rho: mid, ..s };
            match trial.skew(f, t, h) {
                Some(v) if v < target_skew => lo = mid,
                Some(_) => hi = mid,
                None => break,
            }
        }
        s.rho = 0.5 * (lo + hi);
    }

    let ok = s.atm_vol(f, t).is_some_and(|v| (v - target_atm).abs() < 1e-6)
        && s.skew(f, t, h).is_some_and(|v| (v - target_skew).abs() < 1e-7);
    ok.then_some(s)
}

/// The smile's slope in the strike, at the money: the skew.
///
/// Measured with the forward held fixed, which is the whole distinction from
/// [`Sabr::backbone`].
pub fn skew(&self, f: f64, t: f64, h: f64) -> Option<f64> {
    let up = self.implied_vol(f, f * (1.0 + h), t)?;
    let dn = self.implied_vol(f, f * (1.0 - h), t)?;
    Some((up - dn) / (2.0 * h * f))
}