feat: add rough heston

stochastic-rs-core/Cargo.toml

@@ -25,6 +25,7 @@ ndarray-rand = "0.15.0"
 ndrustfft = "0.5.0"
 tikv-jemallocator = { version = "0.6.0", optional = true }
 mimalloc = { version = "0.1.43", optional = true }
+statrs = "0.17.1"
 
 [dev-dependencies]
 

stochastic-rs-core/src/diffusion/cir.rs

@@ -39,17 +39,13 @@ impl Sampling<f64> for Cir {
       "2 * theta * mu < sigma^2"
     );
 
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
-
     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 cir = Array1::<f64>::zeros(self.n + 1);
     cir[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let random = match self.use_sym.unwrap_or(false) {
         true => self.sigma * (cir[i - 1]).abs().sqrt() * gn[i - 1],
         false => self.sigma * (cir[i - 1]).max(0.0).sqrt() * gn[i - 1],

stochastic-rs-core/src/diffusion/fcir.rs

@@ -49,7 +49,7 @@ impl Sampling<f64> for Fcir {
     let mut fcir = Array1::<f64>::zeros(self.n + 1);
     fcir[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let random = match self.use_sym.unwrap_or(false) {
         true => self.sigma * (fcir[i - 1]).abs().sqrt() * fgn[i - 1],
         false => self.sigma * (fcir[i - 1]).max(0.0) * fgn[i - 1],

stochastic-rs-core/src/diffusion/fgbm.rs

@@ -45,7 +45,7 @@ impl Sampling<f64> for Fgbm {
     let mut fgbm = Array1::<f64>::zeros(self.n + 1);
     fgbm[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       fgbm[i] = fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1]
     }
 

stochastic-rs-core/src/diffusion/fjacobi.rs

@@ -51,7 +51,7 @@ impl Sampling<f64> for Fjacobi {
     let mut fjacobi = Array1::<f64>::zeros(self.n + 1);
     fjacobi[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       fjacobi[i] = match fjacobi[i - 1] {
         _ if fjacobi[i - 1] <= 0.0 && i > 0 => 0.0,
         _ if fjacobi[i - 1] >= 1.0 && i > 0 => 1.0,

stochastic-rs-core/src/diffusion/fou.rs

@@ -47,7 +47,7 @@ impl Sampling<f64> for Fou {
     let mut fou = Array1::<f64>::zeros(self.n + 1);
     fou[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       fou[i] = fou[i - 1] + self.theta * (self.mu - fou[i - 1]) * dt + self.sigma * fgn[i - 1]
     }
 

stochastic-rs-core/src/diffusion/gbm.rs

@@ -30,17 +30,13 @@ impl Gbm {
 
 impl Sampling<f64> for Gbm {
   fn sample(&self) -> Array1<f64> {
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
-
     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 gbm = Array1::<f64>::zeros(self.n + 1);
     gbm[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       gbm[i] = gbm[i - 1] + self.mu * gbm[i - 1] * dt + self.sigma * gbm[i - 1] * gn[i - 1]
     }
 

stochastic-rs-core/src/diffusion/jacobi.rs

@@ -37,17 +37,13 @@ impl Sampling<f64> for Jacobi {
     assert!(self.sigma > 0.0, "sigma must be positive");
     assert!(self.alpha < self.beta, "alpha must be less than beta");
 
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
-
     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 jacobi = Array1::<f64>::zeros(self.n + 1);
     jacobi[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       jacobi[i] = match jacobi[i - 1] {
         _ if jacobi[i - 1] <= 0.0 && i > 0 => 0.0,
         _ if jacobi[i - 1] >= 1.0 && i > 0 => 1.0,

stochastic-rs-core/src/diffusion/ou.rs

@@ -32,17 +32,13 @@ impl Ou {
 
 impl Sampling<f64> for Ou {
   fn sample(&self) -> Array1<f64> {
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
-
     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 ou = Array1::<f64>::zeros(self.n + 1);
     ou[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       ou[i] = ou[i - 1] + self.theta * (self.mu - ou[i - 1]) * dt + self.sigma * gn[i - 1]
     }
 

stochastic-rs-core/src/interest/duffie_kan.rs

@@ -69,7 +69,7 @@ impl Sampling2D<f64> for DuffieKan {
     r[0] = self.r0.unwrap_or(0.0);
     x[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       r[i] = r[i - 1]
         + (self.a1 * r[i - 1] + self.b1 * x[i - 1] + self.c1) * dt
         + self.sigma1 * (self.alpha * r[i - 1] + self.beta * x[i - 1] + self.gamma) * cgn1[i - 1];

stochastic-rs-core/src/interest/ho_lee.rs

@@ -48,7 +48,7 @@ impl<'a> Sampling<f64> for HoLee<'a> {
 
     let mut r = Array1::<f64>::zeros(self.n + 1);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let drift = if let Some(r#fn) = self.f_T.as_ref() {
         (r#fn)(i as f64 * dt) + self.sigma.powf(2.0)
       } else {

stochastic-rs-core/src/jump/bates.rs

@@ -91,7 +91,7 @@ impl<D: ProcessDistribution> Sampling2D<f64> for Bates1996<D> {
       _ => self.mu.unwrap(),
     };
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let [.., jumps] = self.cpoisson.sample();
 
       let sqrt_v = self

stochastic-rs-core/src/jump/ig.rs

@@ -30,14 +30,11 @@ impl Ig {
 impl Sampling<f64> for Ig {
   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, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
     let mut ig = Array1::zeros(self.n + 1);
     ig[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       ig[i] = ig[i - 1] + self.gamma * dt + gn[i - 1]
     }
 

stochastic-rs-core/src/jump/jump_fou.rs

@@ -66,7 +66,7 @@ impl<D: ProcessDistribution> Sampling<f64> for JumpFou<D> {
     let mut jump_fou = Array1::<f64>::zeros(self.n + 1);
     jump_fou[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let [.., jumps] = self.cpoisson.sample();
 
       jump_fou[i] = jump_fou[i - 1]

stochastic-rs-core/src/jump/levy_diffusion.rs

@@ -51,12 +51,9 @@ impl<D: ProcessDistribution> Sampling<f64> for LevyDiffusion<D> {
     let dt = self.t.unwrap_or(1.0) / self.n as f64;
     let mut levy = Array1::<f64>::zeros(self.n + 1);
     levy[0] = self.x0.unwrap_or(0.0);
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let [.., jumps] = self.cpoisson.sample();
       levy[i] = levy[i - 1] + self.gamma * dt + self.sigma * gn[i - 1] + jumps.sum();
     }

stochastic-rs-core/src/jump/merton.rs

@@ -53,12 +53,9 @@ impl<D: ProcessDistribution> Sampling<f64> for Merton<D> {
     let dt = self.t.unwrap_or(1.0) / self.n as f64;
     let mut merton = Array1::<f64>::zeros(self.n + 1);
     merton[0] = self.x0.unwrap_or(0.0);
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       let [.., jumps] = self.cpoisson.sample();
       merton[i] = merton[i - 1]
         + (self.alpha * self.sigma.powf(2.0) / 2.0 - self.lambda * self.theta) * dt

stochastic-rs-core/src/jump/nig.rs

@@ -37,14 +37,11 @@ impl Sampling<f64> for Nig {
     let scale = dt.powf(2.0) / self.kappa;
     let mean = dt / scale;
     let ig = Array1::random(self.n, InverseGaussian::new(mean, scale).unwrap());
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
     let mut nig = Array1::zeros(self.n + 1);
     nig[0] = self.x0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       nig[i] = nig[i - 1] + self.theta * ig[i - 1] + self.sigma * ig[i - 1].sqrt() * gn[i - 1]
     }
 

stochastic-rs-core/src/jump/vg.rs

@@ -41,13 +41,10 @@ impl Sampling<f64> for Vg {
     let mut vg = Array1::<f64>::zeros(self.n + 1);
     vg[0] = self.x0.unwrap_or(0.0);
 
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    let gn = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
     let gammas = Array1::random(self.n, Gamma::new(shape, scale).unwrap());
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       vg[i] = vg[i - 1] + self.mu * gammas[i - 1] + self.sigma * gammas[i - 1].sqrt() * gn[i - 1];
     }
 

stochastic-rs-core/src/noise/cfgns.rs

@@ -45,7 +45,7 @@ impl Sampling2D<f64> for Cfgns {
     let fgn1 = self.fgn.sample();
     let fgn2 = self.fgn.sample();
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       cfgns[[0, i]] = fgn1[i - 1];
       cfgns[[1, i]] = self.rho * fgn1[i - 1] + (1.0 - self.rho.powi(2)).sqrt() * fgn2[i - 1];
     }

stochastic-rs-core/src/noise/cgns.rs

@@ -31,17 +31,12 @@ impl Sampling2D<f64> for Cgns {
       "Correlation coefficient must be in [-1, 1]"
     );
 
+    let dt = self.t.unwrap_or(1.0) / self.n as f64;
     let mut cgns = Array2::<f64>::zeros((2, self.n + 1));
-    let gn1 = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
-    let gn2 = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    let gn1 = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
+    let gn2 = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap());
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       cgns[[0, i]] = cgns[[0, i - 1]] + gn1[i - 1];
       cgns[[1, i]] =
         cgns[[1, i - 1]] + self.rho * gn1[i - 1] + (1.0 - self.rho.powi(2)).sqrt() * gn2[i - 1];

stochastic-rs-core/src/process/bm.rs

@@ -24,10 +24,8 @@ impl Bm {
 
 impl Sampling<f64> for Bm {
   fn sample(&self) -> Array1<f64> {
-    let gn = Array1::random(
-      self.n,
-      Normal::new(0.0, (self.t.unwrap_or(1.0) / self.n as f64).sqrt()).unwrap(),
-    );
+    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 bm = Array1::<f64>::zeros(self.n);
     bm.slice_mut(s![1..]).assign(&gn);
 

stochastic-rs-core/src/process/cbms.rs

@@ -41,7 +41,7 @@ impl Sampling2D<f64> for Cbms {
     let mut bms = Array2::<f64>::zeros((2, self.n + 1));
     let [cgn1, cgn2] = self.cgns.sample();
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       bms[[0, i]] = bms[[0, i - 1]] + cgn1[i - 1];
       bms[[1, i]] =
         bms[[1, i - 1]] + self.rho * cgn1[i - 1] + (1.0 - self.rho.powi(2)).sqrt() * cgn2[i - 1];

stochastic-rs-core/src/process/cfbms.rs

@@ -58,7 +58,7 @@ impl Sampling2D<f64> for Cfbms {
     let mut fbms = Array2::<f64>::zeros((2, self.n + 1));
     let [fgn1, fgn2] = self.cfgns.sample();
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       fbms[[0, i]] = fbms[[0, i - 1]] + fgn1[i - 1];
       fbms[[1, i]] =
         fbms[[1, i - 1]] + self.rho * fgn1[i - 1] + (1.0 - self.rho.powi(2)).sqrt() * fgn2[i - 1];

stochastic-rs-core/src/process/fbm.rs

@@ -35,7 +35,7 @@ impl Sampling<f64> for Fbm {
     let mut fbm = Array1::<f64>::zeros(self.n + 1);
     fbm.slice_mut(s![1..]).assign(&fgn);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       fbm[i] += fbm[i - 1];
     }
 

stochastic-rs-core/src/volatility.rs

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