fix: cgmy models

Author: dancixx

src/stochastic/jump/cgmy.rs

@@ -75,28 +75,29 @@ impl Sampling<f64> for CGMY {
 
     let U = Array1::<f64>::random(self.j, Uniform::new(0.0, 1.0));
     let E = Array1::<f64>::random(self.j, Exp::new(1.0).unwrap());
+    let P = Poisson::new(1.0, Some(self.j), None, None);
+    let P = P.sample();
     let tau = Array1::<f64>::random(self.j, Uniform::new(0.0, 1.0));
-    let poisson = Poisson::new(1.0, Some(self.j), None, None);
-    let poisson = poisson.sample();
 
     for i in 1..self.n {
       let mut jump_component = 0.0;
       let t_1 = (i - 1) as f64 * dt;
       let t = i as f64 * dt;
 
       for j in 1..self.j {
-        let v_j = if rng.gen_bool(0.5) {
-          self.lambda_plus
-        } else {
-          -self.lambda_minus
-        };
-
-        let term1 = (self.alpha * poisson[j] / (2.0 * C * t_max)).powf(-1.0 / self.alpha);
-        let term2 = E[j] * U[j].powf(1.0 / self.alpha) / v_j.abs();
-        let jump_size =
-          term1.min(term2) * (v_j / v_j.abs()) * if tau[j] > t_1 && tau[j] < t { 1.0 } else { 0.0 };
-
-        jump_component += jump_size;
+        if tau[j] > t_1 && tau[j] <= t {
+          let v_j = if rng.gen_bool(0.5) {
+            self.lambda_plus
+          } else {
+            -self.lambda_minus
+          };
+
+          let term1 = (self.alpha * P[j] / (2.0 * C * t_max)).powf(-1.0 / self.alpha);
+          let term2 = E[j] * U[j].powf(1.0 / self.alpha) / v_j.abs();
+          let jump_size = term1.min(term2) * (v_j / v_j.abs());
+
+          jump_component += jump_size;
+        }
       }
 
       x[i] = x[i - 1] + jump_component + b_t * dt;
@@ -137,7 +138,7 @@ mod tests {
 
   #[test]
   fn cgmy_plot() {
-    let cgmy = CGMY::new(25.46, 4.604, 0.52, 100, 1024, Some(2.0), Some(1.0), None);
+    let cgmy = CGMY::new(25.46, 4.604, 0.52, 1000, 1024, Some(2.0), Some(1.0), None);
     plot_1d!(cgmy.sample(), "CGMY Process");
   }
 

src/stochastic/jump/cts.rs

@@ -59,18 +59,19 @@ impl Sampling<f64> for CTS {
       let t = i as f64 * dt;
 
       for j in 1..self.j {
-        let v_j = if rng.gen_bool(0.5) {
-          self.lambda_plus
-        } else {
-          -self.lambda_minus
-        };
-
-        let term1 = (self.alpha * poisson[j] / C).powf(-1.0 / self.alpha);
-        let term2 = E[j] * U[j].powf(1.0 / self.alpha) / v_j.abs();
-        let jump_size =
-          term1.min(term2) * (v_j / v_j.abs()) * if tau[j] > t_1 && tau[j] < t { 1.0 } else { 0.0 };
-
-        jump_component += jump_size;
+        if tau[j] > t_1 && tau[j] <= t {
+          let v_j = if rng.gen_bool(0.5) {
+            self.lambda_plus
+          } else {
+            -self.lambda_minus
+          };
+
+          let term1 = (self.alpha * poisson[j] / C).powf(-1.0 / self.alpha);
+          let term2 = E[j] * U[j].powf(1.0 / self.alpha) / v_j.abs();
+          let jump_size = term1.min(term2) * (v_j / v_j.abs());
+
+          jump_component += jump_size;
+        }
       }
 
       x[i] = x[i - 1] + jump_component + b_t * dt;
@@ -111,13 +112,13 @@ mod tests {
 
   #[test]
   fn cts_plot() {
-    let cts = CTS::new(25.46, 4.604, 0.52, 1024, 1024, Some(2.0), Some(1.0), None);
+    let cts = CTS::new(25.46, 4.604, 0.52, N, 1024, Some(2.0), Some(1.0), None);
     plot_1d!(cts.sample(), "CTS Process");
   }
 
   #[test]
   fn cts_plot_multi() {
-    let cts = CTS::new(25.46, 4.604, 0.52, N, 10000, Some(2.0), Some(1.0), Some(10));
+    let cts = CTS::new(25.46, 4.604, 0.52, N, 1024, Some(2.0), Some(1.0), Some(10));
     plot_nd!(cts.sample_par(), "CTS Process");
   }
 }

src/stochastic/jump/rdts.rs

@@ -59,18 +59,19 @@ impl Sampling<f64> for RDTS {
       let t = i as f64 * dt;
 
       for j in 1..self.j {
-        let v_j = if rng.gen_bool(0.5) {
-          self.lambda_plus
-        } else {
-          -self.lambda_minus
-        };
-
-        let term1 = (self.alpha * poisson[j] / (2.0 * C * t_max)).powf(-1.0 / self.alpha);
-        let term2 = 0.5 * E[j].powf(0.5) * U[j].powf(1.0 / self.alpha) / v_j.abs();
-        let jump_size =
-          term1.min(term2) * (v_j / v_j.abs()) * if tau[j] > t_1 && tau[j] < t { 1.0 } else { 0.0 };
-
-        jump_component += jump_size;
+        if tau[j] > t_1 && tau[j] <= t {
+          let v_j = if rng.gen_bool(0.5) {
+            self.lambda_plus
+          } else {
+            -self.lambda_minus
+          };
+
+          let term1 = (self.alpha * poisson[j] / (2.0 * C * t_max)).powf(-1.0 / self.alpha);
+          let term2 = 0.5 * E[j].powf(0.5) * U[j].powf(1.0 / self.alpha) / v_j.abs();
+          let jump_size = term1.min(term2) * (v_j / v_j.abs());
+
+          jump_component += jump_size;
+        }
       }
 
       x[i] = x[i - 1] + jump_component + b_t * dt;

src/stochastic/volatility/svcgmy.rs

@@ -10,6 +10,9 @@ use crate::{
   stochastic::{process::poisson::Poisson, Sampling},
 };
 
+/// CGMY Stochastic Volatility process
+///
+/// https://www.econstor.eu/bitstream/10419/239493/1/175133161X.pdf
 #[derive(ImplNew)]
 pub struct SVCGMY {
   /// Positive jump rate lambda_plus (corresponds to G)
@@ -18,19 +21,19 @@ pub struct SVCGMY {
   pub lambda_minus: f64, // M
   /// Jump activity parameter alpha (corresponds to Y), with 0 < alpha < 2
   pub alpha: f64,
-  ///
+  /// Mean reversion rate
   pub kappa: f64,
-  ///
+  /// Long-term volatility
   pub eta: f64,
-  ///
+  /// Volatility of volatility
   pub zeta: f64,
   ///
   pub rho: f64,
   /// Number of time steps
   pub n: usize,
   /// Jumps
   pub j: usize,
-  ///
+  /// Initial value
   pub x0: Option<f64>,
   /// Initial value
   pub v0: Option<f64>,
@@ -45,80 +48,74 @@ impl Sampling<f64> for SVCGMY {
     let mut rng = rand::thread_rng();
 
     let t_max = self.t.unwrap_or(1.0);
-    let dt = self.t.unwrap_or(1.0) / (self.n - 1) as f64;
-    let mut x = Array1::<f64>::zeros(self.n);
-    x[0] = self.x0.unwrap_or(0.0);
-
-    let C = (gamma(2.0 - self.alpha)
-      * (self.lambda_plus.powf(self.alpha - 2.0) + self.lambda_minus.powf(self.alpha - 2.0)))
-    .powi(-1);
-    let c = (2.0 * self.kappa) / ((1.0 - (-self.kappa * dt).exp()) * self.zeta.powi(2));
+    let dt = t_max / (self.n - 1) as f64;
 
+    let mut x = Array1::<f64>::zeros(self.n);
     let mut v = Array1::<f64>::zeros(self.n);
+    let mut y = Array1::<f64>::zeros(self.n);
+
+    x[0] = self.x0.unwrap_or(0.0);
     v[0] = self.v0.unwrap_or(0.0);
 
+    let C = 1.0
+      / (gamma(2.0 - self.alpha)
+        * (self.lambda_plus.powf(self.alpha - 2.0) + self.lambda_minus.powf(self.alpha - 2.0)));
+    let c = (2.0 * self.kappa) / ((1.0 - (-self.kappa * dt).exp()) * self.zeta.powi(2));
     let df = 4.0 * self.kappa * self.eta / self.zeta.powi(2);
 
+    // Volatilitás folyamat generálása
     for i in 1..self.n {
       let ncp = 2.0 * c * v[i - 1] * (-self.kappa * dt).exp();
       let xi = non_central_chi_squared::sample(df, ncp, &mut rng);
-
       v[i] = xi / (2.0 * c);
     }
 
     let U = Array1::<f64>::random(self.j, Uniform::new(0.0, 1.0));
-    let E = Array1::<f64>::random(self.j, Exp::new(1.0).unwrap());
-    let E_ = Array1::<f64>::random(self.j, Exp::new(1.0).unwrap());
+    let E = Array1::random(self.j, Exp::new(1.0).unwrap());
     let P = Poisson::new(1.0, Some(self.j), None, None);
-    let mut P = P.sample();
-
-    for (idx, p) in P.iter_mut().enumerate() {
-      *p = *p + E_[idx];
-    }
+    let P = P.sample();
+    let tau = Array1::<f64>::random(self.j, Uniform::new(0.0, 1.0)) * t_max;
 
-    let tau = Array1::<f64>::random(self.j, Uniform::new(0.0, t_max));
     let mut c_tau = Array1::<f64>::zeros(self.j);
-
-    for j in 1..self.j {
-      let mut c_tau_ = 0.0;
-      for i in 1..self.n {
-        let t_1 = (i - 1) as f64 * dt;
-        let t = i as f64 * dt;
-
-        if t_1 < tau[j] && tau[j] <= t {
-          c_tau_ += v[i - 1]
-        }
-      }
-
-      c_tau[j] = C * c_tau_;
+    for (idx, tau_j) in tau.iter().enumerate() {
+      let k = ((tau_j / dt).ceil() as usize).min(self.n - 1);
+      let v_k = if k == 0 { v[0] } else { v[k - 1] };
+      c_tau[idx] = C * v_k;
     }
-    println!("{:?}", c_tau);
+
     for i in 1..self.n {
-      let b = -((v[i]
-        * (self.lambda_plus.powf(self.alpha - 1.0) - self.lambda_minus.powf(self.alpha - 1.0)))
-        / ((1.0 - self.alpha)
-          * (self.lambda_plus.powf(self.alpha - 2.0) + self.lambda_minus.powf(self.alpha) - 2.0)));
+      let numerator = v[i - 1]
+        * (self.lambda_plus.powf(self.alpha - 1.0) - self.lambda_minus.powf(self.alpha - 1.0));
+      let denominator = (1.0 - self.alpha)
+        * (self.lambda_plus.powf(self.alpha - 2.0) + self.lambda_minus.powf(self.alpha - 2.0));
+      let b = -numerator / denominator;
 
       let mut jump_component = 0.0;
+
       let t_1 = (i - 1) as f64 * dt;
       let t = i as f64 * dt;
 
-      for j in 1..self.j {
-        let v_j = if rng.gen_bool(0.5) {
-          self.lambda_plus
-        } else {
-          -self.lambda_minus
-        };
-
-        let term1 = ((self.alpha * P[j]) / (2.0 * c_tau[j] * t_max)).powf(-1.0 / self.alpha);
-        let term2 = E[j] * U[j].powf(1.0 / self.alpha) / v_j.abs();
-        let jump_size =
-          term1.min(term2) * (v_j / v_j.abs()) * if tau[j] > t_1 && tau[j] < t { 1.0 } else { 0.0 };
-
-        jump_component += jump_size;
+      for j in 0..self.j {
+        if tau[j] > t_1 && tau[j] <= t {
+          let v_j = if rng.gen_bool(0.5) {
+            self.lambda_plus
+          } else {
+            -self.lambda_minus
+          };
+
+          let term1 = ((self.alpha * P[j]) / (2.0 * c_tau[j] * t_max)).powf(-1.0 / self.alpha);
+          let term2 = E[j] * U[j].powf(1.0 / self.alpha) / v_j.abs();
+          let min_term = term1.min(term2);
+          let jump_size = min_term * (v_j / v_j.abs());
+          jump_component += jump_size;
+        }
       }
 
-      x[i] = x[i - 1] + jump_component + b * dt + self.rho * v[i - 1];
+      y[i] = y[i - 1] + jump_component + b * dt;
+    }
+
+    for i in 0..self.n {
+      x[i] = y[i] + self.rho * v[i];
     }
 
     x
@@ -135,8 +132,10 @@ impl Sampling<f64> for SVCGMY {
 
 #[cfg(test)]
 mod tests {
+  use ndarray::Axis;
+
   use super::*;
-  use crate::{plot_1d, stochastic::N};
+  use crate::{plot_1d, plot_nd, stochastic::N};
 
   #[test]
   fn svcgmy_length_equals_n() {
@@ -188,7 +187,7 @@ mod tests {
       0.0711,
       0.3443,
       -2.0280,
-      100,
+      1000,
       1024,
       Some(-0.25),
       Some(0.0064),
@@ -197,4 +196,24 @@ mod tests {
     );
     plot_1d!(svcgmy.sample(), "SVCGMY Process");
   }
+
+  #[test]
+  fn svcgmy_plot_multi() {
+    let svcgmy = SVCGMY::new(
+      25.46,
+      4.604,
+      0.52,
+      1.003,
+      0.0711,
+      0.3443,
+      -2.0280,
+      1000,
+      1024,
+      Some(-0.25),
+      Some(0.0064),
+      Some(1.0),
+      Some(10),
+    );
+    plot_nd!(svcgmy.sample_par(), "SVCGMY Process");
+  }
 }