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<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);
   }
@@ -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::<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();
   }
 }
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<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 {
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<ValueOrVec<chrono::NaiveDate>>,
   /// Expiration date
   pub expiry: Option<ValueOrVec<chrono::NaiveDate>>,
+  /// Prices of European call and put options
+  pub(crate) prices: Option<ValueOrVec<(f64, f64)>>,
+  /// Partial derivative of the C function with respect to the parameters
+  pub(crate) derivates: Option<ValueOrVec<f64>>,
 }
 
 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::<Vec<f64>>()),
+        });
+
         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<f64> {
+    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<f64> tau
       eval: None,
       expiry: None,
+      prices: None,
+      derivates: None,
     };
 
     let price = heston.price();