feat: fix jumps and gns, cgns

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "stochastic-rs"
-version = "0.6.1"
+version = "0.6.2"
 edition = "2021"
 license = "MIT"
 description = "A Rust library for stochastic processes"
@@ -16,14 +16,14 @@ ndarray = { version = "0.15.6", features = [
     "rayon",
     "matrixmultiply-threading",
 ] }
-num-complex = { version = "0.4.4", features = ["rand"] }
+num-complex = { version = "0.4.6", features = ["rand"] }
 rand = "0.8.5"
 rand_distr = "0.4.3"
-rayon = "1.8.0"
+rayon = "1.10.0"
 indicatif = "0.17.7"
 plotly = "0.8.4"
 ndarray-rand = "0.14.0"
-ndrustfft = "0.4.2"
+ndrustfft = "0.4.5"
 
 [dev-dependencies]
 

README.md

@@ -15,6 +15,7 @@ Documentation is available at [stochastic-rs](https://docs.rs/stochastic-rs/).
 
 # Stochastic processes
 - [x] Gaussian noise
+- [x] Correlated Gaussian noise
 - [x] Brownian motion
 - [x] Correlated Brownian motion
 - [x] Geometric Brownian motion
@@ -43,6 +44,7 @@ Documentation is available at [stochastic-rs](https://docs.rs/stochastic-rs/).
 
 # Fractional Stochastic processes
 - [x] Fractional Gaussian noise
+- [x] Correlated Gaussian noise
 - [x] Fractional Brownian motion
 - [x] Correlated Fractional Brownian motion
 - [x] Fractional Geometric Brownian motion

src/jumps/bates.rs

@@ -1,23 +1,28 @@
 use ndarray::Array1;
 use rand_distr::Distribution;
 
-use crate::prelude::{
-  cbms::{correlated_bms, CorrelatedBms},
-  cpoisson::{compound_poisson, CompoundPoisson},
+use crate::{
+  noises::cgns::{cgns, Cgns},
+  prelude::cpoisson::{compound_poisson, CompoundPoisson},
 };
 
 /// Generates paths for the Bates (1996) model.
 ///
-/// The Bates model combines a stochastic volatility model with jump diffusion, commonly used in financial mathematics to model asset prices.
+/// The Bates model combines a stochastic volatility model with jump diffusion,
+/// commonly used in financial mathematics to model asset prices.
 ///
 /// # Parameters
 ///
 /// - `mu`: Drift parameter of the asset price.
-/// - `kappa`: Rate of mean reversion of the volatility.
-/// - `theta`: Long-term mean level of the volatility.
-/// - `eta`: Volatility of the volatility (vol of vol).
-/// - `rho`: Correlation between the asset price and its volatility.
+/// - `b`: The continuously compounded domestic/foreign interest rate differential.
+/// - `r`: The continuously compounded risk-free interest rate.
+/// - `r_f`: The continuously compounded foreign interest rate.
 /// - `lambda`: Jump intensity.
+/// - `k`: Mean jump size.
+/// - `alpha`: Rate of mean reversion of the volatility.
+/// - `beta`: Long-term mean level of the volatility.
+/// - `sigma`: Volatility of the volatility (vol of vol).
+/// - `rho`: Correlation between the asset price and its volatility.
 /// - `n`: Number of time steps.
 /// - `s0`: Initial value of the asset price (optional, defaults to 0.0).
 /// - `v0`: Initial value of the volatility (optional, defaults to 0.0).
@@ -31,7 +36,23 @@ use crate::prelude::{
 /// # Example
 ///
 /// ```
-/// let paths = bates_1996(0.05, 1.5, 0.04, 0.3, -0.7, 0.1, 1000, Some(100.0), Some(0.04), Some(1.0), Some(false));
+/// let params = Bates1996 {
+///     mu: 0.05,
+///     lambda: 0.1,
+///     k: 0.2,
+///     alpha: 1.5,
+///     beta: 0.04,
+///     sigma: 0.3,
+///     rho: -0.7,
+///     n: 1000,
+///     s0: Some(100.0),
+///     v0: Some(0.04),
+///     t: Some(1.0),
+///     use_sym: Some(false),
+/// };
+///
+/// let jump_distr = Normal::new(0.0, 1.0); // Example jump distribution
+/// let paths = bates_1996(&params, jump_distr);
 /// let asset_prices = paths[0];
 /// let volatilities = paths[1];
 /// ```
@@ -40,14 +61,18 @@ use crate::prelude::{
 ///
 /// This function will panic if the `correlated_bms` or `compound_poisson` functions return invalid lengths or if there are issues with array indexing.
 
-#[derive(Default)]
+#[derive(Default, Debug)]
 pub struct Bates1996 {
-  pub mu: f64,
-  pub kappa: f64,
-  pub theta: f64,
-  pub eta: f64,
-  pub rho: f64,
+  pub mu: Option<f64>,
+  pub b: Option<f64>,
+  pub r: Option<f64>,
+  pub r_f: Option<f64>,
   pub lambda: f64,
+  pub k: f64,
+  pub alpha: f64,
+  pub beta: f64,
+  pub sigma: f64,
+  pub rho: f64,
   pub n: usize,
   pub s0: Option<f64>,
   pub v0: Option<f64>,
@@ -61,19 +86,23 @@ pub fn bates_1996(
 ) -> [Array1<f64>; 2] {
   let Bates1996 {
     mu,
-    kappa,
-    theta,
-    eta,
-    rho,
+    b,
+    r,
+    r_f,
     lambda,
+    k,
+    alpha,
+    beta,
+    sigma,
+    rho,
     n,
     s0,
     v0,
     t,
     use_sym,
   } = *params;
 
-  let correlated_bms = correlated_bms(&CorrelatedBms { rho, n, t });
+  let [cgn1, cgn2] = cgns(&Cgns { rho, n, t });
   let dt = t.unwrap_or(1.0) / n as f64;
 
   let mut s = Array1::<f64>::zeros(n);
@@ -82,6 +111,12 @@ pub fn bates_1996(
   s[0] = s0.unwrap_or(0.0);
   v[0] = v0.unwrap_or(0.0);
 
+  let drift = match (mu, b, r, r_f) {
+    (Some(r), Some(r_f), ..) => r - r_f,
+    (Some(b), ..) => b,
+    _ => mu.unwrap(),
+  };
+
   for i in 1..n {
     let [.., jumps] = compound_poisson(
       &CompoundPoisson {
@@ -92,17 +127,18 @@ pub fn bates_1996(
       jump_distr,
     );
 
+    let sqrt_v = use_sym
+      .unwrap_or(false)
+      .then(|| v[i - 1].abs())
+      .unwrap_or(v[i - 1].max(0.0))
+      .sqrt();
+
     s[i] = s[i - 1]
-      + mu * s[i - 1] * dt
-      + s[i - 1] * v[i - 1].sqrt() * correlated_bms[0][i - 1]
+      + (drift - lambda * k) * s[i - 1] * dt
+      + s[i - 1] * sqrt_v * cgn1[i - 1]
       + jumps.sum();
 
-    let random: f64 = match use_sym.unwrap_or(false) {
-      true => eta * (v[i]).abs().sqrt() * correlated_bms[1][i - 1],
-      false => eta * (v[i]).max(0.0).sqrt() * correlated_bms[1][i - 1],
-    };
-
-    v[i] = v[i - 1] + kappa * (theta - v[i - 1]) * dt + random;
+    v[i] = v[i - 1] + (alpha - beta * v[i - 1]) * dt + sigma * v[i - 1] * cgn2[i - 1];
   }
 
   [s, v]

src/main.rs

@@ -1,30 +1,45 @@
 use std::time::Instant;
 
 use plotly::{common::Line, Plot, Scatter};
-use rand_distr::{Gamma, Normal};
-use stochastic_rs::{
-  jumps::{bates::bates_1996, jump_fou::jump_fou, levy_diffusion::levy_diffusion, merton::merton},
-  processes::{cpoisson::compound_poisson, fbm::Fbm, poisson::poisson},
-  utils::Generator,
-};
+use rand_distr::Normal;
+use stochastic_rs::jumps::bates::{bates_1996, Bates1996};
 
 fn main() {
   let start = Instant::now();
 
   let mut plot = Plot::new();
-  let fbm = Fbm::new(0.9, 1000, Some(1.0), Some(10));
 
-  for i in 0..1 {
-    // let d = poisson(10.0, Some(50), None);
-    // let trace = Scatter::new((0..d.len()).collect::<Vec<_>>(), d.to_vec())
-    //   .mode(plotly::common::Mode::Lines)
-    //   .line(
-    //     Line::new()
-    //       .color("orange")
-    //       .shape(plotly::common::LineShape::Hv),
-    //   )
-    //   .name("Poisson");
-    //plot.add_trace(trace);
+  let jump_distr = Normal::new(0.0, 0.07).unwrap();
+  let params = Bates1996 {
+    mu: Some(0.05),
+    b: Some(0.05),
+    r: None,
+    r_f: None,
+    lambda: 0.01,
+    k: 0.,
+    alpha: 0.0225,
+    beta: 4.0,
+    sigma: 0.1,
+    rho: 0.1,
+    n: 100,
+    s0: Some(40.),
+    v0: Some(0.025),
+    t: Some(25.0),
+    use_sym: Some(false),
+  };
+  println!("{:?}", params);
+  for _ in 0..10 {
+    let [s, _v] = bates_1996(&params, jump_distr);
+
+    let trace = Scatter::new((0..s.len()).collect::<Vec<_>>(), s.to_vec())
+      .mode(plotly::common::Mode::Lines)
+      .line(
+        Line::new()
+          .color("orange")
+          .shape(plotly::common::LineShape::Hv),
+      )
+      .name("Bates");
+    plot.add_trace(trace);
   }
 
   plot.show();

src/models/duffie_kan.rs

@@ -1,6 +1,6 @@
 use ndarray::Array1;
 
-use crate::processes::cbms::{correlated_bms, CorrelatedBms};
+use crate::noises::cgns::{cgns, Cgns};
 
 /// Generates paths for the Duffie-Kan multifactor interest rate model.
 ///
@@ -76,7 +76,7 @@ pub fn duffie_kan(params: &DuffieKan) -> [Array1<f64>; 2] {
     t,
   } = *params;
 
-  let correlated_bms = correlated_bms(&CorrelatedBms { rho, n, t });
+  let [cgn1, cgn2] = cgns(&Cgns { rho, n, t });
   let dt = t.unwrap_or(1.0) / n as f64;
 
   let mut r = Array1::<f64>::zeros(n);
@@ -88,10 +88,10 @@ pub fn duffie_kan(params: &DuffieKan) -> [Array1<f64>; 2] {
   for i in 1..n {
     r[i] = r[i - 1]
       + (a1 * r[i - 1] + b1 * x[i - 1] + c1) * dt
-      + sigma1 * (alpha * r[i - 1] + beta * x[i - 1] + gamma) * correlated_bms[0][i - 1];
+      + sigma1 * (alpha * r[i - 1] + beta * x[i - 1] + gamma) * cgn1[i - 1];
     x[i] = x[i - 1]
       + (a2 * r[i - 1] + b2 * x[i - 1] + c2) * dt
-      + sigma2 * (alpha * r[i - 1] + beta * x[i - 1] + gamma) * correlated_bms[1][i - 1];
+      + sigma2 * (alpha * r[i - 1] + beta * x[i - 1] + gamma) * cgn2[i - 1];
   }
 
   [r, x]

src/models/heston.rs

@@ -1,6 +1,6 @@
 use ndarray::Array1;
 
-use crate::prelude::cbms::{correlated_bms, CorrelatedBms};
+use crate::noises::cgns::{cgns, Cgns};
 
 /// Generates paths for the Heston model.
 ///
@@ -59,7 +59,7 @@ pub fn heston(params: &Heston) -> [Array1<f64>; 2] {
     use_sym,
   } = *params;
 
-  let correlated_bms = correlated_bms(&CorrelatedBms { rho, n, t });
+  let [cgn1, cgn2] = cgns(&Cgns { rho, n, t });
   let dt = t.unwrap_or(1.0) / n as f64;
 
   let mut s = Array1::<f64>::zeros(n);
@@ -69,11 +69,11 @@ pub fn heston(params: &Heston) -> [Array1<f64>; 2] {
   v[0] = v0.unwrap_or(0.0);
 
   for i in 1..n {
-    s[i] = s[i - 1] + mu * s[i - 1] * dt + s[i - 1] * v[i - 1].sqrt() * correlated_bms[0][i - 1];
+    s[i] = s[i - 1] + mu * s[i - 1] * dt + s[i - 1] * v[i - 1].sqrt() * cgn1[i - 1];
 
     let random: f64 = match use_sym.unwrap_or(false) {
-      true => eta * (v[i - 1]).abs().sqrt() * correlated_bms[1][i - 1],
-      false => eta * (v[i - 1]).max(0.0).sqrt() * correlated_bms[1][i - 1],
+      true => eta * (v[i - 1]).abs().sqrt() * cgn2[i - 1],
+      false => eta * (v[i - 1]).max(0.0).sqrt() * cgn2[i - 1],
     };
     v[i] = v[i - 1] + kappa * (theta - v[i - 1]) * dt + random;
   }

src/models/sabr.rs

@@ -1,6 +1,6 @@
 use ndarray::Array1;
 
-use crate::prelude::cbms::{correlated_bms, CorrelatedBms};
+use crate::noises::cgns::{cgns, Cgns};
 
 /// Generates a path of the SABR (Stochastic Alpha, Beta, Rho) model.
 ///
@@ -53,7 +53,7 @@ pub fn sabr(params: &Sabr) -> [Array1<f64>; 2] {
   assert!(-1.0 < rho && rho < 1.0, "Rho parameter must be in (-1, 1)");
   assert!(alpha > 0.0, "Alpha parameter must be positive");
 
-  let correlated_bms = correlated_bms(&CorrelatedBms { rho, n, t });
+  let [cgn1, cgn2] = cgns(&Cgns { rho, n, t });
 
   let mut f = Array1::<f64>::zeros(n);
   let mut v = Array1::<f64>::zeros(n);
@@ -62,8 +62,8 @@ pub fn sabr(params: &Sabr) -> [Array1<f64>; 2] {
   v[0] = v0.unwrap_or(0.0);
 
   for i in 0..n {
-    f[i] = f[i - 1] + v[i - 1] * f[i - 1].powf(beta) * correlated_bms[0][i - 1];
-    v[i] = v[i - 1] + alpha * v[i - 1] * correlated_bms[1][i - 1];
+    f[i] = f[i - 1] + v[i - 1] * f[i - 1].powf(beta) * cgn1[i - 1];
+    v[i] = v[i - 1] + alpha * v[i - 1] * cgn2[i - 1];
   }
 
   [f, v]