market_model/process/
cir.rs

1use crate::process::StochasticProcess;
2use crate::simulatable::Simulatable;
3use rand::Rng;
4use rand_distr::StandardNormal;
5
6/// Cox-Ingersoll-Ross mean-reverting square-root process.
7///
8/// $$dv_t = kappa * (theta - v_t) * dt + sigma * sqrt(v_t) * dW_t$$
9#[derive(Clone)]
10pub struct CoxIngersollRoss {
11    pub kappa: f64,
12    pub theta: f64,
13    pub sigma: f64,
14    state: f64,
15}
16
17impl CoxIngersollRoss {
18    pub fn new(kappa: f64, theta: f64, sigma: f64, initial_value: f64) -> Self {
19        Self {
20            kappa,
21            theta,
22            sigma,
23            state: initial_value,
24        }
25    }
26
27    pub fn value(&self) -> f64 {
28        self.state
29    }
30
31    /// Closed-form expected value.
32    pub fn expected_value(&self, t: f64) -> f64 {
33        let exp = (-self.kappa * t).exp();
34        self.state * exp + self.theta * (1.0 - exp)
35    }
36
37    /// Closed-form variance.
38    pub fn variance(&self, t: f64) -> f64 {
39        let e1 = (-self.kappa * t).exp();
40        let e2 = (-2.0 * self.kappa * t).exp();
41        self.state * (self.sigma.powi(2) / self.kappa) * (e1 - e2)
42            + (self.theta * self.sigma.powi(2) / (2.0 * self.kappa)) * (1.0 - e1).powi(2)
43    }
44}
45
46impl Simulatable for CoxIngersollRoss {
47    type Output = f64;
48
49    fn dim(&self) -> usize {
50        1
51    }
52
53    fn step(&mut self, dt: f64, dw: &[f64], _rng: &mut impl Rng) -> f64 {
54        let v = self.state.max(0.0);
55        let drift = self.kappa * (self.theta - v);
56        let diffusion = if v <= 0.0 { 0.0 } else { self.sigma * v.sqrt() };
57        self.state = v + drift * dt + diffusion * dw[0];
58        if self.state < 0.0 {
59            self.state = 0.0;
60        }
61        self.state
62    }
63
64    fn current(&self) -> f64 {
65        self.state
66    }
67}
68
69impl StochasticProcess for CoxIngersollRoss {
70    type State = f64;
71    type Diffusion = f64;
72    type BrownianIncrement = f64;
73
74    fn drift(&self, x: &f64, _t: f64) -> f64 {
75        self.kappa * (self.theta - x)
76    }
77
78    fn diffusion(&self, x: &f64, _t: f64) -> f64 {
79        if *x < 0.0 { 0.0 } else { self.sigma * x.sqrt() }
80    }
81
82    fn jump_size<R: Rng + ?Sized>(&self, _x: &f64, _t: f64, _rng: &mut R) -> f64 {
83        0.0
84    }
85
86    fn expected_value(&self, x0: &f64, t: f64) -> f64 {
87        let exp_minus_kappa_t = (-self.kappa * t).exp();
88        x0 * exp_minus_kappa_t + self.theta * (1.0 - exp_minus_kappa_t)
89    }
90
91    fn variance(&self, x0: &f64, t: f64) -> f64 {
92        let exp_minus_kappa_t = (-self.kappa * t).exp();
93        let exp_minus_2_kappa_t = (-2.0 * self.kappa * t).exp();
94        let term1 =
95            x0 * (self.sigma.powi(2) / self.kappa) * (exp_minus_kappa_t - exp_minus_2_kappa_t);
96        let term2 = (self.theta * self.sigma.powi(2) / (2.0 * self.kappa))
97            * (1.0 - exp_minus_kappa_t).powi(2);
98        term1 + term2
99    }
100
101    fn step<R: Rng + ?Sized>(&self, x: &f64, t: f64, dt: f64, dw: &f64, _rng: &mut R) -> f64 {
102        let drift = self.drift(x, t);
103        let diffusion = self.diffusion(x, t);
104        let mut next_x = x + drift * dt + diffusion * dw;
105        if next_x < 0.0 {
106            next_x = 0.0;
107        }
108        next_x
109    }
110
111    fn generate_increment<R: Rng + ?Sized>(&self, rng: &mut R, dt: f64) -> f64 {
112        let dw: f64 = rng.sample(StandardNormal);
113        dw * dt.sqrt()
114    }
115
116    fn simulate<R: Rng + ?Sized>(
117        &self,
118        x0: f64,
119        t0: f64,
120        t_end: f64,
121        dt: f64,
122        rng: &mut R,
123    ) -> Vec<(f64, f64)> {
124        let n_steps = ((t_end - t0) / dt).ceil() as usize;
125        let mut path = Vec::with_capacity(n_steps + 1);
126        let mut t = t0;
127        let mut x = x0;
128        path.push((t, x));
129
130        let sqrt_dt = dt.sqrt();
131        for _ in 0..n_steps {
132            let z: f64 = rng.sample(StandardNormal);
133            let dw = z * sqrt_dt;
134            x = self.step(&x, t, dt, &dw, rng);
135            t += dt;
136            path.push((t, x));
137        }
138        path
139    }
140}
141
142#[cfg(test)]
143mod tests {
144    use super::*;
145
146    #[test]
147    fn test_cir_expected_value() {
148        let cir = CoxIngersollRoss::new(3.0, 0.04, 0.3, 0.04);
149        let ev = cir.expected_value(1.0);
150        assert!((ev - 0.04).abs() < 0.01);
151    }
152
153    #[test]
154    fn test_cir_step_nonnegative() {
155        let mut cir = CoxIngersollRoss::new(3.0, 0.04, 0.3, 0.01);
156        for _ in 0..1000 {
157            Simulatable::step(&mut cir, 0.01, &[0.0], &mut rand::rng());
158            assert!(cir.current() >= 0.0);
159        }
160    }
161}
market_model/process/cir.rs

market_model/process/
cir.rs
