feat: stabilze jumps and fix compound poisson

Cargo.toml

@@ -22,7 +22,6 @@ indicatif = "0.17.7"
 plotly = "0.8.4"
 ndarray-rand = "0.14.0"
 ndrustfft = "0.4.2"
-csv = "1.3.0"
 
 [dev-dependencies]
 

README.md

@@ -21,11 +21,21 @@ Minimal documentation is available at [stochastic-rs](https://docs.rs/stochastic
 - [x] Cox-Ingersoll-Ross process
 - [x] Ornstein-Uhlenbeck process
 - [x] Jacobi process
-- [x] Variance Gamma model
+
+# Jumps and Levy processes
+- [x] Poisson process
+- [x] Compound Poisson process
+- [x] Levy jump diffusion
+- [x] Inverse Gaussian
+- [x] Normal Inverse Gaussian
+- [x] Variance Gamma
 
 
 # Stochastic models
 - [x] Heston model
+- [x] Merton model
+- [x] Bates model
+
 
 # Fractional Stochastic processes
 - [x] Fractional Gaussian noise
@@ -40,8 +50,6 @@ Minimal documentation is available at [stochastic-rs](https://docs.rs/stochastic
 - [ ] Rough Heston model
 - [ ] Bergomi model 
 - [ ] Rough Bergomi model
-- [ ] Merton model
-- [ ] Bates model
 - [ ] SABR model
 - [ ] Hull-White model
 - [ ] Barndorff-Nielsen & Shephard model
@@ -57,7 +65,6 @@ Minimal documentation is available at [stochastic-rs](https://docs.rs/stochastic
 - [ ] Vasicek model
 - [ ] Duffie-Kan model
 - [ ] Jump-diffusion model
-- [ ] Normal Inverse Gaussian model
 
 ## Future work
 - [ ] Add more tests

src/jumps/bates.rs

@@ -3,15 +3,14 @@ use ndarray::Array1;
 use crate::prelude::{correlated::correlated_bms, poisson::compound_poisson};
 
 /// Bates (1996) model
-// TODO: under development
 #[allow(clippy::too_many_arguments)]
 pub fn bates_1996(
   mu: f64,
   kappa: f64,
   theta: f64,
   eta: f64,
   rho: f64,
-  lambda: usize,
+  lambda: f64,
   n: usize,
   s0: Option<f64>,
   v0: Option<f64>,
@@ -23,23 +22,17 @@ pub fn bates_1996(
 
   let mut s = Array1::<f64>::zeros(n);
   let mut v = Array1::<f64>::zeros(n);
-
-  let [jump_times, jumps] = compound_poisson(n, lambda, None, Some(t.unwrap_or(1.0)), None, None);
+  let z = compound_poisson(n, lambda, None, t, None, None);
+  println!("{:?}", z.len());
 
   s[0] = s0.unwrap_or(0.0);
   v[0] = v0.unwrap_or(0.0);
 
   for i in 1..n {
-    // find the index of the last jump before time t_i
-    let mut j = 0;
-    while jump_times[j] < i as f64 * dt {
-      j += 1;
-    }
-
     s[i] = s[i - 1]
       + mu * s[i - 1] * dt
       + s[i - 1] * v[i - 1].sqrt() * correlated_bms[0][i - 1]
-      + jumps[j];
+      + z[i - 1];
 
     let random: f64 = match use_sym.unwrap_or(false) {
       true => eta * (v[i]).abs().sqrt() * correlated_bms[1][i - 1],
@@ -50,3 +43,35 @@ pub fn bates_1996(
 
   [s.to_vec(), v.to_vec()]
 }
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+
+  #[test]
+  fn test_bates_1996() {
+    let mu = 0.0;
+    let kappa = 0.0;
+    let theta = 0.0;
+    let eta = 0.0;
+    let rho = 0.0;
+    let lambda = 1.0;
+    let n = 1000;
+    let t = 10.0;
+    let b = bates_1996(
+      mu,
+      kappa,
+      theta,
+      eta,
+      rho,
+      lambda,
+      n,
+      None,
+      None,
+      Some(t),
+      None,
+    );
+    assert_eq!(b[0].len(), n);
+    assert_eq!(b[1].len(), n);
+  }
+}

src/jumps/ig.rs

@@ -0,0 +1,40 @@
+use crate::noises::gn::gn;
+
+use ndarray::Array1;
+use ndarray_rand::{rand_distr::InverseGaussian, RandomExt};
+
+pub fn ig(gamma: f64, n: usize, x0: Option<f64>, t: Option<f64>) -> Vec<f64> {
+  let dt = t.unwrap_or(1.0) / n as f64;
+  let gn = gn(n - 1, t);
+  let mut ig = Array1::zeros(n);
+  ig[0] = x0.unwrap_or(0.0);
+
+  for i in 1..n {
+    ig[i] = ig[i - 1] + gamma * dt + gn[i - 1]
+  }
+
+  ig.to_vec()
+}
+
+pub fn nig(
+  theta: f64,
+  sigma: f64,
+  kappa: f64,
+  n: usize,
+  x0: Option<f64>,
+  t: Option<f64>,
+) -> Vec<f64> {
+  let dt = t.unwrap_or(1.0) / n as f64;
+  let scale = dt.powf(2.0) / kappa;
+  let mean = dt / scale;
+  let ig = Array1::random(n - 1, InverseGaussian::new(mean, scale).unwrap());
+  let gn = gn(n - 1, t);
+  let mut nig = Array1::zeros(n);
+  nig[0] = x0.unwrap_or(0.0);
+
+  for i in 1..n {
+    nig[i] = nig[i - 1] + theta * ig[i - 1] + sigma * ig[i - 1].sqrt() * gn[i - 1]
+  }
+
+  nig.to_vec()
+}

src/jumps/levy_diffusion.rs

@@ -0,0 +1,39 @@
+use crate::{noises::gn::gn, processes::poisson::compound_poisson};
+use ndarray::Array1;
+
+pub fn levy_diffusion(
+  gamma: f64,
+  sigma: f64,
+  lambda: f64,
+  n: usize,
+  x0: Option<f64>,
+  t: Option<f64>,
+) -> Vec<f64> {
+  let dt = t.unwrap_or(1.0) / n as f64;
+  let mut levy = Array1::<f64>::zeros(n);
+  levy[0] = x0.unwrap_or(0.0);
+  let gn = gn(n - 1, t);
+  let z = compound_poisson(n, lambda, None, t, None, None);
+
+  for i in 1..n {
+    levy[i] = levy[i - 1] + gamma * dt + sigma * gn[i - 1] * z[i - 1];
+  }
+
+  levy.to_vec()
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+
+  #[test]
+  fn test_levy_diffusion() {
+    let gamma = 0.0;
+    let sigma = 1.0;
+    let lambda = 10.0;
+    let n = 1000;
+    let t = 10.0;
+    let l = levy_diffusion(gamma, sigma, lambda, n, None, Some(t));
+    assert_eq!(l.len(), n);
+  }
+}

src/jumps/merton.rs

@@ -0,0 +1,45 @@
+use ndarray::Array1;
+
+use crate::{noises::gn::gn, processes::poisson::compound_poisson};
+
+pub fn merton(
+  alpha: f64,
+  sigma: f64,
+  lambda: f64,
+  theta: f64,
+  n: usize,
+  x0: Option<f64>,
+  t: Option<f64>,
+) -> Vec<f64> {
+  let dt = t.unwrap_or(1.0) / n as f64;
+  let mut merton = Array1::<f64>::zeros(n);
+  merton[0] = x0.unwrap_or(0.0);
+  let gn = gn(n - 1, t);
+  let z = compound_poisson(n, lambda, None, t, None, None);
+
+  for i in 1..n {
+    merton[i] = merton[i - 1]
+      + (alpha * sigma.powf(2.0) / 2.0 - lambda * theta) * dt
+      + sigma * gn[i - 1]
+      + z[i - 1];
+  }
+
+  merton.to_vec()
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+
+  #[test]
+  fn test_merton() {
+    let alpha = 0.0;
+    let sigma = 1.0;
+    let lambda = 10.0;
+    let theta = 0.0;
+    let n = 1000;
+    let t = 10.0;
+    let m = merton(alpha, sigma, lambda, theta, n, None, Some(t));
+    assert_eq!(m.len(), n);
+  }
+}

src/jumps/mod.rs

@@ -1,2 +1,5 @@
 pub mod bates;
+pub mod ig;
+pub mod levy_diffusion;
+pub mod merton;
 pub mod vg;

src/jumps/vg.rs

@@ -1,23 +1,19 @@
 use crate::noises::gn;
 use ndarray::Array1;
-use rand::Rng;
+use ndarray_rand::RandomExt;
 use rand_distr::Gamma;
 
 pub fn vg(mu: f64, sigma: f64, nu: f64, n: usize, x0: Option<f64>, t: Option<f64>) -> Vec<f64> {
   let dt = t.unwrap_or(1.0) / n as f64;
-  let rng = rand::thread_rng();
 
   let shape = dt / nu;
   let scale = nu;
 
   let mut vg = Array1::<f64>::zeros(n);
   vg[0] = x0.unwrap_or(0.0);
 
-  let gn = gn::gn(n - 1, Some(t.unwrap_or(1.0)));
-  let gammas = rng
-    .sample_iter(Gamma::new(shape, scale).unwrap())
-    .take(n - 1)
-    .collect::<Vec<f64>>();
+  let gn = gn::gn(n - 1, t);
+  let gammas = Array1::random(n - 1, Gamma::new(shape, scale).unwrap());
 
   for i in 1..n {
     vg[i] = vg[i - 1] + mu * gammas[i - 1] + sigma * gammas[i - 1].sqrt() * gn[i - 1];

src/processes/poisson.rs

@@ -1,6 +1,8 @@
 use ndarray::Array1;
-use rand::{thread_rng, Rng};
-use rand_distr::{Distribution, Exp, Normal};
+use ndarray_rand::rand_distr::{Normal, Poisson};
+use ndarray_rand::RandomExt;
+use rand::thread_rng;
+use rand_distr::{Distribution, Exp};
 
 pub fn poisson(n: usize, lambda: usize, t_max: Option<f64>) -> Vec<f64> {
   if n == 0 || lambda == 0 {
@@ -25,34 +27,36 @@ pub fn poisson(n: usize, lambda: usize, t_max: Option<f64>) -> Vec<f64> {
 
 pub fn compound_poisson(
   n: usize,
-  lambda: usize,
+  lambda: f64,
   jumps: Option<Vec<f64>>,
   t_max: Option<f64>,
   jump_mean: Option<f64>,
   jump_std: Option<f64>,
-) -> [Vec<f64>; 2] {
+) -> Vec<f64> {
   if n == 0 {
     panic!("n must be a positive integer");
   }
 
   let _t_max = t_max.unwrap_or(1.0);
-  let times = poisson(n, lambda, Some(_t_max));
-  let mut cp = Array1::<f64>::zeros(times.len());
-  let jumps = match jumps {
+  let p = Array1::random(n, Poisson::new(lambda).unwrap());
+  let mut cp = Array1::<f64>::zeros(n);
+
+  match jumps {
     Some(jumps) => jumps,
     None => {
-      let _jump_mean = jump_mean.unwrap_or(1.0);
-      let _jump_std = jump_std.unwrap_or(0.0);
-      let norm = Normal::new(_jump_mean, _jump_std).unwrap();
-      thread_rng().sample_iter(norm).take(n).collect()
-    }
-  };
+      let _jump_mean = jump_mean.unwrap_or(0.0);
+      let _jump_std = jump_std.unwrap_or(1.0);
 
-  for i in 1..times.len() {
-    cp[i] = cp[i - 1] + jumps[i - 1];
-  }
+      for i in 0..n {
+        for j in &p {
+          let norm = Array1::random(*j as usize, Normal::new(_jump_mean, _jump_std).unwrap());
+          cp[i] = norm.sum();
+        }
+      }
 
-  [times, cp.to_vec()]
+      cp.to_vec()
+    }
+  }
 }
 
 #[cfg(test)]
@@ -71,10 +75,9 @@ mod tests {
   #[test]
   fn test_compound_poisson() {
     let n = 1000;
-    let lambda = 10;
+    let lambda = 2.0;
     let t = 10.0;
-    let cp = compound_poisson(n, lambda, None, Some(t), Some(1.0), Some(0.0));
-    println!("{:?}", cp[0].len());
-    println!("{:?}", cp[1].len());
+    let cp = compound_poisson(n, lambda, None, Some(t), None, None);
+    assert_eq!(cp.len(), n);
   }
 }