From 6830458f75ecfbe5e9d029cb0ec63267d4cecfc1 Mon Sep 17 00:00:00 2001 From: Daniel Boros Date: Wed, 21 Feb 2024 02:19:18 +0100 Subject: [PATCH] feat: stabilze jumps and fix compound poisson --- Cargo.toml | 1 - README.md | 15 ++++++++---- src/jumps/bates.rs | 47 ++++++++++++++++++++++++++++--------- src/jumps/ig.rs | 40 +++++++++++++++++++++++++++++++ src/jumps/levy_diffusion.rs | 39 ++++++++++++++++++++++++++++++ src/jumps/merton.rs | 45 +++++++++++++++++++++++++++++++++++ src/jumps/mod.rs | 3 +++ src/jumps/vg.rs | 10 +++----- src/processes/poisson.rs | 45 ++++++++++++++++++----------------- 9 files changed, 201 insertions(+), 44 deletions(-) create mode 100644 src/jumps/ig.rs create mode 100644 src/jumps/levy_diffusion.rs create mode 100644 src/jumps/merton.rs diff --git a/Cargo.toml b/Cargo.toml index 460b16f..cb1821c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -22,7 +22,6 @@ indicatif = "0.17.7" plotly = "0.8.4" ndarray-rand = "0.14.0" ndrustfft = "0.4.2" -csv = "1.3.0" [dev-dependencies] diff --git a/README.md b/README.md index c51c3c3..b800318 100644 --- a/README.md +++ b/README.md @@ -21,11 +21,21 @@ Minimal documentation is available at [stochastic-rs](https://docs.rs/stochastic - [x] Cox-Ingersoll-Ross process - [x] Ornstein-Uhlenbeck process - [x] Jacobi process -- [x] Variance Gamma model + +# Jumps and Levy processes +- [x] Poisson process +- [x] Compound Poisson process +- [x] Levy jump diffusion +- [x] Inverse Gaussian +- [x] Normal Inverse Gaussian +- [x] Variance Gamma # Stochastic models - [x] Heston model +- [x] Merton model +- [x] Bates model + # Fractional Stochastic processes - [x] Fractional Gaussian noise @@ -40,8 +50,6 @@ Minimal documentation is available at [stochastic-rs](https://docs.rs/stochastic - [ ] Rough Heston model - [ ] Bergomi model - [ ] Rough Bergomi model -- [ ] Merton model -- [ ] Bates model - [ ] SABR model - [ ] Hull-White model - [ ] Barndorff-Nielsen & Shephard model @@ -57,7 +65,6 @@ Minimal documentation is available at [stochastic-rs](https://docs.rs/stochastic - [ ] Vasicek model - [ ] Duffie-Kan model - [ ] Jump-diffusion model -- [ ] Normal Inverse Gaussian model ## Future work - [ ] Add more tests diff --git a/src/jumps/bates.rs b/src/jumps/bates.rs index 5115f58..f88bd0f 100644 --- a/src/jumps/bates.rs +++ b/src/jumps/bates.rs @@ -3,7 +3,6 @@ use ndarray::Array1; use crate::prelude::{correlated::correlated_bms, poisson::compound_poisson}; /// Bates (1996) model -// TODO: under development #[allow(clippy::too_many_arguments)] pub fn bates_1996( mu: f64, @@ -11,7 +10,7 @@ pub fn bates_1996( theta: f64, eta: f64, rho: f64, - lambda: usize, + lambda: f64, n: usize, s0: Option, v0: Option, @@ -23,23 +22,17 @@ pub fn bates_1996( let mut s = Array1::::zeros(n); let mut v = Array1::::zeros(n); - - let [jump_times, jumps] = compound_poisson(n, lambda, None, Some(t.unwrap_or(1.0)), None, None); + let z = compound_poisson(n, lambda, None, t, None, None); + println!("{:?}", z.len()); s[0] = s0.unwrap_or(0.0); v[0] = v0.unwrap_or(0.0); for i in 1..n { - // find the index of the last jump before time t_i - let mut j = 0; - while jump_times[j] < i as f64 * dt { - j += 1; - } - s[i] = s[i - 1] + mu * s[i - 1] * dt + s[i - 1] * v[i - 1].sqrt() * correlated_bms[0][i - 1] - + jumps[j]; + + z[i - 1]; let random: f64 = match use_sym.unwrap_or(false) { true => eta * (v[i]).abs().sqrt() * correlated_bms[1][i - 1], @@ -50,3 +43,35 @@ pub fn bates_1996( [s.to_vec(), v.to_vec()] } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_bates_1996() { + let mu = 0.0; + let kappa = 0.0; + let theta = 0.0; + let eta = 0.0; + let rho = 0.0; + let lambda = 1.0; + let n = 1000; + let t = 10.0; + let b = bates_1996( + mu, + kappa, + theta, + eta, + rho, + lambda, + n, + None, + None, + Some(t), + None, + ); + assert_eq!(b[0].len(), n); + assert_eq!(b[1].len(), n); + } +} diff --git a/src/jumps/ig.rs b/src/jumps/ig.rs new file mode 100644 index 0000000..deaa5d5 --- /dev/null +++ b/src/jumps/ig.rs @@ -0,0 +1,40 @@ +use crate::noises::gn::gn; + +use ndarray::Array1; +use ndarray_rand::{rand_distr::InverseGaussian, RandomExt}; + +pub fn ig(gamma: f64, n: usize, x0: Option, t: Option) -> Vec { + let dt = t.unwrap_or(1.0) / n as f64; + let gn = gn(n - 1, t); + let mut ig = Array1::zeros(n); + ig[0] = x0.unwrap_or(0.0); + + for i in 1..n { + ig[i] = ig[i - 1] + gamma * dt + gn[i - 1] + } + + ig.to_vec() +} + +pub fn nig( + theta: f64, + sigma: f64, + kappa: f64, + n: usize, + x0: Option, + t: Option, +) -> Vec { + let dt = t.unwrap_or(1.0) / n as f64; + let scale = dt.powf(2.0) / kappa; + let mean = dt / scale; + let ig = Array1::random(n - 1, InverseGaussian::new(mean, scale).unwrap()); + let gn = gn(n - 1, t); + let mut nig = Array1::zeros(n); + nig[0] = x0.unwrap_or(0.0); + + for i in 1..n { + nig[i] = nig[i - 1] + theta * ig[i - 1] + sigma * ig[i - 1].sqrt() * gn[i - 1] + } + + nig.to_vec() +} diff --git a/src/jumps/levy_diffusion.rs b/src/jumps/levy_diffusion.rs new file mode 100644 index 0000000..d810bbd --- /dev/null +++ b/src/jumps/levy_diffusion.rs @@ -0,0 +1,39 @@ +use crate::{noises::gn::gn, processes::poisson::compound_poisson}; +use ndarray::Array1; + +pub fn levy_diffusion( + gamma: f64, + sigma: f64, + lambda: f64, + n: usize, + x0: Option, + t: Option, +) -> Vec { + let dt = t.unwrap_or(1.0) / n as f64; + let mut levy = Array1::::zeros(n); + levy[0] = x0.unwrap_or(0.0); + let gn = gn(n - 1, t); + let z = compound_poisson(n, lambda, None, t, None, None); + + for i in 1..n { + levy[i] = levy[i - 1] + gamma * dt + sigma * gn[i - 1] * z[i - 1]; + } + + levy.to_vec() +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_levy_diffusion() { + let gamma = 0.0; + let sigma = 1.0; + let lambda = 10.0; + let n = 1000; + let t = 10.0; + let l = levy_diffusion(gamma, sigma, lambda, n, None, Some(t)); + assert_eq!(l.len(), n); + } +} diff --git a/src/jumps/merton.rs b/src/jumps/merton.rs new file mode 100644 index 0000000..734dafc --- /dev/null +++ b/src/jumps/merton.rs @@ -0,0 +1,45 @@ +use ndarray::Array1; + +use crate::{noises::gn::gn, processes::poisson::compound_poisson}; + +pub fn merton( + alpha: f64, + sigma: f64, + lambda: f64, + theta: f64, + n: usize, + x0: Option, + t: Option, +) -> Vec { + let dt = t.unwrap_or(1.0) / n as f64; + let mut merton = Array1::::zeros(n); + merton[0] = x0.unwrap_or(0.0); + let gn = gn(n - 1, t); + let z = compound_poisson(n, lambda, None, t, None, None); + + for i in 1..n { + merton[i] = merton[i - 1] + + (alpha * sigma.powf(2.0) / 2.0 - lambda * theta) * dt + + sigma * gn[i - 1] + + z[i - 1]; + } + + merton.to_vec() +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_merton() { + let alpha = 0.0; + let sigma = 1.0; + let lambda = 10.0; + let theta = 0.0; + let n = 1000; + let t = 10.0; + let m = merton(alpha, sigma, lambda, theta, n, None, Some(t)); + assert_eq!(m.len(), n); + } +} diff --git a/src/jumps/mod.rs b/src/jumps/mod.rs index c3c51a5..59bad73 100644 --- a/src/jumps/mod.rs +++ b/src/jumps/mod.rs @@ -1,2 +1,5 @@ pub mod bates; +pub mod ig; +pub mod levy_diffusion; +pub mod merton; pub mod vg; diff --git a/src/jumps/vg.rs b/src/jumps/vg.rs index 76cb041..cd19c4b 100644 --- a/src/jumps/vg.rs +++ b/src/jumps/vg.rs @@ -1,11 +1,10 @@ use crate::noises::gn; use ndarray::Array1; -use rand::Rng; +use ndarray_rand::RandomExt; use rand_distr::Gamma; pub fn vg(mu: f64, sigma: f64, nu: f64, n: usize, x0: Option, t: Option) -> Vec { let dt = t.unwrap_or(1.0) / n as f64; - let rng = rand::thread_rng(); let shape = dt / nu; let scale = nu; @@ -13,11 +12,8 @@ pub fn vg(mu: f64, sigma: f64, nu: f64, n: usize, x0: Option, t: Option::zeros(n); vg[0] = x0.unwrap_or(0.0); - let gn = gn::gn(n - 1, Some(t.unwrap_or(1.0))); - let gammas = rng - .sample_iter(Gamma::new(shape, scale).unwrap()) - .take(n - 1) - .collect::>(); + let gn = gn::gn(n - 1, t); + let gammas = Array1::random(n - 1, Gamma::new(shape, scale).unwrap()); for i in 1..n { vg[i] = vg[i - 1] + mu * gammas[i - 1] + sigma * gammas[i - 1].sqrt() * gn[i - 1]; diff --git a/src/processes/poisson.rs b/src/processes/poisson.rs index f000009..16080d0 100644 --- a/src/processes/poisson.rs +++ b/src/processes/poisson.rs @@ -1,6 +1,8 @@ use ndarray::Array1; -use rand::{thread_rng, Rng}; -use rand_distr::{Distribution, Exp, Normal}; +use ndarray_rand::rand_distr::{Normal, Poisson}; +use ndarray_rand::RandomExt; +use rand::thread_rng; +use rand_distr::{Distribution, Exp}; pub fn poisson(n: usize, lambda: usize, t_max: Option) -> Vec { if n == 0 || lambda == 0 { @@ -25,34 +27,36 @@ pub fn poisson(n: usize, lambda: usize, t_max: Option) -> Vec { pub fn compound_poisson( n: usize, - lambda: usize, + lambda: f64, jumps: Option>, t_max: Option, jump_mean: Option, jump_std: Option, -) -> [Vec; 2] { +) -> Vec { if n == 0 { panic!("n must be a positive integer"); } let _t_max = t_max.unwrap_or(1.0); - let times = poisson(n, lambda, Some(_t_max)); - let mut cp = Array1::::zeros(times.len()); - let jumps = match jumps { + let p = Array1::random(n, Poisson::new(lambda).unwrap()); + let mut cp = Array1::::zeros(n); + + match jumps { Some(jumps) => jumps, None => { - let _jump_mean = jump_mean.unwrap_or(1.0); - let _jump_std = jump_std.unwrap_or(0.0); - let norm = Normal::new(_jump_mean, _jump_std).unwrap(); - thread_rng().sample_iter(norm).take(n).collect() - } - }; + let _jump_mean = jump_mean.unwrap_or(0.0); + let _jump_std = jump_std.unwrap_or(1.0); - for i in 1..times.len() { - cp[i] = cp[i - 1] + jumps[i - 1]; - } + for i in 0..n { + for j in &p { + let norm = Array1::random(*j as usize, Normal::new(_jump_mean, _jump_std).unwrap()); + cp[i] = norm.sum(); + } + } - [times, cp.to_vec()] + cp.to_vec() + } + } } #[cfg(test)] @@ -71,10 +75,9 @@ mod tests { #[test] fn test_compound_poisson() { let n = 1000; - let lambda = 10; + let lambda = 2.0; let t = 10.0; - let cp = compound_poisson(n, lambda, None, Some(t), Some(1.0), Some(0.0)); - println!("{:?}", cp[0].len()); - println!("{:?}", cp[1].len()); + let cp = compound_poisson(n, lambda, None, Some(t), None, None); + assert_eq!(cp.len(), n); } }