market_model/process/
rough_ou.rs

1//! Rough Ornstein-Uhlenbeck process.
2//!
3//! A fractional Ornstein-Uhlenbeck process approximated via a multi-factor
4//! Markovian representation (Abi Jaber, 2019).
5
6use crate::process::StochasticProcess;
7use crate::simulatable::Simulatable;
8use rand::Rng;
9use rand_distr::StandardNormal;
10
11/// State for the multi-factor rough OU process.
12#[derive(Clone, Debug)]
13pub struct RoughOUState {
14    pub value: f64,
15    pub factors: Vec<f64>,
16}
17
18/// Rough Ornstein-Uhlenbeck process with Hurst parameter `H`.
19///
20/// The fractional kernel K(t) = t^{H-1/2} / Gamma(H+1/2) is approximated
21/// by a sum of exponentials via `n_factors` factors.
22#[derive(Clone)]
23pub struct RoughOrnsteinUhlenbeck {
24    pub kappa: f64,
25    pub theta: f64,
26    pub sigma: f64,
27    pub hurst: f64,
28
29    weights: Vec<f64>,
30    mean_reversions: Vec<f64>,
31    state: RoughOUState,
32}
33
34impl RoughOrnsteinUhlenbeck {
35    pub fn new(
36        kappa: f64,
37        theta: f64,
38        sigma: f64,
39        hurst: f64,
40        n_factors: usize,
41        initial_value: f64,
42    ) -> Self {
43        let (weights, mean_reversions) = compute_factors(hurst, n_factors);
44        let factors = vec![0.0; n_factors];
45        Self {
46            kappa,
47            theta,
48            sigma,
49            hurst,
50            weights,
51            mean_reversions,
52            state: RoughOUState {
53                value: initial_value,
54                factors,
55            },
56        }
57    }
58
59    pub fn value(&self) -> f64 {
60        self.state.value
61    }
62
63    /// Closed-form expected value.
64    pub fn expected_value(&self, t: f64) -> f64 {
65        self.theta + (self.state.value - self.theta) * (-self.kappa * t).exp()
66    }
67}
68
69fn compute_factors(hurst: f64, n: usize) -> (Vec<f64>, Vec<f64>) {
70    let alpha = hurst + 0.5;
71    let r_min = 1e-4_f64;
72    let r_max = 100.0_f64;
73
74    if n == 1 {
75        return (vec![1.0], vec![1.0]);
76    }
77
78    let ratio = r_max / r_min;
79    let base = ratio.powf(1.0 / (n as f64 - 1.0));
80
81    let mut x = Vec::with_capacity(n);
82    let mut c = Vec::with_capacity(n);
83
84    for i in 0..n {
85        let xi = r_min * base.powf(i as f64);
86        x.push(xi);
87        c.push(xi.powf(-alpha));
88    }
89
90    // Normalize weights.
91    let mut variance_at_1 = 0.0;
92    for i in 0..n {
93        for j in 0..n {
94            let s = x[i] + x[j];
95            variance_at_1 += c[i] * c[j] * (1.0 - (-s).exp()) / s;
96        }
97    }
98
99    let scale = (1.0 / variance_at_1).sqrt();
100    for val in c.iter_mut() {
101        *val *= scale;
102    }
103
104    (c, x)
105}
106
107impl Simulatable for RoughOrnsteinUhlenbeck {
108    type Output = RoughOUState;
109
110    fn dim(&self) -> usize {
111        1 // Single Brownian driver for all factors
112    }
113
114    fn step(&mut self, dt: f64, dw: &[f64], _rng: &mut impl Rng) -> Self::Output {
115        let sqrt_dt = dt.sqrt();
116        let scaled_dw = dw[0] * sqrt_dt;
117
118        let mut sum_dy = 0.0;
119        let mut new_factors = std::mem::take(&mut self.state.factors);
120
121        for (i, y_i) in new_factors.iter_mut().enumerate() {
122            let xi = self.mean_reversions[i];
123            let ci = self.weights[i];
124            let dy = -xi * *y_i * dt + ci * scaled_dw;
125            *y_i += dy;
126            sum_dy += dy;
127        }
128
129        let dx = self.kappa * (self.theta - self.state.value) * dt + self.sigma * sum_dy;
130
131        self.state = RoughOUState {
132            value: self.state.value + dx,
133            factors: new_factors,
134        };
135
136        self.state.clone()
137    }
138
139    fn current(&self) -> RoughOUState {
140        self.state.clone()
141    }
142}
143
144impl StochasticProcess for RoughOrnsteinUhlenbeck {
145    type State = RoughOUState;
146    type Diffusion = f64;
147    type BrownianIncrement = f64;
148
149    fn drift(&self, x: &Self::State, _t: f64) -> Self::State {
150        let mut drift_factors = Vec::with_capacity(self.mean_reversions.len());
151        let mut sum_factor_drifts = 0.0;
152
153        for (i, &y_i) in x.factors.iter().enumerate() {
154            let d_y = -self.mean_reversions[i] * y_i;
155            drift_factors.push(d_y);
156            sum_factor_drifts += d_y;
157        }
158
159        let drift_val = self.kappa * (self.theta - x.value) + self.sigma * sum_factor_drifts;
160
161        RoughOUState {
162            value: drift_val,
163            factors: drift_factors,
164        }
165    }
166
167    fn diffusion(&self, _x: &Self::State, _t: f64) -> Self::Diffusion {
168        let sum_c: f64 = self.weights.iter().sum();
169        self.sigma * sum_c
170    }
171
172    fn jump_size<R: Rng + ?Sized>(&self, _x: &Self::State, _t: f64, _rng: &mut R) -> Self::State {
173        RoughOUState {
174            value: 0.0,
175            factors: vec![0.0; self.weights.len()],
176        }
177    }
178
179    fn expected_value(&self, x0: &Self::State, t: f64) -> Self::State {
180        RoughOUState {
181            value: self.theta + (x0.value - self.theta) * (-self.kappa * t).exp(),
182            factors: x0.factors.clone(),
183        }
184    }
185
186    fn variance(&self, _x0: &Self::State, t: f64) -> Self::State {
187        RoughOUState {
188            value: self.sigma * self.sigma * t,
189            factors: vec![],
190        }
191    }
192
193    fn step<R: Rng + ?Sized>(
194        &self,
195        x: &Self::State,
196        _t: f64,
197        dt: f64,
198        dw: &Self::BrownianIncrement,
199        _rng: &mut R,
200    ) -> Self::State {
201        let mut new_factors = Vec::with_capacity(x.factors.len());
202        let mut sum_dy = 0.0;
203
204        let scaled_dw = dw * dt.sqrt();
205
206        for (i, &y_i) in x.factors.iter().enumerate() {
207            let xi = self.mean_reversions[i];
208            let ci = self.weights[i];
209            let dy = -xi * y_i * dt + ci * scaled_dw;
210            let new_y = y_i + dy;
211            new_factors.push(new_y);
212            sum_dy += dy;
213        }
214
215        let dx = self.kappa * (self.theta - x.value) * dt + self.sigma * sum_dy;
216        let new_val = x.value + dx;
217
218        RoughOUState {
219            value: new_val,
220            factors: new_factors,
221        }
222    }
223
224    fn generate_increment<R: Rng + ?Sized>(
225        &self,
226        rng: &mut R,
227        _dt: f64,
228    ) -> Self::BrownianIncrement {
229        rng.sample(StandardNormal)
230    }
231
232    fn simulate<R: Rng + ?Sized>(
233        &self,
234        x0: Self::State,
235        t0: f64,
236        t_end: f64,
237        dt: f64,
238        rng: &mut R,
239    ) -> Vec<(f64, Self::State)> {
240        let mut t = t0;
241        let mut x = x0;
242        let mut path = Vec::new();
243        path.push((t, x.clone()));
244
245        while t < t_end {
246            let dw = self.generate_increment(rng, dt);
247            x = self.step(&x, t, dt, &dw, rng);
248            t += dt;
249            path.push((t, x.clone()));
250        }
251        path
252    }
253}
254
255#[cfg(test)]
256mod tests {
257    use super::*;
258
259    #[test]
260    fn test_rough_ou_dim() {
261        let rou = RoughOrnsteinUhlenbeck::new(1.0, 0.05, 0.1, 0.1, 5, 0.03);
262        assert_eq!(rou.dim(), 1);
263    }
264
265    #[test]
266    fn test_rough_ou_initial_state() {
267        let rou = RoughOrnsteinUhlenbeck::new(1.0, 0.05, 0.1, 0.1, 5, 0.03);
268        assert_eq!(rou.value(), 0.03);
269        assert_eq!(rou.state.factors.len(), 5);
270    }
271}
market_model/process/rough_ou.rs

market_model/process/
rough_ou.rs
