option-pricing 0.1.4

//! # Black Scholes pricing
//! This module exposes functions to calculate the price of European call and put options using the [Black Scholes model](https://en.wikipedia.org/wiki/Black%E2%80%93Scholes_model).
//!
//!
//! ## Example
//! ```
//! use crate::black_scholes::{BlackScholes, Input};
//!
//! // define params
//! let spot = 105.6;
//! let strike = 100.0;
//! let mat = 3.2;
//! let vol = 0.25;
//! let rate = 0.03;
//! let div = 0.01;
//!
//! // define inputs
//! let input_call =  input_call: Input {
//!     is_call: true,
//!     spot,
//!     strike,
//!     mat,
//!     vol,
//!     rate,
//!     div,
//! };
//! let input_put =  input_call: Input {
//!     is_call: false,
//!     spot,
//!     strike,
//!     mat,
//!     vol,
//!     rate,
//!     div,
//! };
//!
//! // pricing
//! let bs_call = BlackScholes::new(test.input_call.clone());
//! let bs_put = BlackScholes::new(test.input_put.clone());
//!
//! // check no error
//! assert!(bs_call.output.is_ok());
//! assert!(bs_put.output.is_ok());
//!
//! // extract output
//! let out_call = bs_call.output.unwrap();
//! let out_put = bs_put.output.unwrap();
//!
//! let r = test.input_call.rate;
//! let q = test.input_call.div;
//! let mat = test.input_call.mat;
//! let call = out_call.price;
//! let put = out_put.price;
//! let pv = out_call.pv;
//! let fwd = test.input_call.spot * ((r - q) * mat).exp();
//! let strike = test.input_call.strike;
//!
//! // print results
//! println!("call: {:?}", out_call);
//! println!("put: {:?}", out_put);
//!
//! // check put call parity
//! let epsilon = 1e-5;
//! let parity = (call - put) - pv * (fwd - strike);
//! assert!(parity.abs() < epsilon);
//!

//!```
//!

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

use crate::math;

/// # Black Scholes wrapper
///
/// Contains the input and output of the Black Scholes pricing model.
///
pub struct BlackScholes {
    /// input params
    pub input: Input,
    /// output result
    pub output: Result<Output, PriceError>,
}

/// ## Pricing input
///
/// This is the input to the Black Scholes pricing model.  
/// Valid for both call and put options.  
///
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Input {
    /// `true` for call, `false` for put
    pub is_call: bool,
    /// underlying spot, in currency
    pub spot: f32,
    /// option strike, in currency
    pub strike: f32,
    /// option maturity, in years
    pub mat: f32,
    /// underlying volatility, annualized, e.g. 0.2=20%
    pub vol: f32,
    /// interest rate, annualized, e.g. 0.05=5%
    pub rate: f32,
    /// underlying dividend yield, annualized, e.g. 0.02=2%
    pub div: f32,
}

/// ## Pricing output - full
///
/// This is the output to the Black Scholes calculator.  
/// Valid for both call and put options.  
///
/// **WARNING**: fields cdf_d1 and cdf_d2 do not represent the same value for call and put options.  
///
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct Output {
    /// `true` for call, `false` for put
    pub is_call: bool,
    /// $d_1=\frac{\log(S/K) + (r-q +\sigma^2/2)T}{\sigma \sqrt{T}}$
    pub d1: f32,
    /// $d_2=\frac{\log(S/K) + (r-q -\sigma^2/2)T}{\sigma \sqrt{T}}=d_1 - \sigma \sqrt{T - t}$
    pub d2: f32,
    /// [cumulative distribution function](https://en.wikipedia.org/wiki/Cumulative_distribution_function) for d1  
    /// call: $N(d_1)$  
    /// put:  $N(-d_1)$  
    pub cdf_d1: f32,
    /// [cumulative distribution function](https://en.wikipedia.org/wiki/Cumulative_distribution_function) for d2  
    /// call: $N(d_2)$  
    /// put:  $N(-d_2)$  
    pub cdf_d2: f32,
    /// [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) for d1
    pub pdf_d1: f32,
    /// [probability density function](https://en.wikipedia.org/wiki/Probability_density_function) for d2
    pub pdf_d2: f32,
    /// present value: $e^{-rT}$
    pub pv: f32,
    /// present value of strike: $Ke^{-rT}$
    pub pv_k: f32,
    /// present value impact of dividend yield: $e^{-qT}$
    pub w: f32,
    /// option price
    pub price: f32,
    /// delta: $\Delta=\frac{\partial V}{\partial S}$  
    /// call: $e^{-qT}N(d_1)$  
    /// put:  $-e^{-qT}N(-d_1)$  
    pub delta: f32,
    /// gamma: $\Gamma=\frac{\partial \Delta}{\partial S}=e^{-qT}\frac{N'(d_1)}{S\sigma\sqrt{T}}$  
    pub gamma: f32,
    /// *Note: vega is not a greek letter, so let us use nu instead, it looks "v-" in latin !*  
    /// $\nu=\frac{\partial V}{\partial \sigma}=Se^{-qT}N'(d_1)\sqrt{T}=Ke^{-rT}N'(d_2)\sqrt{T}$
    pub vega: f32,
    /// $\Theta=-\frac{\partial V}{\partial T}$  
    /// call: $-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}-rKe^{-rT}N(d_2)+qSe^{-qT}N(d_1)$  
    /// put: $-e^{qT}\frac{SN'(d_1)\sigma}{2\sqrt{T}}+rKe^{-rT}N(-d_2)-qSe^{-qT}N(-d_1)$  
    pub theta: f32,
    /// $\rho=\frac{\partial V}{\partial r}$  
    /// call: $KTe^{-rT}N(d_2)$  
    /// put: $-KTe^{-rT}N(-d_2)$  
    pub rho: f32,
    /// *Note: voma is not a greek letter, so let us not make up another name*  
    /// $Voma=\frac{\partial \nu}{\partial \sigma}$  
    /// $Se^{-qT}N'(d_1)\sqrt{T}\frac{d_1 d_2}{\sigma}$  
    pub voma: f32,
    /// payoff  
    /// call: $Max(0, S-K)$  
    /// put: $Max(0, K-S)$  
    pub payoff: f32,
    /// present value of payoff: $e^{-rT}$payoff  
    pub pv_payoff: f32,
}

/// # Pricing output - partial
///
/// This is the output of the call pricer restricted to (1) price call and (2) vega.
#[derive(Debug)]
pub struct OutputPriceVega {
    pub price: f32,
    pub vega: f32,
}

/// # Pricing output - partial
///
/// This is the output of the call pricer restricted to (1) price call and (2) vega and (3) voma.
#[derive(Debug)]
pub struct OutputPriceVegaVoma {
    pub price: f32,
    pub vega: f32,
    pub voma: f32,
}

// const EPSILON: f32 = f32::EPSILON;
const TEST_EPSILON: f32 = 1e-3;

/// ## Custom Check
///
/// This is a custom check for the `Output` struct, vebose to help debug.
impl PartialEq for Output {
    fn eq(&self, other: &Self) -> bool {
        if self.is_call != other.is_call {
            println!("is_call: {} != {}", self.is_call, other.is_call);
            return false;
        }

        if (self.d1 - other.d1).abs() > TEST_EPSILON {
            println!("d1: {} != {}", self.d1, other.d1);
            return false;
        }

        if (self.d2 - other.d2).abs() > TEST_EPSILON {
            println!("d2: {} != {}", self.d2, other.d2);
            return false;
        }

        if (self.cdf_d1 - other.cdf_d1).abs() > TEST_EPSILON {
            println!("cdf_d1: {} != {}", self.cdf_d1, other.cdf_d1);
            return false;
        }

        if (self.cdf_d2 - other.cdf_d2).abs() > TEST_EPSILON {
            println!("cdf_d2: {} != {}", self.cdf_d2, other.cdf_d2);
            return false;
        }

        if (self.pdf_d1 - other.pdf_d1).abs() > TEST_EPSILON {
            println!("pdf_d1: {} != {}", self.pdf_d1, other.pdf_d1);
            return false;
        }

        if (self.pdf_d2 - other.pdf_d2).abs() > TEST_EPSILON {
            println!("pdf_d2: {} != {}", self.pdf_d2, other.pdf_d2);
            return false;
        }

        if (self.pv - other.pv).abs() > TEST_EPSILON {
            println!("pv: {} != {}", self.pv, other.pv);
            return false;
        }

        if (self.pv_k - other.pv_k).abs() > TEST_EPSILON {
            println!("pv_k: {} != {}", self.pv_k, other.pv_k);
            return false;
        }

        if (self.w - other.w).abs() > TEST_EPSILON {
            println!("w: {} != {}", self.w, other.w);
            return false;
        }

        if (self.price - other.price).abs() > TEST_EPSILON {
            println!("price: {} != {}", self.price, other.price);
            return false;
        }

        if (self.delta - other.delta).abs() > TEST_EPSILON {
            println!("delta: {} != {}", self.delta, other.delta);
            return false;
        }

        if (self.gamma - other.gamma).abs() > TEST_EPSILON {
            println!("gamma: {} != {}", self.gamma, other.gamma);
            return false;
        }

        if (self.vega - other.vega).abs() > TEST_EPSILON {
            println!("vega: {} != {}", self.vega, other.vega);
            return false;
        }

        if (self.theta - other.theta).abs() > TEST_EPSILON {
            println!("theta: {} != {}", self.theta, other.theta);
            return false;
        }

        if (self.rho - other.rho).abs() > TEST_EPSILON {
            println!("rho: {} != {}", self.rho, other.rho);
            return false;
        }

        if (self.voma - other.voma).abs() > TEST_EPSILON {
            println!("voma: {} != {}", self.voma, other.voma);
            return false;
        }

        if (self.payoff - other.payoff).abs() > TEST_EPSILON {
            println!("payoff: {} != {}", self.payoff, other.payoff);
            return false;
        }

        if (self.pv_payoff - other.pv_payoff).abs() > TEST_EPSILON {
            println!("pv_payoff: {} != {}", self.pv_payoff, other.pv_payoff);
            return false;
        }

        return true;
    }
}

/// # Pricer Input Error
///
/// This covers the various ways pricing params can be invalid.
#[derive(Error, Debug, Serialize, Deserialize)]
pub enum PriceError {
    /// Negative spot
    #[error("negative spot: {0} - must be positive")]
    NegativeSpot(f32),
    /// Negative strike
    #[error("negative strike: {0} - must be positive")]
    NegativeStrike(f32),
    /// Negative volatility
    #[error("negative volatility: {0} - must be positive")]
    NegativeVol(f32),
    /// Negative maturity
    #[error("negative maturity: {0} - must be positive")]
    NegativeMat(f32),
}

impl BlackScholes {
    pub fn new(input: Input) -> BlackScholes {
        let output = price(&input);
        BlackScholes { input, output }
    }
}

/// ## Main pricing function
///
/// This function checks input params are valid and runs the pricing.  
pub fn price(input: &Input) -> Result<Output, PriceError> {
    if input.spot < 0.0 {
        return Err(PriceError::NegativeSpot(input.spot));
    }
    if input.strike < 0.0 {
        return Err(PriceError::NegativeStrike(input.strike));
    }
    if input.vol < 0.0 {
        return Err(PriceError::NegativeVol(input.vol));
    }
    if input.mat < 0.0 {
        return Err(PriceError::NegativeMat(input.mat));
    }

    let output = match input.is_call {
        true => price_call(input),
        false => price_put(input),
    };

    Ok(output)
}

fn price_call(input: &Input) -> Output {
    let s = input.spot;
    let k = input.strike;
    let t = input.mat;
    let v = input.vol;
    let r = input.rate;
    let q = input.div;

    let sqrt_t = t.sqrt();
    let v_sqrt_t = v * sqrt_t;

    let d1 = (1.0 / v_sqrt_t) * ((s / k).ln() + (r - q + v * v / 2.0) * t);
    let d2 = d1 - v_sqrt_t;

    let cdf_d1 = math::cdf(d1);
    let pdf_d1 = math::pdf(d1);

    let cdf_d2 = math::cdf(d2);
    let pdf_d2 = math::pdf(d2);

    let pv = (-r * t).exp();
    let pv_k = pv * k;
    let w: f32 = (-q * t).exp();

    let price = cdf_d1 * s * w - cdf_d2 * pv_k;

    let delta = w * cdf_d1;
    let gamma = w * pdf_d1 / (s * v_sqrt_t);
    let vega = s * w * pdf_d1 * sqrt_t;
    let theta = -s * w * pdf_d1 * v / (2.0 * sqrt_t) - r * pv_k * cdf_d2 + q * s * w * cdf_d1;
    let rho = pv_k * t * cdf_d2;
    let voma = s * w * pdf_d1 * sqrt_t * d1 * d2 / v;

    let payoff = (s - k).max(0.0);
    let pv_payoff = pv * payoff;

    Output {
        is_call: true,
        d1,
        d2,
        cdf_d1,
        cdf_d2,
        pdf_d1,
        pdf_d2,
        pv,
        pv_k,
        w,
        price,
        delta,
        gamma,
        vega,
        theta,
        rho,
        voma,
        payoff,
        pv_payoff,
    }
}

fn price_put(input: &Input) -> Output {
    let s = input.spot;
    let k = input.strike;
    let t = input.mat;
    let v = input.vol;
    let r = input.rate;
    let q = input.div;

    let sqrt_t = t.sqrt();
    let v_sqrt_t = v * sqrt_t;

    let d1 = (1.0 / v_sqrt_t) * ((s / k).ln() + (r - q + v * v / 2.0) * t);
    let d2 = d1 - v_sqrt_t;

    let cdf_minus_d1 = math::cdf(-d1);
    let pdf_d1 = math::pdf(d1);

    let cdf_minus_d2 = math::cdf(-d2);
    let pdf_d2 = math::pdf(d2);

    let pv = (-r * t).exp();
    let pv_k = pv * k;
    let w = (-q * t).exp();

    let price = cdf_minus_d2 * pv_k - cdf_minus_d1 * s * w;

    let delta = -w * cdf_minus_d1;
    let gamma = w * pdf_d1 / (s * v_sqrt_t);
    let vega = s * w * pdf_d1 * sqrt_t;
    let theta =
        -s * w * pdf_d1 * v / (2.0 * sqrt_t) + r * pv_k * cdf_minus_d2 - q * s * w * cdf_minus_d1;
    let rho = -pv_k * t * cdf_minus_d2;
    let voma = s * w * pdf_d1 * sqrt_t * d1 * d2 / v;

    let payoff = (k - s).max(0.0);
    let pv_payoff = pv * payoff;

    Output {
        is_call: false,
        d1,
        d2,
        cdf_d1: cdf_minus_d1,
        cdf_d2: cdf_minus_d2,
        pdf_d1,
        pdf_d2,
        pv,
        pv_k,
        w,
        price,
        delta,
        gamma,
        vega,
        theta,
        rho,
        voma,
        payoff,
        pv_payoff,
    }
}

/// # Custom pricing function
///
/// This function calculates the price and vega of a call option.  
/// It is used by [ImpliedVol](../implied_vol).
pub fn call_price_vega(s: f32, k: f32, r: f32, t: f32, v: f32, q: f32) -> OutputPriceVega {
    let sqrt_t = t.sqrt();
    let v_sqrt_t = v * sqrt_t;

    let d1 = (1.0 / v_sqrt_t) * ((s / k).ln() + (r - q + v * v / 2.0) * t);
    let d2 = d1 - v_sqrt_t;

    let cdf_d1 = math::cdf(d1);
    let pdf_d1 = math::pdf(d1);

    let cdf_d2 = math::cdf(d2);

    let pv = (-r * t).exp();
    let pv_k = pv * k;
    let w = (-q * t).exp();

    let price = cdf_d1 * s * w - cdf_d2 * pv_k;
    let vega = s * w * pdf_d1 * sqrt_t;

    OutputPriceVega { price, vega }
}

/// # Custom pricing function
///
/// This function calculates the price, vega and voma of a call option.
/// It is used by [ImpliedVol](../implied_vol).
pub fn call_price_vega_voma(s: f32, k: f32, r: f32, t: f32, v: f32, q: f32) -> OutputPriceVegaVoma {
    let sqrt_t = t.sqrt();
    let v_sqrt_t = v * sqrt_t;

    let d1 = (1.0 / v_sqrt_t) * ((s / k).ln() + (r - q + v * v / 2.0) * t);
    let d2 = d1 - v_sqrt_t;

    let cdf_d1 = math::cdf(d1);
    let pdf_d1 = math::pdf(d1);

    let cdf_d2 = math::cdf(d2);

    let pv = (-r * t).exp();
    let pv_k = pv * k;
    let w = (-q * t).exp();

    let price = cdf_d1 * s * w - cdf_d2 * pv_k;
    let vega = s * w * pdf_d1 * sqrt_t;
    let voma = s * w * pdf_d1 * sqrt_t * d1 * d2 / v;

    OutputPriceVegaVoma { price, vega, voma }
}

#[cfg(test)]
mod tests {

    use crate::black_scholes::{BlackScholes, Input, Output, TEST_EPSILON};

    struct InputTest {
        input: Input,
        output: Output,
    }
    struct InputTestParity {
        input_call: Input,
        input_put: Input,
    }

    #[test]
    fn sample_prices() {
        let test_data: Vec<InputTest> = vec![
            InputTest {
                input: Input {
                    is_call: true,
                    spot: 105.0,
                    strike: 100.0,
                    mat: 3.0,
                    vol: 0.15,
                    rate: 0.03,
                    div: 0.01,
                },
                output: Output {
                    is_call: true,
                    d1: 0.5486,
                    d2: 0.2888,
                    cdf_d1: 0.7084,
                    cdf_d2: 0.6136,
                    pdf_d1: 0.3432,
                    pdf_d2: 0.3826,
                    pv: 0.9139,
                    pv_k: 91.3931,
                    w: 0.9704,
                    price: 16.098,
                    delta: 0.6874,
                    gamma: 0.0122,
                    vega: 60.5716,
                    theta: -2.4750,
                    rho: 168.2485,
                    voma: 63.9890,
                    payoff: 5.0,
                    pv_payoff: 4.569,
                },
            },
            InputTest {
                input: Input {
                    is_call: false,
                    spot: 90.0,
                    strike: 100.0,
                    mat: 3.0,
                    vol: 0.15,
                    rate: 0.03,
                    div: 0.01,
                },
                output: Output {
                    is_call: false,
                    d1: -0.0447,
                    d2: -0.3045,
                    cdf_d1: 0.5178,
                    cdf_d2: 0.6196,
                    pdf_d1: 0.3985,
                    pdf_d2: 0.3809,
                    pv: 0.9139,
                    pv_k: 91.3931,
                    w: 0.9704,
                    price: 11.402,
                    delta: -0.5025,
                    gamma: 0.0165,
                    vega: 60.2908,
                    theta: -0.2607,
                    rho: -169.8884,
                    voma: 5.4694,
                    payoff: 10.0,
                    pv_payoff: 9.1393,
                },
            },
        ];

        println!("epsilon used in tests: {}", TEST_EPSILON);
        test_data.iter().for_each(|test| {
            let bs = BlackScholes::new(test.input.clone());
            assert!(bs.output.is_ok());
            let out = bs.output.unwrap();
            assert_eq!(out, test.output);
        });
    }

    #[test]
    fn put_call_parity() {
        let test_data: Vec<InputTestParity> = vec![InputTestParity {
            input_call: Input {
                is_call: true,
                spot: 105.6,
                strike: 100.0,
                mat: 3.2,
                vol: 0.25,
                rate: 0.03,
                div: 0.01,
            },
            input_put: Input {
                is_call: false,
                spot: 105.6,
                strike: 100.0,
                mat: 3.2,
                vol: 0.25,
                rate: 0.03,
                div: 0.01,
            },
        }];

        test_data.iter().for_each(|test| {
            let bs_call = BlackScholes::new(test.input_call.clone());
            let bs_put = BlackScholes::new(test.input_put.clone());

            assert!(bs_call.output.is_ok());
            assert!(bs_put.output.is_ok());

            let call_out = bs_call.output.unwrap();
            let put_out = bs_put.output.unwrap();

            let r = test.input_call.rate;
            let q = test.input_call.div;
            let mat = test.input_call.mat;

            let call = call_out.price;
            let put = put_out.price;
            let pv = call_out.pv;
            let fwd = test.input_call.spot * ((r - q) * mat).exp();
            let strike = test.input_call.strike;
            let parity = (call - put) - pv * (fwd - strike);

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