rspp 0.1.7

#![allow(non_camel_case_types)]
#![allow(non_snake_case)]
use serde::{Deserialize, Serialize};
use statrs::distribution::{Continuous, ContinuousCDF, Normal};
use statrs::statistics::Distribution;
use std::f64::consts::E;
use std::f64::consts::PI;
use std::f64::NAN;
///  计算结果都已近和此网站校验过  https://www.barchart.com/options/options-calculator
#[derive(Serialize, Deserialize, Debug, Clone)]
pub struct Option {
    pub is_call: bool,

    pub S: f64, // base asset price
    pub K: f64, // exacuate price
    pub T: f64,
    pub divdend: f64, //divdent rate
    pub r: f64,       //risk free rate
}

pub fn d1(op: &Option, vol: f64) -> f64 {
    let num = (op.S / op.K).log(E) + (op.r - op.divdend + vol.powf(2.0) / 2.0) * op.T;
    let denomin = vol * (op.T.sqrt());
    num / denomin
}
fn d2(op: &Option, vol: f64) -> f64 {
    d1(op, vol) - vol * (op.T.sqrt())
}
fn discount_r(op: &Option) -> f64 {
    (0.0 - op.r * op.T).exp()
}
fn discount_divd(op: &Option) -> f64 {
    (0.0 - op.divdend * op.T).exp()
}
/// # Example
/// ```
/// use rspp::option_price::bsm::Option;
///    let op = Option {
///       is_call: false,
///      S: 217.36, // base asset price
///      K: 219.0,  // exacuate price
///    T: 13.0 / 365.0,
///   divdend: -1.44 / 100.0, //divdent rate
///  r: 3.62 / 100.0,        //risk free rate
///   };

///   let vol = 35.23 / 100.0;
/// let mp = 4.6;
/// let price = op.bsm_price(vol);
/// let theta = op.bsm_theta(vol);
///let vega = op.bsm_vega(vol);
///let rho = op.bsm_rho(vol);
///let gamma = op.bsm_gamma(vol);
///let delta = op.bsm_delta(vol);
///let imply = op.bsm_implied_vol(mp);
/// dbg!("bsm==", imply, price, delta, gamma, theta, vega, rho);
/// ```

impl Option {
    pub fn bsm_call(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let num1 = self.S * discount_divd(&self) * norm.cdf(d1(&self, vol));
        let num2 = self.K * discount_r(&self) * norm.cdf(d2(&self, vol));
        num1 - num2
    }
    pub fn bsm_put(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let num1 = self.S * discount_divd(&self) * norm.cdf(0.0 - d1(self, vol));
        let num2 = self.K * discount_r(&self) * norm.cdf(0.0 - d2(self, vol));
        num2 - num1
    }

    pub fn bsm_price(&self, vol: f64) -> f64 {
        if self.is_call {
            return self.bsm_call(vol);
        }
        self.bsm_put(vol)
    }
    /// 波动率对期权价格
    pub fn bsm_vega(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let x = d1(&self, vol);
        self.S * self.T.sqrt() * norm.pdf(x) * discount_divd(&self) / 100.0
    }
    pub fn bsm_implied_vol(&self, mp: f64) -> f64 {
        let precise = 1e-3;
        let mut l = Vec::new();
        let mut start = 0.5;
        for n in 0..9000 {
            let dif = mp - self.bsm_price(start);
            if dif.abs() < precise {
                return start;
            }
            l.push((dif.abs(), start));
            let mut veg = self.bsm_vega(start);
            if veg == 0. {
                veg = 0.01;
            }
            start = start + precise * 2. * dif / veg;
        }
        l.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap());
        l[0].1
    }

    ///基础价格对期权价格求导
    pub fn bsm_delta(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let nd = norm.cdf(d1(&self, vol));
        if self.is_call {
            return nd * discount_divd(&self);
        }
        (nd - 1.0) * discount_divd(&self)
    }

    /// 基础价格对期权价格二阶导
    pub fn bsm_gamma(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let d = d1(&self, vol);
        norm.pdf(d) * discount_divd(&self) / (self.S * vol * self.T.sqrt())
    }
    ///期限对期权价格         到期时间不是时间减法,所以要加-
    pub fn bsm_theta(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let first = (0.0 - self.S) * norm.pdf(d1(&self, vol)) * vol * discount_divd(&self)
            / (self.T.sqrt() * 2.0);
        if self.is_call {
            let second = self.divdend * self.S * norm.cdf(d1(&self, vol)) * discount_divd(&self);
            let third = (0.0 - self.r) * self.K * discount_r(&self) * norm.cdf(d2(&self, vol));
            return (first + second + third) / 365.0;
        }
        let second =
            -1.0 * self.divdend * self.S * norm.cdf(-1.0 * d1(&self, vol)) * discount_divd(&self);
        let third = self.r * self.K * discount_r(&self) * norm.cdf(-1.0 * d2(&self, vol));
        (first + second + third) / 365.0
    }

    ///期限对利率价格
    pub fn bsm_rho(&self, vol: f64) -> f64 {
        let norm = Normal::standard();
        let n = self.K * self.T * discount_r(&self) / 100.0;
        let x = if self.is_call {
            norm.cdf(d2(&self, vol))
        } else {
            -1.0 * norm.cdf(-1.0 * d2(&self, vol))
        };
        n * x
    }
}
#[test]
#[docify::export]
fn test_check_validate() {
    // let op = Option {
    //     is_call: true,
    //     S: 220.11, // base asset price
    //     K: 220.0,  // exacuate price
    //     T: 10.0 / 365.0,
    //     divdend: 0.44 / 100.0, //divdent rate
    //     r: 4.07 / 100.0,       //risk free rate
    // };

    // let vol = 28.9 / 100.0;
    // let mp = 4.4;
    // let price = op.bsm_price(vol);
    // let theta = op.bsm_theta(vol);
    // let vega = op.bsm_vega(vol);
    // let rho = op.bsm_rho(vol);
    // let gamma = op.bsm_gamma(vol);
    // let delta = op.bsm_delta(vol);
    // let imply = op.bsm_implied_vol(mp);
    // dbg!("bsm==", imply, price, delta, gamma, theta, vega, rho);
    let op2 = Option {
        is_call: false,
        S: 73110.,  // base asset price
        K: 75000.0, // exacuate price
        T: 0.12,
        divdend: 0., //divdent rate
        r: 0.013,    //risk free rate
    };
    let mp = 2600.;
    let iv = op2.bsm_implied_vol(mp);
    dbg!(iv);
}