feat: add bergomi model

stochastic-rs-core/src/volatility.rs

@@ -1,3 +1,4 @@
+pub mod bergomi;
 pub mod fheston;
 pub mod heston;
 pub mod rbergomi;

stochastic-rs-core/src/volatility/bergomi.rs

@@ -0,0 +1,71 @@
+use ndarray::{s, Array1};
+
+use crate::{noise::cgns::Cgns, Sampling2D};
+
+#[derive(Default)]
+pub struct Bergomi {
+  pub nu: f64,
+  pub v0: Option<f64>,
+  pub s0: Option<f64>,
+  pub r: f64,
+  pub rho: f64,
+  pub n: usize,
+  pub t: Option<f64>,
+  pub m: Option<usize>,
+  pub cgns: Cgns,
+}
+
+impl Bergomi {
+  #[must_use]
+  pub fn new(params: &Self) -> Self {
+    let cgns = Cgns::new(&Cgns {
+      rho: params.rho,
+      n: params.n,
+      t: params.t,
+      m: params.m,
+    });
+
+    Self {
+      nu: params.nu,
+      v0: params.v0,
+      s0: params.s0,
+      r: params.r,
+      rho: params.rho,
+      n: params.n,
+      t: params.t,
+      m: params.m,
+      cgns,
+    }
+  }
+}
+
+impl Sampling2D<f64> for Bergomi {
+  fn sample(&self) -> [Array1<f64>; 2] {
+    let dt = self.t.unwrap_or(1.0) / self.n as f64;
+    let [cgn1, z] = self.cgns.sample();
+
+    let mut s = Array1::<f64>::zeros(self.n + 1);
+    let mut v2 = Array1::<f64>::zeros(self.n + 1);
+    s[0] = self.s0.unwrap_or(100.0);
+    v2[0] = self.v0.unwrap_or(1.0).powi(2);
+
+    for i in 0..=self.n {
+      s[i] = s[i - 1] + self.r * s[i - 1] * dt + v2[i - 1].sqrt() * s[i - 1] * cgn1[i - 1];
+
+      let sum_z = z.slice(s![..i]).sum();
+      let t = i as f64 * dt;
+      v2[i] =
+        self.v0.unwrap_or(1.0).powi(2) * (self.nu * t * sum_z - 0.5 * self.nu.powi(2) * t.powi(2))
+    }
+
+    [s, v2]
+  }
+
+  fn n(&self) -> usize {
+    self.n
+  }
+
+  fn m(&self) -> Option<usize> {
+    self.m
+  }
+}

stochastic-rs-core/src/volatility/fheston.rs

@@ -52,7 +52,7 @@ impl Sampling<f64> for RoughHeston {
     for i in 1..=self.n {
       let t = dt * i as f64;
       yt[i] = self.theta + (yt[i - 1] - self.theta) * (-self.kappa * dt).exp();
-      zt[i] = zt[i - 1] + (-self.kappa * dt).exp() * gn[i - 1];
+      zt[i] = zt[i - 1] * (-self.kappa * dt).exp() + (v2[i - 1].powi(2)).sqrt() * gn[i - 1];
 
       let integral = (0..i)
         .map(|j| {

stochastic-rs-core/src/volatility/rbergomi.rs

@@ -7,6 +7,7 @@ pub struct RoughBergomi {
   pub hurst: f64,
   pub nu: f64,
   pub v0: Option<f64>,
+  pub s0: Option<f64>,
   pub r: f64,
   pub rho: f64,
   pub n: usize,
@@ -29,6 +30,7 @@ impl RoughBergomi {
       hurst: params.hurst,
       nu: params.nu,
       v0: params.v0,
+      s0: params.s0,
       r: params.r,
       rho: params.rho,
       n: params.n,
@@ -46,10 +48,11 @@ impl Sampling2D<f64> for RoughBergomi {
 
     let mut s = Array1::<f64>::zeros(self.n + 1);
     let mut v2 = Array1::<f64>::zeros(self.n + 1);
+    s[0] = self.s0.unwrap_or(100.0);
     v2[0] = self.v0.unwrap_or(1.0).powi(2);
 
     for i in 1..=self.n {
-      s[i] = s[i - 1] + self.r * s[i - 1] + v2[i - 1].sqrt() * cgn1[i - 1];
+      s[i] = s[i - 1] + self.r * s[i - 1] * dt + v2[i - 1].sqrt() * s[i - 1] * cgn1[i - 1];
 
       let sum_z = z.slice(s![..i]).sum();
       let t = i as f64 * dt;