From 5f230a4ea45b544b22042525febac718d4fa149b Mon Sep 17 00:00:00 2001 From: Daniel Boros Date: Fri, 27 Sep 2024 22:47:25 +0200 Subject: [PATCH] feat: upgrade heston calibrator --- stochastic-rs-quant/Cargo.toml | 2 + stochastic-rs-quant/src/calibration/heston.rs | 92 ++++--- stochastic-rs-quant/src/lib.rs | 10 + stochastic-rs-quant/src/pricer/heston.rs | 232 ++++++++++++------ 4 files changed, 228 insertions(+), 108 deletions(-) diff --git a/stochastic-rs-quant/Cargo.toml b/stochastic-rs-quant/Cargo.toml index f1bee0e..b18cb22 100644 --- a/stochastic-rs-quant/Cargo.toml +++ b/stochastic-rs-quant/Cargo.toml @@ -11,6 +11,8 @@ mimalloc = { version = "0.1.43", optional = true } nalgebra = "0.33.0" num-complex = "0.4.6" quadrature = "0.1.2" +rand = "0.8.5" +rand_distr = "0.4.3" stochastic-rs = { path = "../stochastic-rs-core" } tikv-jemallocator = { version = "0.6.0", optional = true } diff --git a/stochastic-rs-quant/src/calibration/heston.rs b/stochastic-rs-quant/src/calibration/heston.rs index c54f679..b40d738 100644 --- a/stochastic-rs-quant/src/calibration/heston.rs +++ b/stochastic-rs-quant/src/calibration/heston.rs @@ -1,6 +1,8 @@ use levenberg_marquardt::{LeastSquaresProblem, LevenbergMarquardt}; -use nalgebra::{storage, DMatrix, DVector, Dyn, OMatrix, OVector}; -use stochastic_rs::{volatility::heston::Heston, Sampling2D}; +use nalgebra::{storage, DMatrix, DVector, Dyn, OVector}; +use rand::thread_rng; +use rand_distr::{Distribution, Normal}; +use stochastic_rs::volatility::heston::Heston; use crate::pricer::heston::HestonPricer; @@ -29,50 +31,58 @@ impl<'a> LeastSquaresProblem for Calibrator<'a> { } fn residuals(&self) -> Option> { - let model_prices = vec![100.0]; - let market_prices = vec![89.0]; + let model_prices = self.pricer.prices.as_ref().unwrap(); + let call_prices = unsafe { + model_prices + .v + .clone() + .iter() + .map(|x| x.0) + .collect::>() + }; + // Add some noise to the market prices + let market_prices = call_prices + .iter() + .map(|x| *x + Normal::new(1.0, 0.5).unwrap().sample(&mut thread_rng())) + .collect::>(); - let residuals = model_prices + let residuals = call_prices .iter() .zip(market_prices.iter()) - .map(|(m, p)| m - p) + .map(|(x, y)| x - y) .collect::>(); Some(DVector::from_vec(residuals)) } fn jacobian(&self) -> Option> { - let dC_dv0 = self.pricer.dC_dv0(); - let dC_dtheta = self.pricer.dC_dtheta(); - let dC_drho = self.pricer.dC_drho(); - let dC_dkappa = self.pricer.dC_dkappa(); - let dC_dsigma = self.pricer.dC_dsigma(); + let derivates = self.pricer.derivates.as_ref().unwrap(); + let derivates = unsafe { derivates.v.clone().to_vec() }; - let jacobian = DMatrix::from_vec(1, 5, vec![dC_dv0, dC_dtheta, dC_drho, dC_dkappa, dC_dsigma]); + // Convert flattened vector to a matrix + let jacobian = DMatrix::from_vec(derivates.len() / 5, 5, derivates); Some(jacobian) } } pub struct HestonCalibrator { - //model: Heston, + //model: , pricer: HestonPricer, } impl HestonCalibrator { #[must_use] - pub fn new(model: Heston, pricer: HestonPricer) -> Self { + pub fn new(pricer: HestonPricer) -> Self { Self { pricer } } - pub fn calibrate(&self) { - //let [s, v] = self.model.sample(); - let price = self.pricer.price(); - - let (_result, report) = LevenbergMarquardt::new().minimize(Calibrator::new( + pub fn calibrate(&mut self) { + self.pricer.price(); + let (result, report) = LevenbergMarquardt::new().minimize(Calibrator::new( DVector::from_vec(vec![0.05, 0.05, -0.8, 5.0, 0.5]), &self.pricer, )); - println!("{:?}", report.objective_function); + println!("{:?}", result.p); println!("{:?}", report.number_of_evaluations); println!("{:?}", report.termination); } @@ -80,30 +90,36 @@ impl HestonCalibrator { #[cfg(test)] mod tests { + use std::mem::ManuallyDrop; + use crate::ValueOrVec; use super::*; #[test] fn test_calibrate() { - let calibrator = HestonCalibrator::new( - Heston::default(), - HestonPricer { - s0: 100.0, - v0: 0.05, - k: 100.0, - r: 0.03, - q: 0.02, - rho: -0.8, - kappa: 5.0, - theta: 0.05, - sigma: 0.5, - lambda: Some(0.0), - tau: Some(ValueOrVec { x: 0.5 }), // Single f64 tau value - eval: None, - expiry: None, - }, - ); + let majurities = (0..=100) + .map(|x| 0.5 + 0.1 * x as f64) + .collect::>(); + let mut calibrator = HestonCalibrator::new(HestonPricer { + s0: 100.0, + v0: 0.05, + k: 100.0, + r: 0.03, + q: 0.02, + rho: -0.8, + kappa: 5.0, + theta: 0.05, + sigma: 0.5, + lambda: Some(0.0), + tau: Some(ValueOrVec { + v: ManuallyDrop::new(majurities.clone()), + }), // Single f64 tau value + eval: None, + expiry: None, + prices: None, + derivates: None, + }); calibrator.calibrate(); } } diff --git a/stochastic-rs-quant/src/lib.rs b/stochastic-rs-quant/src/lib.rs index 9e63a56..c7a7ad7 100644 --- a/stochastic-rs-quant/src/lib.rs +++ b/stochastic-rs-quant/src/lib.rs @@ -31,6 +31,16 @@ impl Clone for ValueOrVec { } } +impl Clone for ValueOrVec<(f64, f64)> { + fn clone(&self) -> Self { + unsafe { + Self { + v: ManuallyDrop::new(self.v.clone().to_vec()), + } + } + } +} + /// Implement the `Clone` trait for `ValueOrVec`. impl Clone for ValueOrVec { fn clone(&self) -> Self { diff --git a/stochastic-rs-quant/src/pricer/heston.rs b/stochastic-rs-quant/src/pricer/heston.rs index 9ab1e89..b3fae49 100644 --- a/stochastic-rs-quant/src/pricer/heston.rs +++ b/stochastic-rs-quant/src/pricer/heston.rs @@ -5,7 +5,7 @@ use quadrature::double_exponential; use crate::ValueOrVec; -#[derive(Default)] +#[derive(Default, Clone)] pub struct HestonPricer { /// Initial stock price pub s0: f64, @@ -33,6 +33,10 @@ pub struct HestonPricer { pub eval: Option>, /// Expiration date pub expiry: Option>, + /// Prices of European call and put options + pub(crate) prices: Option>, + /// Partial derivative of the C function with respect to the parameters + pub(crate) derivates: Option>, } impl HestonPricer { @@ -49,93 +53,59 @@ impl HestonPricer { kappa: params.kappa, theta: params.theta, sigma: params.sigma, - lambda: params.lambda, + lambda: Some(params.lambda.unwrap_or(0.0)), tau: params.tau.clone(), eval: params.eval.clone(), expiry: params.expiry.clone(), + prices: None, + derivates: None, } } /// Calculate the price of a European call option using the Heston model /// https://quant.stackexchange.com/a/18686 - pub fn price(&self) -> ValueOrVec<(f64, f64)> { + pub fn price(&mut self) -> ValueOrVec<(f64, f64)> { if self.tau.is_none() && self.eval.is_none() && self.expiry.is_none() { panic!("At least 2 of tau, eval, and expiry must be provided"); } - let lambda = self.lambda.unwrap_or(0.0); - - let u = |j: u8| match j { - 1 => 0.5, - 2 => -0.5, - _ => panic!("Invalid j"), - }; - - let b = |j: u8| match j { - 1 => self.kappa + lambda - self.rho * self.sigma, - 2 => self.kappa + lambda, - _ => panic!("Invalid j"), - }; - - let d = |j: u8, phi: f64| -> Complex64 { - ((b(j) - self.rho * self.sigma * phi * Complex64::i()).powi(2) - - self.sigma.powi(2) * (2.0 * Complex64::i() * u(j) * phi - phi.powi(2))) - .sqrt() - }; - - let g = |j: u8, phi: f64| -> Complex64 { - (b(j) - self.rho * self.sigma * Complex64::i() * phi + d(j, phi)) - / (b(j) - self.rho * self.sigma * Complex64::i() * phi - d(j, phi)) - }; - - let C = |j: u8, phi: f64, tau: f64| -> Complex64 { - (self.r - self.q) * Complex64::i() * phi * tau - + (self.kappa * self.theta / self.sigma.powi(2)) - * ((b(j) - self.rho * self.sigma * Complex64::i() * phi + d(j, phi)) * tau - - 2.0 * ((1.0 - g(j, phi) * (d(j, phi) * tau).exp()) / (1.0 - g(j, phi))).ln()) - }; - - let D = |j: u8, phi: f64, tau: f64| -> Complex64 { - ((b(j) - self.rho * self.sigma * Complex64::i() * phi + d(j, phi)) / self.sigma.powi(2)) - * ((1.0 - (d(j, phi) * tau).exp()) / (1.0 - g(j, phi) * (d(j, phi) * tau).exp())) - }; - - let f = |j: u8, phi: f64, tau: f64| -> Complex64 { - (C(j, phi, tau) + D(j, phi, tau) * self.v0 + Complex64::i() * phi * self.s0.ln()).exp() - }; - - let re = |j: u8, tau: f64| { - move |phi: f64| -> f64 { - (f(j, phi, tau) * (-Complex64::i() * phi * self.k.ln()).exp() / (Complex64::i() * phi)).re - } - }; - - let p = |j: u8, tau: f64| -> f64 { - 0.5 + FRAC_1_PI * double_exponential::integrate(re(j, tau), 0.00001, 50.0, 10e-6).integral - }; - unsafe { let tau = self.tau.as_ref().unwrap(); if tau.v.is_empty() { let tau = tau.x; - let call = - self.s0 * (-self.q * tau).exp() * p(1, tau) - self.k * (-self.r * tau).exp() * p(2, tau); + let call = self.s0 * (-self.q * tau).exp() * self.p(1, tau) + - self.k * (-self.r * tau).exp() * self.p(2, tau); let put = call + self.k * (-self.r * tau).exp() - self.s0 * (-self.q * tau).exp(); + self.prices = Some(ValueOrVec { x: (call, put) }); + self.derivates = Some(ValueOrVec { + v: ManuallyDrop::new(self.derivates(tau)), + }); ValueOrVec { x: (call, put) } } else { let mut prices = Vec::with_capacity(tau.v.len()); + let mut derivatives = Vec::with_capacity(tau.v.len()); for tau in tau.v.iter() { - let call = self.s0 * (-self.q * tau).exp() * p(1, *tau) - - self.k * (-self.r * tau).exp() * p(2, *tau); + let call = self.s0 * (-self.q * tau).exp() * self.p(1, *tau) + - self.k * (-self.r * tau).exp() * self.p(2, *tau); let put = call + self.k * (-self.r * tau).exp() - self.s0 * (-self.q * tau).exp(); prices.push((call, put)); + derivatives.push(self.derivates(*tau)); } + self.prices = Some(ValueOrVec { + v: ManuallyDrop::new(prices.clone()), + }); + + // Flatten the derivatives vector + self.derivates = Some(ValueOrVec { + v: ManuallyDrop::new(derivatives.into_iter().flatten().collect::>()), + }); + ValueOrVec { v: ManuallyDrop::new(prices), } @@ -143,32 +113,150 @@ impl HestonPricer { } } + pub(self) fn u(&self, j: u8) -> f64 { + match j { + 1 => 0.5, + 2 => -0.5, + _ => panic!("Invalid j"), + } + } + + pub(self) fn b(&self, j: u8) -> f64 { + match j { + 1 => self.kappa + self.lambda.unwrap() - self.rho * self.sigma, + 2 => self.kappa + self.lambda.unwrap(), + _ => panic!("Invalid j"), + } + } + + pub(self) fn d(&self, j: u8, phi: f64) -> Complex64 { + ((self.b(j) - self.rho * self.sigma * phi * Complex64::i()).powi(2) + - self.sigma.powi(2) * (2.0 * Complex64::i() * self.u(j) * phi - phi.powi(2))) + .sqrt() + } + + pub(self) fn g(&self, j: u8, phi: f64) -> Complex64 { + (self.b(j) - self.rho * self.sigma * Complex64::i() * phi + self.d(j, phi)) + / (self.b(j) - self.rho * self.sigma * Complex64::i() * phi - self.d(j, phi)) + } + + pub(self) fn C(&self, j: u8, phi: f64, tau: f64) -> Complex64 { + (self.r - self.q) * Complex64::i() * phi * tau + + (self.kappa * self.theta / self.sigma.powi(2)) + * ((self.b(j) - self.rho * self.sigma * Complex64::i() * phi + self.d(j, phi)) * tau + - 2.0 + * ((1.0 - self.g(j, phi) * (self.d(j, phi) * tau).exp()) / (1.0 - self.g(j, phi))).ln()) + } + + pub(self) fn D(&self, j: u8, phi: f64, tau: f64) -> Complex64 { + ((self.b(j) - self.rho * self.sigma * Complex64::i() * phi + self.d(j, phi)) + / self.sigma.powi(2)) + * ((1.0 - (self.d(j, phi) * tau).exp()) + / (1.0 - self.g(j, phi) * (self.d(j, phi) * tau).exp())) + } + + pub(self) fn f(&self, j: u8, phi: f64, tau: f64) -> Complex64 { + (self.C(j, phi, tau) + self.D(j, phi, tau) * self.v0 + Complex64::i() * phi * self.s0.ln()) + .exp() + } + + pub(self) fn re(&self, j: u8, tau: f64) -> impl Fn(f64) -> f64 { + let self_ = self.clone(); + move |phi: f64| -> f64 { + (self_.f(j, phi, tau) * (-Complex64::i() * phi * self_.k.ln()).exp() / (Complex64::i() * phi)) + .re + } + } + + pub(self) fn p(&self, j: u8, tau: f64) -> f64 { + 0.5 + FRAC_1_PI * double_exponential::integrate(self.re(j, tau), 0.00001, 50.0, 10e-6).integral + } + /// Partial derivative of the C function with respect to parameters /// https://www.sciencedirect.com/science/article/abs/pii/S0377221717304460 /// Partial derivative of the C function with respect to the v0 parameter - pub(crate) fn dC_dv0(&self) -> f64 { - 1.0 + pub(crate) fn dC_dv0(&self, tau: f64) -> f64 { + (-self.A(tau) / self.v0).re } /// Partial derivative of the C function with respect to the theta parameter - pub(crate) fn dC_dtheta(&self) -> f64 { - 1.0 + pub(crate) fn dC_dtheta(&self, tau: f64) -> f64 { + ((2.0 * self.kappa / self.sigma.powi(2)) * self.D_(tau) + - self.kappa * self.rho * tau * Complex64::i() * self.u(1) / self.sigma) + .re } /// Partial derivative of the C function with respect to the rho parameter - pub(crate) fn dC_drho(&self) -> f64 { - 1.0 + pub(crate) fn dC_drho(&self, tau: f64) -> f64 { + (-self.kappa * self.theta * tau * Complex64::i() * self.u(1) / self.sigma).re } /// Partial derivative of the C function with respect to the kappa parameter - pub(crate) fn dC_dkappa(&self) -> f64 { - 1.0 + pub(crate) fn dC_dkappa(&self, tau: f64) -> f64 { + (2.0 * self.theta * self.D_(tau) / self.sigma.powi(2) + + ((2.0 * self.kappa * self.theta) / self.sigma.powi(2) * self.B(tau)) * self.dB_dkappa(tau) + - (self.theta * self.rho * tau * Complex64::i() * self.u(1) / self.sigma)) + .re } /// Partial derivative of the C function with respect to the sigma parameter - pub(crate) fn dC_dsigma(&self) -> f64 { - 1.0 + pub(crate) fn dC_dsigma(&self, tau: f64) -> f64 { + ((-4.0 * self.kappa * self.theta / self.sigma.powi(3)) * self.D_(tau) + + ((2.0 * self.kappa * self.theta) / (self.sigma.powi(2) * self.d_())) * self.dd_dsigma() + + self.kappa * self.theta * self.rho * tau * Complex64::i() * self.u(1) / self.sigma.powi(2)) + .re + } + + pub(self) fn xi(&self) -> Complex64 { + self.kappa - self.sigma * self.rho * Complex64::i() * self.u(1) + } + + pub(self) fn d_(&self) -> Complex64 { + (self.xi().powi(2) + self.sigma.powi(2) * (self.u(1).powi(2) + Complex64::i() * self.u(1))) + .sqrt() + } + + pub(self) fn dd_dsigma(&self) -> Complex64 { + (self.sigma * (self.u(1) + Complex64::i() * self.u(1))) / self.d_() + } + + pub(self) fn A1(&self, tau: f64) -> Complex64 { + (self.u(1).powi(2) + Complex64::i() * self.u(1)) * (self.d_() * tau / 2.0).sinh() + } + + pub(self) fn A2(&self, tau: f64) -> Complex64 { + (self.d_() / self.v0) * (self.d_() * tau / 2.0).cosh() + + (self.xi() / self.v0) * (self.d_() * tau / 2.0).sinh() + } + + pub(self) fn A(&self, tau: f64) -> Complex64 { + self.A1(tau) / self.A2(tau) + } + + pub(self) fn D_(&self, tau: f64) -> Complex64 { + (self.d_() / self.v0).ln() + (self.kappa - self.d_() / 2.0) * tau + - (((self.d_() + self.xi()) / (2.0 * self.v0)) + + ((self.d_() - self.xi()) / (2.0 * self.v0)) * (-self.d_() * tau).exp()) + .ln() + } + + pub(self) fn B(&self, tau: f64) -> Complex64 { + (self.d_() * (self.kappa * tau / 2.0).exp()) / (self.v0 * self.A2(tau)) + } + + pub(self) fn dB_dkappa(&self, tau: f64) -> Complex64 { + (self.d_() * tau * (self.kappa * tau / 2.0).exp()) / (2.0 * self.v0 * self.A2(tau)) + } + + pub(crate) fn derivates(&self, tau: f64) -> Vec { + vec![ + self.dC_dv0(tau), + self.dC_dtheta(tau), + self.dC_drho(tau), + self.dC_dkappa(tau), + self.dC_dsigma(tau), + ] } } @@ -178,7 +266,7 @@ mod tests { #[test] fn test_price_single_tau() { - let heston = HestonPricer { + let mut heston = HestonPricer { s0: 100.0, v0: 0.05, k: 100.0, @@ -192,6 +280,8 @@ mod tests { tau: Some(ValueOrVec { x: 0.5 }), // Single f64 tau value eval: None, expiry: None, + prices: None, + derivates: None, }; let price = heston.price(); @@ -207,7 +297,7 @@ mod tests { #[test] fn test_price_vec_tau() { - let heston = HestonPricer { + let mut heston = HestonPricer { s0: 100.0, v0: 0.04, k: 100.0, @@ -223,6 +313,8 @@ mod tests { }), // Vec tau eval: None, expiry: None, + prices: None, + derivates: None, }; let price = heston.price();