pathwise_core/scheme/
sri.rs

1// pathwise-core/src/scheme/sri.rs
2use super::Scheme;
3use crate::state::{Diffusion, Increment};
4use crate::process::markov::Drift;
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 { Self { h } }
44}
45
46impl Scheme<f64> for Sri {
47    type Noise = f64;
48
49    fn step<D, G>(
50        &self,
51        drift: &D,
52        diffusion: &G,
53        x: &f64,
54        t: f64,
55        dt: f64,
56        inc: &Increment<f64>,
57    ) -> f64
58    where
59        D: Drift<f64>,
60        G: Diffusion<f64, f64>,
61    {
62        let dw = inc.dw;
63        let dz = inc.dz;  // I_{(0,1)} = integral_0^dt W(s) ds
64        let h = self.h;
65
66        // f(x,t) and its central finite difference
67        let f = drift(x, t);
68        let f_plus  = drift(&(x + h), t);
69        let f_minus = drift(&(x - h), t);
70        let df_dx   = (f_plus - f_minus) / (2.0 * h);
71
72        // g(x,t) and its finite-difference derivatives (unit noise via blanket impl)
73        let g       = diffusion.apply(x, t, &1.0_f64);
74        let g_plus  = diffusion.apply(&(x + h), t, &1.0_f64);
75        let g_minus = diffusion.apply(&(x - h), t, &1.0_f64);
76        let dg_dx   = (g_plus - g_minus) / (2.0 * h);
77        let d2g_dx2 = (g_plus - 2.0 * g + g_minus) / (h * h);
78
79        // I_{(1,0)} = dt*dW - I_{(0,1)} = dt*dW - dZ
80        let i10 = dt * dw - dz;
81        // I_{(1,1)} = (dW^2 - dt)/2  (Milstein)
82        let i11 = (dw * dw - dt) * 0.5;
83        // I_{(1,1,1)} = (dW^3 - 3*dt*dW)/6  (triple iterated Ito integral, expressed via dW)
84        let i111 = (dw * dw * dw - 3.0 * dt * dw) / 6.0;
85
86        // Milstein correction: L^1 g * I_{(1,1)}
87        let milstein = g * dg_dx * i11;
88
89        // L^1 f * I_{(1,0)}: g * f' * (dt*dW - dZ)
90        let term_l1f = g * df_dx * i10;
91
92        // L^0 g * I_{(0,1)}: (f * g' + 0.5 * g^2 * g'') * dZ
93        let term_l0g = (f * dg_dx + 0.5 * g * g * d2g_dx2) * dz;
94
95        // L^0 f * dt^2/2: f * f' * dt^2/2
96        let term_l0f = 0.5 * f * df_dx * dt * dt;
97
98        // (L^1)^2 g * I_{(1,1,1)}: (g*g'^2 + g^2*g'') * (dW^3 - 3*dt*dW)/6
99        let d111 = g * dg_dx * dg_dx + g * g * d2g_dx2;
100        let term_l11g = d111 * i111;
101
102        x + f * dt + g * dw + milstein + term_l1f + term_l0g + term_l0f + term_l11g
103    }
104}
105
106pub fn sri() -> Sri { Sri::new(1e-4) }
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111    use crate::process::markov::bm;
112    use crate::state::Increment;
113
114    #[test]
115    fn sri_equals_euler_for_constant_diffusion() {
116        // For constant g: g'=0, g''=0, so all derivative-based corrections vanish.
117        // For BM: f=0, so L^0f and L^1f terms also vanish. SRI == Euler for BM.
118        let b = bm();
119        let s = sri();
120        let e = crate::scheme::euler::euler();
121        let x = 1.0_f64;
122        let inc = Increment { dw: 0.3, dz: 0.001 };
123        let x_euler = e.step(&b.drift, &b.diffusion, &x, 0.0, 0.01, &inc);
124        let x_sri = s.step(&b.drift, &b.diffusion, &x, 0.0, 0.01, &inc);
125        assert!((x_euler - x_sri).abs() < 1e-6,
126            "BM: euler={} sri={}", x_euler, x_sri);
127    }
128
129    #[test]
130    fn sri_differs_from_milstein_for_state_dependent_diffusion() {
131        // GBM: f'=mu, g'=sigma, g''=0.
132        // SRI adds L^1 f (drift-noise cross), L^0 g (noise-drift), and I_{(1,1,1)} terms.
133        let gbm = crate::process::markov::gbm(0.05, 0.3);
134        let s = sri();
135        let m = crate::scheme::milstein::milstein();
136        let x = 1.0_f64;
137        let inc = Increment { dw: 0.05, dz: 0.0001 };
138        let dt = 0.01;
139        let x_sri = s.step(&gbm.drift, &gbm.diffusion, &x, 0.0, dt, &inc);
140        let x_mil = m.step(&gbm.drift, &gbm.diffusion, &x, 0.0, dt, &inc);
141        // They differ because SRI includes the I_{(1,0)}, I_{(0,1)}, and I_{(1,1,1)} terms
142        assert!((x_sri - x_mil).abs() > 1e-10, "SRI and Milstein should differ");
143    }
144}
pathwise_core/scheme/sri.rs

pathwise_core/scheme/
sri.rs
