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

pathwise_core/
simulate.rs
