|
| 1 | +use super::vanila_option::{EquityOption}; |
| 2 | +use super::utils::{Engine}; |
| 3 | + |
| 4 | +use crate::core::trade::{OptionType,Transection}; |
| 5 | +use crate::core::utils::{ContractStyle}; |
| 6 | +use ndarray::{Array, Array2,Array1, ArrayBase, Ix1, OwnedRepr, s}; |
| 7 | +//use num_integer::Integer; |
| 8 | +/// finite difference model for European and American options |
| 9 | +pub fn npv(option: &&EquityOption) -> f64 { |
| 10 | + assert!(option.volatility >= 0.0); |
| 11 | + assert!(option.time_to_maturity() >= 0.0); |
| 12 | + assert!(option.underlying_price.value >= 0.0); |
| 13 | + let strike_price = option.strike_price; |
| 14 | + let time_to_maturity = option.time_to_maturity(); |
| 15 | + let underlying_price = option.underlying_price.value; |
| 16 | + let volatility = option.volatility; |
| 17 | + let risk_free_rate = option.risk_free_rate; |
| 18 | + let dividend_yield = option.dividend_yield; |
| 19 | + let time_steps:f64 = 1000.0; |
| 20 | + //let time_steps:f64 = time_to_maturity/0.001 as f64; |
| 21 | + |
| 22 | + let mut spot_steps = (time_steps / 50.0) as usize; //should be even number |
| 23 | + //let spot_steps:usize = 20; |
| 24 | + if spot_steps % 2 != 0{ |
| 25 | + spot_steps = spot_steps + 1; |
| 26 | + }// should be even number |
| 27 | + |
| 28 | + if option.option_type == OptionType::Call { |
| 29 | + return fd(underlying_price,strike_price,risk_free_rate,dividend_yield,volatility, |
| 30 | + time_to_maturity,spot_steps,time_steps); |
| 31 | + } else { |
| 32 | + //TODO implement for put option |
| 33 | + //println!("Not implemented""); |
| 34 | + return 0.0; |
| 35 | + } |
| 36 | + |
| 37 | +} |
| 38 | +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{ |
| 39 | + let ds = 2.0*s0 / spot_steps as f64; |
| 40 | + |
| 41 | + let M:i32 = (spot_steps as f64+(spot_steps as f64)/2.0) as i32; // 40 +20 = 60-20 = 40 |
| 42 | + //let ds = 2.0*s0 / spot_steps as f64; |
| 43 | + //let M:i32 = spot_steps as i32; |
| 44 | + let dt = time_to_mat/(time_steps as f64); |
| 45 | + let time_steps = time_steps as i32; |
| 46 | + // convert float to nearest integer |
| 47 | + |
| 48 | + //println!(" h {:?}",dt / (ds*ds)); |
| 49 | + //let sigma:f64 = 0.3; |
| 50 | + //let r = 0.05; |
| 51 | + //let q = 0.01; |
| 52 | + let r = risk_free_rate-dividend_yield; |
| 53 | + let mut price_grid:Array2<f64> = Array2::zeros((M as usize +1,time_steps as usize+1)); |
| 54 | + // Underlying price Grid |
| 55 | + for j in 0..time_steps+1 { |
| 56 | + for i in 0..M+1{ |
| 57 | + price_grid[[(M-i) as usize,j as usize]] = (i as f64)*ds as f64; |
| 58 | + } |
| 59 | + } |
| 60 | + let mm = M as usize - ii as usize; |
| 61 | + println!("price_ {:?}",price_grid[[mm as usize as usize,0]]); |
| 62 | + let mut v_grid:Array2<f64> = Array2::zeros((M as usize +1,time_steps as usize+1)); |
| 63 | + // Boundary condition |
| 64 | + // for j in 0..time_steps+1 { |
| 65 | + // for i in 0..M+1{ |
| 66 | + // v_grid[[(M-i) as usize,j as usize]] = (price_grid[[(M-i) as usize,j as usize]]-k).max(0.0); |
| 67 | + // } |
| 68 | + // } |
| 69 | + // Boundary condition |
| 70 | + for i in 0..M+1{ |
| 71 | + v_grid[[(M-i) as usize,time_steps as usize]] = (price_grid[[(M-i) as usize,time_steps as usize]]-k).max(0.0); |
| 72 | + } |
| 73 | + |
| 74 | + let mut pd:Vec<f64> = Vec::with_capacity((M + 1) as usize); |
| 75 | + let mut b:Vec<f64> = Vec::with_capacity((M + 1) as usize); |
| 76 | + let mut x_:Vec<f64> = Vec::with_capacity((M + 1) as usize); |
| 77 | + let mut y_:Vec<f64> = Vec::with_capacity((M + 1) as usize); |
| 78 | + let mut pu:Vec<f64> = Vec::with_capacity((M + 1) as usize); |
| 79 | + //let ss = price_grid.slice(s![0..M+1,time_steps]).to_vec(); |
| 80 | + //let mut xx_:Vec<f64> = Vec::with_capacity((M + 1) as usize); |
| 81 | + for j in 0..M + 1 { |
| 82 | + // ssj = (j as f64)*ds; |
| 83 | + //let x = r * ss[j as usize] / ds; |
| 84 | + let x = r*((M-j) as f64); |
| 85 | + //let y = sigma.powi(2) * (ss[j as usize].powi(2)) / ds.powi(2); |
| 86 | + let y = sigma.powi(2) * ((M-j) as f64).powi(2); |
| 87 | + x_.push(x); |
| 88 | + y_.push(y); |
| 89 | + //xx_.push(ss[j as usize] / ds); |
| 90 | + b.push(1.0 + dt * r + y * dt*0.5); //0.5 * dt * (j as f64)*((r-q) - sigma.powi(2)*(j as f64)); |
| 91 | + pu.push(-0.25 * dt * (x + y)); //0.5 * dt * (j as f64)*((r-q) + sigma.powi(2)*(j as f64)) |
| 92 | + pd.push(0.25 * dt * (x - y)); //-0.5*dt*(j as f64)*((r-q) + sigma.powi(2)*(j as f64)) |
| 93 | + } |
| 94 | + |
| 95 | + for i in (1..time_steps+1).rev(){ |
| 96 | + let mut d = v_grid.slice(s![0..M as usize+1,i]).to_vec(); |
| 97 | + d[0] = d[0]*(1.0-pu[0]); |
| 98 | + for j in 1..M{ |
| 99 | + 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]) + |
| 100 | + 0.25*y_[j as usize]*dt*(d[(j-1) as usize]+d[(j+1) as usize] - 2.0*d[j as usize] ); |
| 101 | + } |
| 102 | + let x = thomas_algorithm(&pu[0..M as usize], &b, &pd[1..(M+1) as usize], &d); |
| 103 | + for j in 0..M+1{ |
| 104 | + v_grid[[j as usize,(i-1) as usize]] = x[j as usize]; |
| 105 | + } |
| 106 | + } |
| 107 | + |
| 108 | + return v_grid[[spot_steps as usize,0]]; |
| 109 | + |
| 110 | +} |
| 111 | +pub fn thomas_algorithm(a: &[f64], b: &[f64], c: &[f64], d: &[f64]) -> Vec<f64> { |
| 112 | + ///https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm |
| 113 | + /// Solves Ax = d where A is a tridiagonal matrix consisting of vectors a, b, c |
| 114 | + let n = d.len(); |
| 115 | + let mut c_ = c.to_vec(); |
| 116 | + let mut d_ = d.to_vec(); |
| 117 | + let mut x: Vec<f64> = vec![0.0; n]; |
| 118 | + |
| 119 | + // Adjust for the upper boundary condition |
| 120 | + //d_[0] = d_[0]*(1.0-a[0]); |
| 121 | + |
| 122 | + c_[0] = c_[0] / b[0]; |
| 123 | + d_[0] = d_[0] / b[0]; |
| 124 | + for i in 1..n-1 { |
| 125 | + let id = 1.0 / (b[i] - a[i-1] * c_[i - 1]); |
| 126 | + c_[i] = c_[i] * id; |
| 127 | + d_[i] = (d_[i] - a[i-1] * d_[i - 1]) * id; |
| 128 | + } |
| 129 | + d_[n - 1] = (d_[n - 1] - a[n - 2] * d_[n - 2]) / (b[n - 1] - a[n - 2] * c_[n - 2]); |
| 130 | + |
| 131 | + x[n - 1] = d_[n - 1]; |
| 132 | + for i in (0..n - 1).rev() { |
| 133 | + x[i] = d_[i] - c_[i] * x[i + 1]; |
| 134 | + } |
| 135 | + x |
| 136 | +} |
0 commit comments