pathwise_core/scheme/
milstein.rs

1use super::Scheme;
2use crate::state::{Diffusion, Increment};
3use crate::process::markov::Drift;
4
5pub struct Milstein {
6    h: f64,
7}
8
9impl Milstein {
10    pub fn new(h: f64) -> Self { Self { h } }
11}
12
13/// Scalar Milstein scheme (strong order 1.0 for state-dependent diffusion).
14/// Uses central finite difference to approximate dg/dx.
15/// x_{n+1} = x_n + f*dt + g*dw + 0.5*g*(dg/dx)*(dw^2 - dt)
16impl Scheme<f64> for Milstein {
17    type Noise = f64;
18
19    fn step<D, G>(&self, drift: &D, diffusion: &G, x: &f64, t: f64, dt: f64, inc: &Increment<f64>) -> f64
20    where
21        D: Drift<f64>,
22        G: Diffusion<f64, f64>,
23    {
24        let dw = inc.dw;
25        let f = drift(x, t);
26        // g(x,t) via apply with unit noise
27        let g = diffusion.apply(x, t, &1.0_f64);
28        let g_plus  = diffusion.apply(&(x + self.h), t, &1.0_f64);
29        let g_minus = diffusion.apply(&(x - self.h), t, &1.0_f64);
30        let dg_dx = (g_plus - g_minus) / (2.0 * self.h);
31        x + f * dt + g * dw + 0.5 * g * dg_dx * (dw * dw - dt)
32    }
33}
34
35pub fn milstein() -> Milstein { Milstein::new(1e-5) }
36
37#[cfg(test)]
38mod tests {
39    use super::*;
40    use crate::process::markov::bm;
41    use crate::state::Increment;
42
43    #[test]
44    fn milstein_equals_euler_for_constant_diffusion() {
45        let b = bm();
46        let m = milstein();
47        let e = crate::scheme::euler::euler();
48        let x = 1.0_f64;
49        let inc = Increment { dw: 0.3_f64, dz: 0.0 };
50        let x_euler = e.step(&b.drift, &b.diffusion, &x, 0.0, 0.01, &inc);
51        let x_milstein = m.step(&b.drift, &b.diffusion, &x, 0.0, 0.01, &inc);
52        assert!((x_euler - x_milstein).abs() < 1e-8,
53            "euler={} milstein={}", x_euler, x_milstein);
54    }
55
56    #[test]
57    fn milstein_differs_from_euler_for_state_dependent_diffusion() {
58        let gbm = crate::process::markov::gbm(0.0, 0.2);
59        let m = milstein();
60        let e = crate::scheme::euler();
61        let x = 1.0_f64;
62        let inc = Increment { dw: 0.5_f64, dz: 0.0 };
63        let dt = 0.01;
64        let x_euler = e.step(&gbm.drift, &gbm.diffusion, &x, 0.0, dt, &inc);
65        let x_milstein = m.step(&gbm.drift, &gbm.diffusion, &x, 0.0, dt, &inc);
66        assert!((x_euler - x_milstein).abs() > 1e-4, "should differ by correction term");
67    }
68}
pathwise_core/scheme/milstein.rs

pathwise_core/scheme/
milstein.rs
