diff --git a/stochastic-rs-core/Cargo.toml b/stochastic-rs-core/Cargo.toml
index d090368..2535bbf 100644
--- a/stochastic-rs-core/Cargo.toml
+++ b/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]
 
diff --git a/stochastic-rs-core/src/diffusion/cir.rs b/stochastic-rs-core/src/diffusion/cir.rs
index a232590..c9f0776 100644
--- a/stochastic-rs-core/src/diffusion/cir.rs
+++ b/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],
diff --git a/stochastic-rs-core/src/diffusion/fcir.rs b/stochastic-rs-core/src/diffusion/fcir.rs
index 5a72efb..ca7174b 100644
--- a/stochastic-rs-core/src/diffusion/fcir.rs
+++ b/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],
diff --git a/stochastic-rs-core/src/diffusion/fgbm.rs b/stochastic-rs-core/src/diffusion/fgbm.rs
index f4753dc..67b80b4 100644
--- a/stochastic-rs-core/src/diffusion/fgbm.rs
+++ b/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]
     }
 
diff --git a/stochastic-rs-core/src/diffusion/fjacobi.rs b/stochastic-rs-core/src/diffusion/fjacobi.rs
index dc0303c..340c4b7 100644
--- a/stochastic-rs-core/src/diffusion/fjacobi.rs
+++ b/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,
diff --git a/stochastic-rs-core/src/diffusion/fou.rs b/stochastic-rs-core/src/diffusion/fou.rs
index 5b3a4bb..8b4ac3e 100644
--- a/stochastic-rs-core/src/diffusion/fou.rs
+++ b/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]
     }
 
diff --git a/stochastic-rs-core/src/diffusion/gbm.rs b/stochastic-rs-core/src/diffusion/gbm.rs
index 9b4b757..b43c107 100644
--- a/stochastic-rs-core/src/diffusion/gbm.rs
+++ b/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]
     }
 
diff --git a/stochastic-rs-core/src/diffusion/jacobi.rs b/stochastic-rs-core/src/diffusion/jacobi.rs
index 1fd0962..30abacf 100644
--- a/stochastic-rs-core/src/diffusion/jacobi.rs
+++ b/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,
diff --git a/stochastic-rs-core/src/diffusion/ou.rs b/stochastic-rs-core/src/diffusion/ou.rs
index ed059e5..9c3e6df 100644
--- a/stochastic-rs-core/src/diffusion/ou.rs
+++ b/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]
     }
 
diff --git a/stochastic-rs-core/src/interest/duffie_kan.rs b/stochastic-rs-core/src/interest/duffie_kan.rs
index 225e812..7f1763a 100644
--- a/stochastic-rs-core/src/interest/duffie_kan.rs
+++ b/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];
diff --git a/stochastic-rs-core/src/interest/ho_lee.rs b/stochastic-rs-core/src/interest/ho_lee.rs
index b800f4d..ab676de 100644
--- a/stochastic-rs-core/src/interest/ho_lee.rs
+++ b/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 {
diff --git a/stochastic-rs-core/src/jump/bates.rs b/stochastic-rs-core/src/jump/bates.rs
index cb14f3f..ceff548 100644
--- a/stochastic-rs-core/src/jump/bates.rs
+++ b/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
diff --git a/stochastic-rs-core/src/jump/ig.rs b/stochastic-rs-core/src/jump/ig.rs
index d5f1341..fb658b8 100644
--- a/stochastic-rs-core/src/jump/ig.rs
+++ b/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]
     }
 
diff --git a/stochastic-rs-core/src/jump/jump_fou.rs b/stochastic-rs-core/src/jump/jump_fou.rs
index c060455..e31eb2c 100644
--- a/stochastic-rs-core/src/jump/jump_fou.rs
+++ b/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]
diff --git a/stochastic-rs-core/src/jump/levy_diffusion.rs b/stochastic-rs-core/src/jump/levy_diffusion.rs
index da268e4..5c43db3 100644
--- a/stochastic-rs-core/src/jump/levy_diffusion.rs
+++ b/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();
     }
diff --git a/stochastic-rs-core/src/jump/merton.rs b/stochastic-rs-core/src/jump/merton.rs
index 481f8ea..b33d517 100644
--- a/stochastic-rs-core/src/jump/merton.rs
+++ b/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
diff --git a/stochastic-rs-core/src/jump/nig.rs b/stochastic-rs-core/src/jump/nig.rs
index a0eebbf..7aaab0f 100644
--- a/stochastic-rs-core/src/jump/nig.rs
+++ b/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]
     }
 
diff --git a/stochastic-rs-core/src/jump/vg.rs b/stochastic-rs-core/src/jump/vg.rs
index fef7947..bf01cb3 100644
--- a/stochastic-rs-core/src/jump/vg.rs
+++ b/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];
     }
 
diff --git a/stochastic-rs-core/src/noise/cfgns.rs b/stochastic-rs-core/src/noise/cfgns.rs
index 3831fca..a63947c 100644
--- a/stochastic-rs-core/src/noise/cfgns.rs
+++ b/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];
     }
diff --git a/stochastic-rs-core/src/noise/cgns.rs b/stochastic-rs-core/src/noise/cgns.rs
index d1e8553..5908a21 100644
--- a/stochastic-rs-core/src/noise/cgns.rs
+++ b/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];
diff --git a/stochastic-rs-core/src/process/bm.rs b/stochastic-rs-core/src/process/bm.rs
index cbe0000..19b4bbc 100644
--- a/stochastic-rs-core/src/process/bm.rs
+++ b/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);
 
diff --git a/stochastic-rs-core/src/process/cbms.rs b/stochastic-rs-core/src/process/cbms.rs
index a43de9f..175d377 100644
--- a/stochastic-rs-core/src/process/cbms.rs
+++ b/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];
diff --git a/stochastic-rs-core/src/process/cfbms.rs b/stochastic-rs-core/src/process/cfbms.rs
index b15c50f..82a0c5c 100644
--- a/stochastic-rs-core/src/process/cfbms.rs
+++ b/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];
diff --git a/stochastic-rs-core/src/process/fbm.rs b/stochastic-rs-core/src/process/fbm.rs
index 3444031..83444ce 100644
--- a/stochastic-rs-core/src/process/fbm.rs
+++ b/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];
     }
 
diff --git a/stochastic-rs-core/src/volatility.rs b/stochastic-rs-core/src/volatility.rs
index 7f06953..09dfe61 100644
--- a/stochastic-rs-core/src/volatility.rs
+++ b/stochastic-rs-core/src/volatility.rs
@@ -1,3 +1,4 @@
+pub mod fheston;
 pub mod heston;
 pub mod rbergomi;
 pub mod sabr;
diff --git a/stochastic-rs-core/src/volatility/fheston.rs b/stochastic-rs-core/src/volatility/fheston.rs
new file mode 100644
index 0000000..1e4fd65
--- /dev/null
+++ b/stochastic-rs-core/src/volatility/fheston.rs
@@ -0,0 +1,117 @@
+use ndarray::Array1;
+use ndarray_rand::RandomExt;
+use rand_distr::Normal;
+use statrs::function::gamma::gamma;
+
+use crate::Sampling;
+
+#[derive(Default)]
+pub struct RoughHeston {
+  pub v0: Option<f64>,
+  pub theta: f64,
+  pub kappa: f64,
+  pub nu: f64,
+  pub hurst: f64,
+  pub c1: Option<f64>,
+  pub c2: Option<f64>,
+  pub t: Option<f64>,
+  pub n: usize,
+  pub m: Option<usize>,
+}
+
+impl RoughHeston {
+  #[must_use]
+  pub fn new(params: &Self) -> Self {
+    Self {
+      v0: params.v0,
+      theta: params.theta,
+      kappa: params.kappa,
+      nu: params.nu,
+      hurst: params.hurst,
+      c1: params.c1,
+      c2: params.c2,
+      t: params.t,
+      n: params.n,
+      m: params.m,
+    }
+  }
+}
+
+impl Sampling<f64> for RoughHeston {
+  fn sample(&self) -> ndarray::Array1<f64> {
+    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 yt = Array1::<f64>::zeros(self.n + 1);
+    let mut zt = Array1::<f64>::zeros(self.n + 1);
+    let mut v2 = Array1::zeros(self.n + 1);
+
+    yt[0] = self.theta + (self.v0.unwrap_or(1.0).powi(2) - self.theta) * (-self.kappa * 0.0).exp();
+    zt[0] = 0.0; // Initial condition for Z_t, typically 0 for such integrals.
+    v2[0] = self.v0.unwrap_or(1.0).powi(2);
+
+    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];
+
+      let integral = (0..i)
+        .map(|j| {
+          let tj = j as f64 * dt;
+          ((t - tj).powf(self.hurst - 0.5) * zt[j]) * dt
+        })
+        .sum::<f64>();
+
+      v2[i] = yt[i]
+        + self.c1.unwrap_or(1.0) * self.nu * zt[i]
+        + self.c2.unwrap_or(1.0) * self.nu * integral / gamma(self.hurst + 0.5);
+    }
+
+    v2
+  }
+
+  fn n(&self) -> usize {
+    self.n
+  }
+
+  fn m(&self) -> Option<usize> {
+    self.m
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use plotly::{common::Line, Plot, Scatter};
+
+  use super::*;
+
+  #[test]
+  fn test_rough_heston() {
+    let params = RoughHeston {
+      v0: Some(1.2),
+      theta: 1.0,
+      kappa: 1.0,
+      nu: 0.2,
+      hurst: 0.75,
+      c1: Some(1.0),
+      c2: Some(1.0),
+      t: Some(1.0),
+      n: 1000,
+      m: Some(1000),
+    };
+
+    let rh = RoughHeston::new(&params);
+    let _sample = rh.sample();
+
+    let mut plot = Plot::new();
+    let trace = Scatter::new((0.._sample.len()).collect::<Vec<_>>(), _sample.to_vec())
+      .mode(plotly::common::Mode::Lines)
+      .line(
+        Line::new()
+          .color("orange")
+          .shape(plotly::common::LineShape::Linear),
+      )
+      .name("Fbm");
+    plot.add_trace(trace);
+    plot.show();
+  }
+}
diff --git a/stochastic-rs-core/src/volatility/heston.rs b/stochastic-rs-core/src/volatility/heston.rs
index a4edea1..d839ada 100644
--- a/stochastic-rs-core/src/volatility/heston.rs
+++ b/stochastic-rs-core/src/volatility/heston.rs
@@ -57,7 +57,7 @@ impl Sampling2D<f64> for Heston {
     s[0] = self.s0.unwrap_or(0.0);
     v[0] = self.v0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       s[i] = s[i - 1] + self.mu * s[i - 1] * dt + s[i - 1] * v[i - 1].sqrt() * cgn1[i - 1];
 
       let random: f64 = match self.use_sym.unwrap_or(false) {
diff --git a/stochastic-rs-core/src/volatility/rbergomi.rs b/stochastic-rs-core/src/volatility/rbergomi.rs
index 330533f..36c5f0a 100644
--- a/stochastic-rs-core/src/volatility/rbergomi.rs
+++ b/stochastic-rs-core/src/volatility/rbergomi.rs
@@ -6,7 +6,7 @@ use crate::{noise::cgns::Cgns, Sampling2D};
 pub struct RoughBergomi {
   pub hurst: f64,
   pub nu: f64,
-  pub sigma_0: f64,
+  pub v0: Option<f64>,
   pub r: f64,
   pub rho: f64,
   pub n: usize,
@@ -28,7 +28,7 @@ impl RoughBergomi {
     Self {
       hurst: params.hurst,
       nu: params.nu,
-      sigma_0: params.sigma_0,
+      v0: params.v0,
       r: params.r,
       rho: params.rho,
       n: params.n,
@@ -45,20 +45,21 @@ impl Sampling2D<f64> for RoughBergomi {
     let [cgn1, z] = self.cgns.sample();
 
     let mut s = Array1::<f64>::zeros(self.n + 1);
-    let mut sigma2 = Array1::<f64>::zeros(self.n + 1);
+    let mut v2 = Array1::<f64>::zeros(self.n + 1);
+    v2[0] = self.v0.unwrap_or(1.0).powi(2);
 
-    for i in 1..(self.n + 1) {
-      s[i] = s[i - 1] + self.r * s[i - 1] + sigma2[i - 1] * cgn1[i - 1];
+    for i in 1..=self.n {
+      s[i] = s[i - 1] + self.r * s[i - 1] + v2[i - 1].sqrt() * cgn1[i - 1];
 
       let sum_z = z.slice(s![..i]).sum();
       let t = i as f64 * dt;
-      sigma2[i] = self.sigma_0.powi(2)
+      v2[i] = self.v0.unwrap_or(1.0).powi(2)
         * (self.nu * (2.0 * self.hurst).sqrt() * t.powf(self.hurst - 0.5) * sum_z
           - 0.5 * self.nu.powi(2) * t.powf(2.0 * self.hurst))
         .exp();
     }
 
-    [s, sigma2]
+    [s, v2]
   }
 
   fn n(&self) -> usize {
diff --git a/stochastic-rs-core/src/volatility/sabr.rs b/stochastic-rs-core/src/volatility/sabr.rs
index b30a24f..d5ae988 100644
--- a/stochastic-rs-core/src/volatility/sabr.rs
+++ b/stochastic-rs-core/src/volatility/sabr.rs
@@ -50,7 +50,7 @@ impl Sampling2D<f64> for Sabr {
     f[0] = self.f0.unwrap_or(0.0);
     v[0] = self.v0.unwrap_or(0.0);
 
-    for i in 1..(self.n + 1) {
+    for i in 1..=self.n {
       f[i] = f[i - 1] + v[i - 1] * f[i - 1].powf(self.beta) * cgn1[i - 1];
       v[i] = v[i - 1] + self.alpha * v[i - 1] * cgn2[i - 1];
     }