fix: malliavin derivative calculations

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "stochastic-rs"
-version = "0.10.2"
+version = "0.10.3"
 edition = "2021"
 license = "MIT"
 description = "A Rust library for quant finance and simulating stochastic processes."

src/stochastic.rs

@@ -68,16 +68,6 @@ pub trait Sampling<T: Clone + Send + Sync + Zero>: Send + Sync {
   fn malliavin(&self) -> Array1<T> {
     unimplemented!()
   }
-
-  #[cfg(feature = "malliavin")]
-  fn malliavin_sensitivity(&self) -> Array1<T> {
-    unimplemented!()
-  }
-
-  #[cfg(feature = "malliavin")]
-  fn malliavin_matrix(&self) -> Array2<T> {
-    unimplemented!()
-  }
 }
 
 pub trait Sampling2D<T: Clone + Send + Sync + Zero>: Send + Sync {
@@ -116,16 +106,6 @@ pub trait Sampling2D<T: Clone + Send + Sync + Zero>: Send + Sync {
   fn malliavin(&self) -> [Array1<T>; 2] {
     unimplemented!()
   }
-
-  #[cfg(feature = "malliavin")]
-  fn malliavin_sensitivity(&self) -> Array1<T> {
-    unimplemented!()
-  }
-
-  #[cfg(feature = "malliavin")]
-  fn malliavin_matrix(&self) -> Array2<T> {
-    unimplemented!()
-  }
 }
 
 pub trait Sampling3D<T: Clone + Send + Sync + Zero>: Send + Sync {

src/stochastic/diffusion/cev.rs

@@ -20,8 +20,6 @@ pub struct CEV {
   pub calculate_malliavin: Option<bool>,
   #[cfg(feature = "malliavin")]
   malliavin: Mutex<Option<Array1<f64>>>,
-  #[cfg(feature = "malliavin")]
-  malliavin_sensitivity: Mutex<Option<Array1<f64>>>,
 }
 
 impl Sampling<f64> for CEV {
@@ -44,7 +42,6 @@ impl Sampling<f64> for CEV {
       let mut det_term = Array1::zeros(self.n);
       let mut stochastic_term = Array1::zeros(self.n);
       let mut malliavin = Array1::zeros(self.n);
-      let mut malliavin_sensitivity = Array1::zeros(self.n);
 
       for i in 0..self.n {
         det_term[i] = (self.mu
@@ -53,18 +50,13 @@ impl Sampling<f64> for CEV {
         if i > 0 {
           stochastic_term[i] = self.sigma * self.gamma * cev[i].powf(self.gamma - 1.0) * gn[i - 1];
         }
-        malliavin[i] =
-          self.sigma * cev[i].powf(self.gamma) * (det_term[i] + stochastic_term[i]).exp();
-        if i > 0 {
-          malliavin_sensitivity[i] = malliavin[i] * gn[i - 1];
-        }
+        malliavin[i] = self.sigma
+          * cev[i].powf(self.gamma)
+          * (det_term[i] + stochastic_term[i]).exp()
+          * (i as f64 * dt);
       }
 
       let _ = std::mem::replace(&mut *self.malliavin.lock().unwrap(), Some(malliavin));
-      let _ = std::mem::replace(
-        &mut *self.malliavin_sensitivity.lock().unwrap(),
-        Some(malliavin_sensitivity),
-      );
     }
 
     cev
@@ -90,17 +82,12 @@ impl Sampling<f64> for CEV {
   fn malliavin(&self) -> Array1<f64> {
     self.malliavin.lock().unwrap().clone().unwrap()
   }
-
-  #[cfg(feature = "malliavin")]
-  fn malliavin_sensitivity(&self) -> Array1<f64> {
-    self.malliavin_sensitivity.lock().unwrap().clone().unwrap()
-  }
 }
 
 #[cfg(test)]
 mod tests {
   use crate::{
-    plot_1d, plot_3d,
+    plot_1d, plot_2d,
     stochastic::{N, X0},
   };
 
@@ -133,14 +120,11 @@ mod tests {
     let cev = CEV::new(0.25, 0.5, 0.3, N, Some(X0), Some(1.0), None, Some(true));
     let process = cev.sample();
     let malliavin = cev.malliavin();
-    let malliavin_sensitivity = cev.malliavin_sensitivity();
-    plot_3d!(
+    plot_2d!(
       process,
       "Constant Elasticity of Variance (CEV) process",
       malliavin,
-      "Malliavin derivative of the CEV process",
-      malliavin_sensitivity,
-      "Malliavin sensitivity of the CEV process"
+      "Malliavin derivative of the CEV process"
     );
   }
 }

src/stochastic/diffusion/gbm.rs

@@ -26,10 +26,6 @@ pub struct GBM {
   pub calculate_malliavin: Option<bool>,
   #[cfg(feature = "malliavin")]
   malliavin: Mutex<Option<Array1<f64>>>,
-  #[cfg(feature = "malliavin")]
-  malliavin_sensitivity: Mutex<Option<Array1<f64>>>,
-  #[cfg(feature = "malliavin")]
-  malliavin_matrix: Mutex<Option<Array2<f64>>>,
 }
 
 impl Sampling<f64> for GBM {
@@ -48,33 +44,15 @@ impl Sampling<f64> for GBM {
     #[cfg(feature = "malliavin")]
     if self.calculate_malliavin.is_some() && self.calculate_malliavin.unwrap() {
       let mut malliavin = Array1::zeros(self.n);
-      let mut malliavin_sensitivity = Array1::zeros(self.n);
-      let mut malliavin_matrix = Array2::zeros((self.n, self.n));
-      for i in 0..self.n {
-        malliavin[i] = self.sigma * gbm[i];
-        if i > 0 {
-          malliavin_sensitivity[i] = malliavin[i] * gn[i - 1];
-        }
-      }
 
-      // TODO: this is in very early stage
-      for i in 1..self.n {
-        for j in 1..self.n {
-          malliavin_matrix[(i, j)] = if i < j { 0.0 } else { malliavin[j] * gn[i - 1] }
-        }
+      // reverse due the option pricing
+      for i in 0..self.n {
+        malliavin[i] = self.sigma * gbm.last().unwrap() * (i as f64 * dt);
       }
 
       // This equivalent to the following:
       // self.malliavin.lock().unwrap().replace(Some(malliavin));
       let _ = std::mem::replace(&mut *self.malliavin.lock().unwrap(), Some(malliavin));
-      let _ = std::mem::replace(
-        &mut *self.malliavin_sensitivity.lock().unwrap(),
-        Some(malliavin_sensitivity),
-      );
-      let _ = std::mem::replace(
-        &mut *self.malliavin_matrix.lock().unwrap(),
-        Some(malliavin_matrix),
-      );
     }
 
     gbm
@@ -111,17 +89,6 @@ impl Sampling<f64> for GBM {
   fn malliavin(&self) -> Array1<f64> {
     self.malliavin.lock().unwrap().as_ref().unwrap().clone()
   }
-
-  #[cfg(feature = "malliavin")]
-  fn malliavin_sensitivity(&self) -> Array1<f64> {
-    self.malliavin_sensitivity.lock().unwrap().clone().unwrap()
-  }
-
-  /// Calculate the Malliavin derivative matrix
-  #[cfg(feature = "malliavin")]
-  fn malliavin_matrix(&self) -> Array2<f64> {
-    self.malliavin_matrix.lock().unwrap().clone().unwrap()
-  }
 }
 
 impl Distribution for GBM {
@@ -213,11 +180,8 @@ impl Distribution for GBM {
 
 #[cfg(test)]
 mod tests {
-  use ndarray::Array;
-  use plotly::{Layout, Plot, Surface};
-
   use crate::{
-    plot_1d, plot_3d,
+    plot_1d, plot_2d,
     stochastic::{N, X0},
   };
 
@@ -247,49 +211,11 @@ mod tests {
     let gbm = GBM::new(0.25, 0.5, N, Some(X0), Some(1.0), None, None, Some(true));
     let process = gbm.sample();
     let malliavin = gbm.malliavin();
-    let malliavin_sensitivity = gbm.malliavin_sensitivity();
-    plot_3d!(
+    plot_2d!(
       process,
       "Geometric Brownian Motion (GBM) process",
       malliavin,
-      "Malliavin derivative of the GBM process",
-      malliavin_sensitivity,
-      "Malliavin sensitivity of the GBM process"
+      "Malliavin derivative of the GBM process"
     );
   }
-
-  #[test]
-  #[cfg(feature = "malliavin")]
-  fn gbm_plot_malliavin_matrix() {
-    let gbm = GBM::new(0.25, 0.5, N, Some(X0), Some(1.0), None, None, Some(true));
-    let _ = gbm.sample();
-    let malliavin_matrix: ndarray::ArrayBase<ndarray::OwnedRepr<f64>, ndarray::Dim<[usize; 2]>> =
-      gbm.malliavin_matrix();
-
-    let x = Array::linspace(-1., 1., N).into_raw_vec_and_offset().0;
-    let y = Array::linspace(-1., 1., N).into_raw_vec_and_offset().0;
-
-    let z = malliavin_matrix
-      .outer_iter() // Iterálás a sorokon
-      .map(|row| row.to_vec()) // A sorokat átalakítjuk vektorrá
-      .collect();
-
-    let trace = Surface::new(z)
-      .x(x.clone())
-      .y(y.clone())
-      .connect_gaps(true)
-      .cauto(true);
-    let mut plot = Plot::new();
-    plot.add_trace(trace);
-
-    plot.set_layout(
-      Layout::new()
-        .title("Malliavin derivative matrix of the GBM process")
-        .auto_size(false)
-        .height(800)
-        .width(800),
-    );
-
-    plot.show();
-  }
 }

src/stochastic/process/fbm.rs

@@ -36,7 +36,8 @@ impl Sampling<f64> for FBM {
       let mut malliavin = Array1::zeros(self.n);
       let dt = self.t.unwrap_or(1.0) / (self.n) as f64;
       for i in 0..self.n {
-        malliavin[i] = 1.0 / (gamma::gamma(self.hurst + 0.5)) * dt.powf(self.hurst - 0.5);
+        malliavin[i] =
+          1.0 / (gamma::gamma(self.hurst + 0.5)) * (i as f64 * dt).powf(self.hurst - 0.5);
       }
 
       let _ = std::mem::replace(&mut *self.malliavin.lock().unwrap(), Some(malliavin));

src/stochastic/volatility/heston.rs

@@ -45,9 +45,6 @@ pub struct Heston {
   /// Malliavin derivative of the volatility
   #[cfg(feature = "malliavin")]
   malliavin_of_vol: Mutex<Option<Array1<f64>>>,
-  /// Sensitivity of the Malliavin derivative of the volatility
-  #[cfg(feature = "malliavin")]
-  malliavin_sensitivity_of_vol: Mutex<Option<Array1<f64>>>,
   /// Malliavin derivative of the price
   #[cfg(feature = "malliavin")]
   malliavin_of_price: Mutex<Option<Array1<f64>>>,
@@ -85,42 +82,31 @@ impl Sampling2D<f64> for Heston {
     if self.calculate_malliavin.is_some() && self.calculate_malliavin.unwrap() {
       let mut det_term = Array1::zeros(self.n);
       let mut malliavin = Array1::zeros(self.n);
-      let mut malliavin_sensitivity = Array1::zeros(self.n);
 
       for i in 0..self.n {
         match self.pow {
           HestonPow::Sqrt => {
             det_term[i] = ((-(self.kappa * self.theta / 2.0 - self.sigma.powi(2) / 8.0)
-              * (1.0 / v[i - 1])
+              * (1.0 / v.last().unwrap())
               - self.kappa / 2.0)
-              * dt)
+              * ((self.n - i) as f64 * dt))
               .exp();
-            malliavin[i] = (self.sigma * v[i - 1].sqrt() / 2.0) * det_term[i];
-
-            if i > 0 {
-              malliavin_sensitivity[i] = malliavin[i] * cgn2[i - 1];
-            }
+            malliavin[i] =
+              (self.sigma * v.last().unwrap().sqrt() / 2.0) * det_term[i] * (i as f64 * dt);
           }
           HestonPow::ThreeHalves => {
             det_term[i] = ((-(self.kappa * self.theta / 2.0 + 3.0 * self.sigma.powi(2) / 8.0)
-              * v[i - 1]
+              * v.last().unwrap()
               - (self.kappa * self.theta) / 2.0)
-              * dt)
+              * ((self.n - i) as f64 * dt))
               .exp();
-            malliavin[i] = (self.sigma * v[i - 1].powf(1.5) / 2.0) * det_term[i];
-
-            if i > 0 {
-              malliavin_sensitivity[i] = malliavin[i] * cgn2[i - 1];
-            }
+            malliavin[i] =
+              (self.sigma * v.last().unwrap().powf(1.5) / 2.0) * det_term[i] * (i as f64 * dt);
           }
         };
       }
 
       let _ = std::mem::replace(&mut *self.malliavin_of_vol.lock().unwrap(), Some(malliavin));
-      let _ = std::mem::replace(
-        &mut *self.malliavin_sensitivity_of_vol.lock().unwrap(),
-        Some(malliavin_sensitivity),
-      );
     }
 
     [s, v]
@@ -157,3 +143,42 @@ impl Sampling2D<f64> for Heston {
     ]
   }
 }
+
+#[cfg(test)]
+mod tests {
+  use crate::{
+    plot_2d,
+    stochastic::{N, S0, X0},
+  };
+
+  use super::*;
+
+  #[test]
+  #[cfg(feature = "malliavin")]
+  fn heston_malliavin() {
+    let heston = Heston::new(
+      Some(S0),
+      Some(X0),
+      0.5,
+      1.0,
+      1.0,
+      1.0,
+      1.0,
+      N,
+      Some(1.0),
+      HestonPow::Sqrt,
+      None,
+      None,
+      CGNS::new(0.7, N, None, None),
+      Some(true),
+    );
+    let process = heston.sample();
+    let malliavin = heston.malliavin();
+    plot_2d!(
+      process[1],
+      "Heston volatility process",
+      malliavin[1],
+      "Malliavin derivative of the Heston volatility process"
+    );
+  }
+}