feat: add fft support for fractional processes

Cargo.toml

@@ -1,9 +1,9 @@
 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ## RustQuant: A Rust library for quantitative finance tools.
 ## Copyright (C) 2023 https://github.com/avhz
-## Dual licensed under Apache 2.0 and MIT. 
+## Dual licensed under Apache 2.0 and MIT.
 ## See:
-##      - LICENSE-APACHE.md 
+##      - LICENSE-APACHE.md
 ##      - LICENSE-MIT.md
 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -50,6 +50,7 @@ rustdoc-args = ["--html-in-header", "katex_header.html", "--cfg", "docsrs"]
 [dependencies]
 nalgebra = "0.32.2"   # https://docs.rs/nalgebra/latest/nalgebra/
 ndarray = "0.15.6"    # https://docs.rs/ndarray/latest/ndarray/
+ndrustfft = "0.4.1"   # https://docs.rs/ndrustfft/latest/ndrustfft/
 num = "0.4.1"         # https://docs.rs/num/latest/num/
 num-complex = "0.4.2" # https://docs.rs/num-complex/latest/num_complex/
 num-traits = "0.2.16" # https://docs.rs/num-traits/latest/num_traits/

examples/stochastic_processes.rs

@@ -21,7 +21,7 @@ fn main() {
     let hl = HoLee::new(0.2, theta_t);
     let hw = HullWhite::new(0.1, 0.2, theta_t);
     let ou = OrnsteinUhlenbeck::new(0.05, 0.9, 0.1);
-    let fbm = FractionalBrownianMotion::new(0.7);
+    let fbm = FractionalBrownianMotion::new(0.7, FractionalProcessGeneratorMethod::FFT);
 
     // Generate path using Euler-Maruyama scheme.
     // Parameters: x_0, t_0, t_n, n, sims, parallel.

src/stochastics/fractional_brownian_motion.rs

@@ -9,31 +9,43 @@
 
 use crate::stochastics::*;
 use nalgebra::{DMatrix, DVector, Dim, Dyn, RowDVector};
+use ndarray::{concatenate, prelude::*};
+use ndrustfft::{ndfft, FftHandler};
+use num_complex::Complex;
 use rand::Rng;
-#[cfg(feature = "seedable")]
-use rand::{rngs::StdRng, SeedableRng};
 use rand_distr::StandardNormal;
 use rayon::prelude::*;
 
+/// Method used to generate the Fractional Brownian Motion.
+#[derive(Debug)]
+pub enum FractionalProcessGeneratorMethod {
+    /// Chooses the Cholesky decomposition method.
+    CHOLESKY,
+    /// Chooses the Davies-Harte method.
+    FFT,
+}
+
 /// Struct containing the Fractional Brownian Motion parameters.
 #[derive(Debug)]
 pub struct FractionalBrownianMotion {
     /// Hurst parameter of the process.
     pub hurst: f64,
+    /// Method used to generate the process.
+    pub method: FractionalProcessGeneratorMethod,
 }
 
 impl Default for FractionalBrownianMotion {
     fn default() -> Self {
-        Self::new(0.5)
+        Self::new(0.5, FractionalProcessGeneratorMethod::FFT)
     }
 }
 
 impl FractionalBrownianMotion {
     /// Create a new Fractional Brownian Motion process.
-    pub fn new(hurst: f64) -> Self {
+    pub fn new(hurst: f64, method: FractionalProcessGeneratorMethod) -> Self {
         assert!((0.0..=1.0).contains(&hurst));
 
-        Self { hurst }
+        Self { hurst, method }
     }
 
     /// Autocovariance function (ACF).
@@ -73,28 +85,68 @@
     }
 
     /// Fractional Gaussian noise.
-    pub fn fgn_cholesky(&self, n: usize, t_n: f64) -> RowDVector<f64> {
+    pub fn fgn_cholesky(&self, n: usize, t_n: f64) -> Vec<f64> {
         let acf_sqrt = self.acf_matrix_sqrt(n);
         let noise = rand::thread_rng()
             .sample_iter::<f64, StandardNormal>(StandardNormal)
             .take(n)
             .collect();
         let noise = DVector::<f64>::from_vec(noise);
+        let noise = (acf_sqrt * noise).transpose() * (1.0 * t_n / n as f64).powf(self.hurst);
 
-        (acf_sqrt * noise).transpose() * (1.0 * t_n / n as f64).powf(self.hurst)
+        noise.data.as_vec().clone()
     }
 
-    #[cfg(feature = "seedable")]
-    fn seedable_fgn_cholesky(&self, n: usize, t_n: f64, seed: u64) -> RowDVector<f64> {
-        let acf_sqrt = self.acf_matrix_sqrt(n);
-        let rng = StdRng::seed_from_u64(seed);
-        let noise = rng
-            .sample_iter::<f64, StandardNormal>(StandardNormal)
-            .take(n)
-            .collect();
-        let noise = DVector::<f64>::from_vec(noise);
+    /// Fractional Gaussian noise via FFT.
+    pub fn fgn_fft(&self, n: usize, t_n: f64) -> Vec<f64> {
+        if !(0.0..=1.0).contains(&self.hurst) {
+            panic!("Hurst parameter must be between 0 and 1");
+        }
 
-        (acf_sqrt * noise).transpose() * (1.0 * t_n / n as f64).powf(self.hurst)
+        let r = Array::from_shape_fn((n + 1,), |i| {
+            if i == 0 {
+                1.0
+            } else {
+                0.5 * ((i as f64 + 1.0).powf(2.0 * self.hurst)
+                    - 2.0 * (i as f64).powf(2.0 * self.hurst)
+                    + (i as f64 - 1.0).powf(2.0 * self.hurst))
+            }
+        });
+
+        let r = concatenate(
+            Axis(0),
+            &[r.view(), r.slice(s![..;-1]).slice(s![1..-1]).view()],
+        )
+        .unwrap();
+
+        let mut data = Array1::<Complex<f64>>::zeros(r.len());
+        for (i, v) in r.iter().enumerate() {
+            data[i] = Complex::new(*v, 0.0);
+        }
+        let mut r_fft = FftHandler::new(r.len());
+        let mut lambda = Array1::<Complex<f64>>::zeros(r.len());
+        ndfft(&data, &mut lambda, &mut r_fft, 0);
+        let lambda = lambda.mapv(|x| x.re / (2.0 * n as f64)).mapv(f64::sqrt);
+
+        let mut rng = rand::thread_rng();
+        let mut rnd = Array1::<Complex<f64>>::zeros(2 * n);
+        for (_, v) in rnd.iter_mut().enumerate() {
+            let real: f64 = rng.sample(StandardNormal);
+            let imag: f64 = rng.sample(StandardNormal);
+            *v = Complex::new(real, imag);
+        }
+
+        let mut fgn = Array1::<Complex<f64>>::zeros(2 * n);
+        for (i, v) in rnd.iter().enumerate() {
+            fgn[i] = lambda[i] * v;
+        }
+        let mut fgn_fft_handler = FftHandler::new(2 * n);
+        let mut fgn_fft = Array1::<Complex<f64>>::zeros(2 * n);
+        ndfft(&fgn, &mut fgn_fft, &mut fgn_fft_handler, 0);
+
+        let fgn = fgn_fft.slice(s![1..n + 1]).mapv(|x| x.re);
+        fgn.mapv(|x| (x * (n as f64).powf(-self.hurst)) * t_n.powf(self.hurst))
+            .to_vec()
     }
 }
 
@@ -129,7 +181,10 @@
         let times: Vec<f64> = (0..=n_steps).map(|t| t_0 + dt * (t as f64)).collect();
 
         let path_generator = |path: &mut Vec<f64>| {
-            let fgn = self.fgn_cholesky(n_steps, t_n);
+            let fgn = match self.method {
+                FractionalProcessGeneratorMethod::FFT => self.fgn_fft(n_steps, t_n),
+                FractionalProcessGeneratorMethod::CHOLESKY => self.fgn_cholesky(n_steps, t_n),
+            };
 
             for t in 0..n_steps {
                 path[t + 1] = path[t]
@@ -195,49 +250,69 @@
     // use std::time::Instant;
 
     use super::*;
-    use crate::{assert_approx_equal, statistics::*};
+    use crate::ml::{Decomposition, LinearRegressionInput};
 
-    #[test]
-    fn test_chol() {
-        let fbm = FractionalBrownianMotion::new(0.7);
-        let n = 3;
-        let hurst = 0.7;
+    pub(self) fn higuchi_fd(x: &Vec<f64>, kmax: usize) -> f64 {
+        let n_times = x.len();
 
-        let acf_vector = fbm.acf_vector(n);
-        let acf_matrix = fbm.acf_matrix_sqrt(n);
+        let mut lk = vec![0.0; kmax];
+        let mut x_reg = vec![0.0; kmax];
+        let mut y_reg = vec![0.0; kmax];
 
-        println!("ACF vector = {:?}", acf_vector);
-        println!("ACF matrix = {:?}", acf_matrix);
+        for k in 1..=kmax {
+            let mut lm = vec![0.0; k];
 
-        let noise = rand::thread_rng()
-            .sample_iter::<f64, StandardNormal>(StandardNormal)
-            .take(n)
-            .collect();
-        let noise = DVector::<f64>::from_vec(noise);
+            for m in 0..k {
+                let mut ll = 0.0;
+                let n_max = ((n_times - m - 1) as f64 / k as f64).floor() as usize;
+
+                for j in 1..n_max {
+                    ll += num::abs(x[m + j * k] - x[m + (j - 1) * k]);
+                }
+
+                ll /= k as f64;
+                ll *= (n_times - 1) as f64 / (k * n_max) as f64;
+                lm[m] = ll;
+            }
 
-        let fgn = (acf_matrix * noise).transpose() * (n as f64).powf(-hurst);
+            lk[k - 1] = lm.iter().sum::<f64>() / k as f64;
+            x_reg[k - 1] = (1.0 / k as f64).ln();
+            y_reg[k - 1] = lk[k - 1].ln();
+        }
+
+        let x_reg = DMatrix::from_vec(kmax, 1, x_reg);
+        let y_reg = DVector::from_vec(y_reg);
+        let linear_regression = LinearRegressionInput::new(x_reg, y_reg);
+        let regression_result = linear_regression.fit(Decomposition::None).unwrap();
 
-        println!("{:?}", fgn);
+        regression_result.coefficients[0]
     }
 
     #[test]
-    fn test_brownian_motion() -> Result<(), Box<dyn std::error::Error>> {
-        let fbm = FractionalBrownianMotion::new(0.7);
-        let output_serial = fbm.euler_maruyama(0.0, 0.0, 0.5, 100, 1000, false);
-        // let output_parallel = (&bm).euler_maruyama(10.0, 0.0, 0.5, 100, 10, true);
-
-        // Test the distribution of the final values.
-        let X_T: Vec<f64> = output_serial
-            .paths
-            .iter()
-            .filter_map(|v| v.last().cloned())
-            .collect();
-
-        // E[X_T] = 0
-        assert_approx_equal!(X_T.mean(), 0.0, 0.5);
-        // V[X_T] = T
-        assert_approx_equal!(X_T.variance(), 0.5, 0.5);
+    fn test_chol() {
+        let fbm = FractionalBrownianMotion::new(0.7, FractionalProcessGeneratorMethod::FFT);
+        let hursts = vec![0.1, 0.3, 0.5, 0.7, 0.9];
+
+        for hurst in hursts {
+            let fbm = FractionalBrownianMotion::fgn_cholesky(&fbm, 2000, 1.0);
+            let higuchi_fd = higuchi_fd(&fbm.to_vec(), 10);
+            let est_hurst = 2.0 - higuchi_fd;
+            print!("hurst: {}, higuchi_fd: {}\n", hurst, est_hurst);
+            assert!(est_hurst - hurst < 0.05);
+        }
+    }
 
-        std::result::Result::Ok(())
+    #[test]
+    fn test_fft() {
+        let fbm = FractionalBrownianMotion::new(0.7, FractionalProcessGeneratorMethod::FFT);
+        let hursts = vec![0.1, 0.3, 0.5, 0.7, 0.9];
+
+        for hurst in hursts {
+            let fbm = FractionalBrownianMotion::fgn_fft(&fbm, 2000, 1.0);
+            let higuchi_fd = higuchi_fd(&fbm.to_vec(), 10);
+            let est_hurst = 2.0 - higuchi_fd;
+            print!("hurst: {}, higuchi_fd: {}\n", hurst, est_hurst);
+            assert!(est_hurst - hurst < 0.05);
+        }
     }
 }