market_model/process/
jump_diffusion.rs

1use crate::process::StochasticProcess;
2use crate::simulatable::Simulatable;
3use rand::Rng;
4use rand_distr::StandardNormal;
5
6/// Jump diffusion process adding compound Poisson jumps to GBM.
7///
8/// $$dS_t = mu * S_t * dt + sigma * S_t * dW_t + dJ_t$$
9/// $$dJ_t = (e^Z - 1) dN_t,  Z ~ N(mu_jump, sigma_jump^2), N_t ~ Poisson(lambda * t)$$
10#[derive(Clone)]
11pub struct JumpDiffusion {
12    pub mu: f64,
13    pub sigma: f64,
14    pub lambda: f64,
15    pub mu_jump: f64,
16    pub sigma_jump: f64,
17    state: f64,
18}
19
20impl JumpDiffusion {
21    pub fn new(
22        mu: f64,
23        sigma: f64,
24        lambda: f64,
25        mu_jump: f64,
26        sigma_jump: f64,
27        initial_price: f64,
28    ) -> Self {
29        Self {
30            mu,
31            sigma,
32            lambda,
33            mu_jump,
34            sigma_jump,
35            state: initial_price,
36        }
37    }
38
39    pub fn value(&self) -> f64 {
40        self.state
41    }
42}
43
44impl Simulatable for JumpDiffusion {
45    type Output = f64;
46
47    fn dim(&self) -> usize {
48        1
49    }
50
51    fn step(&mut self, dt: f64, dw: &[f64], rng: &mut impl Rng) -> f64 {
52        // GBM step: S_{t+dt} = S_t * exp((mu - sigma^2/2) * dt + sigma * dW)
53        let drift_log = (self.mu - 0.5 * self.sigma.powi(2)) * dt;
54        let diff_log = self.sigma * dw[0];
55        self.state *= (drift_log + diff_log).exp();
56
57        // Apply compound Poisson jumps.
58        if self.lambda > 0.0 {
59            let p_jump = self.lambda * dt;
60            if rng.random::<f64>() < p_jump {
61                let z: f64 = rng.sample(StandardNormal);
62                let log_jump = self.mu_jump + self.sigma_jump * z;
63                self.state *= log_jump.exp();
64            }
65        }
66        self.state
67    }
68
69    fn current(&self) -> f64 {
70        self.state
71    }
72}
73
74impl StochasticProcess for JumpDiffusion {
75    type State = f64;
76    type Diffusion = f64;
77    type BrownianIncrement = f64;
78
79    fn drift(&self, x: &f64, _t: f64) -> f64 {
80        self.mu * x
81    }
82
83    fn diffusion(&self, x: &f64, _t: f64) -> f64 {
84        self.sigma * x
85    }
86
87    fn jump_intensity(&self, _x: &f64, _t: f64) -> f64 {
88        self.lambda
89    }
90
91    fn jump_size<R: Rng + ?Sized>(&self, _x: &f64, _t: f64, rng: &mut R) -> f64 {
92        let z: f64 = rng.sample(StandardNormal);
93        self.mu_jump + self.sigma_jump * z
94    }
95
96    fn expected_value(&self, x0: &f64, t: f64) -> f64 {
97        let drift_gbm = x0 * (self.mu * t).exp();
98        let jump_drift = (self.mu_jump.exp() - 1.0) * self.lambda * t;
99        drift_gbm * (1.0 + jump_drift)
100    }
101
102    fn variance(&self, x0: &f64, t: f64) -> f64 {
103        let sigma_sq_t = self.sigma.powi(2) * t;
104        let jump_var = self.lambda * t * (self.sigma_jump.powi(2) + self.mu_jump.powi(2));
105        x0.powi(2) * (sigma_sq_t + jump_var)
106    }
107
108    fn step<R: Rng + ?Sized>(&self, x: &f64, t: f64, dt: f64, dw: &f64, rng: &mut R) -> f64 {
109        let drift_log = (self.mu - 0.5 * self.sigma.powi(2)) * dt;
110        let diff_log = self.sigma * dw;
111        let mut next_x = x * (drift_log + diff_log).exp();
112
113        if self.lambda > 0.0 {
114            let p = self.lambda * dt;
115            if rng.random::<f64>() < p {
116                let jump = self.jump_size(x, t, rng);
117                next_x *= jump.exp();
118            }
119        }
120        next_x
121    }
122
123    fn generate_increment<R: Rng + ?Sized>(&self, rng: &mut R, dt: f64) -> f64 {
124        let dw: f64 = rng.sample(StandardNormal);
125        dw * dt.sqrt()
126    }
127
128    fn simulate<R: Rng + ?Sized>(
129        &self,
130        x0: f64,
131        t0: f64,
132        t_end: f64,
133        dt: f64,
134        rng: &mut R,
135    ) -> Vec<(f64, f64)> {
136        let n_steps = ((t_end - t0) / dt).ceil() as usize;
137        let mut path = Vec::with_capacity(n_steps + 1);
138        let mut t = t0;
139        let mut x = x0;
140        path.push((t, x));
141
142        for _ in 0..n_steps {
143            let dw = self.generate_increment(rng, dt);
144            x = self.step(&x, t, dt, &dw, rng);
145            t += dt;
146            path.push((t, x));
147        }
148        path
149    }
150}
151
152#[cfg(test)]
153mod tests {
154    use super::*;
155
156    #[test]
157    fn test_jump_diffusion_no_jumps() {
158        let jd = JumpDiffusion::new(0.0, 0.2, 0.0, 0.0, 0.0, 100.0);
159        assert_eq!(jd.dim(), 1);
160        assert_eq!(jd.current(), 100.0);
161    }
162
163    #[test]
164    fn test_jump_diffusion_negative_jumps_reduce_price() {
165        let mut jd = JumpDiffusion::new(
166            0.0,   // mu: zero drift
167            0.0,   // sigma: zero diffusion
168            10.0,  // lambda: 10 jumps per year
169            -0.5,  // mu_jump: negative average log-jump
170            0.1,   // sigma_jump
171            100.0, // initial price
172        );
173        let initial = jd.current();
174        for _ in 0..200 {
175            Simulatable::step(&mut jd, 0.01, &[0.0], &mut rand::rng());
176        }
177        assert!(jd.current() < initial);
178    }
179}
market_model/process/jump_diffusion.rs

market_model/process/
jump_diffusion.rs
