diff --git a/src/stochastic.rs b/src/stochastic.rs
index 258bc1f..3d6f02c 100644
--- a/src/stochastic.rs
+++ b/src/stochastic.rs
@@ -17,7 +17,10 @@ use rand_distr::Distribution as RandDistribution;
 pub trait ProcessDistribution: RandDistribution<f64> + Copy + Send + Sync + Default {}
 
 pub trait Sampling<T: Clone + Send + Sync + Zero>: Send + Sync {
+  /// Sample the process
   fn sample(&self) -> Array1<T>;
+
+  /// Parallel sampling
   fn sample_par(&self) -> Array2<T> {
     if self.m().is_none() {
       panic!("m must be specified for parallel sampling");
@@ -31,13 +34,27 @@ pub trait Sampling<T: Clone + Send + Sync + Zero>: Send + Sync {
 
     xs
   }
+
+  /// Number of time steps
   fn n(&self) -> usize;
+
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize>;
+
+  /// Distribution of the process
   fn distribution(&mut self) {}
+
+  /// Malliavin derivative of the process
+  fn malliavin(&self) -> Array1<T> {
+    unimplemented!()
+  }
 }
 
 pub trait Sampling2D<T: Clone + Send + Sync + Zero>: Send + Sync {
+  /// Sample the process
   fn sample(&self) -> [Array1<T>; 2];
+
+  /// Parallel sampling
   fn sample_par(&self) -> [Array2<T>; 2] {
     if self.m().is_none() {
       panic!("m must be specified for parallel sampling");
@@ -57,16 +74,32 @@ pub trait Sampling2D<T: Clone + Send + Sync + Zero>: Send + Sync {
     let xs2 = xs2.lock().unwrap().clone();
     [xs1, xs2]
   }
+
+  /// Number of time steps
   fn n(&self) -> usize;
+
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize>;
+
+  /// Malliavin derivative of the process
+  fn malliavin(&self) -> [Array1<T>; 2] {
+    unimplemented!()
+  }
 }
 
 pub trait Sampling3D<T: Clone + Send + Sync + Zero>: Send + Sync {
+  /// Sample the process
   fn sample(&self) -> [Array1<T>; 3];
+
+  /// Parallel sampling
   fn sample_par(&self) -> [Array2<T>; 3] {
     unimplemented!()
   }
+
+  /// Number of time steps
   fn n(&self) -> usize;
+
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize>;
 }
 
diff --git a/src/stochastic/diffusion.rs b/src/stochastic/diffusion.rs
index ff334ca..8e622cf 100644
--- a/src/stochastic/diffusion.rs
+++ b/src/stochastic/diffusion.rs
@@ -1,3 +1,4 @@
+pub mod cev;
 pub mod cir;
 pub mod fcir;
 pub mod fgbm;
diff --git a/src/stochastic/diffusion/cev.rs b/src/stochastic/diffusion/cev.rs
new file mode 100644
index 0000000..ed4f8da
--- /dev/null
+++ b/src/stochastic/diffusion/cev.rs
@@ -0,0 +1,95 @@
+use std::sync::Mutex;
+
+use ndarray::Array1;
+use ndarray_rand::RandomExt;
+use rand_distr::Normal;
+
+use crate::stochastic::Sampling;
+
+#[derive(Default)]
+pub struct CEV {
+  pub mu: f64,
+  pub sigma: f64,
+  pub gamma: f64,
+  pub n: usize,
+  pub x0: Option<f64>,
+  pub t: Option<f64>,
+  pub m: Option<usize>,
+  pub calculate_malliavin: Option<bool>,
+  malliavin: Mutex<Option<Array1<f64>>>,
+}
+
+impl CEV {
+  #[must_use]
+  pub fn new(params: &Self) -> Self {
+    Self {
+      mu: params.mu,
+      sigma: params.sigma,
+      gamma: params.gamma,
+      n: params.n,
+      x0: params.x0,
+      t: params.t,
+      m: params.m,
+      calculate_malliavin: Some(false),
+      malliavin: Mutex::new(None),
+    }
+  }
+}
+
+impl Sampling<f64> for CEV {
+  /// Sample the CEV process
+  fn sample(&self) -> Array1<f64> {
+    let dt = self.t.unwrap_or(1.0) / self.n as f64;
+    let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
+
+    let mut cev = Array1::<f64>::zeros(self.n + 1);
+    cev[0] = self.x0.unwrap_or(0.0);
+
+    for i in 1..=self.n {
+      cev[i] = cev[i - 1]
+        + self.mu * cev[i - 1] * dt
+        + self.sigma * cev[i - 1].powf(self.gamma) * gn[i - 1]
+    }
+
+    if self.calculate_malliavin.is_some() && self.calculate_malliavin.unwrap() {
+      let mut det_term = Array1::zeros(self.n + 1);
+      let mut stochastic_term = Array1::zeros(self.n + 1);
+      let mut malliavin = Array1::zeros(self.n + 1);
+
+      for i in 1..=self.n {
+        det_term[i] = (self.mu
+          - (self.gamma.powi(2) * self.sigma.powi(2) * cev[i - 1].powf(2.0 * self.gamma - 2.0)
+            / 2.0))
+          * dt;
+        stochastic_term[i] =
+          self.sigma * self.gamma * cev[i - 1].powf(self.gamma - 1.0) * gn[i - 1];
+        malliavin[i] =
+          self.sigma * cev[i - 1].powf(self.gamma) * (det_term[i] + stochastic_term[i]).exp();
+      }
+
+      let _ = std::mem::replace(&mut *self.malliavin.lock().unwrap(), Some(malliavin));
+    }
+
+    cev
+  }
+
+  /// Number of time steps
+  fn n(&self) -> usize {
+    self.n
+  }
+
+  /// Number of samples for parallel sampling
+  fn m(&self) -> Option<usize> {
+    self.m
+  }
+
+  /// Calculate the Malliavin derivative of the CEV process
+  ///
+  /// The Malliavin derivative of the CEV process is given by
+  /// D_r S_t = \sigma S_t^{\gamma} * 1_{[0, r]}(r) exp(\int_0^r (\mu - \frac{\gamma^2 \sigma^2 S_u^{2\gamma - 2}}{2}) du + \int_0^r \gamma \sigma S_u^{\gamma - 1} dW_u)
+  ///
+  /// The Malliavin derivative of the CEV process shows the sensitivity of the stock price with respect to the Wiener process.
+  fn malliavin(&self) -> Array1<f64> {
+    self.malliavin.lock().unwrap().clone().unwrap()
+  }
+}
diff --git a/src/stochastic/diffusion/cir.rs b/src/stochastic/diffusion/cir.rs
index 98d9718..d89bfe9 100644
--- a/src/stochastic/diffusion/cir.rs
+++ b/src/stochastic/diffusion/cir.rs
@@ -36,6 +36,7 @@ impl CIR {
 }
 
 impl Sampling<f64> for CIR {
+  /// Sample the Cox-Ingersoll-Ross (CIR) process
   fn sample(&self) -> Array1<f64> {
     assert!(
       2.0 * self.theta * self.mu < self.sigma.powi(2),
@@ -61,10 +62,12 @@ impl Sampling<f64> for CIR {
     cir
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/diffusion/fcir.rs b/src/stochastic/diffusion/fcir.rs
index 612e671..0de2cab 100644
--- a/src/stochastic/diffusion/fcir.rs
+++ b/src/stochastic/diffusion/fcir.rs
@@ -37,6 +37,7 @@ impl FCIR {
 }
 
 impl Sampling<f64> for FCIR {
+  /// Sample the Fractional Cox-Ingersoll-Ross (FCIR) process
   fn sample(&self) -> Array1<f64> {
     assert!(
       2.0 * self.theta * self.mu < self.sigma.powi(2),
@@ -62,10 +63,12 @@ impl Sampling<f64> for FCIR {
     fcir.slice(s![..self.n()]).to_owned()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/diffusion/fgbm.rs b/src/stochastic/diffusion/fgbm.rs
index 1548a2b..68638a4 100644
--- a/src/stochastic/diffusion/fgbm.rs
+++ b/src/stochastic/diffusion/fgbm.rs
@@ -33,6 +33,7 @@ impl FGBM {
 }
 
 impl Sampling<f64> for FGBM {
+  /// Sample the Fractional Geometric Brownian Motion (FGBM) process
   fn sample(&self) -> Array1<f64> {
     assert!(
       self.hurst > 0.0 && self.hurst < 1.0,
@@ -52,10 +53,13 @@ impl Sampling<f64> for FGBM {
     fgbm.slice(s![..self.n()]).to_owned()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of paths
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/diffusion/fjacobi.rs b/src/stochastic/diffusion/fjacobi.rs
index 3fc3e04..c71f1db 100644
--- a/src/stochastic/diffusion/fjacobi.rs
+++ b/src/stochastic/diffusion/fjacobi.rs
@@ -35,6 +35,7 @@ impl FJacobi {
 }
 
 impl Sampling<f64> for FJacobi {
+  /// Sample the Fractional Jacobi process
   fn sample(&self) -> Array1<f64> {
     assert!(
       self.hurst > 0.0 && self.hurst < 1.0,
@@ -66,10 +67,13 @@ impl Sampling<f64> for FJacobi {
     fjacobi.slice(s![..self.n()]).to_owned()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of paths
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/diffusion/fou.rs b/src/stochastic/diffusion/fou.rs
index 63c642d..cd33326 100644
--- a/src/stochastic/diffusion/fou.rs
+++ b/src/stochastic/diffusion/fou.rs
@@ -35,6 +35,7 @@ impl FOU {
 }
 
 impl Sampling<f64> for FOU {
+  /// Sample the Fractional Ornstein-Uhlenbeck (FOU) process
   fn sample(&self) -> Array1<f64> {
     assert!(
       self.hurst > 0.0 && self.hurst < 1.0,
@@ -54,10 +55,13 @@ impl Sampling<f64> for FOU {
     fou.slice(s![..self.n()]).to_owned()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of paths
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/diffusion/gbm.rs b/src/stochastic/diffusion/gbm.rs
index 94f0d33..ee3add7 100644
--- a/src/stochastic/diffusion/gbm.rs
+++ b/src/stochastic/diffusion/gbm.rs
@@ -1,3 +1,5 @@
+use std::sync::Mutex;
+
 use ndarray::Array1;
 use ndarray_rand::RandomExt;
 use num_complex::Complex64;
@@ -18,6 +20,8 @@ pub struct GBM {
   pub t: Option<f64>,
   pub m: Option<usize>,
   pub distribution: Option<LogNormal>,
+  pub calculate_malliavin: Option<bool>,
+  malliavin: Mutex<Option<Array1<f64>>>,
 }
 
 impl GBM {
@@ -31,11 +35,14 @@ impl GBM {
       t: params.t,
       m: params.m,
       distribution: None,
+      calculate_malliavin: Some(false),
+      malliavin: Mutex::new(None),
     }
   }
 }
 
 impl Sampling<f64> for GBM {
+  /// Sample the GBM process
   fn sample(&self) -> Array1<f64> {
     let dt = self.t.unwrap_or(1.0) / self.n as f64;
     let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
@@ -47,17 +54,31 @@ impl Sampling<f64> for GBM {
       gbm[i] = gbm[i - 1] + self.mu * gbm[i - 1] * dt + self.sigma * gbm[i - 1] * gn[i - 1]
     }
 
+    if self.calculate_malliavin.is_some() && self.calculate_malliavin.unwrap() {
+      let mut malliavin = Array1::zeros(self.n + 1);
+      for i in 1..=self.n {
+        malliavin[i] = self.sigma * gbm[i - 1];
+      }
+
+      // This equivalent to the following:
+      // self.malliavin.lock().unwrap().replace(Some(malliavin));
+      let _ = std::mem::replace(&mut *self.malliavin.lock().unwrap(), Some(malliavin));
+    }
+
     gbm
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
 
+  /// Distribution of the GBM process
   fn distribution(&mut self) {
     let mu = self.x0.unwrap() * (self.mu * self.t.unwrap()).exp();
     let sigma = (self.x0.unwrap().powi(2)
@@ -67,6 +88,16 @@ impl Sampling<f64> for GBM {
 
     self.distribution = Some(LogNormal::new(mu, sigma).unwrap());
   }
+
+  /// Mallaivin derivative of the GBM process
+  ///
+  /// The Malliavin derivative of the CEV process is given by
+  /// D_r S_t = \sigma S_t * 1_[0, r](r)
+  ///
+  /// The Malliavin derivate of the GBM shows the sensitivity of the stock price with respect to the Wiener process.
+  fn malliavin(&self) -> Array1<f64> {
+    self.malliavin.lock().unwrap().as_ref().unwrap().clone()
+  }
 }
 
 impl Distribution for GBM {
diff --git a/src/stochastic/diffusion/jacobi.rs b/src/stochastic/diffusion/jacobi.rs
index a23faf6..2fd0b02 100644
--- a/src/stochastic/diffusion/jacobi.rs
+++ b/src/stochastic/diffusion/jacobi.rs
@@ -31,6 +31,7 @@ impl Jacobi {
 }
 
 impl Sampling<f64> for Jacobi {
+  /// Sample the Jacobi process
   fn sample(&self) -> Array1<f64> {
     assert!(self.alpha > 0.0, "alpha must be positive");
     assert!(self.beta > 0.0, "beta must be positive");
@@ -58,10 +59,12 @@ impl Sampling<f64> for Jacobi {
     jacobi
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/diffusion/ou.rs b/src/stochastic/diffusion/ou.rs
index 8ee85b8..1c48b01 100644
--- a/src/stochastic/diffusion/ou.rs
+++ b/src/stochastic/diffusion/ou.rs
@@ -31,6 +31,7 @@ impl OU {
 }
 
 impl Sampling<f64> for OU {
+  /// Sample the Ornstein-Uhlenbeck (OU) process
   fn sample(&self) -> Array1<f64> {
     let dt = self.t.unwrap_or(1.0) / self.n as f64;
     let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
@@ -45,10 +46,12 @@ impl Sampling<f64> for OU {
     ou
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/interest/duffie_kan.rs b/src/stochastic/interest/duffie_kan.rs
index d9202e9..8557a7b 100644
--- a/src/stochastic/interest/duffie_kan.rs
+++ b/src/stochastic/interest/duffie_kan.rs
@@ -59,6 +59,7 @@ impl DuffieKan {
 }
 
 impl Sampling2D<f64> for DuffieKan {
+  /// Sample the Duffie-Kan process
   fn sample(&self) -> [Array1<f64>; 2] {
     let [cgn1, cgn2] = self.cgns.sample();
     let dt = self.t.unwrap_or(1.0) / self.n as f64;
@@ -81,10 +82,12 @@ impl Sampling2D<f64> for DuffieKan {
     [r, x]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/interest/fvasicek.rs b/src/stochastic/interest/fvasicek.rs
index b5f3f58..12dd6de 100644
--- a/src/stochastic/interest/fvasicek.rs
+++ b/src/stochastic/interest/fvasicek.rs
@@ -49,10 +49,12 @@ impl Sampling<f64> for FVasicek {
     self.fou.sample()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/interest/ho_lee.rs b/src/stochastic/interest/ho_lee.rs
index 8838978..34140bf 100644
--- a/src/stochastic/interest/ho_lee.rs
+++ b/src/stochastic/interest/ho_lee.rs
@@ -61,10 +61,12 @@ impl<'a> Sampling<f64> for HoLee<'a> {
     r
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/interest/hull_white.rs b/src/stochastic/interest/hull_white.rs
index c4bc256..1280f1d 100644
--- a/src/stochastic/interest/hull_white.rs
+++ b/src/stochastic/interest/hull_white.rs
@@ -49,10 +49,12 @@ impl Sampling<f64> for HullWhite {
     hw
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/interest/hull_white_2f.rs b/src/stochastic/interest/hull_white_2f.rs
index 1e15250..3bb930c 100644
--- a/src/stochastic/interest/hull_white_2f.rs
+++ b/src/stochastic/interest/hull_white_2f.rs
@@ -66,10 +66,12 @@ impl Sampling2D<f64> for HullWhite2F {
     [x, u]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/interest/vasicek.rs b/src/stochastic/interest/vasicek.rs
index c504392..e79f706 100644
--- a/src/stochastic/interest/vasicek.rs
+++ b/src/stochastic/interest/vasicek.rs
@@ -45,10 +45,12 @@ impl Sampling<f64> for Vasicek {
     self.ou.sample()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/bates.rs b/src/stochastic/jump/bates.rs
index 9611e3e..63e9564 100644
--- a/src/stochastic/jump/bates.rs
+++ b/src/stochastic/jump/bates.rs
@@ -110,10 +110,12 @@ impl<D: ProcessDistribution> Sampling2D<f64> for Bates1996<D> {
     [s, v]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/ig.rs b/src/stochastic/jump/ig.rs
index 4f68fdd..adf8dbb 100644
--- a/src/stochastic/jump/ig.rs
+++ b/src/stochastic/jump/ig.rs
@@ -41,10 +41,12 @@ impl Sampling<f64> for IG {
     ig
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/jump_fou.rs b/src/stochastic/jump/jump_fou.rs
index cf75876..4a70f6f 100644
--- a/src/stochastic/jump/jump_fou.rs
+++ b/src/stochastic/jump/jump_fou.rs
@@ -78,10 +78,12 @@ impl<D: ProcessDistribution> Sampling<f64> for JumpFOU<D> {
     jump_fou.slice(s![..self.n()]).to_owned()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/levy_diffusion.rs b/src/stochastic/jump/levy_diffusion.rs
index 0627185..f56803b 100644
--- a/src/stochastic/jump/levy_diffusion.rs
+++ b/src/stochastic/jump/levy_diffusion.rs
@@ -63,10 +63,12 @@ impl<D: ProcessDistribution> Sampling<f64> for LevyDiffusion<D> {
     levy
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/merton.rs b/src/stochastic/jump/merton.rs
index cd8f7bf..d56e5b3 100644
--- a/src/stochastic/jump/merton.rs
+++ b/src/stochastic/jump/merton.rs
@@ -68,10 +68,12 @@ impl<D: ProcessDistribution> Sampling<f64> for Merton<D> {
     merton
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/nig.rs b/src/stochastic/jump/nig.rs
index 5d6ebef..ed71ce7 100644
--- a/src/stochastic/jump/nig.rs
+++ b/src/stochastic/jump/nig.rs
@@ -48,10 +48,12 @@ impl Sampling<f64> for NIG {
     nig
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/jump/vg.rs b/src/stochastic/jump/vg.rs
index bf1c315..0818c58 100644
--- a/src/stochastic/jump/vg.rs
+++ b/src/stochastic/jump/vg.rs
@@ -51,10 +51,12 @@ impl Sampling<f64> for VG {
     vg
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/noise/cfgns.rs b/src/stochastic/noise/cfgns.rs
index e379ebe..c99ed99 100644
--- a/src/stochastic/noise/cfgns.rs
+++ b/src/stochastic/noise/cfgns.rs
@@ -56,10 +56,12 @@ impl Sampling2D<f64> for CFGNS {
     ]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/noise/cgns.rs b/src/stochastic/noise/cgns.rs
index da0cb95..23003d0 100644
--- a/src/stochastic/noise/cgns.rs
+++ b/src/stochastic/noise/cgns.rs
@@ -47,10 +47,12 @@ impl Sampling2D<f64> for CGNS {
     ]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/noise/fgn.rs b/src/stochastic/noise/fgn.rs
index 40cb488..13a6c6a 100644
--- a/src/stochastic/noise/fgn.rs
+++ b/src/stochastic/noise/fgn.rs
@@ -95,10 +95,12 @@ impl Sampling<f64> for FGN {
     fgn
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n - self.offset
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/bm.rs b/src/stochastic/process/bm.rs
index 36b0f07..d8fd708 100644
--- a/src/stochastic/process/bm.rs
+++ b/src/stochastic/process/bm.rs
@@ -36,10 +36,12 @@ impl Sampling<f64> for BM {
     bm
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/cbms.rs b/src/stochastic/process/cbms.rs
index afac6c3..8a605d8 100644
--- a/src/stochastic/process/cbms.rs
+++ b/src/stochastic/process/cbms.rs
@@ -50,10 +50,12 @@ impl Sampling2D<f64> for CBMS {
     [bms.row(0).into_owned(), bms.row(1).into_owned()]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/ccustom.rs b/src/stochastic/process/ccustom.rs
index e56eb2a..2faf894 100644
--- a/src/stochastic/process/ccustom.rs
+++ b/src/stochastic/process/ccustom.rs
@@ -66,10 +66,12 @@ where
     [p, cum_jupms, jumps]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n.unwrap_or(0)
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/cfbms.rs b/src/stochastic/process/cfbms.rs
index ef31d23..e4e4be8 100644
--- a/src/stochastic/process/cfbms.rs
+++ b/src/stochastic/process/cfbms.rs
@@ -67,10 +67,12 @@ impl Sampling2D<f64> for Cfbms {
     [fbms.row(0).to_owned(), fbms.row(1).to_owned()]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/cpoisson.rs b/src/stochastic/process/cpoisson.rs
index 5bad913..f4c7f20 100644
--- a/src/stochastic/process/cpoisson.rs
+++ b/src/stochastic/process/cpoisson.rs
@@ -57,10 +57,12 @@ impl<D: ProcessDistribution> Sampling3D<f64> for CompoundPoisson<D> {
     [poisson, cum_jupms, jumps]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n.unwrap_or(0)
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/customjt.rs b/src/stochastic/process/customjt.rs
index 48b4a6d..0e16c5d 100644
--- a/src/stochastic/process/customjt.rs
+++ b/src/stochastic/process/customjt.rs
@@ -53,10 +53,12 @@ impl<D: ProcessDistribution> Sampling<f64> for CustomJt<D> {
     }
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n.unwrap_or(0)
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/fbm.rs b/src/stochastic/process/fbm.rs
index b84ffdf..0d32848 100644
--- a/src/stochastic/process/fbm.rs
+++ b/src/stochastic/process/fbm.rs
@@ -42,10 +42,12 @@ impl Sampling<f64> for Fbm {
     fbm.slice(s![..self.n()]).to_owned()
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/process/poisson.rs b/src/stochastic/process/poisson.rs
index 516b980..10bec57 100644
--- a/src/stochastic/process/poisson.rs
+++ b/src/stochastic/process/poisson.rs
@@ -54,10 +54,12 @@ impl Sampling<f64> for Poisson {
     }
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n.unwrap_or(0)
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/volatility/bergomi.rs b/src/stochastic/volatility/bergomi.rs
index d171ef6..d976ceb 100644
--- a/src/stochastic/volatility/bergomi.rs
+++ b/src/stochastic/volatility/bergomi.rs
@@ -61,10 +61,12 @@ impl Sampling2D<f64> for Bergomi {
     [s, v2]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/volatility/fheston.rs b/src/stochastic/volatility/fheston.rs
index 4b3b605..6039bdb 100644
--- a/src/stochastic/volatility/fheston.rs
+++ b/src/stochastic/volatility/fheston.rs
@@ -69,10 +69,12 @@ impl Sampling<f64> for RoughHeston {
     v2
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/volatility/heston.rs b/src/stochastic/volatility/heston.rs
index 01f6ae2..304f301 100644
--- a/src/stochastic/volatility/heston.rs
+++ b/src/stochastic/volatility/heston.rs
@@ -96,10 +96,12 @@ impl Sampling2D<f64> for Heston {
     [s, v]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/volatility/rbergomi.rs b/src/stochastic/volatility/rbergomi.rs
index b127bdd..6251c44 100644
--- a/src/stochastic/volatility/rbergomi.rs
+++ b/src/stochastic/volatility/rbergomi.rs
@@ -65,10 +65,12 @@ impl Sampling2D<f64> for RoughBergomi {
     [s, v2]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
diff --git a/src/stochastic/volatility/sabr.rs b/src/stochastic/volatility/sabr.rs
index 291a646..e35092b 100644
--- a/src/stochastic/volatility/sabr.rs
+++ b/src/stochastic/volatility/sabr.rs
@@ -1,3 +1,5 @@
+use std::sync::Mutex;
+
 use ndarray::Array1;
 
 use crate::stochastic::{noise::cgns::CGNS, Sampling2D};
@@ -14,6 +16,9 @@ pub struct Sabr {
   pub t: Option<f64>,
   pub m: Option<usize>,
   pub cgns: CGNS,
+  pub calculate_malliavin: Option<bool>,
+  malliavin_of_vol: Mutex<Option<Array1<f64>>>,
+  malliavin_of_price: Mutex<Option<Array1<f64>>>,
 }
 
 impl Sabr {
@@ -36,6 +41,9 @@ impl Sabr {
       t: params.t,
       m: params.m,
       cgns,
+      calculate_malliavin: Some(false),
+      malliavin_of_vol: Mutex::new(None),
+      malliavin_of_price: Mutex::new(None),
     }
   }
 }
@@ -55,14 +63,46 @@ impl Sampling2D<f64> for Sabr {
       v[i] = v[i - 1] + self.alpha * v[i - 1] * cgn2[i - 1];
     }
 
+    if self.calculate_malliavin.is_some() && self.calculate_malliavin.unwrap() {
+      // Only volatility Malliavin derivative is supported
+      let mut malliavin_of_vol = Array1::<f64>::zeros(self.n + 1);
+      for i in 1..=self.n {
+        malliavin_of_vol[i] = self.alpha * v[i - 1];
+      }
+
+      let _ = std::mem::replace(
+        &mut *self.malliavin_of_vol.lock().unwrap(),
+        Some(malliavin_of_vol),
+      );
+    }
+
     [f, v]
   }
 
+  /// Number of time steps
   fn n(&self) -> usize {
     self.n
   }
 
+  /// Number of samples for parallel sampling
   fn m(&self) -> Option<usize> {
     self.m
   }
+
+  /// Calculate the Malliavin derivative of the SABR model
+  ///
+  /// The Malliavin derivative of the volaility process in the SABR model is given by:
+  /// D_r \sigma_t = \alpha \sigma_t 1_{[0, T]}(r)
+  fn malliavin(&self) -> [Array1<f64>; 2] {
+    [
+      Array1::zeros(self.n + 1),
+      self
+        .malliavin_of_vol
+        .lock()
+        .unwrap()
+        .as_ref()
+        .unwrap()
+        .clone(),
+    ]
+  }
 }