Skip to main content

pathwise_core/scheme/
euler.rs

1use super::Scheme;
2use crate::process::markov::Drift;
3use crate::state::{Diffusion, Increment, NoiseIncrement, State};
4
5pub struct EulerMaruyama;
6
7/// Generic Euler-Maruyama: works for any State S with Noise = S.
8/// x_{n+1} = x_n + f(x_n, t)*dt + g(x_n, t)*dW
9impl<S: State + NoiseIncrement> Scheme<S> for EulerMaruyama {
10    type Noise = S;
11
12    fn step<D, G>(&self, drift: &D, diffusion: &G, x: &S, t: f64, dt: f64, inc: &Increment<S>) -> S
13    where
14        D: Drift<S>,
15        G: Diffusion<S, S>,
16    {
17        let f_dt = drift(x, t) * dt;
18        let g_dw = diffusion.apply(x, t, &inc.dw);
19        x.clone() + f_dt + g_dw
20    }
21}
22
23pub fn euler() -> EulerMaruyama {
24    EulerMaruyama
25}
26
27#[cfg(test)]
28mod tests {
29    use super::*;
30    use crate::process::markov::{bm, ou};
31    use crate::state::Increment;
32
33    #[test]
34    fn euler_step_on_ou_is_deterministic_with_zero_noise() {
35        let scheme = euler();
36        let o = ou(1.0, 0.0, 0.5);
37        // x=1.0, dt=0.01, dw=0 => x_new = 1.0 + (-1.0)*0.01 + 0.5*0 = 0.99
38        let inc = Increment {
39            dw: 0.0_f64,
40            dz: 0.0_f64,
41        };
42        let x_new = scheme.step(&o.drift, &o.diffusion, &1.0_f64, 0.0, 0.01, &inc);
43        assert!((x_new - 0.99).abs() < 1e-12);
44    }
45
46    #[test]
47    fn euler_step_adds_diffusion_term_with_positive_noise() {
48        let scheme = euler();
49        let o = ou(1.0, 0.0, 0.5);
50        // dw = 1.0 => x_new = 1.0 + (-1.0)*0.01 + 0.5*1.0 = 1.49
51        let inc = Increment {
52            dw: 1.0_f64,
53            dz: 0.0_f64,
54        };
55        let x_new = scheme.step(&o.drift, &o.diffusion, &1.0_f64, 0.0, 0.01, &inc);
56        assert!((x_new - 1.49).abs() < 1e-12);
57    }
58
59    #[test]
60    fn euler_step_bm_is_x_plus_dw() {
61        let b = bm();
62        let e = euler();
63        let inc = Increment {
64            dw: 0.3_f64,
65            dz: 0.0,
66        };
67        let x_new = e.step(&b.drift, &b.diffusion, &1.0_f64, 0.0, 0.01, &inc);
68        // BM drift=0, diffusion=1: x_new = 1.0 + 0.0*0.01 + 1.0*0.3 = 1.3
69        assert!((x_new - 1.3).abs() < 1e-12, "got {}", x_new);
70    }
71
72    #[test]
73    fn euler_step_via_increment() {
74        let e = euler();
75        let b = bm();
76        let inc = Increment {
77            dw: 0.3_f64,
78            dz: 0.0_f64,
79        };
80        let x_new: f64 = e.step(&b.drift, &b.diffusion, &0.0_f64, 0.0, 0.01, &inc);
81        assert!(x_new.is_finite());
82    }
83}