chore: update ferro-ta version to 1.1.3

- Bumped version numbers across Cargo.toml, Cargo.lock, pyproject.toml, and conda/meta.yaml to 1.1.3.
- Added new features including American option pricing, digital options, extended Greeks, and historical volatility estimators.
- Enhanced documentation and tests for new functionalities.
- Updated CHANGELOG.md to reflect changes for version 1.1.3.

crates/ferro_ta_core/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "ferro_ta_core"
-version = "1.1.2"
+version = "1.1.3"
 edition = "2021"
 description = "Pure Rust core indicator library — no PyO3, no numpy dependency"
 license = "MIT"

crates/ferro_ta_core/README.md

@@ -13,7 +13,7 @@ PyO3, NumPy, or Python runtime dependency, which makes it a good fit for:
 
 ```toml
 [dependencies]
-ferro_ta_core = "1.1.2"
+ferro_ta_core = "1.1.3"
 ```
 
 ## Design

crates/ferro_ta_core/src/options/american.rs

@@ -0,0 +1,410 @@
+//! American option pricing via the Barone-Adesi-Whaley (1987) quadratic approximation.
+
+use super::normal::cdf;
+use super::pricing::black_scholes_price;
+use super::OptionKind;
+
+fn invalid_inputs(spot: f64, strike: f64, time_to_expiry: f64, volatility: f64) -> bool {
+    !spot.is_finite()
+        || !strike.is_finite()
+        || !time_to_expiry.is_finite()
+        || !volatility.is_finite()
+        || spot <= 0.0
+        || strike <= 0.0
+        || time_to_expiry < 0.0
+        || volatility < 0.0
+}
+
+/// Compute d1 for BSM given spot S* (used inside the Newton-Raphson loop).
+fn d1_fn(s: f64, strike: f64, rate: f64, carry: f64, time_to_expiry: f64, volatility: f64) -> f64 {
+    let sigma_sqrt_t = volatility * time_to_expiry.sqrt();
+    ((s / strike).ln() + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry)
+        / sigma_sqrt_t
+}
+
+/// Find the critical spot price S* for American call early exercise using Newton-Raphson.
+///
+/// S* satisfies: C(S*) - (S* - K) = (S*/q2) * (1 - e^{-q*T} * N(d1(S*)))
+/// Rearranged as F(S*) = 0:
+/// F(x) = C(x) - (x - K) - (x/q2) * (1 - carry_discount * N(d1(x))) = 0
+fn find_critical_call(
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    q2: f64,
+) -> f64 {
+    let carry_discount = (-carry * time_to_expiry).exp();
+
+    // Initial guess: S* ≈ K * q2 / (q2 - 1), clamped to be above strike
+    let mut s = if q2 > 1.0 {
+        strike * q2 / (q2 - 1.0)
+    } else {
+        // q2 <= 1 means the denominator is small/negative; fall back to a safe value
+        strike * 2.0
+    };
+    // Ensure starting guess is positive
+    if s <= 0.0 {
+        s = strike * 1.5;
+    }
+
+    for _ in 0..50 {
+        let c = black_scholes_price(
+            s,
+            strike,
+            rate,
+            carry,
+            time_to_expiry,
+            volatility,
+            OptionKind::Call,
+        );
+        let d1 = d1_fn(s, strike, rate, carry, time_to_expiry, volatility);
+        let nd1 = cdf(d1);
+        let lhs = c - (s - strike);
+        let rhs = (s / q2) * (1.0 - carry_discount * nd1);
+        let f = lhs - rhs;
+
+        // Derivative of F with respect to s:
+        // dC/ds = e^{-q*T} * N(d1) (BSM delta for call)
+        // d(s - K)/ds = 1
+        // d(rhs)/ds = (1/q2) * (1 - carry_discount * N(d1))
+        //            + (s/q2) * (-carry_discount * phi(d1) / (s * vol * sqrt(T)))
+        //           = (1/q2) * (1 - carry_discount * N(d1)) - carry_discount * phi(d1) / (q2 * vol * sqrt(T))
+        let sigma_sqrt_t = volatility * time_to_expiry.sqrt();
+        let phi_d1 = super::normal::pdf(d1);
+        let d_lhs_ds = carry_discount * nd1 - 1.0;
+        let d_rhs_ds = (1.0 / q2) * (1.0 - carry_discount * nd1)
+            - carry_discount * phi_d1 / (q2 * sigma_sqrt_t);
+        let df = d_lhs_ds - d_rhs_ds;
+
+        if df.abs() < 1e-14 {
+            break;
+        }
+        let step = f / df;
+        s -= step;
+        // Keep s positive
+        if s <= 0.0 {
+            s = strike * 0.1;
+        }
+        if step.abs() < 1e-8 {
+            break;
+        }
+    }
+    s
+}
+
+/// Find the critical spot price S** for American put early exercise using Newton-Raphson.
+///
+/// S** satisfies: P(S**) - (K - S**) = -(S**/q1) * (1 - e^{-q*T} * N(-d1(S**)))
+/// F(x) = P(x) - (K - x) + (x/q1) * (1 - carry_discount * N(-d1(x))) = 0
+fn find_critical_put(
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    q1: f64,
+) -> f64 {
+    let carry_discount = (-carry * time_to_expiry).exp();
+
+    // Initial guess for put: S** ≈ K * q1 / (q1 - 1)
+    // q1 is negative, so q1 - 1 < 0, and the guess should be below strike.
+    let mut s = if (q1 - 1.0).abs() > 1e-10 {
+        strike * q1 / (q1 - 1.0)
+    } else {
+        strike * 0.5
+    };
+    if s <= 0.0 || s >= strike {
+        s = strike * 0.5;
+    }
+
+    for _ in 0..50 {
+        let p = black_scholes_price(
+            s,
+            strike,
+            rate,
+            carry,
+            time_to_expiry,
+            volatility,
+            OptionKind::Put,
+        );
+        let d1 = d1_fn(s, strike, rate, carry, time_to_expiry, volatility);
+        let n_neg_d1 = cdf(-d1);
+        let lhs = p - (strike - s);
+        // rhs = -(s/q1) * (1 - carry_discount * N(-d1))
+        let rhs = -(s / q1) * (1.0 - carry_discount * n_neg_d1);
+        let f = lhs - rhs;
+
+        // Derivative:
+        // dP/ds = -e^{-q*T} * N(-d1)  (BSM delta for put = e^{-q*T}*(N(d1)-1))
+        // d(K - s)/ds = -1  so  d(lhs)/ds = dP/ds - (-1) = dP/ds + 1
+        // d(rhs)/ds = -(1/q1)*(1 - carry_discount*N(-d1))
+        //           + -(s/q1)*carry_discount*phi(d1)/(s*vol*sqrt(T))   [since d(N(-d1))/ds = -phi(d1)*dd1/ds]
+        //           = -(1/q1)*(1 - carry_discount*N(-d1))
+        //           - carry_discount*phi(d1)/(q1*vol*sqrt(T))
+        let sigma_sqrt_t = volatility * time_to_expiry.sqrt();
+        let phi_d1 = super::normal::pdf(d1);
+        let d_lhs_ds = -carry_discount * n_neg_d1 + 1.0;
+        let d_rhs_ds = -(1.0 / q1) * (1.0 - carry_discount * n_neg_d1)
+            - carry_discount * phi_d1 / (q1 * sigma_sqrt_t);
+        let df = d_lhs_ds - d_rhs_ds;
+
+        if df.abs() < 1e-14 {
+            break;
+        }
+        let step = f / df;
+        s -= step;
+        if s <= 0.0 {
+            s = strike * 0.01;
+        }
+        if s >= strike {
+            s = strike * 0.99;
+        }
+        if step.abs() < 1e-8 {
+            break;
+        }
+    }
+    s
+}
+
+/// American option price using the Barone-Adesi-Whaley (1987) quadratic approximation.
+///
+/// # Parameters
+/// - `spot`: current underlying price
+/// - `strike`: option strike price
+/// - `rate`: risk-free rate (annualized, decimal)
+/// - `carry`: continuous dividend yield / carry rate
+/// - `time_to_expiry`: time to expiry in years
+/// - `volatility`: implied vol (annualized, decimal)
+/// - `kind`: call or put
+pub fn american_price_baw(
+    spot: f64,
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    kind: OptionKind,
+) -> f64 {
+    if invalid_inputs(spot, strike, time_to_expiry, volatility)
+        || !rate.is_finite()
+        || !carry.is_finite()
+    {
+        return f64::NAN;
+    }
+
+    // At expiry: immediate exercise value
+    if time_to_expiry == 0.0 {
+        return match kind {
+            OptionKind::Call => (spot - strike).max(0.0),
+            OptionKind::Put => (strike - spot).max(0.0),
+        };
+    }
+
+    // At zero vol: deterministic — exercise if ITM
+    if volatility == 0.0 {
+        return match kind {
+            OptionKind::Call => (spot - strike).max(0.0),
+            OptionKind::Put => (strike - spot).max(0.0),
+        };
+    }
+
+    let european = black_scholes_price(spot, strike, rate, carry, time_to_expiry, volatility, kind);
+
+    match kind {
+        OptionKind::Call => {
+            // No early exercise premium when there are no dividends (carry == 0 means q==0
+            // in BSM parameterisation where carry = q).
+            if carry <= 0.0 {
+                return european;
+            }
+
+            let sigma2 = volatility * volatility;
+            let m = 2.0 * rate / sigma2;
+            let n = 2.0 * (rate - carry) / sigma2;
+            let h = 1.0 - (-rate * time_to_expiry).exp();
+
+            if h.abs() < 1e-14 {
+                return european;
+            }
+
+            let discriminant = (n - 1.0) * (n - 1.0) + 4.0 * m / h;
+            if discriminant < 0.0 {
+                return european;
+            }
+
+            let q2 = (-(n - 1.0) + discriminant.sqrt()) / 2.0;
+
+            // Find critical price S*
+            let s_star = find_critical_call(strike, rate, carry, time_to_expiry, volatility, q2);
+
+            if s_star <= strike {
+                // Degenerate critical price; fall back to European
+                return european;
+            }
+
+            // A2 = (S*/q2) * (1 - e^{-q*T} * N(d1(S*)))
+            let carry_discount = (-carry * time_to_expiry).exp();
+            let d1_star = d1_fn(s_star, strike, rate, carry, time_to_expiry, volatility);
+            let a2 = (s_star / q2) * (1.0 - carry_discount * cdf(d1_star));
+
+            if spot >= s_star {
+                // Immediate exercise is optimal
+                (spot - strike).max(0.0)
+            } else {
+                (european + a2 * (spot / s_star).powf(q2)).max(european)
+            }
+        }
+
+        OptionKind::Put => {
+            // No early exercise when rate == 0 (no time value of money)
+            if rate <= 0.0 {
+                return european;
+            }
+
+            let sigma2 = volatility * volatility;
+            let m = 2.0 * rate / sigma2;
+            let n = 2.0 * (rate - carry) / sigma2;
+            let h = 1.0 - (-rate * time_to_expiry).exp();
+
+            if h.abs() < 1e-14 {
+                return european;
+            }
+
+            let discriminant = (n - 1.0) * (n - 1.0) + 4.0 * m / h;
+            if discriminant < 0.0 {
+                return european;
+            }
+
+            let q1 = (-(n - 1.0) - discriminant.sqrt()) / 2.0;
+
+            // Find critical price S**
+            let s_star_star =
+                find_critical_put(strike, rate, carry, time_to_expiry, volatility, q1);
+
+            if s_star_star <= 0.0 || s_star_star >= strike {
+                return european;
+            }
+
+            // A1 = -(S**/q1) * (1 - e^{-q*T} * N(-d1(S**)))
+            let carry_discount = (-carry * time_to_expiry).exp();
+            let d1_star = d1_fn(s_star_star, strike, rate, carry, time_to_expiry, volatility);
+            let a1 = -(s_star_star / q1) * (1.0 - carry_discount * cdf(-d1_star));
+
+            if spot <= s_star_star {
+                // Immediate exercise is optimal
+                (strike - spot).max(0.0)
+            } else {
+                (european + a1 * (spot / s_star_star).powf(q1)).max(european)
+            }
+        }
+    }
+}
+
+/// Early exercise premium = american_price - european_bsm_price.
+///
+/// Always non-negative for valid inputs.
+pub fn early_exercise_premium(
+    spot: f64,
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    kind: OptionKind,
+) -> f64 {
+    let american = american_price_baw(spot, strike, rate, carry, time_to_expiry, volatility, kind);
+    let european = black_scholes_price(spot, strike, rate, carry, time_to_expiry, volatility, kind);
+    if american.is_nan() || european.is_nan() {
+        return f64::NAN;
+    }
+    (american - european).max(0.0)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::options::OptionKind;
+
+    #[test]
+    fn american_call_gte_european_call() {
+        let european = crate::options::pricing::black_scholes_price(
+            100.0,
+            100.0,
+            0.05,
+            0.03,
+            1.0,
+            0.2,
+            OptionKind::Call,
+        );
+        let american = american_price_baw(100.0, 100.0, 0.05, 0.03, 1.0, 0.2, OptionKind::Call);
+        assert!(american >= european - 1e-10);
+    }
+
+    #[test]
+    fn american_put_gte_european_put() {
+        let european = crate::options::pricing::black_scholes_price(
+            100.0,
+            100.0,
+            0.05,
+            0.0,
+            1.0,
+            0.2,
+            OptionKind::Put,
+        );
+        let american = american_price_baw(100.0, 100.0, 0.05, 0.0, 1.0, 0.2, OptionKind::Put);
+        assert!(american >= european - 1e-10);
+    }
+
+    #[test]
+    fn early_exercise_premium_nonneg() {
+        let prem = early_exercise_premium(100.0, 100.0, 0.05, 0.03, 1.0, 0.2, OptionKind::Call);
+        assert!(prem >= 0.0);
+    }
+
+    #[test]
+    fn american_call_no_dividends_equals_european() {
+        // With no dividends (carry == 0), no early exercise is optimal for calls
+        let european = crate::options::pricing::black_scholes_price(
+            100.0,
+            100.0,
+            0.05,
+            0.0,
+            1.0,
+            0.2,
+            OptionKind::Call,
+        );
+        let american = american_price_baw(100.0, 100.0, 0.05, 0.0, 1.0, 0.2, OptionKind::Call);
+        assert!((american - european).abs() < 1e-10);
+    }
+
+    #[test]
+    fn american_price_returns_nan_for_invalid() {
+        let price = american_price_baw(-1.0, 100.0, 0.05, 0.0, 1.0, 0.2, OptionKind::Call);
+        assert!(price.is_nan());
+    }
+
+    #[test]
+    fn american_price_at_expiry_is_intrinsic() {
+        let call = american_price_baw(110.0, 100.0, 0.05, 0.03, 0.0, 0.2, OptionKind::Call);
+        assert!((call - 10.0).abs() < 1e-10);
+        let put = american_price_baw(90.0, 100.0, 0.05, 0.0, 0.0, 0.2, OptionKind::Put);
+        assert!((put - 10.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn american_put_itm_has_positive_premium() {
+        // Deep ITM put with high rate should have meaningful early exercise premium
+        let prem = early_exercise_premium(80.0, 100.0, 0.10, 0.0, 1.0, 0.2, OptionKind::Put);
+        assert!(prem >= 0.0);
+    }
+
+    #[test]
+    fn american_prices_are_finite_for_valid_inputs() {
+        let call = american_price_baw(100.0, 100.0, 0.05, 0.02, 1.0, 0.25, OptionKind::Call);
+        let put = american_price_baw(100.0, 100.0, 0.05, 0.02, 1.0, 0.25, OptionKind::Put);
+        assert!(call.is_finite());
+        assert!(put.is_finite());
+    }
+}

crates/ferro_ta_core/src/options/digital.rs

@@ -0,0 +1,382 @@
+//! Digital (binary) option pricing.
+
+use super::normal::cdf;
+use super::OptionKind;
+
+/// Type of digital option payoff.
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+pub enum DigitalKind {
+    /// Pays 1 unit of cash if option expires in the money.
+    CashOrNothing,
+    /// Pays the underlying asset if option expires in the money.
+    AssetOrNothing,
+}
+
+fn invalid_inputs(spot: f64, strike: f64, time_to_expiry: f64, volatility: f64) -> bool {
+    !spot.is_finite()
+        || !strike.is_finite()
+        || !time_to_expiry.is_finite()
+        || !volatility.is_finite()
+        || spot <= 0.0
+        || strike <= 0.0
+        || time_to_expiry < 0.0
+        || volatility < 0.0
+}
+
+/// Price a digital (binary) option under BSM.
+///
+/// # Parameters
+/// - `spot`: current underlying price
+/// - `strike`: option strike price
+/// - `rate`: risk-free rate (annualized, decimal)
+/// - `carry`: continuous dividend yield / carry rate
+/// - `time_to_expiry`: time to expiry in years
+/// - `volatility`: implied vol (annualized, decimal)
+/// - `option_kind`: call or put
+/// - `digital_kind`: cash-or-nothing or asset-or-nothing
+#[allow(clippy::too_many_arguments)]
+pub fn digital_price(
+    spot: f64,
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    option_kind: OptionKind,
+    digital_kind: DigitalKind,
+) -> f64 {
+    if invalid_inputs(spot, strike, time_to_expiry, volatility)
+        || !rate.is_finite()
+        || !carry.is_finite()
+    {
+        return f64::NAN;
+    }
+
+    // At expiry: pay intrinsic based on ITM status
+    if time_to_expiry == 0.0 {
+        let itm = match option_kind {
+            OptionKind::Call => spot > strike,
+            OptionKind::Put => spot < strike,
+        };
+        return if itm {
+            match digital_kind {
+                DigitalKind::CashOrNothing => 1.0,
+                DigitalKind::AssetOrNothing => spot,
+            }
+        } else {
+            0.0
+        };
+    }
+
+    let discount = (-rate * time_to_expiry).exp();
+    let carry_discount = (-carry * time_to_expiry).exp();
+
+    // At zero vol: deterministic payoff
+    if volatility == 0.0 {
+        let forward = spot * (carry_discount / discount); // S * e^{(r-q)*T} equivalent: S*e^{-q*T}/e^{-r*T}
+                                                          // forward = S * e^{(r-q)*T}; ITM if forward > K for call
+        let itm = match option_kind {
+            OptionKind::Call => spot * carry_discount > strike * discount,
+            OptionKind::Put => spot * carry_discount < strike * discount,
+        };
+        let _ = forward; // suppress unused warning
+        return if itm {
+            match digital_kind {
+                DigitalKind::CashOrNothing => discount,
+                DigitalKind::AssetOrNothing => spot * carry_discount,
+            }
+        } else {
+            0.0
+        };
+    }
+
+    let sqrt_t = time_to_expiry.sqrt();
+    let sigma_sqrt_t = volatility * sqrt_t;
+    let d1 = ((spot / strike).ln()
+        + (rate - carry + 0.5 * volatility * volatility) * time_to_expiry)
+        / sigma_sqrt_t;
+    let d2 = d1 - sigma_sqrt_t;
+
+    match digital_kind {
+        DigitalKind::CashOrNothing => match option_kind {
+            OptionKind::Call => discount * cdf(d2),
+            OptionKind::Put => discount * cdf(-d2),
+        },
+        DigitalKind::AssetOrNothing => match option_kind {
+            OptionKind::Call => spot * carry_discount * cdf(d1),
+            OptionKind::Put => spot * carry_discount * cdf(-d1),
+        },
+    }
+}
+
+/// Compute numerical delta, gamma, and vega for a digital option.
+///
+/// Uses central finite differences:
+/// - delta/gamma: bump spot by ε = spot * 1e-3
+/// - vega: bump volatility by 1e-3
+///
+/// Returns `(delta, gamma, vega)`.
+#[allow(clippy::too_many_arguments)]
+pub fn digital_greeks(
+    spot: f64,
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    option_kind: OptionKind,
+    digital_kind: DigitalKind,
+) -> (f64, f64, f64) {
+    let eps = spot * 1e-3;
+    if eps <= 0.0 {
+        return (f64::NAN, f64::NAN, f64::NAN);
+    }
+
+    let price_mid = digital_price(
+        spot,
+        strike,
+        rate,
+        carry,
+        time_to_expiry,
+        volatility,
+        option_kind,
+        digital_kind,
+    );
+    let price_up = digital_price(
+        spot + eps,
+        strike,
+        rate,
+        carry,
+        time_to_expiry,
+        volatility,
+        option_kind,
+        digital_kind,
+    );
+    let price_dn = digital_price(
+        spot - eps,
+        strike,
+        rate,
+        carry,
+        time_to_expiry,
+        volatility,
+        option_kind,
+        digital_kind,
+    );
+
+    let delta = (price_up - price_dn) / (2.0 * eps);
+    let gamma = (price_up - 2.0 * price_mid + price_dn) / (eps * eps);
+
+    let vol_bump = 1e-3;
+    let vega = if volatility + vol_bump > 0.0 && volatility - vol_bump > 0.0 {
+        let price_vup = digital_price(
+            spot,
+            strike,
+            rate,
+            carry,
+            time_to_expiry,
+            volatility + vol_bump,
+            option_kind,
+            digital_kind,
+        );
+        let price_vdn = digital_price(
+            spot,
+            strike,
+            rate,
+            carry,
+            time_to_expiry,
+            volatility - vol_bump,
+            option_kind,
+            digital_kind,
+        );
+        (price_vup - price_vdn) / (2.0 * vol_bump)
+    } else {
+        // vol too close to zero; one-sided bump
+        let price_vup = digital_price(
+            spot,
+            strike,
+            rate,
+            carry,
+            time_to_expiry,
+            volatility + vol_bump,
+            option_kind,
+            digital_kind,
+        );
+        (price_vup - price_mid) / vol_bump
+    };
+
+    (delta, gamma, vega)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::options::OptionKind;
+
+    #[test]
+    fn cash_or_nothing_call_atm() {
+        // ATM cash-or-nothing call: price = e^{-rT} * N(d2)
+        // At S=K=100, r=0.05, q=0, T=1, σ=0.2:
+        //   d1 = (0 + 0.07) / 0.2 = 0.35,  d2 = 0.15  →  N(0.15) ≈ 0.5596
+        //   price ≈ e^{-0.05} * 0.5596 ≈ 0.532
+        let price = digital_price(
+            100.0,
+            100.0,
+            0.05,
+            0.0,
+            1.0,
+            0.2,
+            OptionKind::Call,
+            DigitalKind::CashOrNothing,
+        );
+        assert!(
+            price > 0.0 && price < 1.0,
+            "price should be between 0 and 1"
+        );
+        assert!((price - 0.532).abs() < 0.01, "price ≈ 0.532, got {price}");
+    }
+
+    #[test]
+    fn asset_or_nothing_call_at_zero_vol() {
+        // At zero vol, ITM asset-or-nothing call should equal S * e^{-q*T}
+        let price = digital_price(
+            110.0,
+            100.0,
+            0.05,
+            0.0,
+            1.0,
+            0.0,
+            OptionKind::Call,
+            DigitalKind::AssetOrNothing,
+        );
+        assert!((price - 110.0).abs() < 1e-6);
+    }
+
+    #[test]
+    fn digital_price_returns_nan_for_invalid() {
+        let price = digital_price(
+            -1.0,
+            100.0,
+            0.05,
+            0.0,
+            1.0,
+            0.2,
+            OptionKind::Call,
+            DigitalKind::CashOrNothing,
+        );
+        assert!(price.is_nan());
+    }
+
+    #[test]
+    fn cash_or_nothing_put_call_parity() {
+        // Cash-or-nothing call + cash-or-nothing put = e^{-rT}
+        let call = digital_price(
+            100.0,
+            100.0,
+            0.05,
+            0.02,
+            1.0,
+            0.25,
+            OptionKind::Call,
+            DigitalKind::CashOrNothing,
+        );
+        let put = digital_price(
+            100.0,
+            100.0,
+            0.05,
+            0.02,
+            1.0,
+            0.25,
+            OptionKind::Put,
+            DigitalKind::CashOrNothing,
+        );
+        let discount = (-0.05_f64).exp();
+        assert!((call + put - discount).abs() < 1e-10);
+    }
+
+    #[test]
+    fn asset_or_nothing_put_call_parity() {
+        // Asset-or-nothing call + asset-or-nothing put = S * e^{-q*T}
+        let s = 100.0_f64;
+        let q = 0.02_f64;
+        let call = digital_price(
+            s,
+            100.0,
+            0.05,
+            q,
+            1.0,
+            0.25,
+            OptionKind::Call,
+            DigitalKind::AssetOrNothing,
+        );
+        let put = digital_price(
+            s,
+            100.0,
+            0.05,
+            q,
+            1.0,
+            0.25,
+            OptionKind::Put,
+            DigitalKind::AssetOrNothing,
+        );
+        let expected = s * (-q).exp();
+        assert!((call + put - expected).abs() < 1e-10);
+    }
+
+    #[test]
+    fn digital_greeks_are_finite_for_valid_inputs() {
+        let (delta, gamma, vega) = digital_greeks(
+            100.0,
+            100.0,
+            0.05,
+            0.0,
+            1.0,
+            0.2,
+            OptionKind::Call,
+            DigitalKind::CashOrNothing,
+        );
+        assert!(delta.is_finite());
+        assert!(gamma.is_finite());
+        assert!(vega.is_finite());
+    }
+
+    #[test]
+    fn digital_at_expiry_itm_returns_intrinsic() {
+        let price = digital_price(
+            110.0,
+            100.0,
+            0.05,
+            0.0,
+            0.0,
+            0.2,
+            OptionKind::Call,
+            DigitalKind::CashOrNothing,
+        );
+        assert!((price - 1.0).abs() < 1e-10);
+        let price2 = digital_price(
+            110.0,
+            100.0,
+            0.05,
+            0.0,
+            0.0,
+            0.2,
+            OptionKind::Call,
+            DigitalKind::AssetOrNothing,
+        );
+        assert!((price2 - 110.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn digital_at_expiry_otm_returns_zero() {
+        let price = digital_price(
+            90.0,
+            100.0,
+            0.05,
+            0.0,
+            0.0,
+            0.2,
+            OptionKind::Call,
+            DigitalKind::CashOrNothing,
+        );
+        assert!((price - 0.0).abs() < 1e-10);
+    }
+}

crates/ferro_ta_core/src/options/greeks.rs

@@ -2,7 +2,7 @@
 
 use super::normal::{cdf, pdf};
 use super::pricing::{black_76_price, black_scholes_price};
-use super::{Greeks, OptionEvaluation, OptionKind, PricingModel};
+use super::{ExtendedGreeks, Greeks, OptionEvaluation, OptionKind, PricingModel};
 
 fn bs_inputs_valid(
     underlying: f64,
@@ -203,9 +203,94 @@ pub fn model_theta(input: OptionEvaluation) -> f64 {
     })
 }
 
+/// Extended Greeks under Black-Scholes-Merton (closed-form).
+///
+/// All inputs must be positive finite; returns NaN fields for invalid inputs.
+pub fn black_scholes_extended_greeks(
+    spot: f64,
+    strike: f64,
+    rate: f64,
+    dividend_yield: f64,
+    time_to_expiry: f64,
+    volatility: f64,
+    _kind: OptionKind,
+) -> ExtendedGreeks {
+    if !bs_inputs_valid(
+        spot,
+        strike,
+        rate,
+        dividend_yield,
+        time_to_expiry,
+        volatility,
+    ) {
+        return ExtendedGreeks {
+            vanna: f64::NAN,
+            volga: f64::NAN,
+            charm: f64::NAN,
+            speed: f64::NAN,
+            color: f64::NAN,
+        };
+    }
+
+    let sqrt_t = time_to_expiry.sqrt();
+    let sigma_sqrt_t = volatility * sqrt_t;
+    let carry_discount = (-dividend_yield * time_to_expiry).exp();
+    let d1 = ((spot / strike).ln()
+        + (rate - dividend_yield + 0.5 * volatility * volatility) * time_to_expiry)
+        / sigma_sqrt_t;
+    let d2 = d1 - sigma_sqrt_t;
+    let pdf_d1 = pdf(d1);
+
+    let gamma = carry_discount * pdf_d1 / (spot * sigma_sqrt_t);
+
+    let vanna = -carry_discount * pdf_d1 * d2 / volatility;
+    let volga = spot * carry_discount * pdf_d1 * sqrt_t * d1 * d2 / volatility;
+    let charm = -carry_discount
+        * pdf_d1
+        * (2.0 * (rate - dividend_yield) * time_to_expiry - d2 * sigma_sqrt_t)
+        / (2.0 * time_to_expiry * sigma_sqrt_t);
+    let speed = -gamma / spot * (d1 / sigma_sqrt_t + 1.0);
+    let color = -carry_discount * pdf_d1 / (2.0 * spot * time_to_expiry * sigma_sqrt_t)
+        * (2.0 * (rate - dividend_yield) * time_to_expiry + 1.0
+            - d1 * (2.0 * (rate - dividend_yield) * time_to_expiry - d2 * sigma_sqrt_t)
+                / sigma_sqrt_t);
+
+    ExtendedGreeks {
+        vanna,
+        volga,
+        charm,
+        speed,
+        color,
+    }
+}
+
+/// Model-dispatched extended Greeks.
+/// Only BSM is supported with closed-form; Black-76 is not yet supported (returns NaN).
+pub fn model_extended_greeks(input: OptionEvaluation) -> ExtendedGreeks {
+    let contract = input.contract;
+    match contract.model {
+        PricingModel::BlackScholes => black_scholes_extended_greeks(
+            contract.underlying,
+            contract.strike,
+            contract.rate,
+            contract.carry,
+            contract.time_to_expiry,
+            input.volatility,
+            contract.kind,
+        ),
+        PricingModel::Black76 => ExtendedGreeks {
+            vanna: f64::NAN,
+            volga: f64::NAN,
+            charm: f64::NAN,
+            speed: f64::NAN,
+            color: f64::NAN,
+        },
+    }
+}
+
 #[cfg(test)]
 mod tests {
-    use super::{black_76_greeks, black_scholes_greeks};
+    use super::{black_76_greeks, black_scholes_extended_greeks, black_scholes_greeks};
     use crate::options::OptionKind;
 
     #[test]
@@ -227,4 +312,16 @@ mod tests {
         assert!(g.theta.is_finite());
         assert!(g.rho.is_finite());
     }
+
+    #[test]
+    fn extended_greeks_finite_for_valid_inputs() {
+        let eg = black_scholes_extended_greeks(100.0, 100.0, 0.05, 0.0, 1.0, 0.2, OptionKind::Call);
+        assert!(eg.vanna.is_finite());
+        assert!(eg.volga.is_finite());
+        assert!(eg.charm.is_finite());
+        assert!(eg.speed.is_finite());
+        assert!(eg.color.is_finite());
+        // Volga must be positive (convex in vol)
+        assert!(eg.volga >= 0.0);
+    }
 }

crates/ferro_ta_core/src/options/mod.rs

@@ -4,11 +4,15 @@
 //! IV-series helpers, and smile/chain utilities. The public API is scalar-first
 //! and is used by the PyO3 bridge to build vectorized batch functions.
 
+pub mod american;
 pub mod chain;
+pub mod digital;
 pub mod greeks;
 pub mod iv;
 pub mod normal;
+pub mod payoff;
 pub mod pricing;
+pub mod realized_vol;
 pub mod surface;
 
 /// Option side.
@@ -49,6 +53,16 @@ pub struct Greeks {
     pub rho: f64,
 }
 
+/// Second-order and cross Greeks.
+#[derive(Clone, Copy, Debug, PartialEq)]
+pub struct ExtendedGreeks {
+    pub vanna: f64, // ∂Δ/∂σ
+    pub volga: f64, // ∂²V/∂σ² (vomma)
+    pub charm: f64, // ∂Δ/∂t
+    pub speed: f64, // ∂Γ/∂S
+    pub color: f64, // ∂Γ/∂t
+}
+
 /// Shared contract fields for model-based option analytics.
 #[derive(Clone, Copy, Debug, PartialEq)]
 pub struct OptionContract {

crates/ferro_ta_core/src/options/payoff.rs

@@ -0,0 +1,392 @@
+//! Pure-Rust (no PyO3, no numpy) strategy payoff and value functions.
+//!
+//! NOTE: `crates/ferro_ta_core/src/options/mod.rs` must declare `pub mod payoff;`
+//! for this module to be reachable from the rest of the crate and from the PyO3 bridge.
+
+use super::pricing::black_scholes_price;
+use super::OptionKind;
+
+// ---------------------------------------------------------------------------
+// Internal helpers
+// ---------------------------------------------------------------------------
+
+/// Instrument codes: 0=option, 1=future, 2=stock.
+const INSTRUMENT_OPTION: i64 = 0;
+const INSTRUMENT_FUTURE: i64 = 1;
+const INSTRUMENT_STOCK: i64 = 2;
+
+/// Side sign from encoded value: 1=long (+1.0), -1=short (-1.0).
+#[inline]
+fn side_sign(v: i64) -> f64 {
+    if v == 1 {
+        1.0
+    } else if v == -1 {
+        -1.0
+    } else {
+        f64::NAN
+    }
+}
+
+/// Option kind from encoded value: 1=call, -1=put.
+#[inline]
+fn option_kind(v: i64) -> Option<OptionKind> {
+    match v {
+        1 => Some(OptionKind::Call),
+        -1 => Some(OptionKind::Put),
+        _ => None,
+    }
+}
+
+// ---------------------------------------------------------------------------
+// strategy_payoff_dense
+// ---------------------------------------------------------------------------
+
+/// Aggregate strategy payoff over a spot grid.
+///
+/// Parameters (all slices of length n_legs):
+/// - `instruments`: 0=option, 1=future, 2=stock
+/// - `sides`: 1=long, -1=short
+/// - `option_types`: 1=call, -1=put (ignored for futures/stocks)
+/// - `strikes`: strike for options
+/// - `premiums`: premium for options
+/// - `entry_prices`: entry price for futures/stocks
+/// - `quantities`, `multipliers`: applied to all instruments
+///
+/// Returns a Vec<f64> of length spot_grid.len() with aggregate P&L per spot point.
+#[allow(clippy::too_many_arguments)]
+pub fn strategy_payoff_dense(
+    spot_grid: &[f64],
+    instruments: &[i64],
+    sides: &[i64],
+    option_types: &[i64],
+    strikes: &[f64],
+    premiums: &[f64],
+    entry_prices: &[f64],
+    quantities: &[f64],
+    multipliers: &[f64],
+) -> Vec<f64> {
+    let n_legs = instruments.len();
+    // Validate that all leg slices are the same length; return zeros if not.
+    if sides.len() != n_legs
+        || option_types.len() != n_legs
+        || strikes.len() != n_legs
+        || premiums.len() != n_legs
+        || entry_prices.len() != n_legs
+        || quantities.len() != n_legs
+        || multipliers.len() != n_legs
+    {
+        return vec![0.0; spot_grid.len()];
+    }
+
+    let mut total = vec![0.0_f64; spot_grid.len()];
+
+    for leg_idx in 0..n_legs {
+        let inst = instruments[leg_idx];
+        let sign = side_sign(sides[leg_idx]);
+        if sign.is_nan() {
+            // Invalid side — skip leg (treat as zero contribution).
+            continue;
+        }
+        let leg_scale = sign * quantities[leg_idx] * multipliers[leg_idx];
+
+        match inst {
+            INSTRUMENT_OPTION => {
+                let kind = match option_kind(option_types[leg_idx]) {
+                    Some(k) => k,
+                    None => continue, // Invalid option type — skip.
+                };
+                let k = strikes[leg_idx];
+                let p = premiums[leg_idx];
+                for (i, &s) in spot_grid.iter().enumerate() {
+                    let intrinsic = match kind {
+                        OptionKind::Call => (s - k).max(0.0),
+                        OptionKind::Put => (k - s).max(0.0),
+                    };
+                    total[i] += leg_scale * (intrinsic - p);
+                }
+            }
+            INSTRUMENT_FUTURE | INSTRUMENT_STOCK => {
+                let e = entry_prices[leg_idx];
+                for (i, &s) in spot_grid.iter().enumerate() {
+                    total[i] += leg_scale * (s - e);
+                }
+            }
+            _ => {
+                // Unknown instrument code — skip leg (NaN would propagate; zeros are safer).
+            }
+        }
+    }
+
+    total
+}
+
+// ---------------------------------------------------------------------------
+// strategy_value_dense / strategy_value_grid
+// ---------------------------------------------------------------------------
+
+/// Current BSM value of a strategy at a single spot (pre-expiry).
+///
+/// Unlike `strategy_payoff_dense`, this uses BSM pricing for option legs rather
+/// than intrinsic value.
+///
+/// Parameters: same as `strategy_payoff_dense` plus per-leg BSM inputs:
+/// - `time_to_expiries`: TTE for each option leg (ignored for futures/stocks)
+/// - `volatilities`: vol for each option leg (ignored for futures/stocks)
+/// - `rates`: risk-free rate for each leg
+/// - `carries`: carry/dividend yield for each option leg
+///
+/// Returns a scalar f64 (strategy P&L at the given spot).
+#[allow(clippy::too_many_arguments)]
+pub fn strategy_value_dense(
+    spot: f64,
+    instruments: &[i64],
+    sides: &[i64],
+    option_types: &[i64],
+    strikes: &[f64],
+    premiums: &[f64],
+    entry_prices: &[f64],
+    quantities: &[f64],
+    multipliers: &[f64],
+    time_to_expiries: &[f64],
+    volatilities: &[f64],
+    rates: &[f64],
+    carries: &[f64],
+) -> f64 {
+    let n_legs = instruments.len();
+    // Validate that all leg slices are the same length; return NaN if not.
+    if sides.len() != n_legs
+        || option_types.len() != n_legs
+        || strikes.len() != n_legs
+        || premiums.len() != n_legs
+        || entry_prices.len() != n_legs
+        || quantities.len() != n_legs
+        || multipliers.len() != n_legs
+        || time_to_expiries.len() != n_legs
+        || volatilities.len() != n_legs
+        || rates.len() != n_legs
+        || carries.len() != n_legs
+    {
+        return f64::NAN;
+    }
+
+    let mut total = 0.0_f64;
+
+    for leg_idx in 0..n_legs {
+        let inst = instruments[leg_idx];
+        let sign = side_sign(sides[leg_idx]);
+        if sign.is_nan() {
+            continue;
+        }
+        let leg_scale = sign * quantities[leg_idx] * multipliers[leg_idx];
+
+        match inst {
+            INSTRUMENT_OPTION => {
+                let kind = match option_kind(option_types[leg_idx]) {
+                    Some(k) => k,
+                    None => continue,
+                };
+                let bsm = black_scholes_price(
+                    spot,
+                    strikes[leg_idx],
+                    rates[leg_idx],
+                    carries[leg_idx],
+                    time_to_expiries[leg_idx],
+                    volatilities[leg_idx],
+                    kind,
+                );
+                total += leg_scale * (bsm - premiums[leg_idx]);
+            }
+            INSTRUMENT_FUTURE | INSTRUMENT_STOCK => {
+                total += leg_scale * (spot - entry_prices[leg_idx]);
+            }
+            _ => {}
+        }
+    }
+
+    total
+}
+
+// ---------------------------------------------------------------------------
+// aggregate_greeks_dense
+// ---------------------------------------------------------------------------
+
+/// Aggregate BSM Greeks for a multi-leg strategy at a single spot.
+///
+/// Parameters (all slices of length n_legs):
+/// - `instruments`: 0=option, 1=future, 2=stock
+/// - `sides`: 1=long, -1=short
+/// - `option_types`: 1=call, -1=put (ignored for futures/stocks)
+/// - `strikes`: strike price for option legs
+/// - `volatilities`: implied vol for option legs
+/// - `time_to_expiries`: TTE in years for option legs
+/// - `rates`: risk-free rate for each leg
+/// - `carries`: carry/dividend yield for option legs
+/// - `quantities`, `multipliers`: applied to all instruments
+///
+/// Returns `(delta, gamma, vega, theta, rho)` aggregate across all legs.
+/// Future/stock legs contribute `leg_scale` to delta only (all other Greeks = 0).
+#[allow(clippy::too_many_arguments)]
+pub fn aggregate_greeks_dense(
+    spot: f64,
+    instruments: &[i64],
+    sides: &[i64],
+    option_types: &[i64],
+    strikes: &[f64],
+    volatilities: &[f64],
+    time_to_expiries: &[f64],
+    rates: &[f64],
+    carries: &[f64],
+    quantities: &[f64],
+    multipliers: &[f64],
+) -> (f64, f64, f64, f64, f64) {
+    use super::greeks::model_greeks;
+    use super::{OptionContract, OptionEvaluation, PricingModel};
+
+    let n_legs = instruments.len();
+    if sides.len() != n_legs
+        || option_types.len() != n_legs
+        || strikes.len() != n_legs
+        || volatilities.len() != n_legs
+        || time_to_expiries.len() != n_legs
+        || rates.len() != n_legs
+        || carries.len() != n_legs
+        || quantities.len() != n_legs
+        || multipliers.len() != n_legs
+    {
+        return (f64::NAN, f64::NAN, f64::NAN, f64::NAN, f64::NAN);
+    }
+
+    let mut delta = 0.0_f64;
+    let mut gamma = 0.0_f64;
+    let mut vega = 0.0_f64;
+    let mut theta = 0.0_f64;
+    let mut rho = 0.0_f64;
+
+    for i in 0..n_legs {
+        let sign = side_sign(sides[i]);
+        if sign.is_nan() {
+            continue;
+        }
+        let leg_scale = sign * quantities[i] * multipliers[i];
+
+        match instruments[i] {
+            INSTRUMENT_FUTURE | INSTRUMENT_STOCK => {
+                delta += leg_scale;
+            }
+            INSTRUMENT_OPTION => {
+                let kind = match option_kind(option_types[i]) {
+                    Some(k) => k,
+                    None => continue,
+                };
+                let greeks = model_greeks(OptionEvaluation {
+                    contract: OptionContract {
+                        model: PricingModel::BlackScholes,
+                        underlying: spot,
+                        strike: strikes[i],
+                        rate: rates[i],
+                        carry: carries[i],
+                        time_to_expiry: time_to_expiries[i],
+                        kind,
+                    },
+                    volatility: volatilities[i],
+                });
+                delta += leg_scale * greeks.delta;
+                gamma += leg_scale * greeks.gamma;
+                vega += leg_scale * greeks.vega;
+                theta += leg_scale * greeks.theta;
+                rho += leg_scale * greeks.rho;
+            }
+            _ => {}
+        }
+    }
+
+    (delta, gamma, vega, theta, rho)
+}
+
+/// Evaluate `strategy_value_dense` for each point in `spot_grid`.
+///
+/// Returns a `Vec<f64>` of length `spot_grid.len()`.
+#[allow(clippy::too_many_arguments)]
+pub fn strategy_value_grid(
+    spot_grid: &[f64],
+    instruments: &[i64],
+    sides: &[i64],
+    option_types: &[i64],
+    strikes: &[f64],
+    premiums: &[f64],
+    entry_prices: &[f64],
+    quantities: &[f64],
+    multipliers: &[f64],
+    time_to_expiries: &[f64],
+    volatilities: &[f64],
+    rates: &[f64],
+    carries: &[f64],
+) -> Vec<f64> {
+    spot_grid
+        .iter()
+        .map(|&s| {
+            strategy_value_dense(
+                s,
+                instruments,
+                sides,
+                option_types,
+                strikes,
+                premiums,
+                entry_prices,
+                quantities,
+                multipliers,
+                time_to_expiries,
+                volatilities,
+                rates,
+                carries,
+            )
+        })
+        .collect()
+}
+
+// ---------------------------------------------------------------------------
+// Tests
+// ---------------------------------------------------------------------------
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn payoff_single_call() {
+        let grid = vec![90.0, 100.0, 110.0, 120.0];
+        let out = strategy_payoff_dense(
+            &grid,
+            &[0],
+            &[1],
+            &[1],
+            &[100.0],
+            &[5.0],
+            &[0.0],
+            &[1.0],
+            &[1.0],
+        );
+        assert!(out[0] < 0.0); // below strike, loss = premium
+        assert!((out[0] - (-5.0)).abs() < 1e-10);
+        assert!((out[2] - 5.0).abs() < 1e-10); // at 110, intrinsic=10, net=10-5=5
+    }
+
+    #[test]
+    fn stock_leg_linear() {
+        let grid = vec![90.0, 100.0, 110.0];
+        let out = strategy_payoff_dense(
+            &grid,
+            &[2],
+            &[1],
+            &[0],
+            &[0.0],
+            &[0.0],
+            &[100.0],
+            &[1.0],
+            &[1.0],
+        );
+        assert!((out[0] - (-10.0)).abs() < 1e-10);
+        assert!((out[1] - 0.0).abs() < 1e-10);
+        assert!((out[2] - 10.0).abs() < 1e-10);
+    }
+}

crates/ferro_ta_core/src/options/pricing.rs

@@ -118,6 +118,37 @@ pub fn model_price(input: OptionEvaluation) -> f64 {
     }
 }
 
+/// Put-call parity deviation: `C - P - (S·e^{-q·T} - K·e^{-r·T})`.
+///
+/// Returns 0.0 when no arbitrage exists. A non-zero value indicates the
+/// magnitude of mispricing or data error.
+pub fn put_call_parity_deviation(
+    call_price: f64,
+    put_price: f64,
+    spot: f64,
+    strike: f64,
+    rate: f64,
+    carry: f64,
+    time_to_expiry: f64,
+) -> f64 {
+    if !call_price.is_finite()
+        || !put_price.is_finite()
+        || !spot.is_finite()
+        || !strike.is_finite()
+        || !rate.is_finite()
+        || !carry.is_finite()
+        || !time_to_expiry.is_finite()
+        || spot <= 0.0
+        || strike <= 0.0
+        || time_to_expiry < 0.0
+    {
+        return f64::NAN;
+    }
+    let pv_forward = spot * (-carry * time_to_expiry).exp();
+    let pv_strike = strike * (-rate * time_to_expiry).exp();
+    call_price - put_price - (pv_forward - pv_strike)
+}
+
 /// Lower no-arbitrage bound for the option price.
 pub fn price_lower_bound(contract: OptionContract) -> f64 {
     match contract.model {

crates/ferro_ta_core/src/options/realized_vol.rs

@@ -0,0 +1,445 @@
+//! Historical (realized) volatility estimators and volatility cone.
+
+/// Rolling close-to-close realized volatility.
+///
+/// Returns a `Vec<f64>` of the same length as `close`. The first `window` values
+/// are NaN (we need `window` log-returns, which require `window+1` prices, so the
+/// first valid output sits at index `window`).
+///
+/// Annualization: `sqrt(sum(r²) / window * trading_days)`.
+pub fn close_to_close_vol(close: &[f64], window: usize, trading_days: f64) -> Vec<f64> {
+    let n = close.len();
+    let mut out = vec![f64::NAN; n];
+    if window == 0 || n <= window {
+        return out;
+    }
+
+    // Precompute log-returns; returns[i] = ln(close[i+1] / close[i])
+    let mut returns = vec![f64::NAN; n - 1];
+    for i in 0..(n - 1) {
+        if close[i] > 0.0 && close[i + 1] > 0.0 {
+            returns[i] = (close[i + 1] / close[i]).ln();
+        }
+    }
+
+    // Rolling sum of squared returns over `window` bars.
+    // The output at position `end` (in the original close array) uses
+    // returns[end-window .. end-1], i.e. `window` returns.
+    for end in window..n {
+        let slice = &returns[(end - window)..end];
+        let sum_sq: f64 = slice.iter().map(|&r| r * r).sum();
+        let var = sum_sq / window as f64 * trading_days;
+        out[end] = if var >= 0.0 { var.sqrt() } else { f64::NAN };
+    }
+    out
+}
+
+/// Rolling Parkinson high-low realized volatility estimator.
+///
+/// Returns a `Vec<f64>` of the same length as `high`. The first `window-1` values
+/// are NaN.
+#[allow(clippy::needless_range_loop)]
+pub fn parkinson_vol(high: &[f64], low: &[f64], window: usize, trading_days: f64) -> Vec<f64> {
+    let n = high.len();
+    let mut out = vec![f64::NAN; n];
+    if window == 0 || n < window || low.len() != n {
+        return out;
+    }
+
+    let factor = 1.0 / (4.0 * 2_f64.ln());
+
+    for end in (window - 1)..n {
+        let start = end + 1 - window;
+        let mut sum_sq = 0.0;
+        let mut valid = true;
+        for i in start..=end {
+            if high[i] <= 0.0 || low[i] <= 0.0 || !high[i].is_finite() || !low[i].is_finite() {
+                valid = false;
+                break;
+            }
+            let u = (high[i] / low[i]).ln();
+            sum_sq += u * u;
+        }
+        if valid {
+            let var = factor * sum_sq / window as f64 * trading_days;
+            out[end] = if var >= 0.0 { var.sqrt() } else { f64::NAN };
+        }
+    }
+    out
+}
+
+/// Rolling Garman-Klass OHLC realized volatility estimator.
+///
+/// Returns a `Vec<f64>` of the same length as the inputs. The first `window-1`
+/// values are NaN. All four slices must have the same length.
+pub fn garman_klass_vol(
+    open: &[f64],
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    window: usize,
+    trading_days: f64,
+) -> Vec<f64> {
+    let n = open.len();
+    let mut out = vec![f64::NAN; n];
+    if window == 0 || n < window || high.len() != n || low.len() != n || close.len() != n {
+        return out;
+    }
+
+    let ln2 = 2_f64.ln();
+
+    // Precompute per-bar GK contributions.
+    let mut gk = vec![f64::NAN; n];
+    for i in 0..n {
+        let o = open[i];
+        let h = high[i];
+        let l = low[i];
+        let c = close[i];
+        if o > 0.0
+            && h > 0.0
+            && l > 0.0
+            && c > 0.0
+            && o.is_finite()
+            && h.is_finite()
+            && l.is_finite()
+            && c.is_finite()
+        {
+            let u = (h / o).ln();
+            let d = (l / o).ln();
+            let ci = (c / o).ln();
+            gk[i] = 0.5 * (u - d).powi(2) - (2.0 * ln2 - 1.0) * ci * ci;
+        }
+    }
+
+    for end in (window - 1)..n {
+        let start = end + 1 - window;
+        let slice = &gk[start..=end];
+        if slice.iter().all(|v| v.is_finite()) {
+            let sum: f64 = slice.iter().sum();
+            let var = sum / window as f64 * trading_days;
+            out[end] = if var >= 0.0 { var.sqrt() } else { f64::NAN };
+        }
+    }
+    out
+}
+
+/// Compute the Rogers-Satchell per-bar variance contribution.
+fn rs_bar(open: f64, high: f64, low: f64, close: f64) -> f64 {
+    let u = (high / close).ln();
+    let d = (low / close).ln();
+    let uo = (high / open).ln();
+    let do_ = (low / open).ln();
+    u * uo + d * do_
+}
+
+/// Rolling Rogers-Satchell OHLC realized volatility estimator.
+///
+/// Returns a `Vec<f64>` of the same length as the inputs. The first `window-1`
+/// values are NaN. All four slices must have the same length.
+pub fn rogers_satchell_vol(
+    open: &[f64],
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    window: usize,
+    trading_days: f64,
+) -> Vec<f64> {
+    let n = open.len();
+    let mut out = vec![f64::NAN; n];
+    if window == 0 || n < window || high.len() != n || low.len() != n || close.len() != n {
+        return out;
+    }
+
+    // Precompute per-bar RS contributions.
+    let mut rs = vec![f64::NAN; n];
+    for i in 0..n {
+        let o = open[i];
+        let h = high[i];
+        let l = low[i];
+        let c = close[i];
+        if o > 0.0
+            && h > 0.0
+            && l > 0.0
+            && c > 0.0
+            && o.is_finite()
+            && h.is_finite()
+            && l.is_finite()
+            && c.is_finite()
+        {
+            rs[i] = rs_bar(o, h, l, c);
+        }
+    }
+
+    for end in (window - 1)..n {
+        let start = end + 1 - window;
+        let slice = &rs[start..=end];
+        if slice.iter().all(|v| v.is_finite()) {
+            let sum: f64 = slice.iter().sum();
+            let var = sum / window as f64 * trading_days;
+            out[end] = if var >= 0.0 { var.sqrt() } else { f64::NAN };
+        }
+    }
+    out
+}
+
+/// Rolling Yang-Zhang OHLC realized volatility estimator.
+///
+/// Handles overnight gaps. Returns a `Vec<f64>` of the same length as the inputs.
+/// The first `window` values are NaN (we need `window` bars plus the prior close
+/// for overnight returns, so valid output starts at index `window`).
+/// All four slices must have the same length.
+pub fn yang_zhang_vol(
+    open: &[f64],
+    high: &[f64],
+    low: &[f64],
+    close: &[f64],
+    window: usize,
+    trading_days: f64,
+) -> Vec<f64> {
+    let n = open.len();
+    let mut out = vec![f64::NAN; n];
+    if window == 0 || n <= window || high.len() != n || low.len() != n || close.len() != n {
+        return out;
+    }
+
+    let k = 0.34 / (1.34 + (window as f64 + 1.0) / (window as f64 - 1.0).max(1e-10));
+
+    // Precompute per-bar components; index 0 has no overnight return.
+    // overnight[i] = ln(O_i / C_{i-1}), valid for i >= 1
+    // openclose[i] = ln(C_i / O_i)
+    // rs[i] = Rogers-Satchell for bar i
+    let mut overnight = vec![f64::NAN; n];
+    let mut openclose = vec![f64::NAN; n];
+    let mut rs = vec![f64::NAN; n];
+
+    for i in 0..n {
+        let o = open[i];
+        let h = high[i];
+        let l = low[i];
+        let c = close[i];
+        if o > 0.0
+            && h > 0.0
+            && l > 0.0
+            && c > 0.0
+            && o.is_finite()
+            && h.is_finite()
+            && l.is_finite()
+            && c.is_finite()
+        {
+            openclose[i] = (c / o).ln();
+            rs[i] = rs_bar(o, h, l, c);
+
+            if i > 0 {
+                let prev_c = close[i - 1];
+                if prev_c > 0.0 && prev_c.is_finite() {
+                    overnight[i] = (o / prev_c).ln();
+                }
+            }
+        }
+    }
+
+    // Valid windows start at index `window` (using bars [end-window+1 .. end],
+    // all of which have valid overnight returns since they start at index >= 1).
+    for end in window..n {
+        let start = end + 1 - window; // start >= 1 because end >= window
+
+        let o_slice = &overnight[start..=end];
+        let c_slice = &openclose[start..=end];
+        let r_slice = &rs[start..=end];
+
+        if !o_slice.iter().all(|v| v.is_finite())
+            || !c_slice.iter().all(|v| v.is_finite())
+            || !r_slice.iter().all(|v| v.is_finite())
+        {
+            continue;
+        }
+
+        let w = window as f64;
+
+        let o_sum: f64 = o_slice.iter().sum();
+        let o_sum_sq: f64 = o_slice.iter().map(|&x| x * x).sum();
+        let overnight_var = o_sum_sq / (w - 1.0) - (o_sum / w).powi(2) * w / (w - 1.0);
+
+        let c_sum: f64 = c_slice.iter().sum();
+        let c_sum_sq: f64 = c_slice.iter().map(|&x| x * x).sum();
+        let openclose_var = c_sum_sq / (w - 1.0) - (c_sum / w).powi(2) * w / (w - 1.0);
+
+        let rs_sum: f64 = r_slice.iter().sum();
+        let rs_var = rs_sum / w;
+
+        let yz_var = overnight_var + k * openclose_var + (1.0 - k) * rs_var;
+        let annualized = yz_var * trading_days;
+        out[end] = if annualized >= 0.0 {
+            annualized.sqrt()
+        } else {
+            f64::NAN
+        };
+    }
+    out
+}
+
+/// Summary statistics of realized vol distribution for one window length.
+#[derive(Clone, Copy, Debug)]
+pub struct VolConeSlice {
+    pub window: usize,
+    pub min: f64,
+    pub p25: f64,
+    pub median: f64,
+    pub p75: f64,
+    pub max: f64,
+}
+
+/// Compute a percentile via linear interpolation on a sorted slice.
+///
+/// `sorted` must be non-empty and already sorted ascending.
+fn percentile_sorted(sorted: &[f64], p: f64) -> f64 {
+    let n = sorted.len();
+    if n == 1 {
+        return sorted[0];
+    }
+    let idx = (n - 1) as f64 * p;
+    let lo = idx.floor() as usize;
+    let hi = idx.ceil() as usize;
+    let frac = idx - lo as f64;
+    sorted[lo] + frac * (sorted[hi] - sorted[lo])
+}
+
+/// Compute vol cone: distribution of realized vols across multiple window lengths.
+///
+/// For each window in `windows`, the close-to-close rolling vol is computed,
+/// NaN values are filtered out, and the distribution statistics (min, p25,
+/// median, p75, max) are derived via linear interpolation.
+pub fn vol_cone(close: &[f64], windows: &[usize], trading_days: f64) -> Vec<VolConeSlice> {
+    windows
+        .iter()
+        .map(|&w| {
+            let vols = close_to_close_vol(close, w, trading_days);
+            let mut valid: Vec<f64> = vols.into_iter().filter(|v| v.is_finite()).collect();
+            valid.sort_by(|a, b| a.partial_cmp(b).unwrap());
+
+            if valid.is_empty() {
+                return VolConeSlice {
+                    window: w,
+                    min: f64::NAN,
+                    p25: f64::NAN,
+                    median: f64::NAN,
+                    p75: f64::NAN,
+                    max: f64::NAN,
+                };
+            }
+
+            VolConeSlice {
+                window: w,
+                min: valid[0],
+                p25: percentile_sorted(&valid, 0.25),
+                median: percentile_sorted(&valid, 0.5),
+                p75: percentile_sorted(&valid, 0.75),
+                max: *valid.last().unwrap(),
+            }
+        })
+        .collect()
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn fake_prices(n: usize) -> Vec<f64> {
+        // simple synthetic price series
+        let mut prices = vec![100.0_f64; n];
+        for i in 1..n {
+            prices[i] = prices[i - 1] * (1.0 + 0.01 * (i as f64 % 7_f64 - 3.0) * 0.01);
+        }
+        prices
+    }
+
+    #[test]
+    fn close_to_close_returns_nans_for_warmup() {
+        let close = fake_prices(100);
+        let result = close_to_close_vol(&close, 20, 252.0);
+        assert_eq!(result.len(), 100);
+        // first 20 values should be NaN (window-1 of returns warmup + 1 for diff)
+        for i in 0..20 {
+            assert!(result[i].is_nan(), "result[{i}] should be NaN");
+        }
+        assert!(result[20].is_finite());
+    }
+
+    #[test]
+    fn parkinson_vol_is_positive() {
+        let close = fake_prices(100);
+        let high: Vec<f64> = close.iter().map(|&c| c * 1.01).collect();
+        let low: Vec<f64> = close.iter().map(|&c| c * 0.99).collect();
+        let result = parkinson_vol(&high, &low, 20, 252.0);
+        for v in result.iter().skip(19) {
+            assert!(v.is_finite() && *v >= 0.0);
+        }
+    }
+
+    #[test]
+    fn vol_cone_is_ordered() {
+        let close = fake_prices(300);
+        let cones = vol_cone(&close, &[20, 60], 252.0);
+        assert_eq!(cones.len(), 2);
+        for cone in &cones {
+            assert!(cone.min <= cone.p25);
+            assert!(cone.p25 <= cone.median);
+            assert!(cone.median <= cone.p75);
+            assert!(cone.p75 <= cone.max);
+        }
+    }
+
+    #[test]
+    fn garman_klass_returns_nans_for_warmup() {
+        let close = fake_prices(50);
+        let high: Vec<f64> = close.iter().map(|&c| c * 1.01).collect();
+        let low: Vec<f64> = close.iter().map(|&c| c * 0.99).collect();
+        let result = garman_klass_vol(&close, &high, &low, &close, 10, 252.0);
+        assert_eq!(result.len(), 50);
+        for i in 0..9 {
+            assert!(result[i].is_nan(), "result[{i}] should be NaN");
+        }
+        assert!(result[9].is_finite());
+    }
+
+    #[test]
+    fn rogers_satchell_returns_nans_for_warmup() {
+        let close = fake_prices(50);
+        let high: Vec<f64> = close.iter().map(|&c| c * 1.01).collect();
+        let low: Vec<f64> = close.iter().map(|&c| c * 0.99).collect();
+        let result = rogers_satchell_vol(&close, &high, &low, &close, 10, 252.0);
+        assert_eq!(result.len(), 50);
+        for i in 0..9 {
+            assert!(result[i].is_nan(), "result[{i}] should be NaN");
+        }
+        assert!(result[9].is_finite());
+    }
+
+    #[test]
+    fn yang_zhang_returns_nans_for_warmup() {
+        let close = fake_prices(50);
+        let high: Vec<f64> = close.iter().map(|&c| c * 1.01).collect();
+        let low: Vec<f64> = close.iter().map(|&c| c * 0.99).collect();
+        let result = yang_zhang_vol(&close, &high, &low, &close, 10, 252.0);
+        assert_eq!(result.len(), 50);
+        for i in 0..10 {
+            assert!(result[i].is_nan(), "result[{i}] should be NaN");
+        }
+        assert!(result[10].is_finite());
+    }
+
+    #[test]
+    fn mismatched_lengths_return_all_nan() {
+        let a = vec![100.0_f64; 20];
+        let b = vec![101.0_f64; 15]; // wrong length
+        let result = parkinson_vol(&a, &b, 5, 252.0);
+        assert!(result.iter().all(|v| v.is_nan()));
+    }
+
+    #[test]
+    fn window_larger_than_data_returns_all_nan() {
+        let close = fake_prices(10);
+        let result = close_to_close_vol(&close, 20, 252.0);
+        assert!(result.iter().all(|v| v.is_nan()));
+    }
+}

crates/ferro_ta_core/src/options/surface.rs

@@ -202,6 +202,35 @@ pub fn term_structure_slope(tenors: &[f64], atm_ivs: &[f64]) -> f64 {
     regression_slope(tenors, atm_ivs)
 }
 
+/// Expected ±1σ move over `days_to_expiry` calendar days.
+///
+/// Returns `(lower_move, upper_move)` as absolute changes from `spot`.
+/// Example: if spot=100 and upper_move=5.0 then the 1σ upper bound is 105.
+///
+/// Uses the log-normal approximation: `spot × e^{±σ√(days/trading_days)} − spot`.
+pub fn expected_move(
+    spot: f64,
+    iv: f64,
+    days_to_expiry: f64,
+    trading_days_per_year: f64,
+) -> (f64, f64) {
+    if !spot.is_finite()
+        || !iv.is_finite()
+        || !days_to_expiry.is_finite()
+        || !trading_days_per_year.is_finite()
+        || spot <= 0.0
+        || iv < 0.0
+        || days_to_expiry < 0.0
+        || trading_days_per_year <= 0.0
+    {
+        return (f64::NAN, f64::NAN);
+    }
+    let sigma_sqrt_t = iv * (days_to_expiry / trading_days_per_year).sqrt();
+    let upper = spot * sigma_sqrt_t.exp() - spot;
+    let lower = spot * (-sigma_sqrt_t).exp() - spot;
+    (lower, upper)
+}
+
 #[cfg(test)]
 mod tests {
     use super::{atm_iv, smile_metrics, term_structure_slope};