tsorder_core/
sim.rs

1//! Simulation utilities for tests and benchmarks (feature `sim`).
2//!
3//! A tiny xorshift64 + Box–Muller generator keeps the crate dependency-free
4//! and the draws reproducible across platforms: the same seed gives the same
5//! series, bit for bit. [`arma`] is the plain Gaussian ARMA generator the
6//! integration tests use; [`ArmaSim`] adds Student-t innovations, additive
7//! outliers and observation noise for benchmark scenarios.
8
9use crate::{Error, Result};
10
11/// Xorshift64 with Box–Muller normals. Deterministic and platform-independent.
12#[derive(Debug, Clone)]
13pub struct Xorshift(u64);
14
15impl Xorshift {
16    /// A zero seed would make xorshift stick at zero, so it is remapped.
17    pub fn new(seed: u64) -> Self {
18        Xorshift(if seed == 0 { 0x9E3779B97F4A7C15 } else { seed })
19    }
20
21    pub fn next_u64(&mut self) -> u64 {
22        let mut x = self.0;
23        x ^= x << 13;
24        x ^= x >> 7;
25        x ^= x << 17;
26        self.0 = x;
27        x
28    }
29
30    /// Uniform on `[0, 1)` with 53 random bits.
31    pub fn uniform(&mut self) -> f64 {
32        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
33    }
34
35    /// Standard normal via Box–Muller (one draw per call, cosine branch).
36    pub fn normal(&mut self) -> f64 {
37        let u1 = self.uniform().max(1e-300);
38        let u2 = self.uniform();
39        (-2.0 * u1.ln()).sqrt() * (2.0 * core::f64::consts::PI * u2).cos()
40    }
41}
42
43/// Innovation distribution.
44#[derive(Debug, Clone, Copy, PartialEq)]
45pub enum Innovations {
46    /// Standard normal.
47    Gaussian,
48    /// Student t with `df` degrees of freedom (`df >= 1`), sampled as
49    /// `Z / sqrt(χ²_df / df)`. Not rescaled to unit variance (which is
50    /// `df/(df-2)` for `df > 2`); every statistic and criterion in the crate
51    /// is scale-invariant, so only the tail weight matters.
52    StudentT { df: u32 },
53}
54
55/// Additive-outlier contamination: with probability `prob`, add `±magnitude`
56/// (random sign) to the observed value.
57#[derive(Debug, Clone, Copy, PartialEq)]
58pub struct Contamination {
59    pub prob: f64,
60    pub magnitude: f64,
61}
62
63/// A simulation design: ARMA(φ, θ) plus the perturbations the benchmark needs.
64#[derive(Debug, Clone, Copy)]
65pub struct ArmaSim<'a> {
66    pub phi: &'a [f64],
67    pub theta: &'a [f64],
68    pub innovations: Innovations,
69    /// Applied to the observed series, after generation.
70    pub additive_outliers: Option<Contamination>,
71    /// Standard deviation of Gaussian measurement noise added to the observed
72    /// series (0 = none).
73    pub obs_noise_sd: f64,
74    /// Transient discarded at the start.
75    pub burn_in: usize,
76}
77
78impl Default for ArmaSim<'static> {
79    fn default() -> Self {
80        ArmaSim {
81            phi: &[],
82            theta: &[],
83            innovations: Innovations::Gaussian,
84            additive_outliers: None,
85            obs_noise_sd: 0.0,
86            burn_in: 500,
87        }
88    }
89}
90
91impl ArmaSim<'_> {
92    /// Simulate `n` observations. Draw order is fixed (innovations first, then
93    /// one uniform per observation for outliers plus a sign when one fires,
94    /// then one normal per observation if `obs_noise_sd > 0`), so results are
95    /// reproducible for a given seed whatever the configuration.
96    pub fn simulate(&self, n: usize, seed: u64) -> Result<Vec<f64>> {
97        if n == 0 {
98            return Err(Error::Invalid("n must be >= 1"));
99        }
100        if let Innovations::StudentT { df } = self.innovations {
101            if df == 0 {
102                return Err(Error::Invalid("StudentT df must be >= 1"));
103            }
104        }
105        if let Some(c) = self.additive_outliers {
106            if !(0.0..=1.0).contains(&c.prob) || !c.magnitude.is_finite() {
107                return Err(Error::Invalid("contamination prob must be in [0, 1]"));
108            }
109        }
110        if !self.obs_noise_sd.is_finite() || self.obs_noise_sd < 0.0 {
111            return Err(Error::Invalid("obs_noise_sd must be finite and >= 0"));
112        }
113        let mut rng = Xorshift::new(seed);
114        let total = n + self.burn_in;
115        let e: Vec<f64> = (0..total)
116            .map(|_| draw(&mut rng, self.innovations))
117            .collect();
118        let start = self.phi.len().max(self.theta.len());
119        let mut x = vec![0.0; total];
120        for t in start..total {
121            let mut v = e[t];
122            for (j, p) in self.phi.iter().enumerate() {
123                v += p * x[t - j - 1];
124            }
125            for (j, th) in self.theta.iter().enumerate() {
126                v += th * e[t - j - 1];
127            }
128            x[t] = v;
129        }
130        let mut out = x[self.burn_in..].to_vec();
131        if let Some(c) = self.additive_outliers {
132            for v in out.iter_mut() {
133                if rng.uniform() < c.prob {
134                    let sign = if rng.uniform() < 0.5 { -1.0 } else { 1.0 };
135                    *v += sign * c.magnitude;
136                }
137            }
138        }
139        if self.obs_noise_sd > 0.0 {
140            for v in out.iter_mut() {
141                *v += self.obs_noise_sd * rng.normal();
142            }
143        }
144        Ok(out)
145    }
146}
147
148fn draw(rng: &mut Xorshift, innovations: Innovations) -> f64 {
149    match innovations {
150        Innovations::Gaussian => rng.normal(),
151        Innovations::StudentT { df } => {
152            let z = rng.normal();
153            let chi2: f64 = (0..df).map(|_| rng.normal().powi(2)).sum();
154            z / (chi2 / df as f64).sqrt()
155        }
156    }
157}
158
159/// Gaussian ARMA with a 500-observation burn-in: the generator the integration
160/// tests (and the landing-page example) use.
161pub fn arma(phi: &[f64], theta: &[f64], n: usize, seed: u64) -> Vec<f64> {
162    ArmaSim {
163        phi,
164        theta,
165        ..Default::default()
166    }
167    .simulate(n, seed)
168    .expect("default ArmaSim configuration is valid")
169}
170
171#[cfg(test)]
172mod tests {
173    use super::*;
174
175    #[test]
176    fn same_seed_same_series() {
177        let a = arma(&[0.5], &[], 100, 42);
178        let b = arma(&[0.5], &[], 100, 42);
179        assert_eq!(a, b);
180        let c = arma(&[0.5], &[], 100, 43);
181        assert_ne!(a, c);
182    }
183
184    #[test]
185    fn student_t_has_heavier_tails_than_gaussian() {
186        let base = ArmaSim::default();
187        let t3 = ArmaSim {
188            innovations: Innovations::StudentT { df: 3 },
189            ..Default::default()
190        };
191        let g = base.simulate(20000, 7).unwrap();
192        let t = t3.simulate(20000, 7).unwrap();
193        let kurt = |x: &[f64]| {
194            let n = x.len() as f64;
195            let m2 = x.iter().map(|v| v * v).sum::<f64>() / n;
196            let m4 = x.iter().map(|v| v.powi(4)).sum::<f64>() / n;
197            m4 / (m2 * m2)
198        };
199        assert!(kurt(&g) < 3.5, "gaussian kurtosis {}", kurt(&g));
200        assert!(kurt(&t) > 4.5, "t3 kurtosis {}", kurt(&t));
201    }
202
203    #[test]
204    fn invalid_configs_are_errors() {
205        assert!(matches!(
206            ArmaSim {
207                innovations: Innovations::StudentT { df: 0 },
208                ..Default::default()
209            }
210            .simulate(100, 1),
211            Err(Error::Invalid(_))
212        ));
213        assert!(matches!(
214            ArmaSim {
215                additive_outliers: Some(Contamination {
216                    prob: 1.5,
217                    magnitude: 9.0
218                }),
219                ..Default::default()
220            }
221            .simulate(100, 1),
222            Err(Error::Invalid(_))
223        ));
224    }
225
226    #[test]
227    fn outliers_and_noise_change_only_the_observed_series() {
228        let clean = ArmaSim {
229            phi: &[0.7],
230            ..Default::default()
231        };
232        let dirty = ArmaSim {
233            phi: &[0.7],
234            additive_outliers: Some(Contamination {
235                prob: 0.05,
236                magnitude: 10.0,
237            }),
238            obs_noise_sd: 0.5,
239            ..Default::default()
240        };
241        let a = clean.simulate(500, 11).unwrap();
242        let b = dirty.simulate(500, 11).unwrap();
243        // Same innovations: the two series differ only by the perturbations.
244        let n_big = a
245            .iter()
246            .zip(&b)
247            .filter(|(u, v)| (*u - *v).abs() > 5.0)
248            .count();
249        assert!(n_big > 5, "expected some outliers, got {n_big}");
250    }
251}
tsorder_core/sim.rs

tsorder_core/
sim.rs
