pathwise_core/scheme/
milstein.rs

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

pathwise_core/scheme/
milstein.rs
