feat: add sabr and duffie-lan

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "stochastic-rs"
-version = "0.5.1"
+version = "0.5.2"
 edition = "2021"
 license = "MIT"
 description = "A Rust library for stochastic processes"

README.md

@@ -37,6 +37,8 @@ Documentation is available at [stochastic-rs](https://docs.rs/stochastic-rs/).
 - [x] Merton model
 - [x] Bates model
 - [x] Vasicek model
+- [x] SABR model (unstable)
+- [x] Duffie-Kan model (unstable)
 
 
 # Fractional Stochastic processes
@@ -52,7 +54,6 @@ Documentation is available at [stochastic-rs](https://docs.rs/stochastic-rs/).
 - [ ] Rough Heston model
 - [ ] Bergomi model
 - [ ] Rough Bergomi model
-- [ ] SABR model
 - [ ] Hull-White model
 - [ ] Barndorff-Nielsen & Shephard model
 - [ ] Alpha-stable models
@@ -63,7 +64,6 @@ Documentation is available at [stochastic-rs](https://docs.rs/stochastic-rs/).
 - [ ] Wu-Zhang model
 - [ ] Affine model
 - [ ] Heath-Jarrow-Morton model & Multi-factor Heath-Jarrow-Morton model
-- [ ] Duffie-Kan model
 
 ## Future work
 - [x] Add more tests

src/models/duffie_kan.rs

@@ -0,0 +1,76 @@
+use ndarray::Array1;
+
+use crate::processes::correlated::correlated_bms;
+
+/// Generates paths for the Duffie-Kan multifactor interest rate model.
+///
+/// The Duffie-Kan model is a multifactor interest rate model incorporating correlated Brownian motions,
+/// used in financial mathematics for modeling interest rates.
+///
+/// # Parameters
+///
+/// - `alpha`: The drift term coefficient for the Brownian motion.
+/// - `beta`: The drift term coefficient for the Brownian motion.
+/// - `gamma`: The drift term coefficient for the Brownian motion.
+/// - `rho`: The correlation between the two Brownian motions.
+/// - `a1`: The coefficient for the `r` term in the drift of `r`.
+/// - `b1`: The coefficient for the `x` term in the drift of `r`.
+/// - `c1`: The constant term in the drift of `r`.
+/// - `sigma1`: The diffusion coefficient for the `r` process.
+/// - `a2`: The coefficient for the `r` term in the drift of `x`.
+/// - `b2`: The coefficient for the `x` term in the drift of `x`.
+/// - `c2`: The constant term in the drift of `x`.
+/// - `sigma2`: The diffusion coefficient for the `x` process.
+/// - `n`: Number of time steps.
+/// - `r0`: Initial value of the `r` process (optional, defaults to 0.0).
+/// - `x0`: Initial value of the `x` process (optional, defaults to 0.0).
+/// - `t`: Total time (optional, defaults to 1.0).
+///
+/// # Returns
+///
+/// A tuple of two `Vec<f64>` representing the generated paths for the `r` and `x` processes.
+///
+/// # Example
+///
+/// ```
+/// let (r_path, x_path) = duffie_kan(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1000, Some(0.05), Some(0.02), Some(1.0));
+/// ```
+#[allow(clippy::many_single_char_names)]
+pub fn duffie_kan(
+  alpha: f64,
+  beta: f64,
+  gamma: f64,
+  rho: f64,
+  a1: f64,
+  b1: f64,
+  c1: f64,
+  sigma1: f64,
+  a2: f64,
+  b2: f64,
+  c2: f64,
+  sigma2: f64,
+  n: usize,
+  r0: Option<f64>,
+  x0: Option<f64>,
+  t: Option<f64>,
+) -> (Vec<f64>, Vec<f64>) {
+  let correlated_bms = correlated_bms(rho, n, t);
+  let dt = t.unwrap_or(1.0) / n as f64;
+
+  let mut r = Array1::<f64>::zeros(n);
+  let mut x = Array1::<f64>::zeros(n);
+
+  r[0] = r0.unwrap_or(0.0);
+  x[0] = x0.unwrap_or(0.0);
+
+  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];
+    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];
+  }
+
+  (r.to_vec(), x.to_vec())
+}

src/models/mod.rs

@@ -2,11 +2,23 @@
 //!
 //! The following diffusion processes are implemented:
 //!
+//! - Duffie-Kan Model
+//!   - The Duffie-Kan model is a multifactor interest rate model incorporating correlated Brownian motions.
+//!   - SDE: dr(t) = (a1 * r(t) + b1 * x(t) + c1) * dt + sigma1 * (alpha * r(t) + beta * x(t) + gamma) * dW_r(t)
+//!   - SDE: dx(t) = (a2 * r(t) + b2 * x(t) + c2) * dt + sigma2 * (alpha * r(t) + beta * x(t) + gamma) * dW_x(t)
+//!   - where Corr(W_r(t), W_x(t)) = rho
+//!
 //! - **Heston Model**
 //!   - A stochastic volatility model used to describe the evolution of the volatility of an underlying asset.
 //!   - SDE: `dS(t) = mu * S(t) * dt + S(t) * sqrt(V(t)) * dW_1(t)`
 //!   - SDE: `dV(t) = kappa * (theta - V(t)) * dt + eta * sqrt(V(t)) * dW_2(t)`
 //!
+//! - SABR (Stochastic Alpha, Beta, Rho) Model
+//! - Widely used in financial mathematics for modeling stochastic volatility.
+//!   - SDE: dF(t) = V(t) * F(t)^beta * dW_F(t)
+//!   - SDE: dV(t) = alpha * V(t) * dW_V(t)
+//!   - where Corr(W_F(t), W_V(t)) = rho
+//!
 //! - **Vasicek Model**
 //!   - An Ornstein-Uhlenbeck process used to model interest rates.
 //!   - SDE: `dX(t) = theta * (mu - X(t)) * dt + sigma * dW(t)`
@@ -17,5 +29,7 @@
 //!
 //! Each process has its own module and functions to generate sample paths.
 
+pub mod duffie_kan;
 pub mod heston;
+pub mod sabr;
 pub mod vasicek;

src/models/sabr.rs

@@ -0,0 +1,83 @@
+use ndarray::Array1;
+
+use crate::prelude::correlated::correlated_bms;
+
+/// Generates a path of the SABR (Stochastic Alpha, Beta, Rho) model.
+///
+/// The SABR model is widely used in financial mathematics for modeling stochastic volatility.
+/// It incorporates correlated Brownian motions to simulate the underlying asset price and volatility.
+///
+/// # Parameters
+///
+/// - `alpha`: The volatility of volatility.
+/// - `beta`: The elasticity parameter, must be in the range (0, 1).
+/// - `rho`: The correlation between the asset price and volatility, must be in the range (-1, 1).
+/// - `n`: Number of time steps.
+/// - `f0`: Initial value of the forward rate (optional, defaults to 0.0).
+/// - `v0`: Initial value of the volatility (optional, defaults to 0.0).
+/// - `t`: Total time (optional, defaults to 1.0).
+///
+/// # Returns
+///
+/// A tuple of two `Vec<f64>` representing the generated paths for the forward rate and volatility.
+///
+/// # Example
+///
+/// ```
+/// let (forward_rate_path, volatility_path) = sabr(0.2, 0.5, -0.3, 1000, Some(0.04), Some(0.2), Some(1.0));
+/// ```
+pub fn sabr(
+  alpha: f64,
+  beta: f64,
+  rho: f64,
+  n: usize,
+  f0: Option<f64>,
+  v0: Option<f64>,
+  t: Option<f64>,
+) -> (Vec<f64>, Vec<f64>) {
+  if !(0.0..1.0).contains(&beta) {
+    panic!("Beta parameter must be in (0, 1)")
+  }
+
+  if !(-1.0..1.0).contains(&rho) {
+    panic!("Rho parameter must be in (-1, 1)")
+  }
+
+  if alpha < 0.0 {
+    panic!("Alpha parameter must be positive")
+  }
+
+  let correlated_bms = correlated_bms(rho, n, t);
+
+  let mut f = Array1::<f64>::zeros(n);
+  let mut v = Array1::<f64>::zeros(n);
+
+  f[0] = f0.unwrap_or(0.0);
+  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.to_vec(), v.to_vec())
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+
+  #[test]
+  fn test_sabr() {
+    let alpha = 0.2;
+    let beta = 0.5;
+    let rho = -0.3;
+    let n = 1000;
+    let f0 = Some(0.04);
+    let v0 = Some(0.2);
+    let t = Some(1.0);
+    let (f, v) = sabr(alpha, beta, rho, n, f0, v0, t);
+    assert_eq!(f.len(), n);
+    assert_eq!(v.len(), n);
+  }
+}

src/noises/fgn.rs

@@ -69,9 +69,9 @@ impl FgnFft {
     )
     .unwrap();
     let data = r.mapv(|v| Complex::new(v, 0.0));
-    let mut r_fft = FftHandler::new(r.len());
+    let r_fft = FftHandler::new(r.len());
     let mut sqrt_eigenvalues = Array1::<Complex<f64>>::zeros(r.len());
-    ndfft_par(&data, &mut sqrt_eigenvalues, &mut r_fft, 0);
+    ndfft_par(&data, &mut sqrt_eigenvalues, &r_fft, 0);
     sqrt_eigenvalues.par_mapv_inplace(|x| Complex::new((x.re / (2.0 * n as f64)).sqrt(), x.im));
 
     Self {
@@ -106,9 +106,9 @@ impl Generator for FgnFft {
       ComplexDistribution::new(StandardNormal, StandardNormal),
     );
     let fgn = &self.sqrt_eigenvalues * &rnd;
-    let mut fft_handler = self.fft_handler.clone();
+    let fft_handler = self.fft_handler.clone();
     let mut fgn_fft = self.fft_fgn.clone();
-    ndfft_par(&fgn, &mut fgn_fft, &mut fft_handler, 0);
+    ndfft_par(&fgn, &mut fgn_fft, &fft_handler, 0);
     let fgn = fgn_fft
       .slice(s![1..self.n - self.offset + 1])
       .mapv(|x: Complex<f64>| (x.re * (self.n as f64).powf(-self.hurst)) * self.t.powf(self.hurst));

src/processes/poisson.rs

@@ -1,8 +1,6 @@
-use std::ops::Add;
-
 use ndarray::{Array0, Array1, Axis, Dim};
+use ndarray_rand::rand_distr::Normal;
 use ndarray_rand::rand_distr::{Distribution, Exp};
-use ndarray_rand::rand_distr::{Normal, Poisson};
 use ndarray_rand::RandomExt;
 use rand::thread_rng;
 
@@ -132,7 +130,7 @@ mod tests {
     let n = 1000;
     let lambda = 2.0;
     let t = 10.0;
-    let cp = compound_poisson(n, lambda, None, Some(t), None, None);
+    let (.., cp) = compound_poisson(n, lambda, None, Some(t), None);
     assert_eq!(cp.len(), n);
   }
 }