fix: heston nmle

stochastic-rs-stats/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "stochastic-rs-stats"
-version = "0.1.2"
+version = "0.1.4"
 edition = "2021"
 license = "MIT"
 description = "A Rust library for the statistical analysis of stochastic processes"

stochastic-rs-stats/src/mle.rs

@@ -1,4 +1,5 @@
 /// Maximum likelihood estimation for Heston model
+/// http://scis.scichina.com/en/2018/042202.pdf
 ///
 /// # Arguments
 /// s: Vec<f64> - stock prices
@@ -7,43 +8,42 @@
 ///
 /// # Returns
 /// Vec<f64> - estimated parameters
-pub fn mle_heston(s: Vec<f64>, v: Vec<f64>, r: f64) -> Vec<f64> {
+pub fn nmle_heston(s: Vec<f64>, v: Vec<f64>, r: f64) -> Vec<f64> {
   let n = v.len();
   let delta = 1.0 / n as f64;
-  let mut sum = [0.0; 6];
+  let mut sum = [0.0; 4];
 
   for i in 1..n {
-    // sum of sqrt(V_i * V_{i-1})
-    sum[0] += (v[i] * v[i - 1]).sqrt();
-
-    // sum of sqrt(V_i / V_{i-1})
-    sum[1] += (v[i] / v[i - 1]).sqrt();
+    // sum of (V_i / V_{i-1}
+    sum[0] += v[i] / v[i - 1];
 
     // sum of V_i
-    sum[2] += v[i];
+    sum[1] += v[i];
 
     // sum of V_{i-1}
-    sum[3] += v[i - 1];
-
-    // sum of sqrt(V_i)
-    sum[4] += v[i].sqrt();
+    sum[2] += v[i - 1];
 
-    // sum of sqrt(V_{i-1})
-    sum[5] += v[i - 1].sqrt();
+    // sum of 1 / V_{i-1}
+    sum[3] += 1.0 / v[i - 1];
   }
 
-  let P_hat = ((1.0 / n as f64) * sum[0] - (1.0 / n as f64).powi(2) * sum[1] * sum[3])
-    / ((delta / 2.0) - (delta / 2.0) * (1.0 / n as f64).powi(2) * (1.0 / sum[3]) * sum[3]);
-
-  let kappa_hat = (2.0 / delta)
-    * (1.0 + (P_hat * delta / 2.0) * (1.0 / n as f64) * (1.0 / sum[3]) - (1.0 / n as f64) * sum[1]);
-
-  let sigma_hat = ((4.0 / delta)
-    * (1.0 / n as f64)
-    * (sum[4] - sum[5] - (delta / (2.0 * sum[5])) * (P_hat - kappa_hat * sum[3])).powi(2))
-  .sqrt();
+  let beta_hat1 = ((n as f64).powi(-2) * sum[1] * sum[3] - (n as f64).powi(-1) * sum[0])
+    / ((n as f64).powi(-2) * sum[2] * sum[3] - 1.0);
+  let beta_hat2 =
+    ((n as f64).powi(-1) * sum[0] - beta_hat1) / ((1.0 - beta_hat1) * (n as f64).powi(-1) * sum[3]);
 
-  let theta_hat = (P_hat + 0.25 * sigma_hat.powi(2)) / kappa_hat;
+  let sum_beta_hat3 = {
+    let mut sum = 0.0;
+    for i in 1..n {
+      sum +=
+        (v[i] - v[i - 1] * beta_hat1 - beta_hat2 * (1.0 - beta_hat1).powi(2)) * (1.0 / v[i - 1])
+    }
+    sum
+  };
+  let beta_hat3 = (n as f64).powi(-1) * sum_beta_hat3;
+  let kappa_hat = -(1.0 / delta) * beta_hat1.ln();
+  let theta_hat = beta_hat2;
+  let sigma_hat = (2.0 * kappa_hat * beta_hat3) / (1.0 - beta_hat1.powi(2)).sqrt();
 
   let mut sum_dw1dw2 = 0.0;
   for i in 1..n {