diff --git a/Cargo.toml b/Cargo.toml index 39caf906..025b9898 100644 --- a/Cargo.toml +++ b/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 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -51,8 +51,10 @@ rustdoc-args = ["--html-in-header", "katex_header.html", "--cfg", "docsrs"] csv = "1.3.0" # https://docs.rs/csv/latest/csv/ 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/ +ndarray-rand = "0.14.0" num = "0.4.1" # https://docs.rs/num/latest/num/ -num-complex = "0.4.2" # https://docs.rs/num-complex/latest/num_complex/ +num-complex = {version = "0.4.2", features = ["rand"]} # https://docs.rs/num-complex/latest/num_complex/ num-traits = "0.2.16" # https://docs.rs/num-traits/latest/num_traits/ plotters = "0.3.4" # https://docs.rs/plotters/latest/plotters/ rand = "0.8.5" # https://docs.rs/rand/latest/rand/ diff --git a/examples/stochastic_processes.rs b/examples/stochastic_processes.rs index e8a497dd..bb66db54 100644 --- a/examples/stochastic_processes.rs +++ b/examples/stochastic_processes.rs @@ -28,7 +28,7 @@ fn main() { let hl = HoLee::new(0.2, 0.1); let hw = HullWhite::new(0.1, 0.2, 0.1); let ou = OrnsteinUhlenbeck::new(0.05, 0.9, 0.1); - let fbm = FractionalBrownianMotion::new(0.7); + let fbm = FractionalBrownianMotion::new(0.7, FractionalProcessGeneratorMethod::FFT); let mjd = MertonJumpDiffusion::new(0.05, 0.5, 30.0, 0.0, 5.0); let gbb = GeometricBrownianBridge::new(0.05, 0.9, INITIAL_VALUE, END_TIME); let cev = ConstantElasticityOfVariance::new(0.05, 0.9, f64::sin); diff --git a/src/stochastics/fractional_brownian_motion.rs b/src/stochastics/fractional_brownian_motion.rs index 9b8253a2..60aec005 100644 --- a/src/stochastics/fractional_brownian_motion.rs +++ b/src/stochastics/fractional_brownian_motion.rs @@ -9,22 +9,37 @@ use crate::stochastics::{StochasticProcess, Trajectories}; use nalgebra::{DMatrix, DVector, Dim, Dyn, RowDVector}; +use ndarray::{concatenate, prelude::*}; +use ndarray_rand::RandomExt; +use ndrustfft::{ndfft_par, FftHandler}; +use num_complex::{Complex, ComplexDistribution}; 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) } } @@ -35,10 +50,10 @@ impl FractionalBrownianMotion { /// /// Will panic if Hurst parameter is not in [0, 1]. #[must_use] - 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). @@ -78,29 +93,73 @@ impl FractionalBrownianMotion { } /// Fractional Gaussian noise. - #[must_use] - pub fn fgn_cholesky(&self, n: usize, t_n: f64) -> RowDVector { + pub fn fgn_cholesky(&self, n: usize, t_n: f64) -> Vec { let acf_sqrt = self.acf_matrix_sqrt(n); let noise = rand::thread_rng() .sample_iter::(StandardNormal) .take(n) .collect(); let noise = DVector::::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 { + /// Seedable Fractional Gaussian noise. + pub fn seedable_fgn_cholesky(&self, n: usize, t_n: f64, seed: u64) -> Vec { let acf_sqrt = self.acf_matrix_sqrt(n); - let rng = StdRng::seed_from_u64(seed); - let noise = rng + let noise = StdRng::seed_from_u64(seed) .sample_iter::(StandardNormal) .take(n) .collect(); let noise = DVector::::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() + } + + /// Fractional Gaussian noise via FFT. + pub fn fgn_fft(&self, n: usize, t_n: f64) -> Vec { + if !(0.0..=1.0).contains(&self.hurst) { + panic!("Hurst parameter must be between 0 and 1"); + } + let mut r = Array1::linspace(0.0, n as f64, n + 1); + r.par_mapv_inplace(|x| { + if x == 0.0 { + 1.0 + } else { + 0.5 * ((x as f64 + 1.0).powf(2.0 * self.hurst) + - 2.0 * (x as f64).powf(2.0 * self.hurst) + + (x as f64 - 1.0).powf(2.0 * self.hurst)) + } + }); + let r = concatenate( + Axis(0), + #[allow(clippy::reversed_empty_ranges)] + &[r.view(), r.slice(s![..;-1]).slice(s![1..-1]).view()], + ) + .unwrap(); + let mut data = Array1::>::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 sqrt_eigenvalues = Array1::>::zeros(r.len()); + ndfft_par(&data, &mut sqrt_eigenvalues, &mut r_fft, 0); + sqrt_eigenvalues.par_mapv_inplace(|x| Complex::new((x.re / (2.0 * n as f64)).sqrt(), x.im)); + let rnd = Array1::>::random( + 2 * n, + ComplexDistribution::new(StandardNormal, StandardNormal), + ); + let fgn = &sqrt_eigenvalues * &rnd; + let mut fft_handler = FftHandler::new(2 * n); + let mut fgn_fft = Array1::>::zeros(2 * n); + ndfft_par(&fgn, &mut fgn_fft, &mut fft_handler, 0); + let fgn = fgn_fft + .slice(s![1..n + 1]) + .mapv(|x: Complex| (x.re * (n as f64).powf(-self.hurst)) * t_n.powf(self.hurst)); + fgn.to_vec() } } @@ -135,7 +194,10 @@ impl StochasticProcess for FractionalBrownianMotion { let times: Vec = (0..=n_steps).map(|t| t_0 + dt * (t as f64)).collect(); let path_generator = |path: &mut Vec| { - 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] @@ -201,34 +263,64 @@ mod test_fractional_brownian_motion { // use std::time::Instant; use super::*; - use crate::{assert_approx_equal, statistics::*}; + use crate::{ + ml::{Decomposition, LinearRegressionInput}, + statistics::Statistic, + }; - #[test] - fn test_chol() { - let fbm = FractionalBrownianMotion::new(0.7); - let n = 3; - let hurst = 0.7; + fn higuchi_fd(x: &Vec, 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::(StandardNormal) - .take(n) - .collect(); - let noise = DVector::::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 += (x[m + j * k] - x[m + (j - 1) * k]).abs(); + } - let fgn = (acf_matrix * noise).transpose() * (n as f64).powf(-hurst); + ll /= k as f64; + ll *= (n_times - 1) as f64 / (k * n_max) as f64; + lm[m] = ll; + } + + lk[k - 1] = lm.iter().sum::() / k as f64; + x_reg[k - 1] = (1.0 / k as f64).ln(); + y_reg[k - 1] = lk[k - 1].ln(); + } - println!("{:?}", fgn); + 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(); + + regression_result.coefficients[0] + } + + #[test] + 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); + } } #[test] fn test_brownian_motion() { - let fbm = FractionalBrownianMotion::new(0.7); + let fbm = FractionalBrownianMotion::new(0.7, FractionalProcessGeneratorMethod::FFT); 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); @@ -240,8 +332,8 @@ mod test_fractional_brownian_motion { .collect(); // E[X_T] = 0 - assert_approx_equal!(X_T.mean(), 0.0, 0.5); + assert_approx_equal!(X_T.clone().mean(), 0.0, 0.5); // V[X_T] = T - assert_approx_equal!(X_T.variance(), 0.5, 0.5); + assert_approx_equal!(X_T.clone().variance(), 0.5, 0.5); } } diff --git a/src/stochastics/fractional_cox_ingersoll_ross.rs b/src/stochastics/fractional_cox_ingersoll_ross.rs new file mode 100644 index 00000000..adf5252c --- /dev/null +++ b/src/stochastics/fractional_cox_ingersoll_ross.rs @@ -0,0 +1,140 @@ +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// RustQuant: A Rust library for quantitative finance tools. +// Copyright (C) 2023 https://github.com/avhz +// Dual licensed under Apache 2.0 and MIT. +// See: +// - LICENSE-APACHE.md +// - LICENSE-MIT.md +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +use crate::stochastics::*; +use rayon::prelude::*; + +/// Struct containing the Ornstein-Uhlenbeck process parameters. +pub struct FractionalCoxIngersollRoss { + /// The long-run mean ($\mu$). + pub mu: TimeDependent, + + /// The diffusion, or instantaneous volatility ($\sigma$). + pub sigma: TimeDependent, + + /// Mean reversion parameter ($\theta$). + /// Defines the speed at which the process reverts to the long-run mean. + pub theta: TimeDependent, + + /// Hurst parameter of the process. + /// The Hurst parameter is a measure of the long-term memory of the process. + pub hurst: f64, + + /// Method to generate Fractional Gaussian Noise. + pub method: FractionalProcessGeneratorMethod, +} + +impl FractionalCoxIngersollRoss { + /// Create a new Ornstein-Uhlenbeck process. + pub fn new( + mu: impl Into, + sigma: impl Into, + theta: impl Into, + hurst: f64, + method: FractionalProcessGeneratorMethod, + ) -> Self { + assert!((0.0..=1.0).contains(&hurst)); + Self { + mu: mu.into(), + sigma: sigma.into(), + theta: theta.into(), + hurst, + method, + } + } +} + +impl StochasticProcess for FractionalCoxIngersollRoss { + fn drift(&self, x: f64, t: f64) -> f64 { + self.theta.0(t) * (self.mu.0(t) - x) + } + + fn diffusion(&self, x: f64, t: f64) -> f64 { + self.sigma.0(t) * x.sqrt() + } + + fn jump(&self, _x: f64, _t: f64) -> Option { + Some(0.0) + } + + fn euler_maruyama( + &self, + x_0: f64, + t_0: f64, + t_n: f64, + n_steps: usize, + m_paths: usize, + parallel: bool, + ) -> Trajectories { + let fgn = match self.method { + FractionalProcessGeneratorMethod::CHOLESKY => { + let fbm = FractionalBrownianMotion::new( + self.hurst, + FractionalProcessGeneratorMethod::CHOLESKY, + ); + FractionalBrownianMotion::fgn_cholesky(&fbm, n_steps, t_n) + } + FractionalProcessGeneratorMethod::FFT => { + let fbm = FractionalBrownianMotion::new( + self.hurst, + FractionalProcessGeneratorMethod::FFT, + ); + FractionalBrownianMotion::fgn_fft(&fbm, n_steps, t_n) + } + }; + + let dt: f64 = (t_n - t_0) / (n_steps as f64); + + // Initialise empty paths and fill in the time points. + let mut paths = vec![vec![x_0; n_steps + 1]; m_paths]; + let times: Vec = (0..=n_steps).map(|t| t_0 + dt * (t as f64)).collect(); + + let path_generator = |path: &mut Vec| { + for t in 0..n_steps { + path[t + 1] = path[t] + + self.drift(path[t], times[t]) * dt + + self.diffusion(path[t], times[t]) * fgn[t]; + } + }; + + if parallel { + paths.par_iter_mut().for_each(path_generator); + } else { + paths.iter_mut().for_each(path_generator); + } + + Trajectories { times, paths } + } +} + +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// TESTS +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +#[cfg(test)] +mod test_fractional_cir { + use super::*; + + #[test] + #[ignore = "Hard to test."] + fn test_fractional_cir() -> Result<(), Box> { + let fou = FractionalOrnsteinUhlenbeck::new( + 0.15, + 0.45, + 0.01, + 0.7, + FractionalProcessGeneratorMethod::FFT, + ); + + #[allow(dead_code)] + let _output = fou.euler_maruyama(10.0, 0.0, 0.5, 100, 100, false); + + std::result::Result::Ok(()) + } +} diff --git a/src/stochastics/fractional_ornstein_uhlenbeck.rs b/src/stochastics/fractional_ornstein_uhlenbeck.rs index 73a50ba0..269330e8 100644 --- a/src/stochastics/fractional_ornstein_uhlenbeck.rs +++ b/src/stochastics/fractional_ornstein_uhlenbeck.rs @@ -12,6 +12,8 @@ use crate::stochastics::{ }; use rayon::prelude::*; +use super::FractionalProcessGeneratorMethod; + /// Struct containing the Ornstein-Uhlenbeck process parameters. pub struct FractionalOrnsteinUhlenbeck { /// The long-run mean ($\mu$). @@ -27,6 +29,9 @@ pub struct FractionalOrnsteinUhlenbeck { /// Hurst parameter of the process. /// The Hurst parameter is a measure of the long-term memory of the process. pub hurst: f64, + + /// Method to generate Fractional Gaussian Noise. + pub method: FractionalProcessGeneratorMethod, } impl FractionalOrnsteinUhlenbeck { @@ -36,6 +41,7 @@ impl FractionalOrnsteinUhlenbeck { sigma: impl Into, theta: impl Into, hurst: f64, + method: FractionalProcessGeneratorMethod, ) -> Self { assert!((0.0..=1.0).contains(&hurst)); Self { @@ -43,6 +49,7 @@ impl FractionalOrnsteinUhlenbeck { sigma: sigma.into(), theta: theta.into(), hurst, + method, } } } @@ -70,8 +77,22 @@ impl StochasticProcess for FractionalOrnsteinUhlenbeck { m_paths: usize, parallel: bool, ) -> Trajectories { - let fbm = FractionalBrownianMotion::new(self.hurst); - let fgn = FractionalBrownianMotion::fgn_cholesky(&fbm, n_steps, t_n); + let fgn = match self.method { + FractionalProcessGeneratorMethod::CHOLESKY => { + let fbm = FractionalBrownianMotion::new( + self.hurst, + FractionalProcessGeneratorMethod::CHOLESKY, + ); + FractionalBrownianMotion::fgn_cholesky(&fbm, n_steps, t_n) + } + FractionalProcessGeneratorMethod::FFT => { + let fbm = FractionalBrownianMotion::new( + self.hurst, + FractionalProcessGeneratorMethod::FFT, + ); + FractionalBrownianMotion::fgn_fft(&fbm, n_steps, t_n) + } + }; let dt: f64 = (t_n - t_0) / (n_steps as f64); @@ -102,13 +123,19 @@ impl StochasticProcess for FractionalOrnsteinUhlenbeck { // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #[cfg(test)] -mod tests_ornstein_uhlenbeck { +mod tests_fractional_ornstein_uhlenbeck { use super::*; #[test] #[ignore = "Hard to test."] fn test_fractional_ornstein_uhlenbeck() { - let fou = FractionalOrnsteinUhlenbeck::new(0.15, 0.45, 0.01, 0.7); + let fou = FractionalOrnsteinUhlenbeck::new( + 0.15, + 0.45, + 0.01, + 0.7, + FractionalProcessGeneratorMethod::FFT, + ); #[allow(dead_code)] let _output = fou.euler_maruyama(10.0, 0.0, 0.5, 100, 100, false); diff --git a/src/stochastics/mod.rs b/src/stochastics/mod.rs index 26e09223..51c85781 100644 --- a/src/stochastics/mod.rs +++ b/src/stochastics/mod.rs @@ -47,6 +47,7 @@ pub use constant_elasticity_of_variance::*; pub use cox_ingersoll_ross::*; pub use extended_vasicek::*; pub use fractional_brownian_motion::*; +pub use fractional_cox_ingersoll_ross::*; pub use fractional_ornstein_uhlenbeck::*; pub use geometric_brownian_bridge::*; pub use geometric_brownian_motion::*; @@ -70,6 +71,8 @@ pub mod cox_ingersoll_ross; pub mod extended_vasicek; /// Fractional Brownian Motion. pub mod fractional_brownian_motion; +/// Fractional Cox-Ingersoll-Ross process. +pub mod fractional_cox_ingersoll_ross; /// Fractional Ornstein-Uhlenbeck process. pub mod fractional_ornstein_uhlenbeck; /// Geometric brownian bridge process.