pathwise_core/
simulate.rs

1use ndarray::{Array2, Array3};
2use rand::SeedableRng;
3use rayon::prelude::*;
4
5use crate::error::PathwiseError;
6use crate::process::markov::Drift;
7use crate::rng::splitmix64;
8use crate::scheme::Scheme;
9use crate::state::{Diffusion, NoiseIncrement};
10
11/// Simulate `n_paths` paths of a scalar SDE from `t0` to `t1` with `n_steps` steps.
12///
13/// Returns Array2<f64> of shape `(n_paths, n_steps + 1)`.
14///
15/// Generates both `dw` and `dz` per step; schemes that do not use `dz` ignore it.
16/// Non-finite values are stored as NaN.
17#[allow(clippy::too_many_arguments)]
18pub fn simulate<D, G, SC>(
19    drift: &D,
20    diffusion: &G,
21    scheme: &SC,
22    x0: f64,
23    t0: f64,
24    t1: f64,
25    n_paths: usize,
26    n_steps: usize,
27    seed: u64,
28) -> Result<Array2<f64>, PathwiseError>
29where
30    D: Drift<f64> + Sync,
31    G: Diffusion<f64, f64> + Sync,
32    SC: Scheme<f64, Noise = f64>,
33{
34    if n_paths == 0 || n_steps == 0 {
35        return Err(PathwiseError::InvalidParameters(
36            "n_paths and n_steps must be > 0".into(),
37        ));
38    }
39    if t1 <= t0 {
40        return Err(PathwiseError::InvalidParameters("t1 must be > t0".into()));
41    }
42
43    let dt = (t1 - t0) / n_steps as f64;
44    let base_seed = splitmix64(seed);
45
46    let rows: Vec<Vec<f64>> = (0..n_paths)
47        .into_par_iter()
48        .map(|i| {
49            let path_seed = splitmix64(base_seed.wrapping_add(i as u64));
50            let mut rng = rand::rngs::SmallRng::seed_from_u64(path_seed);
51            let mut path = Vec::with_capacity(n_steps + 1);
52            let mut x = x0;
53            path.push(x);
54            for step in 0..n_steps {
55                let t = t0 + step as f64 * dt;
56                let inc = <f64 as NoiseIncrement>::sample(&mut rng, dt);
57                x = scheme.step(drift, diffusion, &x, t, dt, &inc);
58                if !x.is_finite() {
59                    x = f64::NAN;
60                }
61                path.push(x);
62            }
63            path
64        })
65        .collect();
66
67    let mut out = Array2::zeros((n_paths, n_steps + 1));
68    for (i, row) in rows.iter().enumerate() {
69        for (j, &v) in row.iter().enumerate() {
70            out[[i, j]] = v;
71        }
72    }
73    Ok(out)
74}
75
76/// Simulate `n_paths` paths of an N-dimensional SDE from `t0` to `t1`.
77///
78/// Returns Array3<f64> of shape `(n_paths, n_steps + 1, N)`.
79/// Negative values in each component are NOT clipped unless the process does so
80/// in its diffusion (e.g. CIR, Heston full-truncation).
81#[allow(clippy::too_many_arguments)]
82pub fn simulate_nd<const N: usize, D, G, SC>(
83    drift: &D,
84    diffusion: &G,
85    scheme: &SC,
86    x0: nalgebra::SVector<f64, N>,
87    t0: f64,
88    t1: f64,
89    n_paths: usize,
90    n_steps: usize,
91    seed: u64,
92) -> Result<Array3<f64>, PathwiseError>
93where
94    D: Fn(&nalgebra::SVector<f64, N>, f64) -> nalgebra::SVector<f64, N> + Send + Sync,
95    G: Diffusion<nalgebra::SVector<f64, N>, nalgebra::SVector<f64, N>> + Sync,
96    SC: Scheme<nalgebra::SVector<f64, N>, Noise = nalgebra::SVector<f64, N>>,
97{
98    if n_paths == 0 || n_steps == 0 {
99        return Err(PathwiseError::InvalidParameters(
100            "n_paths and n_steps must be > 0".into(),
101        ));
102    }
103    if t1 <= t0 {
104        return Err(PathwiseError::InvalidParameters("t1 must be > t0".into()));
105    }
106
107    let dt = (t1 - t0) / n_steps as f64;
108    let base_seed = splitmix64(seed);
109
110    let rows: Vec<Vec<nalgebra::SVector<f64, N>>> = (0..n_paths)
111        .into_par_iter()
112        .map(|i| {
113            let path_seed = splitmix64(base_seed.wrapping_add(i as u64));
114            let mut rng = rand::rngs::SmallRng::seed_from_u64(path_seed);
115            let mut path = Vec::with_capacity(n_steps + 1);
116            let mut x = x0;
117            path.push(x);
118            for step in 0..n_steps {
119                let t = t0 + step as f64 * dt;
120                let inc = <nalgebra::SVector<f64, N> as NoiseIncrement>::sample(&mut rng, dt);
121                x = scheme.step(drift, diffusion, &x, t, dt, &inc);
122                // Check for non-finite in any component; freeze path
123                if x.iter().any(|v| !v.is_finite()) {
124                    x = nalgebra::SVector::from_fn(|_, _| f64::NAN);
125                }
126                path.push(x);
127            }
128            path
129        })
130        .collect();
131
132    let mut out = Array3::zeros((n_paths, n_steps + 1, N));
133    for (i, path) in rows.iter().enumerate() {
134        for (j, state) in path.iter().enumerate() {
135            for k in 0..N {
136                out[[i, j, k]] = state[k];
137            }
138        }
139    }
140    Ok(out)
141}
142
143#[cfg(test)]
144mod tests {
145    use super::*;
146    use crate::process::markov::{bm, ou};
147    use crate::scheme::euler;
148
149    #[test]
150    fn simulate_returns_correct_shape() {
151        let b = bm();
152        let out = simulate(&b.drift, &b.diffusion, &euler(), 0.0, 0.0, 1.0, 10, 100, 0).unwrap();
153        assert_eq!(out.shape(), &[10, 101]);
154    }
155
156    #[test]
157    fn simulate_first_column_is_x0() {
158        let b = bm();
159        let out = simulate(&b.drift, &b.diffusion, &euler(), 2.5, 0.0, 1.0, 50, 200, 0).unwrap();
160        for i in 0..50 {
161            assert!((out[[i, 0]] - 2.5).abs() < 1e-12, "path {} x0 wrong", i);
162        }
163    }
164
165    #[test]
166    fn simulate_errors_on_zero_paths() {
167        let b = bm();
168        let result = simulate(&b.drift, &b.diffusion, &euler(), 0.0, 0.0, 1.0, 0, 100, 0);
169        assert!(result.is_err());
170    }
171
172    #[test]
173    fn simulate_ou_mean_reverts_on_average() {
174        let o = ou(5.0, 3.0, 0.1);
175        let out = simulate(
176            &o.drift,
177            &o.diffusion,
178            &euler(),
179            0.0,
180            0.0,
181            1.0,
182            2000,
183            500,
184            0,
185        )
186        .unwrap();
187        let last_col = out.column(500);
188        let mean: f64 = last_col.iter().sum::<f64>() / 2000.0;
189        assert!((mean - 3.0).abs() < 0.1, "OU mean={} expected ~3.0", mean);
190    }
191
192    #[test]
193    fn simulate_is_reproducible_with_same_seed() {
194        let b = bm();
195        let r1 = simulate(&b.drift, &b.diffusion, &euler(), 0.0, 0.0, 1.0, 10, 50, 42).unwrap();
196        let r2 = simulate(&b.drift, &b.diffusion, &euler(), 0.0, 0.0, 1.0, 10, 50, 42).unwrap();
197        assert_eq!(r1, r2, "same seed must produce identical output");
198    }
199
200    #[test]
201    fn simulate_different_seeds_produce_different_paths() {
202        let b = bm();
203        let r1 = simulate(&b.drift, &b.diffusion, &euler(), 0.0, 0.0, 1.0, 10, 50, 0).unwrap();
204        let r2 = simulate(&b.drift, &b.diffusion, &euler(), 0.0, 0.0, 1.0, 10, 50, 1).unwrap();
205        assert_ne!(r1, r2, "different seeds must produce different paths");
206    }
207}
pathwise_core/simulate.rs

pathwise_core/
simulate.rs
