feat: upgrade heston calibrator

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 }
 

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,81 +31,95 @@ impl<'a> LeastSquaresProblem<f64, Dyn, Dyn> for Calibrator<'a> {
   }
 
   fn residuals(&self) -> Option<DVector<f64>> {
-    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::<Vec<f64>>()
+    };
+    // 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::<Vec<f64>>();
 
-    let residuals = model_prices
+    let residuals = call_prices
       .iter()
       .zip(market_prices.iter())
-      .map(|(m, p)| m - p)
+      .map(|(x, y)| x - y)
       .collect::<Vec<f64>>();
 
     Some(DVector::from_vec(residuals))
   }
 
   fn jacobian(&self) -> Option<DMatrix<f64>> {
-    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);
   }
 }
 
 #[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::<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();
   }
 }

stochastic-rs-quant/src/lib.rs

@@ -31,6 +31,16 @@ impl Clone for ValueOrVec<f64> {
   }
 }
 
+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<T>`.
 impl Clone for ValueOrVec<chrono::NaiveDate> {
   fn clone(&self) -> Self {