feat: add gbm, cev, sabr mallaivin derivates

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>;
 }
 

src/stochastic/diffusion.rs

@@ -1,3 +1,4 @@
+pub mod cev;
 pub mod cir;
 pub mod fcir;
 pub mod fgbm;

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()
+  }
+}

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
   }

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
   }

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
   }

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
   }

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
   }

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 {

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
   }