option-pricing 0.1.4

//! # Implied Volatility
//! This module provides the implied volatility calculation for a given option price.
//!
//! ## Example
//! ```
//! use crate::black_scholes::BlackScholes;
//! use crate::black_scholes::Input as BSInput;
//! use crate::implied_vol::Input as IVInput;
//! use crate::implied_vol::{ImpliedVol, Method};
//!
//! let input = BSInput {
//!     is_call: true,
//!     spot: 135.6,
//!     strike: 100.0,
//!     mat: 3.2,
//!     vol: 0.25,
//!     rate: 0.03,
//!     div: 0.01,
//! };
//!
//! let iv_input = IVInput {
//!     price: bs_call.output.unwrap().price,
//!     spot: input.spot,
//!     strike: input.strike,
//!     mat: input.mat,
//!     rate: input.rate,
//!     div: input.div,
//!     iter: 10,
//!     prec: 1e-5,
//! };
//! let iv = ImpliedVol::new(iv_input, Method::Halley);
//! let iv_vol = iv.output.unwrap().vol;
//!
//! let bs_call = BlackScholes::new(input.clone());
//! let bs_vol = bs_call.input.vol;
//!
//! println!("bs_vol: {}, iv_vol: {}", bs_vol, iv_vol);
//!
//! let epsilon = 1e-5;
//! assert!((bs_vol - iv_vol).abs() < epsilon);
//! ```
//!

use serde::{Deserialize, Serialize};
use thiserror::Error;

use crate::black_scholes;

/// # Implied Volatility Scholes wrapper
///
/// Contains the input and output for the implied volatility calculation.
pub struct ImpliedVol {
    /// input parameters
    pub input: Input,
    /// output result
    pub output: Result<Output, ImpliedVolError>,
}

/// # Implied volatility input
///
/// This is the input to the implied volatility function reversing the Black Scholes model.  
///
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Input {
    /// call option price observed in the market
    pub price: f32,
    /// underlying spot, in currency
    pub spot: f32,
    /// option strike, in currency
    pub strike: f32,
    /// option maturity, in years
    pub mat: f32,
    /// interest rate, annualized, e.g. 0.05=5%
    pub rate: f32,
    /// underlying dividend yield, annualized, e.g. 0.02=2%
    pub div: f32,
    /// maximum number of iterations - optional
    pub iter: u32,
    /// target precision - optional
    pub prec: f32,
}

/// # Implied volatility output
///
/// This is the output of the implied volatility function reversing the Black Scholes model.
///
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Output {
    /// volatility calculated
    pub vol: f32,
    /// iterations run
    pub iter: u32,
    /// precision reached
    pub prec: f32,
}

/// # Implied Vol Input Error
///
/// This covers the various ways implied volatility input params can be invalid.
#[derive(Error, Debug, Serialize, Deserialize)]
pub enum ImpliedVolError {
    /// Negative spot
    #[error("negative spot: {0} - must be positive")]
    NegativeSpot(f32),
    /// Negative strike
    #[error("negative strike: {0} - must be positive")]
    NegativeStrike(f32),
    /// Negative maturity
    #[error("negative maturity: {0} - must be positive")]
    NegativeMat(f32),
    /// Out of bound price
    #[error("out of bound price: {0} - must be between max(0, S-PV(K)) and spot")]
    OutOfBoundPrice(f32),
}

/// # Convergence Method
///
/// This enum defines the method used to search for the implied volatility.  
#[derive(Debug, Clone, Serialize, Deserialize)]
pub enum Method {
    /// [Newton-Raphson](https://en.wikipedia.org/wiki/Newton%27s_method) method
    Newton,
    /// [Halley's](https://en.wikipedia.org/wiki/Halley%27s_method) method
    Halley,
}

impl ImpliedVol {
    pub fn new(input: Input, method: Method) -> ImpliedVol {
        let output = find_vol(&input, method);
        ImpliedVol { input, output }
    }
}

fn find_vol(input: &Input, method: Method) -> Result<Output, ImpliedVolError> {
    if input.spot < 0.0 {
        return Err(ImpliedVolError::NegativeSpot(input.spot));
    }
    if input.strike < 0.0 {
        return Err(ImpliedVolError::NegativeStrike(input.strike));
    }
    if input.mat < 0.0 {
        return Err(ImpliedVolError::NegativeMat(input.mat));
    }

    let pv_r = (-input.rate * input.mat).exp();
    let pv_q = (-input.div * input.mat).exp();
    let min_price = (input.spot * pv_q - input.strike * pv_r).max(0.0);
    let max_price = input.spot;

    if input.price > max_price || input.price < min_price {
        return Err(ImpliedVolError::OutOfBoundPrice(input.price));
    }

    let output = match method {
        Method::Newton => newton(input),
        Method::Halley => halley(input),
    };

    Ok(output)
}

const VOL_START_MIN: f32 = 1e-10;

fn newton(input: &Input) -> Output {
    let pv_r = (-input.rate * input.mat).exp();
    let pv_q = (-input.div * input.mat).exp();
    let moneyness = (input.spot * pv_q) / (input.strike * pv_r);
    let vol_start = ((2.0 * moneyness.ln().abs() / input.mat).sqrt()).max(VOL_START_MIN);

    let mut vol = vol_start;
    let mut iter = 0;
    let prec = input.prec;

    while iter < input.iter {
        let calc = black_scholes::call_price_vega(
            input.spot,
            input.strike,
            input.rate,
            input.mat,
            vol,
            input.div,
        );

        let diff = calc.price - input.price;
        if diff.abs() < input.prec {
            break;
        }

        // newton method
        let f = diff;
        let f_prime = calc.vega;
        let step = -f / f_prime;
        vol += step;

        iter += 1;
    }

    Output { vol, iter, prec }
}

fn halley(input: &Input) -> Output {
    let pv_r = (-input.rate * input.mat).exp();
    let pv_q = (-input.div * input.mat).exp();
    let moneyness = (input.spot * pv_q) / (input.strike * pv_r);
    let vol_start = ((2.0 * moneyness.ln().abs() / input.mat).sqrt()).max(VOL_START_MIN);

    let mut vol = vol_start;
    let mut iter = 0;
    let mut prec = 0.0;

    while iter < input.iter {
        let calc = black_scholes::call_price_vega_voma(
            input.spot,
            input.strike,
            input.rate,
            input.mat,
            vol,
            input.div,
        );

        let diff = calc.price - input.price;
        if diff.abs() < input.prec {
            break;
        }

        // halley method
        let f = diff;
        let f_prime = calc.vega;
        let f_second = calc.voma;
        let step = -(2.0 * f * f_prime) / (2.0 * f_prime.powi(2) - f * f_second);
        vol += step;

        iter += 1;
        prec = diff;
    }

    Output { vol, iter, prec }
}

#[cfg(test)]
mod tests {

    use crate::black_scholes::BlackScholes;
    use crate::black_scholes::Input as BSInput;
    use crate::implied_vol::Input as IVInput;
    use crate::implied_vol::{ImpliedVol, Method};

    #[test]
    fn test_newton() {
        let test_data: Vec<BSInput> = vec![
            BSInput {
                is_call: true,
                spot: 135.6,
                strike: 100.0,
                mat: 3.2,
                vol: 0.25,
                rate: 0.03,
                div: 0.01,
            },
            BSInput {
                is_call: true,
                spot: 60.6,
                strike: 100.0,
                mat: 3.2,
                vol: 0.25,
                rate: 0.03,
                div: 0.01,
            },
            BSInput {
                is_call: true,
                spot: 100.0,
                strike: 100.0,
                mat: 3.2,
                vol: 0.35,
                rate: 0.03,
                div: 0.00,
            },
            BSInput {
                is_call: true,
                spot: 100.0,
                strike: 100.0,
                mat: 1.0,
                vol: 0.10,
                rate: 0.00,
                div: 0.00,
            },
        ];

        test_data.iter().for_each(|input| {
            let bs_call = BlackScholes::new(input.clone());
            let bs_vol = bs_call.input.vol;

            let iv_input = IVInput {
                price: bs_call.output.unwrap().price,
                spot: input.spot,
                strike: input.strike,
                mat: input.mat,
                rate: input.rate,
                div: input.div,
                iter: 10,
                prec: 1e-6,
            };

            let iv = ImpliedVol::new(iv_input, Method::Newton);
            let iv_vol = iv.output.unwrap().vol;

            println!("bs_vol: {}, iv_vol: {}", bs_vol, iv_vol);

            let epsilon = 1e-5;
            assert!((bs_vol - iv_vol).abs() < epsilon);
        })
    }

    #[test]
    fn test_halley() {
        let test_data: Vec<BSInput> = vec![
            BSInput {
                is_call: true,
                spot: 135.6,
                strike: 100.0,
                mat: 3.2,
                vol: 0.25,
                rate: 0.03,
                div: 0.01,
            },
            BSInput {
                is_call: true,
                spot: 60.6,
                strike: 100.0,
                mat: 3.2,
                vol: 0.25,
                rate: 0.03,
                div: 0.01,
            },
            BSInput {
                is_call: true,
                spot: 100.0,
                strike: 100.0,
                mat: 3.2,
                vol: 0.35,
                rate: 0.03,
                div: 0.00,
            },
            BSInput {
                is_call: true,
                spot: 100.0,
                strike: 100.0,
                mat: 1.0,
                vol: 0.10,
                rate: 0.00,
                div: 0.00,
            },
        ];

        test_data.iter().for_each(|input| {
            let bs_call = BlackScholes::new(input.clone());
            let bs_vol = bs_call.input.vol;

            let iv_input = IVInput {
                price: bs_call.output.unwrap().price,
                spot: input.spot,
                strike: input.strike,
                mat: input.mat,
                rate: input.rate,
                div: input.div,
                iter: 10,
                prec: 1e-6,
            };

            let iv = ImpliedVol::new(iv_input, Method::Halley);
            let iv_vol = iv.output.unwrap().vol;

            println!("bs_vol: {}, iv_vol: {}", bs_vol, iv_vol);

            let epsilon = 1e-5;
            assert!((bs_vol - iv_vol).abs() < epsilon);
        })
    }
}