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
106
107
108
109
110
111
112

//! # Cir 2f
//!
//! $$
//! r_t=x_t+y_t,\ dx_t=\kappa_1(\theta_1-x_t)dt+\sigma_1\sqrt{x_t}dW_t^1,\ dy_t=\kappa_2(\theta_2-y_t)dt+\sigma_2\sqrt{y_t}dW_t^2
//! $$
//!
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 super::cir::Cir;
use crate::traits::FloatExt;
use crate::traits::Fn1D;
use crate::traits::ProcessExt;

pub struct Cir2F<T: FloatExt, S: SeedExt = Unseeded> {
  /// Model parameter controlling process dynamics.
  pub x: Cir<T, S>,
  /// Model parameter controlling process dynamics.
  pub y: Cir<T, S>,
  /// Autoregressive coefficient vector.
  pub phi: Fn1D<T>,
  /// Seed strategy (compile-time: [`Unseeded`] or [`Deterministic`]).
  pub seed: S,
}

impl<T: FloatExt> Cir2F<T> {
  pub fn new(x: Cir<T>, y: Cir<T>, phi: impl Into<Fn1D<T>>) -> Self {
    assert_eq!(x.n, y.n, "x and y Cir factors must use the same n");
    if let (Some(tx), Some(ty)) = (x.t, y.t) {
      assert!(
        (tx - ty).abs() <= T::from_f64_fast(1e-12),
        "x and y Cir factors must use the same time horizon"
      );
    }
    Self {
      x,
      y,
      phi: phi.into(),
      seed: Unseeded,
    }
  }
}

impl<T: FloatExt> Cir2F<T, Deterministic> {
  pub fn seeded(
    x: Cir<T, Deterministic>,
    y: Cir<T, Deterministic>,
    phi: impl Into<Fn1D<T>>,
    seed: u64,
  ) -> Self {
    assert_eq!(x.n, y.n, "x and y Cir factors must use the same n");
    if let (Some(tx), Some(ty)) = (x.t, y.t) {
      assert!(
        (tx - ty).abs() <= T::from_f64_fast(1e-12),
        "x and y Cir factors must use the same time horizon"
      );
    }
    Self {
      x,
      y,
      phi: phi.into(),
      seed: Deterministic::new(seed),
    }
  }
}

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

  fn sample(&self) -> Self::Output {
    let x = self.x.sample();
    let y = self.y.sample();

    let n = x.len();

    let dt = self.x.t.unwrap_or(T::one()) / T::from_usize_(n - 1);
    let phi = Array1::<T>::from_shape_fn(n, |i| self.phi.call(T::from_usize_(i) * dt));

    x + y + phi
  }
}

#[cfg(test)]
mod tests {
  use super::*;

  fn phi_fn(t: f64) -> f64 {
    t
  }

  #[test]
  fn default_time_horizon_is_one() {
    let x = Cir::new(0.0_f64, 0.0, 0.0, 3, Some(0.0), None, Some(false));
    let y = Cir::new(0.0_f64, 0.0, 0.0, 3, Some(0.0), None, Some(false));
    let model = Cir2F::new(x, y, phi_fn as fn(f64) -> f64);

    let out = model.sample();
    assert!((out[0] - 0.0).abs() < 1e-12);
    assert!((out[1] - 0.5).abs() < 1e-12);
    assert!((out[2] - 1.0).abs() < 1e-12);
  }

  #[test]
  #[should_panic(expected = "x and y Cir factors must use the same n")]
  fn mismatched_lengths_panic() {
    let x = Cir::new(0.0_f64, 0.0, 0.0, 3, Some(0.0), Some(1.0), Some(false));
    let y = Cir::new(0.0_f64, 0.0, 0.0, 4, Some(0.0), Some(1.0), Some(false));
    let _ = Cir2F::new(x, y, phi_fn as fn(f64) -> f64);
  }
}