diff --git a/src/stochastic/process/fbm.rs b/src/stochastic/process/fbm.rs
index 0d32848..9fec3aa 100644
--- a/src/stochastic/process/fbm.rs
+++ b/src/stochastic/process/fbm.rs
@@ -1,4 +1,7 @@
+use std::sync::Mutex;
+
 use ndarray::{s, Array1};
+use statrs::function::gamma;
 
 use crate::stochastic::{noise::fgn::FGN, Sampling};
 
@@ -9,6 +12,8 @@ pub struct Fbm {
   pub t: Option<f64>,
   pub m: Option<usize>,
   pub fgn: FGN,
+  pub calculate_malliavin: Option<bool>,
+  malliavin: Mutex<Option<Array1<f64>>>,
 }
 
 impl Fbm {
@@ -25,6 +30,8 @@ impl Fbm {
       n: params.n,
       m: params.m,
       fgn,
+      calculate_malliavin: Some(false),
+      malliavin: Mutex::new(None),
     }
   }
 }
@@ -39,6 +46,15 @@ impl Sampling<f64> for Fbm {
       fbm[i] += fbm[i - 1];
     }
 
+    if self.calculate_malliavin.is_some() && self.calculate_malliavin.unwrap() {
+      let mut malliavin = Array1::zeros(self.n + 1);
+      let dt = self.t.unwrap_or(1.0) / self.n as f64;
+      for i in 1..=self.n {
+        malliavin[i] = 1.0 / (gamma::gamma(self.hurst + 0.5)) * dt.powf(self.hurst - 0.5);
+      }
+
+      let _ = std::mem::replace(&mut *self.malliavin.lock().unwrap(), Some(malliavin));
+    }
     fbm.slice(s![..self.n()]).to_owned()
   }
 
@@ -51,6 +67,18 @@ impl Sampling<f64> for Fbm {
   fn m(&self) -> Option<usize> {
     self.m
   }
+
+  /// Calculate the Malliavin derivative
+  ///
+  /// The Malliavin derivative of the fractional Brownian motion is given by:
+  /// D_s B^H_t = 1 / Γ(H + 1/2) (t - s)^{H - 1/2}
+  ///
+  /// where B^H_t is the fractional Brownian motion with Hurst parameter H in Mandelbrot-Van Ness representation as
+  /// B^H_t = 1 / Γ(H + 1/2) ∫_0^t (t - s)^{H - 1/2} dW_s
+  /// which is a truncated Wiener integral.
+  fn malliavin(&self) -> Array1<f64> {
+    self.malliavin.lock().unwrap().clone().unwrap()
+  }
 }
 
 #[cfg(test)]