market_model/process/
gbm.rs

1use crate::process::StochasticProcess;
2use crate::simulatable::Simulatable;
3use rand::Rng;
4use rand_distr::StandardNormal;
5
6/// Geometric Brownian Motion.
7///
8/// $$dS_t = mu * S_t * dt + sigma * S_t * dW_t$$
9#[derive(Clone)]
10pub struct GeometricBrownianMotion {
11    pub mu: f64,
12    pub sigma: f64,
13    pub(crate) state: f64,
14    /// Precomputed mu - sigma^2 / 2 (used in analytical step).
15    drift_log: f64,
16}
17
18impl GeometricBrownianMotion {
19    pub fn new(mu: f64, sigma: f64, initial_price: f64) -> Self {
20        let drift_log = mu - 0.5 * sigma * sigma;
21        Self {
22            mu,
23            sigma,
24            state: initial_price,
25            drift_log,
26        }
27    }
28
29    /// Current price.
30    pub fn price(&self) -> f64 {
31        self.state
32    }
33
34    /// Closed-form expected value at horizon `t`.
35    pub fn expected_value(&self, t: f64) -> f64 {
36        self.state * (self.mu * t).exp()
37    }
38
39    /// Closed-form variance at horizon `t`.
40    pub fn variance(&self, t: f64) -> f64 {
41        self.state.powi(2) * (2.0 * self.mu * t).exp() * ((self.sigma.powi(2) * t).exp() - 1.0)
42    }
43}
44
45impl Simulatable for GeometricBrownianMotion {
46    type Output = f64;
47
48    fn dim(&self) -> usize {
49        1
50    }
51
52    fn step(&mut self, dt: f64, dw: &[f64], _rng: &mut impl Rng) -> f64 {
53        // Analytical solution: S_{t+dt} = S_t * exp(drift_log * dt + sigma * dW)
54        let drift_log = self.drift_log * dt;
55        let diff_log = self.sigma * dw[0];
56        self.state *= (drift_log + diff_log).exp();
57        self.state
58    }
59
60    fn current(&self) -> f64 {
61        self.state
62    }
63}
64
65impl StochasticProcess for GeometricBrownianMotion {
66    type State = f64;
67    type Diffusion = f64;
68    type BrownianIncrement = f64;
69
70    fn drift(&self, x: &f64, _t: f64) -> f64 {
71        self.mu * x
72    }
73
74    fn diffusion(&self, x: &f64, _t: f64) -> f64 {
75        self.sigma * x
76    }
77
78    fn jump_size<R: Rng + ?Sized>(&self, _x: &f64, _t: f64, _rng: &mut R) -> f64 {
79        0.0
80    }
81
82    fn expected_value(&self, x0: &f64, t: f64) -> f64 {
83        x0 * (self.mu * t).exp()
84    }
85
86    fn variance(&self, x0: &f64, t: f64) -> f64 {
87        x0.powi(2) * (2.0 * self.mu * t).exp() * ((self.sigma.powi(2) * t).exp() - 1.0)
88    }
89
90    fn step<R: Rng + ?Sized>(&self, x: &f64, t: f64, dt: f64, dw: &f64, _rng: &mut R) -> f64 {
91        let drift = self.drift(x, t);
92        let diffusion = self.diffusion(x, t);
93        x + drift * dt + diffusion * dw
94    }
95
96    fn generate_increment<R: Rng + ?Sized>(&self, rng: &mut R, dt: f64) -> f64 {
97        let dw: f64 = rng.sample(StandardNormal);
98        dw * dt.sqrt()
99    }
100
101    fn simulate<R: Rng + ?Sized>(
102        &self,
103        x0: f64,
104        t0: f64,
105        t_end: f64,
106        dt: f64,
107        rng: &mut R,
108    ) -> Vec<(f64, f64)> {
109        let n_steps = ((t_end - t0) / dt).ceil() as usize;
110        let mut path = Vec::with_capacity(n_steps + 1);
111        let mut t = t0;
112        let mut x = x0;
113        path.push((t, x));
114
115        let sqrt_dt = dt.sqrt();
116        for _ in 0..n_steps {
117            let z: f64 = rng.sample(StandardNormal);
118            let dw = z * sqrt_dt;
119            x = self.step(&x, t, dt, &dw, rng);
120            t += dt;
121            path.push((t, x));
122        }
123        path
124    }
125}
126
127#[cfg(test)]
128mod tests {
129    use super::*;
130
131    #[test]
132    fn test_gbm_expected_value() {
133        let gbm = GeometricBrownianMotion::new(0.1, 0.2, 100.0);
134        let expected = gbm.expected_value(1.0);
135        let computed = 100.0 * 0.1_f64.exp();
136        assert!((expected - computed).abs() < 1e-10);
137    }
138
139    #[test]
140    fn test_gbm_dim() {
141        let gbm = GeometricBrownianMotion::new(0.0, 0.2, 100.0);
142        assert_eq!(gbm.dim(), 1);
143    }
144
145    #[test]
146    fn test_gbm_step_updates_state() {
147        let mut gbm = GeometricBrownianMotion::new(0.0, 0.0, 100.0);
148        let initial = gbm.current();
149        assert_eq!(initial, 100.0);
150        let s1 = Simulatable::step(&mut gbm, 1.0 / 252.0, &[0.0], &mut rand::rng());
151        // With zero drift and zero vol, price stays the same.
152        assert!((s1 - 100.0_f64).abs() < 1e-10);
153        assert_eq!(gbm.current(), s1);
154    }
155
156    #[test]
157    fn test_gbm_monte_carlo_moments() {
158        use rand::SeedableRng;
159        use rand::rngs::StdRng;
160        use rand_distr::StandardNormal;
161        let base = GeometricBrownianMotion::new(0.05, 0.2, 100.0);
162        let dt: f64 = 1.0 / 252.0;
163        let n_steps = 252;
164        let n_paths = 10_000;
165
166        let mut terminals = Vec::with_capacity(n_paths);
167        let mut rng = StdRng::seed_from_u64(42);
168
169        for _ in 0..n_paths {
170            let mut model = base.clone();
171            for _ in 0..n_steps {
172                let z: f64 = rng.sample::<f64, _>(StandardNormal);
173                let dw = z * dt.sqrt();
174                Simulatable::step(&mut model, dt, &[dw], &mut rng);
175            }
176            terminals.push(model.current());
177        }
178
179        let mean: f64 = terminals.iter().sum::<f64>() / n_paths as f64;
180        let var: f64 =
181            terminals.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / (n_paths - 1) as f64;
182
183        let expected_mean = 100.0 * (0.05_f64).exp();
184        let expected_var = base.variance(1.0);
185
186        assert!((mean - expected_mean).abs() < 1.5);
187        assert!((var - expected_var).abs() < 50.0);
188    }
189}
market_model/process/gbm.rs

market_model/process/
gbm.rs
