market_model/process/

heston.rs

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

/// Heston stochastic volatility model.
///
/// $$dS_t = mu * S_t * dt + sqrt(v_t) * S_t * dW_t^S$$
/// $$dv_t = kappa * (theta - v_t) * dt + sigma * sqrt(v_t) * dW_t^v$$
///
/// with correlation `rho` between the two Brownian motions.
#[derive(Clone)]
pub struct HestonProcess {
    pub mu: f64,
    pub kappa: f64,
    pub theta: f64,
    pub sigma: f64,
    pub rho: f64,
    pub price: f64,
    pub variance: f64,
    /// Precomputed sqrt(1 - rho^2) for Cholesky.
    sqrt_one_minus_rho2: f64,
}

impl HestonProcess {
    #[allow(clippy::too_many_arguments)]
    pub fn new(
        mu: f64,
        kappa: f64,
        theta: f64,
        sigma: f64,
        rho: f64,
        initial_price: f64,
        initial_variance: f64,
    ) -> Self {
        let sqrt_one_minus_rho2 = if rho.abs() < 1.0 {
            (1.0 - rho * rho).sqrt()
        } else {
            0.0
        };
        Self {
            mu,
            kappa,
            theta,
            sigma,
            rho,
            price: initial_price,
            variance: initial_variance,
            sqrt_one_minus_rho2,
        }
    }

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

    pub fn volatility(&self) -> f64 {
        self.variance.max(0.0).sqrt()
    }

    /// Closed-form expected variance.
    pub fn expected_variance(&self, t: f64) -> f64 {
        let exp = (-self.kappa * t).exp();
        self.variance * exp + self.theta * (1.0 - exp)
    }

    /// Closed-form variance of variance.
    pub fn variance_of_variance(&self, t: f64) -> f64 {
        let e1 = (-self.kappa * t).exp();
        let e2 = (-2.0 * self.kappa * t).exp();
        self.variance * (self.sigma.powi(2) / self.kappa) * (e1 - e2)
            + (self.theta * self.sigma.powi(2) / (2.0 * self.kappa)) * (1.0 - e1).powi(2)
    }
}

impl Simulatable for HestonProcess {
    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) {
        // Cholesky decomposition of the correlation matrix:
        // dW^S = dW_1
        // dW^v = rho * dW_1 + sqrt(1 - rho^2) * dW_2
        let z1 = dw[0];
        let z2 = if self.rho.abs() < 1.0 {
            self.rho * dw[0] + self.sqrt_one_minus_rho2 * dw[1]
        } else {
            self.rho.signum() * dw[0]
        };

        let v = self.variance.max(0.0);
        let sqrt_v = v.sqrt();

        // Price SDE
        self.price *= (self.mu * dt + sqrt_v * z1).exp();

        // Variance SDE (full truncation)
        self.variance += self.kappa * (self.theta - v) * dt + self.sigma * sqrt_v * z2;
        if self.variance < 0.0 {
            self.variance = 0.0;
        }

        (self.price, self.variance)
    }

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

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

    fn drift(&self, x: &Array1<f64>, _t: f64) -> Array1<f64> {
        let s = x[0];
        let v = x[1];
        array![self.mu * s, self.kappa * (self.theta - v)]
    }

    fn diffusion(&self, x: &Array1<f64>, _t: f64) -> Array2<f64> {
        let s = x[0];
        let v = if x[1] < 0.0 { 0.0 } else { x[1] };
        let sqrt_v = v.sqrt();
        let rho = self.rho;
        let sqrt_one_minus_rho2 = (1.0 - rho * rho).sqrt();

        array![
            [sqrt_v * s * sqrt_one_minus_rho2, sqrt_v * s * rho],
            [0.0, self.sigma * sqrt_v]
        ]
    }

    fn jump_size<R: Rng + ?Sized>(&self, x: &Array1<f64>, _t: f64, _rng: &mut R) -> Array1<f64> {
        Array1::zeros(x.len())
    }

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

    fn variance(&self, x0: &Array1<f64>, t: f64) -> Array1<f64> {
        let v0 = x0[1];
        let var_s = 0.0;
        let exp_minus_kappa_t = (-self.kappa * t).exp();
        let exp_minus_2_kappa_t = (-2.0 * self.kappa * t).exp();
        let var_v =
            v0 * (self.sigma.powi(2) / self.kappa) * (exp_minus_kappa_t - exp_minus_2_kappa_t)
                + (self.theta * self.sigma.powi(2) / (2.0 * self.kappa))
                    * (1.0 - exp_minus_kappa_t).powi(2);
        array![var_s, var_v]
    }

    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;
        }
        next_x
    }

    fn generate_increment<R: Rng + ?Sized>(&self, rng: &mut R, dt: f64) -> Array1<f64> {
        let dw1: f64 = rng.sample(StandardNormal);
        let dw2: f64 = rng.sample(StandardNormal);
        let dw = array![dw1, dw2];
        dw * dt.sqrt()
    }

    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);
        let mut t = t0;
        let mut x = x0;
        path.push((t, x.clone()));
        let sqrt_dt = dt.sqrt();
        for _ in 0..n_steps {
            let dw1: f64 = rng.sample(StandardNormal);
            let dw2: f64 = rng.sample(StandardNormal);
            let dw = array![dw1, dw2] * sqrt_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_heston_dim() {
        let h = HestonProcess::new(0.05, 2.0, 0.04, 0.3, -0.7, 100.0, 0.04);
        assert_eq!(h.dim(), 2);
    }

    #[test]
    fn test_heston_initial_state() {
        let h = HestonProcess::new(0.05, 2.0, 0.04, 0.3, -0.7, 100.0, 0.04);
        let (p, v) = h.current();
        assert_eq!(p, 100.0);
        assert_eq!(v, 0.04);
    }

    #[test]
    fn test_heston_variance_stays_nonnegative() {
        let mut h = HestonProcess::new(0.05, 1.0, 0.01, 0.5, -0.7, 100.0, 0.01);
        for _ in 0..1000 {
            Simulatable::step(&mut h, 0.001, &[0.0, 0.0], &mut rand::rng());
            let (_, v) = h.current();
            assert!(v >= 0.0);
        }
    }
}