Skip to main content

pathwise_core/scheme/
euler.rs

1use super::Scheme;
2use crate::state::{Diffusion, Increment, NoiseIncrement, State};
3use crate::process::markov::Drift;
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 { EulerMaruyama }
24
25#[cfg(test)]
26mod tests {
27    use super::*;
28    use crate::state::Increment;
29    use crate::process::markov::{bm, ou};
30
31    #[test]
32    fn euler_step_on_ou_is_deterministic_with_zero_noise() {
33        let scheme = euler();
34        let o = ou(1.0, 0.0, 0.5);
35        // x=1.0, dt=0.01, dw=0 => x_new = 1.0 + (-1.0)*0.01 + 0.5*0 = 0.99
36        let inc = Increment { dw: 0.0_f64, dz: 0.0_f64 };
37        let x_new = scheme.step(&o.drift, &o.diffusion, &1.0_f64, 0.0, 0.01, &inc);
38        assert!((x_new - 0.99).abs() < 1e-12);
39    }
40
41    #[test]
42    fn euler_step_adds_diffusion_term_with_positive_noise() {
43        let scheme = euler();
44        let o = ou(1.0, 0.0, 0.5);
45        // dw = 1.0 => x_new = 1.0 + (-1.0)*0.01 + 0.5*1.0 = 1.49
46        let inc = Increment { dw: 1.0_f64, dz: 0.0_f64 };
47        let x_new = scheme.step(&o.drift, &o.diffusion, &1.0_f64, 0.0, 0.01, &inc);
48        assert!((x_new - 1.49).abs() < 1e-12);
49    }
50
51    #[test]
52    fn euler_step_bm_is_x_plus_dw() {
53        let b = bm();
54        let e = euler();
55        let inc = Increment { dw: 0.3_f64, dz: 0.0 };
56        let x_new = e.step(&b.drift, &b.diffusion, &1.0_f64, 0.0, 0.01, &inc);
57        // BM drift=0, diffusion=1: x_new = 1.0 + 0.0*0.01 + 1.0*0.3 = 1.3
58        assert!((x_new - 1.3).abs() < 1e-12, "got {}", x_new);
59    }
60
61    #[test]
62    fn euler_step_via_increment() {
63        let e = euler();
64        let b = bm();
65        let inc = Increment { dw: 0.3_f64, dz: 0.0_f64 };
66        let x_new: f64 = e.step(&b.drift, &b.diffusion, &0.0_f64, 0.0, 0.01, &inc);
67        assert!(x_new.is_finite());
68    }
69}