From f837623db7bda533775f579a8dbcd62d0393fb8f Mon Sep 17 00:00:00 2001 From: Siddharth Date: Fri, 29 Dec 2023 04:25:57 -0500 Subject: [PATCH] added finite difference for call option --- README.md | 1 + src/equity/finite_difference.rs | 136 ++++++++++++++++++++++++++++++++ src/equity/mod.rs | 1 + src/equity/utils.rs | 1 + src/equity/vanila_option.rs | 26 +++--- src/utils/parse_json.rs | 7 +- 6 files changed, 157 insertions(+), 15 deletions(-) create mode 100644 src/equity/finite_difference.rs diff --git a/README.md b/README.md index c3831b6..5a19ae9 100644 --- a/README.md +++ b/README.md @@ -2,6 +2,7 @@ [![License: MIT](https://img.shields.io/badge/License-MIT-blue.svg)](https://opensource.org/licenses/MIT) ![Crates.io](https://img.shields.io/crates/dr/rustyqlib) ![Crates.io](https://img.shields.io/crates/v/rustyqlib) +[![codecov](https://codecov.io/gh/siddharthqs/RustyQLib/graph/badge.svg?token=879K6LTTR4)](https://codecov.io/gh/siddharthqs/RustyQLib) # RUSTYQLib :Pricing Options with Confidence using JSON RustyQLib is a lightweight yet robust quantitative finance library designed for pricing options. Built entirely in Rust, it offers a unique combination of safety, performance, and expressiveness that is crucial diff --git a/src/equity/finite_difference.rs b/src/equity/finite_difference.rs new file mode 100644 index 0000000..907233a --- /dev/null +++ b/src/equity/finite_difference.rs @@ -0,0 +1,136 @@ +use super::vanila_option::{EquityOption}; +use super::utils::{Engine}; + +use crate::core::trade::{OptionType,Transection}; +use crate::core::utils::{ContractStyle}; +use ndarray::{Array, Array2,Array1, ArrayBase, Ix1, OwnedRepr, s}; +//use num_integer::Integer; +/// finite difference model for European and American options +pub fn npv(option: &&EquityOption) -> f64 { + assert!(option.volatility >= 0.0); + assert!(option.time_to_maturity() >= 0.0); + assert!(option.underlying_price.value >= 0.0); + let strike_price = option.strike_price; + let time_to_maturity = option.time_to_maturity(); + let underlying_price = option.underlying_price.value; + let volatility = option.volatility; + let risk_free_rate = option.risk_free_rate; + let dividend_yield = option.dividend_yield; + let time_steps:f64 = 1000.0; + //let time_steps:f64 = time_to_maturity/0.001 as f64; + + let mut spot_steps = (time_steps / 50.0) as usize; //should be even number + //let spot_steps:usize = 20; + if spot_steps % 2 != 0{ + spot_steps = spot_steps + 1; + }// should be even number + + if option.option_type == OptionType::Call { + return fd(underlying_price,strike_price,risk_free_rate,dividend_yield,volatility, + time_to_maturity,spot_steps,time_steps); + } else { + //TODO implement for put option + //println!("Not implemented""); + return 0.0; + } + +} +fn fd(s0:f64,k:f64,risk_free_rate:f64,dividend_yield:f64,sigma:f64,time_to_mat:f64,spot_steps:usize,time_steps:f64)->f64{ + let ds = 2.0*s0 / spot_steps as f64; + + let M:i32 = (spot_steps as f64+(spot_steps as f64)/2.0) as i32; // 40 +20 = 60-20 = 40 + //let ds = 2.0*s0 / spot_steps as f64; + //let M:i32 = spot_steps as i32; + let dt = time_to_mat/(time_steps as f64); + let time_steps = time_steps as i32; + // convert float to nearest integer + + //println!(" h {:?}",dt / (ds*ds)); + //let sigma:f64 = 0.3; + //let r = 0.05; + //let q = 0.01; + let r = risk_free_rate-dividend_yield; + let mut price_grid:Array2 = Array2::zeros((M as usize +1,time_steps as usize+1)); + // Underlying price Grid + for j in 0..time_steps+1 { + for i in 0..M+1{ + price_grid[[(M-i) as usize,j as usize]] = (i as f64)*ds as f64; + } + } + let mm = M as usize - ii as usize; + println!("price_ {:?}",price_grid[[mm as usize as usize,0]]); + let mut v_grid:Array2 = Array2::zeros((M as usize +1,time_steps as usize+1)); + // Boundary condition + // for j in 0..time_steps+1 { + // for i in 0..M+1{ + // v_grid[[(M-i) as usize,j as usize]] = (price_grid[[(M-i) as usize,j as usize]]-k).max(0.0); + // } + // } + // Boundary condition + for i in 0..M+1{ + v_grid[[(M-i) as usize,time_steps as usize]] = (price_grid[[(M-i) as usize,time_steps as usize]]-k).max(0.0); + } + + let mut pd:Vec = Vec::with_capacity((M + 1) as usize); + let mut b:Vec = Vec::with_capacity((M + 1) as usize); + let mut x_:Vec = Vec::with_capacity((M + 1) as usize); + let mut y_:Vec = Vec::with_capacity((M + 1) as usize); + let mut pu:Vec = Vec::with_capacity((M + 1) as usize); + //let ss = price_grid.slice(s![0..M+1,time_steps]).to_vec(); + //let mut xx_:Vec = Vec::with_capacity((M + 1) as usize); + for j in 0..M + 1 { + // ssj = (j as f64)*ds; + //let x = r * ss[j as usize] / ds; + let x = r*((M-j) as f64); + //let y = sigma.powi(2) * (ss[j as usize].powi(2)) / ds.powi(2); + let y = sigma.powi(2) * ((M-j) as f64).powi(2); + x_.push(x); + y_.push(y); + //xx_.push(ss[j as usize] / ds); + b.push(1.0 + dt * r + y * dt*0.5); //0.5 * dt * (j as f64)*((r-q) - sigma.powi(2)*(j as f64)); + pu.push(-0.25 * dt * (x + y)); //0.5 * dt * (j as f64)*((r-q) + sigma.powi(2)*(j as f64)) + pd.push(0.25 * dt * (x - y)); //-0.5*dt*(j as f64)*((r-q) + sigma.powi(2)*(j as f64)) + } + + for i in (1..time_steps+1).rev(){ + let mut d = v_grid.slice(s![0..M as usize+1,i]).to_vec(); + d[0] = d[0]*(1.0-pu[0]); + for j in 1..M{ + d[j as usize] = d[j as usize] +0.25*x_[j as usize]*dt*(d[(j-1) as usize]-d[(j+1) as usize]) + + 0.25*y_[j as usize]*dt*(d[(j-1) as usize]+d[(j+1) as usize] - 2.0*d[j as usize] ); + } + let x = thomas_algorithm(&pu[0..M as usize], &b, &pd[1..(M+1) as usize], &d); + for j in 0..M+1{ + v_grid[[j as usize,(i-1) as usize]] = x[j as usize]; + } + } + + return v_grid[[spot_steps as usize,0]]; + +} +pub fn thomas_algorithm(a: &[f64], b: &[f64], c: &[f64], d: &[f64]) -> Vec { + ///https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm + /// Solves Ax = d where A is a tridiagonal matrix consisting of vectors a, b, c + let n = d.len(); + let mut c_ = c.to_vec(); + let mut d_ = d.to_vec(); + let mut x: Vec = vec![0.0; n]; + + // Adjust for the upper boundary condition + //d_[0] = d_[0]*(1.0-a[0]); + + c_[0] = c_[0] / b[0]; + d_[0] = d_[0] / b[0]; + for i in 1..n-1 { + let id = 1.0 / (b[i] - a[i-1] * c_[i - 1]); + c_[i] = c_[i] * id; + d_[i] = (d_[i] - a[i-1] * d_[i - 1]) * id; + } + d_[n - 1] = (d_[n - 1] - a[n - 2] * d_[n - 2]) / (b[n - 1] - a[n - 2] * c_[n - 2]); + + x[n - 1] = d_[n - 1]; + for i in (0..n - 1).rev() { + x[i] = d_[i] - c_[i] * x[i + 1]; + } + x +} diff --git a/src/equity/mod.rs b/src/equity/mod.rs index 0ec5813..300c749 100644 --- a/src/equity/mod.rs +++ b/src/equity/mod.rs @@ -6,3 +6,4 @@ pub mod utils; pub mod binomial; pub mod build_contracts; pub mod vol_surface; +pub mod finite_difference; diff --git a/src/equity/utils.rs b/src/equity/utils.rs index 09ac36c..e4b2f86 100644 --- a/src/equity/utils.rs +++ b/src/equity/utils.rs @@ -4,4 +4,5 @@ pub enum Engine{ BlackScholes, MonteCarlo, Binomial, + FiniteDifference } \ No newline at end of file diff --git a/src/equity/vanila_option.rs b/src/equity/vanila_option.rs index b4ab5ce..08c9503 100644 --- a/src/equity/vanila_option.rs +++ b/src/equity/vanila_option.rs @@ -1,6 +1,5 @@ use chrono::{Datelike, Local, NaiveDate}; -use crate::equity::montecarlo; -use crate::equity::binomial; +use crate::equity::{binomial,finite_difference,montecarlo}; use super::super::core::termstructure::YieldTermStructure; use super::super::core::quotes::Quote; use super::super::core::traits::{Instrument,Greeks}; @@ -27,9 +26,11 @@ impl Instrument for EquityOption { let value = binomial::npv(&self); value } - _ => { - 0.0 + Engine::FiniteDifference => { + let value = finite_difference::npv(&self); + value } + } } } @@ -67,10 +68,10 @@ impl EquityOption { let rates = vec![0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05]; let ts = YieldTermStructure::new(date, rates); let option_type = &market_data.option_type; - let side: trade::OptionType; + let side: OptionType; match option_type.trim() { - "C" | "c" | "Call" | "call" => side = trade::OptionType::Call, - "P" | "p" | "Put" | "put" => side = trade::OptionType::Put, + "C" | "c" | "Call" | "call" => side = OptionType::Call, + "P" | "p" | "Put" | "put" => side = OptionType::Put, _ => panic!("Invalide side argument! Side has to be either 'C' or 'P'."), } let maturity_date = &market_data.maturity; @@ -79,7 +80,7 @@ impl EquityOption { let risk_free_rate = Some(market_data.risk_free_rate).unwrap(); let dividend = Some(market_data.dividend).unwrap(); - let mut op = 0.0; + //let mut op = 0.0; let option_price = Quote::new(match market_data.option_price { Some(x) => x, @@ -110,15 +111,18 @@ impl EquityOption { valuation_date: today.naive_utc(), }; match data.pricer.trim() { - "Analytical" | "analytical" => { + "Analytical" | "analytical"|"bs" => { option.engine = Engine::BlackScholes; } - "MonteCarlo" | "montecarlo" | "MC" => { + "MonteCarlo" | "montecarlo" | "MC"|"mc" => { option.engine = Engine::MonteCarlo; } - "Binomial" | "binomial" => { + "Binomial" | "binomial"|"bino" => { option.engine = Engine::Binomial; } + "FiniteDifference" | "finitdifference" |"FD" |"fd" => { + option.engine = Engine::FiniteDifference; + } _ => { panic!("Invalid pricer"); } diff --git a/src/utils/parse_json.rs b/src/utils/parse_json.rs index e398ca5..f2f997c 100644 --- a/src/utils/parse_json.rs +++ b/src/utils/parse_json.rs @@ -91,7 +91,7 @@ pub fn parse_contract(mut file: &mut File,output_filename: &str) { println!("No contracts found in JSON file"); return; } - // parallel processing of each contract + // parallel processing of each contract using rayon let mut output_vec: Vec<_> = list_contracts.contracts.par_iter().enumerate() .map(|(index,data)| (index,process_contract(data))) .collect(); @@ -104,7 +104,7 @@ pub fn parse_contract(mut file: &mut File,output_filename: &str) { file.write_all(output_str.as_bytes()).expect("Failed to write to file"); } pub fn process_contract(data: &Contract) -> String { - //println!("Processing {:?}",data); + let date = vec![0.01,0.02,0.05,0.1,0.5,1.0,2.0,3.0]; let rates = vec![0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05]; let ts = YieldTermStructure::new(date,rates); @@ -112,7 +112,6 @@ pub fn process_contract(data: &Contract) -> String { if data.action=="PV" && data.asset=="EQ"{ //let market_data = data.market_data.clone().unwrap(); let option = EquityOption::from_json(&data); - let contract_output = ContractOutput{pv:option.npv(),delta:option.delta(),gamma:option.gamma(), vega:option.vega(),theta:option.theta(),rho:option.rho(), error: None }; println!("Theoretical Price ${}", contract_output.pv); @@ -180,7 +179,7 @@ pub fn process_contract(data: &Contract) -> String { let maturity_date = rates::utils::convert_mm_to_date(maturity_date_str); let start_date = rates::utils::convert_mm_to_date(start_date_str); println!("Maturity Date {:?}",maturity_date); - let mut deposit = rates::deposits::Deposit { + let mut deposit = Deposit { start_date: start_date, maturity_date: maturity_date, valuation_date: current_date.naive_utc(),