fix: heston calibration derivates

Author: dancixx

src/quant/calibration/heston.rs

@@ -222,7 +222,7 @@ mod tests {
       30.75, 25.88, 21.00, 16.50, 11.88, 7.69, 4.44, 2.10, 0.78, 0.25, 0.10, 0.10,
     ];
 
-    let v0 = Array1::linspace(0.0, 0.01, 1);
+    let v0 = Array1::linspace(0.0, 1.0, 10);
 
     for v in v0.iter() {
       let calibrator = HestonCalibrator::new(
@@ -249,6 +249,58 @@ mod tests {
     Ok(())
   }
 
+  #[test]
+  fn test_heston_calibrate_v2() -> Result<()> {
+    let tau = 24.0 / 365.0;
+    println!("Time to maturity: {}", tau);
+
+    let s = vec![
+      8700.53, 8700.53, 8700.53, 8700.53, 8700.53, 8700.53, 8700.53, 8700.53, 8700.53, 8700.53,
+      8700.53, 8700.53,
+    ];
+
+    let k = vec![
+      5220.318, 6090.371, 6960.424, 7830.477, 8265.5035, 8483.01675, 8700.53, 8918.04325,
+      9135.5565, 9570.583, 10440.636, 11310.689,
+    ];
+
+    let c_market = vec![
+      3499.564, 2632.74, 1765.933, 901.846, 479.055, 279.313, 118.848, 28.79, 4.23, 0.143,
+      1.799e-5, 1.259e-9,
+    ];
+
+    let r = 0.04;
+    let q = 0.06;
+    let t = 0.083;
+
+    let v0 = Array1::linspace(1e-5, 2.0, 5);
+
+    for v in v0.iter() {
+      let calibrator = HestonCalibrator::new(
+        Some(HestonParams {
+          v0: *v,
+          theta: 6.47e-5,
+          rho: -1.98e-3,
+          kappa: 6.57e-3,
+          sigma: 5.09e-4,
+        }),
+        c_market.clone().into(),
+        s.clone().into(),
+        k.clone().into(),
+        t,
+        r,
+        Some(q),
+        OptionType::Call,
+      );
+
+      let data = calibrator.calibrate()?;
+      let heston_params = data.iter().map(|d| d.params.clone()).collect::<Vec<_>>();
+      println!("Calibration data: {:?}", heston_params);
+    }
+
+    Ok(())
+  }
+
   #[test]
   fn test_heston_calibrate_guess_params() -> Result<()> {
     let tau = 24.0 / 365.0;
@@ -283,7 +335,8 @@ mod tests {
       calibrator.set_initial_params(s.clone().into(), Array1::from_elem(s.len(), *v), 6.40e-4);
 
       let data = calibrator.calibrate()?;
-      println!("Calibration data: {:?}", data);
+      let heston_params = data.iter().map(|d| d.params.clone()).collect::<Vec<_>>();
+      println!("Calibration data: {:?}", heston_params);
     }
 
     Ok(())

src/quant/pricing/heston.rs

@@ -2,8 +2,10 @@ use std::f64::consts::FRAC_1_PI;
 
 use impl_new_derive::ImplNew;
 use implied_vol::implied_black_volatility;
+use nalgebra::Complex;
 use num_complex::Complex64;
 use quadrature::double_exponential;
+use rand_distr::num_traits::ConstOne;
 
 use crate::quant::{
   r#trait::{PricerExt, TimeExt},
@@ -58,11 +60,11 @@ impl PricerExt for HestonPricer {
     let tau = self.tau().unwrap_or(1.0);
 
     vec![
-      self.dC_dv0(tau),
-      self.dC_dtheta(tau),
-      self.dC_drho(tau),
-      self.dC_dkappa(tau),
-      self.dC_dsigma(tau),
+      self.h1(tau),
+      self.h2(tau),
+      self.h3(tau),
+      self.h4(tau),
+      self.h5(tau),
     ]
   }
 
@@ -154,38 +156,108 @@ impl HestonPricer {
   /// 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, tau: f64) -> f64 {
-    (-self.A(tau) / self.v0).re
+  pub(crate) fn h1(&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, tau: f64) -> f64 {
+  pub(crate) fn h2(&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, tau: f64) -> f64 {
-    (-self.kappa * self.theta * tau * Complex64::i() * self.u(1) / self.sigma).re
+  pub(crate) fn h3(&self, tau: f64) -> f64 {
+    (-self.dA_drho(tau)
+      + ((2.0 * self.kappa * self.theta) / (self.sigma.powi(2) * self.d_()))
+        * (self.dd__drho() - (self.d_() / self.A2(tau)) * self.dA2_drho(tau))
+      - (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, 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))
+  pub(crate) fn h4(&self, tau: f64) -> f64 {
+    ((1.0 / (self.sigma * Complex64::i() * self.u(1))) * self.dA_drho(tau)
+      + ((2.0 * self.theta) / self.sigma.powi(2)) * self.D_(tau)
+      + ((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, 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))
+  pub(crate) fn h5(&self, tau: f64) -> f64 {
+    (-self.dA_dsigma(tau) - ((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.d_() / self.A2(tau)) * self.dA2_dsigma(tau))
+      + (self.kappa * self.theta * self.rho * tau * Complex64::i() * self.u(1))
+        / self.sigma.powi(2))
     .re
   }
 
+  // helpers derivates
+  pub(self) fn dxi_drho(&self) -> Complex64 {
+    -self.sigma * Complex64::i() * self.u(1)
+  }
+
+  pub(self) fn dd__drho(&self) -> Complex64 {
+    -(self.xi() * self.sigma * Complex64::i() * self.u(1)) / self.d_()
+  }
+
+  pub(self) fn dA1_drho(&self, tau: f64) -> Complex64 {
+    -((Complex64::i()
+      * self.u(1)
+      * (self.u(1).powi(2) + Complex64::i() * self.u(1))
+      * tau
+      * self.xi()
+      * self.sigma)
+      / (2.0 * self.d_())
+      * (self.d_() * tau / 2.0).cosh())
+  }
+
+  pub(self) fn dA2_drho(&self, tau: f64) -> Complex64 {
+    ((self.dxi_drho() * (2.0 + self.xi() * tau)) / (2.0 * self.d_() * self.v0))
+      * (self.xi() * (self.d_() * tau / 2.0).cosh() + self.d_() * (self.d_() * tau / 2.0).sinh())
+  }
+
+  pub(self) fn dA_drho(&self, tau: f64) -> Complex64 {
+    (1.0 / self.A2(tau)) * self.dA1_drho(tau) - (self.A(tau) / self.A2(tau)) * self.dA2_drho(tau)
+  }
+
+  pub(self) fn dd__dsigma(&self) -> Complex64 {
+    (self.rho / self.sigma - 1.0 / self.xi()) * self.dd__drho()
+      + (self.sigma * self.u(1).powi(2) / self.d_())
+  }
+
+  pub(self) fn dA1_dsigma(&self, tau: f64) -> Complex64 {
+    (((self.u(1).powi(2) + Complex64::i() * self.u(1)) * tau) / 2.0)
+      * self.dd__dsigma()
+      * (self.d_() * tau / 2.0).cosh()
+  }
+
+  pub(self) fn dA2_dsigma(&self, tau: f64) -> Complex64 {
+    (self.rho / self.sigma) * self.dA2_drho(tau)
+      - ((2.0 + tau * self.xi()) / (self.v0 * tau * self.xi() * Complex64::i() * self.u(1)))
+        * self.dA1_drho(tau)
+      + (self.sigma * tau * self.A1(tau)) / (2.0 * self.v0)
+  }
+
+  pub(self) fn dA_dsigma(&self, tau: f64) -> Complex64 {
+    (1.0 / self.A2(tau)) * self.dA1_dsigma(tau)
+      - (self.A(tau) / self.A2(tau)) * self.dA2_dsigma(tau)
+  }
+
+  pub(self) fn dB_drho(&self, tau: f64) -> Complex64 {
+    ((self.kappa * tau / 2.0).exp() / self.v0)
+      * ((1.0 / self.A2(tau)) * self.dd__drho()
+        - (self.d_() / self.A2(tau).powi(2)) * self.dA2_drho(tau))
+  }
+
+  pub(self) fn dB_dkappa(&self, tau: f64) -> Complex64 {
+    (Complex64::i() / (self.sigma * self.u(1))) * self.dB_drho(tau) + (self.B(tau) * tau) / 2.0
+  }
+
   pub(self) fn xi(&self) -> Complex64 {
     self.kappa - self.sigma * self.rho * Complex64::i() * self.u(1)
   }
@@ -195,10 +267,6 @@ impl HestonPricer {
       .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()
   }
@@ -212,19 +280,12 @@ impl HestonPricer {
     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(self) fn D_(&self, tau: f64) -> Complex64 {
+    self.B(tau).ln()
   }
 }