stochastic-rs-stochastic 2.0.0

Stochastic process simulation.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

//! # Arch
//!
//! $$
//! \sigma_t^2=\omega+\sum_{i=1}^m\alpha_iX_{t-i}^2,\qquad X_t=\sigma_t z_t,\ z_t\sim\mathcal N(0,1)
//! $$
//!
use ndarray::Array1;
use stochastic_rs_core::simd_rng::Deterministic;
use stochastic_rs_core::simd_rng::SeedExt;
use stochastic_rs_core::simd_rng::Unseeded;
use stochastic_rs_distributions::normal::SimdNormal;

use crate::traits::FloatExt;
use crate::traits::ProcessExt;

/// Implements an Arch(m) model:
///
/// \[
///   \sigma_t^2 = \omega + \sum_{i=1}^m \alpha_i X_{t-i}^2,
///   \quad X_t = \sigma_t \cdot z_t, \quad z_t \sim \mathcal{N}(0,1).
/// \]
///
/// # Fields
/// - `omega`: Constant term.
/// - `alpha`: Array of Arch coefficients.
/// - `n`: Number of observations.
/// - `m`: Optional batch size.
#[derive(Debug, Clone)]
pub struct Arch<T: FloatExt, S: SeedExt = Unseeded> {
  /// Omega (constant term in variance)
  pub omega: T,
  /// Coefficients alpha_i
  pub alpha: Array1<T>,
  /// Length of series
  pub n: usize,
  /// Seed strategy (compile-time: [`Unseeded`] or [`Deterministic`]).
  pub seed: S,
}

impl<T: FloatExt> Arch<T> {
  /// Create a new Arch model.
  pub fn new(omega: T, alpha: Array1<T>, n: usize) -> Self {
    assert!(omega > T::zero(), "Arch requires omega > 0");
    Self {
      omega,
      alpha,
      n,
      seed: Unseeded,
    }
  }
}

impl<T: FloatExt> Arch<T, Deterministic> {
  /// Create a new Arch model with a deterministic seed for reproducible output.
  pub fn seeded(omega: T, alpha: Array1<T>, n: usize, seed: u64) -> Self {
    assert!(omega > T::zero(), "Arch requires omega > 0");
    Self {
      omega,
      alpha,
      n,
      seed: Deterministic::new(seed),
    }
  }
}

impl<T: FloatExt, S: SeedExt> ProcessExt<T> for Arch<T, S> {
  type Output = Array1<T>;

  fn sample(&self) -> Self::Output {
    let m = self.alpha.len();
    let mut z = Array1::<T>::zeros(self.n);
    if self.n > 0 {
      let slice = z.as_slice_mut().expect("contiguous");
      let normal = SimdNormal::<T>::from_seed_source(T::zero(), T::one(), &self.seed);
      normal.fill_slice_fast(slice);
    }
    let mut x = Array1::<T>::zeros(self.n);
    let var_floor = T::from_f64_fast(1e-12);

    for t in 0..self.n {
      // compute sigma_t^2
      let mut var_t = self.omega;
      for i in 1..=m {
        if t >= i {
          let x_lag = x[t - i];
          var_t += self.alpha[i - 1] * x_lag.powi(2);
        }
      }
      assert!(
        var_t.is_finite() && var_t > T::zero(),
        "Arch produced non-positive or non-finite conditional variance at t={}",
        t
      );
      let sigma_t = var_t.max(var_floor).sqrt();
      x[t] = sigma_t * z[t];
    }

    x
  }
}

py_process_1d!(PyArch, Arch,
  sig: (omega, alpha, n, seed=None, dtype=None),
  params: (omega: f64, alpha: Vec<f64>, n: usize)
);