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; 3];
for i in 1..n {
sum[0] += v[i] / v[i - 1];
sum[1] += v[i];
sum[2] += v[i - 1];
}
let beta_hat1 = ((n as f64).powi(-2) * sum[1] * sum[2].powi(-1) - (n as f64).powi(-1) * sum[0])
/ ((n as f64).powi(-2) - 1.0);
let beta_hat2 =
((n as f64).powi(-1) * sum[0] - beta_hat1) / ((1.0 - beta_hat1) * (n as f64).powi(-1) * sum[2]);
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 {
let dw1_i = (s[i].ln() - s[i - 1].ln() - (r - 0.5 * v[i - 1]) * delta) / v[i - 1].sqrt();
let dw2_i = (v[i] - v[i - 1] - kappa_hat * (theta_hat - v[i - 1]) * delta)
/ (sigma_hat * v[i - 1].sqrt());
sum_dw1dw2 += dw1_i * dw2_i;
}
let rho_hat = sum_dw1dw2 / (n as f64 * delta);
vec![
v.first().unwrap().clone(),
theta_hat,
rho_hat,
kappa_hat,
sigma_hat,
]
}