market_model/process/

bates.rs

use crate::process::StochasticProcess;
use crate::process::heston::HestonProcess;
use crate::simulatable::Simulatable;
use ndarray::{Array1, Array2, array};
use rand::Rng;
use rand_distr::{Poisson, StandardNormal};

/// Bates model: Heston stochastic volatility + Merton jump diffusion.
///
/// Same as Heston, but with compound Poisson jumps on the price.
#[derive(Clone)]
pub struct BatesProcess {
    heston: HestonProcess,
    pub lambda: f64,
    pub mu_jump: f64,
    pub sigma_jump: f64,
}

impl BatesProcess {
    #[allow(clippy::too_many_arguments)]
    pub fn new(
        mu: f64,
        kappa: f64,
        theta: f64,
        sigma: f64,
        rho: f64,
        lambda: f64,
        mu_jump: f64,
        sigma_jump: f64,
        initial_price: f64,
        initial_variance: f64,
    ) -> Self {
        Self {
            heston: HestonProcess::new(
                mu,
                kappa,
                theta,
                sigma,
                rho,
                initial_price,
                initial_variance,
            ),
            lambda,
            mu_jump,
            sigma_jump,
        }
    }

    pub fn price(&self) -> f64 {
        self.heston.price()
    }

    pub fn volatility(&self) -> f64 {
        self.heston.volatility()
    }

    pub fn mu(&self) -> f64 {
        self.heston.mu
    }

    pub fn kappa(&self) -> f64 {
        self.heston.kappa
    }

    pub fn theta(&self) -> f64 {
        self.heston.theta
    }

    pub fn sigma(&self) -> f64 {
        self.heston.sigma
    }

    pub fn rho(&self) -> f64 {
        self.heston.rho
    }
}

impl Simulatable for BatesProcess {
    type Output = (f64, f64); // (price, variance)

    fn dim(&self) -> usize {
        2
    }

    fn step(&mut self, dt: f64, dw: &[f64], rng: &mut impl Rng) -> (f64, f64) {
        // 1. Heston step (disambiguate via UFCS — both Simulatable and StochasticProcess have step)
        let (_price, _variance) = Simulatable::step(&mut self.heston, dt, dw, rng);

        // 2. Apply jumps to the price
        if self.lambda > 0.0 {
            let p_jump = self.lambda * dt;
            if rng.random::<f64>() < p_jump {
                let z: f64 = rng.sample(StandardNormal);
                let log_jump = self.mu_jump + self.sigma_jump * z;
                let multiplier = log_jump.exp();
                self.heston.price *= multiplier;
            }
        }

        self.heston.current()
    }

    fn current(&self) -> (f64, f64) {
        self.heston.current()
    }
}

impl StochasticProcess for BatesProcess {
    type State = Array1<f64>;
    type Diffusion = Array2<f64>;
    type BrownianIncrement = Array1<f64>;

    fn drift(&self, x: &Array1<f64>, t: f64) -> Array1<f64> {
        self.heston.drift(x, t)
    }

    fn diffusion(&self, x: &Array1<f64>, t: f64) -> Array2<f64> {
        self.heston.diffusion(x, t)
    }

    fn jump_intensity(&self, _x: &Array1<f64>, _t: f64) -> f64 {
        self.lambda
    }

    fn jump_size<R: Rng + ?Sized>(&self, x: &Array1<f64>, _t: f64, rng: &mut R) -> Array1<f64> {
        if self.lambda <= 0.0 {
            return Array1::zeros(2);
        }
        let z: f64 = rng.sample(StandardNormal);
        let log_jump = self.mu_jump + self.sigma_jump * z;
        let jump_multiplier = log_jump.exp();
        let s = x[0];
        let jump_amount = s * (jump_multiplier - 1.0);
        array![jump_amount, 0.0]
    }

    fn expected_value(&self, x0: &Array1<f64>, t: f64) -> Array1<f64> {
        let s0 = x0[0];
        let v0 = x0[1];
        let k = (self.mu_jump + 0.5 * self.sigma_jump.powi(2)).exp() - 1.0;
        let total_drift = self.heston.mu + self.lambda * k;
        let expected_s = s0 * (total_drift * t).exp();
        let exp_minus_kappa_t = (-self.heston.kappa * t).exp();
        let expected_v = v0 * exp_minus_kappa_t + self.heston.theta * (1.0 - exp_minus_kappa_t);
        array![expected_s, expected_v]
    }

    fn variance(&self, x0: &Array1<f64>, t: f64) -> Array1<f64> {
        self.heston.variance(x0, t)
    }

    fn step<R: Rng + ?Sized>(
        &self,
        x: &Array1<f64>,
        t: f64,
        dt: f64,
        dw: &Array1<f64>,
        rng: &mut R,
    ) -> Array1<f64> {
        let drift = self.drift(x, t);
        let diffusion = self.diffusion(x, t);
        let diffusion_term = diffusion.dot(dw);
        let mut next_x = x + &(&drift * dt) + &diffusion_term;
        if next_x[1] < 0.0 {
            next_x[1] = 0.0;
        }
        if self.lambda > 0.0 {
            let poisson = Poisson::new(self.lambda * dt).unwrap();
            let num_jumps = rng.sample(poisson) as u64;
            for _ in 0..num_jumps {
                let jump = self.jump_size(&next_x, t, rng);
                next_x = next_x + jump;
            }
        }
        next_x
    }

    fn generate_increment<R: Rng + ?Sized>(&self, rng: &mut R, dt: f64) -> Array1<f64> {
        self.heston.generate_increment(rng, dt)
    }

    fn simulate<R: Rng + ?Sized>(
        &self,
        x0: Array1<f64>,
        t0: f64,
        t_end: f64,
        dt: f64,
        rng: &mut R,
    ) -> Vec<(f64, Array1<f64>)> {
        let n_steps = ((t_end - t0) / dt).ceil() as usize;
        let mut path = Vec::with_capacity(n_steps + 1);
        path.push((t0, x0.clone()));
        let mut x = x0;
        let mut t = t0;
        for _ in 0..n_steps {
            let dw = self.generate_increment(rng, dt);
            x = self.step(&x, t, dt, &dw, rng);
            t += dt;
            path.push((t, x.clone()));
        }
        path
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_bates_no_jumps() {
        let b = BatesProcess::new(0.05, 2.0, 0.04, 0.3, -0.7, 0.0, 0.0, 0.0, 100.0, 0.04);
        assert_eq!(b.dim(), 2);
        let (p, v) = b.current();
        assert_eq!(p, 100.0);
        assert_eq!(v, 0.04);
    }
}