feat: improve malliavin derivate

src/stochastic/malliavin.rs

@@ -4,43 +4,49 @@ use super::Sampling;
 
 pub trait Malliavin {
   /// Compute the Malliavin derivative of the stochastic process using perturbation.
+  ///
+  /// - `f`: the function that maps the path to a scalar value.
+  /// - `epsilon`: the perturbation value.
+  ///
+  /// The Malliavin derivative is defined as the derivative of the function `f` with respect to the path.
   fn malliavin_derivate<F>(&self, f: F, epsilon: f64) -> Array1<f64>
   where
     F: Fn(&Array1<f64>) -> f64,
   {
-    let path = self.path();
+    let mut path = self.path();
     let mut derivates = Array1::zeros(path.len() + 1);
-
     let f_original = f(&path);
 
     for i in 1..self.path().len() {
-      let mut perturbed_path = path.clone();
-      perturbed_path[i] += epsilon;
-
-      let f_perturbed = f(&perturbed_path);
-
+      let original_value = path[i];
+      path[i] += epsilon;
+      let f_perturbed = f(&path);
       derivates[i] = (f_perturbed - f_original) / epsilon;
+      path[i] = original_value;
     }
 
     derivates
   }
 
   /// Compute the Malliavin derivative of the stochastic process using perturbation for latest value.
+  ///
+  ///  - `epsilon`: the perturbation value.
+  ///
+  /// The Malliavin derivative is defined as the derivative of the function `f` with respect to the path.
   /// For example we want to know how the option price changes if the stock price changes.
   fn malliavin_derivate_latest(&self, epsilon: f64) -> Array1<f64> {
-    let path = self.path();
+    let mut path = self.path();
     let mut derivates = Array1::zeros(path.len() + 1);
 
     let final_value = |path: &Array1<f64>| -> f64 { *path.last().unwrap() };
     let f_original = final_value(&path);
 
     for i in 1..path.len() {
-      let mut perturbed_path = path.clone();
-      perturbed_path[i] += epsilon; // Kis perturbáció a folyamat adott pontján
-
-      let f_perturbed = final_value(&perturbed_path);
-
+      let original_value = path[i];
+      path[i] += epsilon;
+      let f_perturbed = final_value(&path);
       derivates[i] = (f_perturbed - f_original) / epsilon;
+      path[i] = original_value;
     }
 
     derivates