pathwise_core/scheme/
sri.rs

1// pathwise-core/src/scheme/sri.rs
2use super::Scheme;
3use crate::process::markov::Drift;
4use crate::state::{Diffusion, Increment};
5
6/// Kloeden-Platen strong order 1.5 Taylor scheme for scalar autonomous SDEs.
7///
8/// For the Ito SDE `dX = f(X) dt + g(X) dW`, the step is:
9///
10///   x_{n+1} = x + f*dt + g*dW
11///             + 0.5*g*g'*(dW^2 - dt)                         [Milstein: L^1 g * I_{(1,1)}]
12///             + g*f' * (dt*dW - dZ)                           [L^1 f * I_{(1,0)} term]
13///             + (f*g' + 0.5*g^2*g'') * dZ                    [L^0 g * I_{(0,1)} term]
14///             + 0.5*f*f' * dt^2                               [L^0 f * dt^2/2 term]
15///             + (g*g'^2 + g^2*g'') * (dW^3 - 3*dt*dW)/6     [(L^1)^2 g * I_{(1,1,1)}]
16///
17/// where:
18///   `dZ = integral_0^dt W(s) ds = I_{(0,1)}`
19///   `I_{(1,0)} = dt*dW - dZ`
20///   `I_{(1,1)} = (dW^2 - dt)/2`
21///   `I_{(1,1,1)} = (dW^3 - 3*dt*dW)/6`
22///
23/// All spatial derivatives are approximated by central finite differences with step `h`.
24///
25/// # Reference
26/// Kloeden & Platen, "Numerical Solution of Stochastic Differential Equations", 1992,
27/// Chapter 5.5, the strong Taylor 1.5 scheme (Theorem 5.5.1 / eq. 5.5.4).
28///
29/// # Notes
30/// - The `I_{(1,1,1)}` triple iterated integral is expressible in terms of `dW` alone:
31///   `I_{(1,1,1)} = (dW^3 - 3*dt*dW)/6`. This is the key 1.5-order correction beyond Milstein.
32/// - For BM (constant diffusion), all correction terms vanish, recovering Euler.
33/// - For GBM (`g = sigma*x`, `g' = sigma`, `g'' = 0`):
34///   - `d111 = g*g'^2 + g^2*g'' = sigma^2*x*sigma = sigma^3*x`
35///   - `I_{(1,1,1)} = (dW^3 - 3*dt*dW)/6`
36///   - Combined with the other 1.5 terms, local error reduces to O(dt^2), giving O(dt^{1.5}) globally.
37/// - Requires scalar (1D) SDE with diagonal noise.
38pub struct Sri {
39    h: f64,
40}
41
42impl Sri {
43    pub fn new(h: f64) -> Self {
44        Self { h }
45    }
46}
47
48impl Scheme<f64> for Sri {
49    type Noise = f64;
50
51    fn step<D, G>(
52        &self,
53        drift: &D,
54        diffusion: &G,
55        x: &f64,
56        t: f64,
57        dt: f64,
58        inc: &Increment<f64>,
59    ) -> f64
60    where
61        D: Drift<f64>,
62        G: Diffusion<f64, f64>,
63    {
64        let dw = inc.dw;
65        let dz = inc.dz; // I_{(0,1)} = integral_0^dt W(s) ds
66        let h = self.h;
67
68        // f(x,t) and its central finite difference
69        let f = drift(x, t);
70        let f_plus = drift(&(x + h), t);
71        let f_minus = drift(&(x - h), t);
72        let df_dx = (f_plus - f_minus) / (2.0 * h);
73
74        // g(x,t) and its finite-difference derivatives (unit noise via blanket impl)
75        let g = diffusion.apply(x, t, &1.0_f64);
76        let g_plus = diffusion.apply(&(x + h), t, &1.0_f64);
77        let g_minus = diffusion.apply(&(x - h), t, &1.0_f64);
78        let dg_dx = (g_plus - g_minus) / (2.0 * h);
79        let d2g_dx2 = (g_plus - 2.0 * g + g_minus) / (h * h);
80
81        // I_{(1,0)} = dt*dW - I_{(0,1)} = dt*dW - dZ
82        let i10 = dt * dw - dz;
83        // I_{(1,1)} = (dW^2 - dt)/2  (Milstein)
84        let i11 = (dw * dw - dt) * 0.5;
85        // I_{(1,1,1)} = (dW^3 - 3*dt*dW)/6  (triple iterated Ito integral, expressed via dW)
86        let i111 = (dw * dw * dw - 3.0 * dt * dw) / 6.0;
87
88        // Milstein correction: L^1 g * I_{(1,1)}
89        let milstein = g * dg_dx * i11;
90
91        // L^1 f * I_{(1,0)}: g * f' * (dt*dW - dZ)
92        let term_l1f = g * df_dx * i10;
93
94        // L^0 g * I_{(0,1)}: (f * g' + 0.5 * g^2 * g'') * dZ
95        let term_l0g = (f * dg_dx + 0.5 * g * g * d2g_dx2) * dz;
96
97        // L^0 f * dt^2/2: f * f' * dt^2/2
98        let term_l0f = 0.5 * f * df_dx * dt * dt;
99
100        // (L^1)^2 g * I_{(1,1,1)}: (g*g'^2 + g^2*g'') * (dW^3 - 3*dt*dW)/6
101        let d111 = g * dg_dx * dg_dx + g * g * d2g_dx2;
102        let term_l11g = d111 * i111;
103
104        x + f * dt + g * dw + milstein + term_l1f + term_l0g + term_l0f + term_l11g
105    }
106}
107
108pub fn sri() -> Sri {
109    Sri::new(1e-4)
110}
111
112#[cfg(test)]
113mod tests {
114    use super::*;
115    use crate::process::markov::bm;
116    use crate::state::Increment;
117
118    #[test]
119    fn sri_equals_euler_for_constant_diffusion() {
120        // For constant g: g'=0, g''=0, so all derivative-based corrections vanish.
121        // For BM: f=0, so L^0f and L^1f terms also vanish. SRI == Euler for BM.
122        let b = bm();
123        let s = sri();
124        let e = crate::scheme::euler::euler();
125        let x = 1.0_f64;
126        let inc = Increment { dw: 0.3, dz: 0.001 };
127        let x_euler = e.step(&b.drift, &b.diffusion, &x, 0.0, 0.01, &inc);
128        let x_sri = s.step(&b.drift, &b.diffusion, &x, 0.0, 0.01, &inc);
129        assert!(
130            (x_euler - x_sri).abs() < 1e-6,
131            "BM: euler={} sri={}",
132            x_euler,
133            x_sri
134        );
135    }
136
137    #[test]
138    fn sri_differs_from_milstein_for_state_dependent_diffusion() {
139        // GBM: f'=mu, g'=sigma, g''=0.
140        // SRI adds L^1 f (drift-noise cross), L^0 g (noise-drift), and I_{(1,1,1)} terms.
141        let gbm = crate::process::markov::gbm(0.05, 0.3);
142        let s = sri();
143        let m = crate::scheme::milstein::milstein();
144        let x = 1.0_f64;
145        let inc = Increment {
146            dw: 0.05,
147            dz: 0.0001,
148        };
149        let dt = 0.01;
150        let x_sri = s.step(&gbm.drift, &gbm.diffusion, &x, 0.0, dt, &inc);
151        let x_mil = m.step(&gbm.drift, &gbm.diffusion, &x, 0.0, dt, &inc);
152        // They differ because SRI includes the I_{(1,0)}, I_{(0,1)}, and I_{(1,1,1)} terms
153        assert!(
154            (x_sri - x_mil).abs() > 1e-10,
155            "SRI and Milstein should differ"
156        );
157    }
158}
pathwise_core/scheme/sri.rs

pathwise_core/scheme/
sri.rs
