feat: add heston calibration

stochastic-rs-quant/Cargo.toml

@@ -7,17 +7,19 @@ edition = "2021"
 chrono = "0.4.38"
 derive_builder = "0.20.1"
 levenberg-marquardt = "0.14.0"
-mathru = "0.15.4"
 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 }
 
 [features]
 mimalloc = ["dep:mimalloc"]
 jemalloc = ["dep:tikv-jemallocator"]
-default = ["jemalloc"]
+default = []
 
 [lib]
 name = "stochastic_rs_quant"

stochastic-rs-quant/src/calibrator.rs

@@ -0,0 +1,107 @@
+use levenberg_marquardt::LeastSquaresProblem;
+use nalgebra::{DMatrix, DVector, Dyn, Owned};
+use rand::thread_rng;
+use rand_distr::{Distribution, Normal};
+
+use crate::pricer::Pricer;
+
+/// A calibrator.
+pub(crate) struct Calibrator<'a, P>
+where
+  P: Pricer,
+{
+  pub p: DVector<f64>,
+  pub market_prices: Option<DVector<f64>>,
+  pricer: &'a P,
+}
+
+impl<'a, P> Calibrator<'a, P>
+where
+  P: Pricer,
+{
+  #[must_use]
+  pub(crate) fn new(p: DVector<f64>, market_prices: Option<DVector<f64>>, pricer: &'a P) -> Self {
+    Self {
+      p,
+      market_prices,
+      pricer,
+    }
+  }
+}
+
+impl<'a, P> Calibrator<'a, P>
+where
+  P: Pricer,
+{
+  fn generate_market_prices(&self) -> DVector<f64> {
+    let model_prices = self.pricer.prices().unwrap();
+    let call_prices = unsafe {
+      model_prices
+        .v
+        .clone()
+        .iter()
+        .map(|x| x.0)
+        .collect::<Vec<f64>>()
+    };
+
+    // Add some noise to the market prices
+    let market_prices = call_prices
+      .iter()
+      .map(|x| *x + Normal::new(20.0, 4.5).unwrap().sample(&mut thread_rng()))
+      .collect::<Vec<f64>>();
+
+    DVector::from_vec(market_prices)
+  }
+}
+
+impl<'a, P> LeastSquaresProblem<f64, Dyn, Dyn> for Calibrator<'a, P>
+where
+  P: Pricer,
+{
+  type JacobianStorage = Owned<f64, Dyn, Dyn>;
+  type ParameterStorage = Owned<f64, Dyn>;
+  type ResidualStorage = Owned<f64, Dyn>;
+
+  fn set_params(&mut self, p: &DVector<f64>) {
+    self.p.copy_from(p);
+  }
+
+  fn params(&self) -> DVector<f64> {
+    self.p.clone()
+  }
+
+  fn residuals(&self) -> Option<DVector<f64>> {
+    let model_prices = self.pricer.prices().unwrap();
+    let call_prices = unsafe {
+      model_prices
+        .v
+        .clone()
+        .iter()
+        .map(|x| x.0)
+        .collect::<Vec<f64>>()
+    };
+
+    let market_prices = match &self.market_prices {
+      Some(x) => x.clone(),
+      None => self.generate_market_prices(),
+    };
+
+    let residuals = call_prices
+      .iter()
+      .zip(market_prices.iter())
+      .map(|(x, y)| x - y)
+      .collect::<Vec<f64>>();
+
+    Some(DVector::from_vec(residuals))
+  }
+
+  fn jacobian(&self) -> Option<DMatrix<f64>> {
+    let derivates = self.pricer.derivates().unwrap();
+    let derivates = unsafe { derivates.v.clone().to_vec() };
+
+    // The Jacobian matrix is a matrix of partial derivatives
+    // of the residuals with respect to the parameters.
+    let jacobian = DMatrix::from_vec(derivates.len() / self.p.len(), 5, derivates);
+    Some(jacobian)
+  }
+}

stochastic-rs-quant/src/heston.rs

@@ -0,0 +1,2 @@
+pub mod calibrator;
+pub mod pricer;

stochastic-rs-quant/src/heston/calibrator.rs

@@ -0,0 +1,63 @@
+use levenberg_marquardt::LevenbergMarquardt;
+use nalgebra::DVector;
+
+use crate::calibrator::Calibrator;
+
+use super::pricer::HestonPricer;
+
+pub struct HestonCalibrator {
+  pricer: HestonPricer,
+}
+
+impl HestonCalibrator {
+  #[must_use]
+  pub fn new(pricer: HestonPricer) -> Self {
+    Self { pricer }
+  }
+
+  pub fn calibrate(&mut self) {
+    self.pricer.price();
+    let (result, ..) = LevenbergMarquardt::new().minimize(Calibrator::new(
+      DVector::from_vec(vec![0.05, 0.05, -0.8, 5.0, 0.5]),
+      None,
+      &self.pricer,
+    ));
+    println!("{:?}", result.p);
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use std::mem::ManuallyDrop;
+
+  use crate::ValueOrVec;
+
+  use super::*;
+
+  #[test]
+  fn test_calibrate() {
+    let majurities = (0..=100)
+      .map(|x| 0.5 + 0.1 * x as f64)
+      .collect::<Vec<f64>>();
+    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();
+  }
+}