From 2ded4113dc17cda5e142f8371e1bb341667c89fe Mon Sep 17 00:00:00 2001 From: Daniel Boros Date: Fri, 8 Nov 2024 09:15:58 +0100 Subject: [PATCH] feat: add multiple fou estimartors --- src/stats.rs | 1 + src/stats/fou_estimator.rs | 328 +++++++++++++++++++++++++++++++++++-- 2 files changed, 313 insertions(+), 16 deletions(-) diff --git a/src/stats.rs b/src/stats.rs index 75e445c..6386a0f 100644 --- a/src/stats.rs +++ b/src/stats.rs @@ -2,5 +2,6 @@ pub mod cir; pub mod double_exp; pub mod fd; pub mod fou_estimator; +pub mod fou_estimator_v2; pub mod mle; pub mod non_central_chi_squared; diff --git a/src/stats/fou_estimator.rs b/src/stats/fou_estimator.rs index eae28bf..b056362 100644 --- a/src/stats/fou_estimator.rs +++ b/src/stats/fou_estimator.rs @@ -1,12 +1,14 @@ use impl_new_derive::ImplNew; -use ndarray::{array, Array1}; +use ndarray::{array, s, Array1}; use statrs::function::gamma::gamma; use std::f64::consts::SQRT_2; +use crate::stochastic::{noise::fgn::FGN, Sampling}; + +// Version 1: FOUParameterEstimationV1 with linear filter methods #[derive(ImplNew)] -pub struct FOUParameterEstimation { +pub struct FOUParameterEstimationV1 { pub path: Array1, - pub T: f64, pub filter_type: FilterType, // Estimated parameters hurst: Option, @@ -26,7 +28,7 @@ pub enum FilterType { Classical, } -impl FOUParameterEstimation { +impl FOUParameterEstimationV1 { pub fn estimate_parameters(&mut self) -> (f64, f64, f64, f64) { self.linear_filter(); self.hurst_estimator(); @@ -86,7 +88,7 @@ impl FOUParameterEstimation { let denominator = sigma.powi(2) * gamma(2.0 * hurst + 1.0); let theta = (numerator / denominator).powf(-1.0 / (2.0 * hurst)); - self.theta = Some(theta / self.T); + self.theta = Some(theta); } fn linear_filter(&mut self) { @@ -138,7 +140,6 @@ impl FOUParameterEstimation { } fn lfilter(&self, b: &Array1, a: &Array1, x: &Array1) -> Array1 { - // Implementing the difference equation: y[n] = b[0]*x[n] + b[1]*x[n-1] + ... - a[1]*y[n-1] - ... let n = x.len(); let mut y = Array1::::zeros(n); @@ -161,20 +162,268 @@ impl FOUParameterEstimation { } } +// Version 2: FOUParameterEstimationV2 without linear filters +#[derive(ImplNew)] +pub struct FOUParameterEstimationV2 { + pub path: Array1, + pub delta: f64, + pub series_length: usize, + // Estimated parameters + hurst: Option, + sigma: Option, + mu: Option, + theta: Option, +} + +impl FOUParameterEstimationV2 { + pub fn estimate_parameters(&mut self) -> (f64, f64, f64, f64) { + self.hurst_estimator(); + self.sigma_estimator(); + self.mu_estimator(); + self.theta_estimator(); + + ( + self.hurst.unwrap(), + self.sigma.unwrap(), + self.mu.unwrap(), + self.theta.unwrap(), + ) + } + + fn hurst_estimator(&mut self) { + let X = &self.path; + let N = self.series_length; + + let sum1: f64 = (0..(N - 4)) + .map(|i| { + let diff = X[i + 4] - 2.0 * X[i + 2] + X[i]; + diff * diff + }) + .sum(); + + let sum2: f64 = (0..(N - 2)) + .map(|i| { + let diff = X[i + 2] - 2.0 * X[i + 1] + X[i]; + diff * diff + }) + .sum(); + + let estimated_hurst = 0.5 * (sum1 / sum2).log2(); + self.hurst = Some(estimated_hurst); + } + + fn sigma_estimator(&mut self) { + let H = self.hurst.unwrap(); + let X = &self.path; + let N = self.series_length as f64; + let delta = self.delta; + + let numerator: f64 = (0..(self.series_length - 2)) + .map(|i| { + let diff = X[i + 2] - 2.0 * X[i + 1] + X[i]; + diff * diff + }) + .sum(); + + let denominator = N * (4.0 - 2.0_f64.powf(2.0 * H)) * delta.powf(2.0 * H); + let estimated_sigma = (numerator / denominator).sqrt(); + self.sigma = Some(estimated_sigma); + } + + fn mu_estimator(&mut self) { + let mean = self.path.mean().unwrap(); + self.mu = Some(mean); + } + + fn theta_estimator(&mut self) { + let X = &self.path; + let H = self.hurst.unwrap(); + let N = self.series_length as f64; + let sigma = self.sigma.unwrap(); + + let sum_X_squared = X.mapv(|x| x * x).sum(); + let sum_X = X.sum(); + let numerator = N * sum_X_squared - sum_X.powi(2); + let denominator = N.powi(2) * sigma.powi(2) * H * gamma(2.0 * H); + + let estimated_theta = (numerator / denominator).powf(-1.0 / (2.0 * H)); + self.theta = Some(estimated_theta); + } +} + +// Version 3: FOUParameterEstimationV3 with get_path method +pub struct FOUParameterEstimationV3 { + alpha: f64, + mu: f64, + sigma: f64, + initial_value: f64, + T: f64, + delta: f64, + series_length: usize, + hurst: f64, + path: Option>, + // Estimated parameters + estimated_hurst: Option, + estimated_sigma: Option, + estimated_mu: Option, + estimated_alpha: Option, +} + +impl FOUParameterEstimationV3 { + pub fn new( + series_length: usize, + hurst: f64, + sigma: f64, + alpha: f64, + mu: f64, + initial_value: f64, + T: f64, + delta: f64, + ) -> Self { + FOUParameterEstimationV3 { + alpha, + mu, + sigma, + initial_value, + T, + delta, + series_length, + hurst, + path: None, + estimated_hurst: None, + estimated_sigma: None, + estimated_mu: None, + estimated_alpha: None, + } + } + + pub fn estimate_parameters(&mut self) -> (f64, f64, f64, f64) { + self.get_path(); + self.hurst_estimator(); + self.sigma_estimator(); + self.mu_estimator(); + self.alpha_estimator(); + + ( + self.estimated_hurst.unwrap(), + self.estimated_sigma.unwrap(), + self.estimated_mu.unwrap(), + self.estimated_alpha.unwrap(), + ) + } + + fn get_path(&mut self) { + let M = 8; + let gamma = self.delta / M as f64; + + let fgn_length = self.series_length * M; + + // Generate fGN sample of length fgn_length + let fgn = FGN::new(self.hurst, fgn_length - 1, Some(self.T), None); + let fgn_sample = fgn.sample(); + + // Initialize full_fou array + let mut full_fou = Array1::::zeros(fgn_length); + full_fou[0] = self.initial_value; + + for i in 1..fgn_length { + full_fou[i] = full_fou[i - 1] + + self.alpha * (self.mu - full_fou[i - 1]) * gamma + + self.sigma * fgn_sample[i - 1]; + } + + // Initialize fou array + let mut fou = Array1::::zeros(self.series_length); + fou[0] = self.initial_value; + + for i in 1..self.series_length { + let start = (i - 1) * M; + let end = i * M; + + let sum_sub_series: f64 = full_fou.slice(s![start..end]).sum() * gamma / M as f64; + fou[i] = full_fou[end - 1] + self.alpha * sum_sub_series; + } + + // Store the path + self.path = Some(fou); + } + + fn hurst_estimator(&mut self) { + let X = self.path.as_ref().unwrap(); + let N = self.series_length; + + let sum1: f64 = (0..(N - 4)) + .map(|i| { + let diff = X[i + 4] - 2.0 * X[i + 2] + X[i]; + diff * diff + }) + .sum(); + + let sum2: f64 = (0..(N - 2)) + .map(|i| { + let diff = X[i + 2] - 2.0 * X[i + 1] + X[i]; + diff * diff + }) + .sum(); + + let estimated_hurst = 0.5 * (sum1 / sum2).log2(); + self.estimated_hurst = Some(estimated_hurst); + } + + fn sigma_estimator(&mut self) { + let H = self.estimated_hurst.unwrap(); + let X = self.path.as_ref().unwrap(); + let N = self.series_length as f64; + let delta = self.delta; + + let numerator: f64 = (0..(self.series_length - 2)) + .map(|i| { + let diff = X[i + 2] - 2.0 * X[i + 1] + X[i]; + diff * diff + }) + .sum(); + + let denominator = N * (4.0 - 2.0_f64.powf(2.0 * H)) * delta.powf(2.0 * H); + let estimated_sigma = (numerator / denominator).sqrt(); + self.estimated_sigma = Some(estimated_sigma); + } + + fn mu_estimator(&mut self) { + let X = self.path.as_ref().unwrap(); + let mean = X.mean().unwrap(); + self.estimated_mu = Some(mean); + } + + fn alpha_estimator(&mut self) { + let X = self.path.as_ref().unwrap(); + let H = self.estimated_hurst.unwrap(); + let N = self.series_length as f64; + let sigma = self.estimated_sigma.unwrap(); + + let sum_X_squared = X.mapv(|x| x * x).sum(); + let sum_X = X.sum(); + let numerator = N * sum_X_squared - sum_X.powi(2); + let denominator = N.powi(2) * sigma.powi(2) * H * gamma(2.0 * H); + + let estimated_alpha = (numerator / denominator).powf(-1.0 / (2.0 * H)); + self.estimated_alpha = Some(estimated_alpha); + } +} + #[cfg(test)] mod tests { use super::*; use crate::stochastic::{diffusion::fou::FOU, noise::fgn::FGN, Sampling}; #[test] - fn test_fou_parameter_estimation() { + fn test_fou_parameter_estimation_v1() { const N: usize = 10000; const X0: f64 = 0.0; - let fgn = FGN::new(0.75, 1599, Some(1.0), None); - let fou = FOU::new(3.0, 5.0, 3.0, 1600, Some(X0), Some(16.0), None, fgn); + let fgn = FGN::new(0.70, 4095, Some(1.0), None); + let fou = FOU::new(5.0, 2.8, 1.0, 4096, Some(X0), Some(16.0), None, fgn); let path = fou.sample(); - let mut estimator = FOUParameterEstimation::new(path, 1.0, FilterType::Daubechies); + let mut estimator = FOUParameterEstimationV1::new(path, FilterType::Daubechies); // Estimate the parameters let (estimated_hurst, estimated_sigma, estimated_mu, estimated_theta) = @@ -185,13 +434,60 @@ mod tests { println!("Estimated sigma: {}", estimated_sigma); println!("Estimated mu: {}", estimated_mu); println!("Estimated theta: {}", estimated_theta); + } - // Assert that the estimated parameters are close to the original ones - let tolerance = 0.1; // Adjust tolerance as needed + #[test] + fn test_fou_parameter_estimation_v2() { + const N: usize = 4096; + const X0: f64 = 0.0; + let delta = 1.0 / 256.0; + + let fgn = FGN::new(0.70, N - 1, Some(1.0), None); + let fou = FOU::new(5.0, 2.8, 2.0, N, Some(X0), Some(16.0), None, fgn); + let path = fou.sample(); + let mut estimator = FOUParameterEstimationV2::new(path, delta, N); - assert!((estimated_hurst - 0.75).abs() < tolerance); - assert!((estimated_sigma - 2.0).abs() < tolerance); - assert!((estimated_mu - 3.0).abs() < tolerance); - assert!((estimated_theta - 2.0).abs() < tolerance); + // Estimate the parameters + let (estimated_hurst, estimated_sigma, estimated_mu, estimated_theta) = + estimator.estimate_parameters(); + + // Print the estimated parameters + println!("Estimated Hurst exponent: {}", estimated_hurst); + println!("Estimated sigma: {}", estimated_sigma); + println!("Estimated mu: {}", estimated_mu); + println!("Estimated theta: {}", estimated_theta); + } + + #[test] + fn test_fou_parameter_estimation_v3() { + let series_length = 4096; + let hurst = 0.70; + let sigma = 2.0; + let alpha = 5.0; + let mu = 2.8; + let initial_value = 0.0; + let T = 16.0; + let delta = 1.0 / 256.0; + + let mut estimator = FOUParameterEstimationV3::new( + series_length, + hurst, + sigma, + alpha, + mu, + initial_value, + T, + delta, + ); + + // Estimate the parameters + let (estimated_hurst, estimated_sigma, estimated_mu, estimated_alpha) = + estimator.estimate_parameters(); + + // Print the estimated parameters + println!("Estimated Hurst exponent: {}", estimated_hurst); + println!("Estimated sigma: {}", estimated_sigma); + println!("Estimated mu: {}", estimated_mu); + println!("Estimated alpha: {}", estimated_alpha); } }