Skip to content

Commit 0fa9896

Browse files
dancixxdancixx
committed
feat: add sample_simd for SamplingExt<f32,f64>
1 parent b183cd5 commit 0fa9896

File tree

33 files changed

+1696
-20
lines changed

33 files changed

+1696
-20
lines changed

‎Cargo.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
[package]
22
name = "stochastic-rs"
3-
version = "0.15.0"
3+
version = "0.15.1"
44
edition = "2021"
55
license = "MIT"
66
description = "A Rust library for quant finance and simulating stochastic processes."

‎src/stochastic/autoregressive/agrach.rs‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,54 @@ impl SamplingExt<f64> for AGARCH<f64> {
8787
x
8888
}
8989

90+
fn sample_simd(&self) -> Array1<f64> {
91+
use crate::stats::distr::normal::SimdNormal;
92+
93+
let p = self.alpha.len();
94+
let q = self.beta.len();
95+
96+
// Generate white noise
97+
let z = Array1::random(self.n, SimdNormal::new(0.0, 1.0));
98+
99+
// Arrays for X_t and sigma_t^2
100+
let mut x = Array1::<f64>::zeros(self.n);
101+
let mut sigma2 = Array1::<f64>::zeros(self.n);
102+
103+
// Summation for unconditional variance init
104+
let sum_alpha: f64 = self.alpha.iter().sum();
105+
let sum_delta_half: f64 = self.delta.iter().sum::<f64>() * 0.5;
106+
let sum_beta: f64 = self.beta.iter().sum();
107+
let denom = (1.0 - sum_alpha - sum_delta_half - sum_beta).max(1e-8);
108+
109+
for t in 0..self.n {
110+
if t == 0 {
111+
sigma2[t] = self.omega / denom;
112+
} else {
113+
let mut var_t = self.omega;
114+
// p-lag terms
115+
for i in 1..=p {
116+
if t >= i {
117+
let x_lag = x[t - i];
118+
let indicator = if x_lag < 0.0 { 1.0 } else { 0.0 };
119+
120+
var_t +=
121+
self.alpha[i - 1] * x_lag.powi(2) + self.delta[i - 1] * x_lag.powi(2) * indicator;
122+
}
123+
}
124+
// q-lag terms
125+
for j in 1..=q {
126+
if t >= j {
127+
var_t += self.beta[j - 1] * sigma2[t - j];
128+
}
129+
}
130+
sigma2[t] = var_t;
131+
}
132+
x[t] = sigma2[t].sqrt() * z[t] as f64;
133+
}
134+
135+
x
136+
}
137+
90138
fn n(&self) -> usize {
91139
self.n
92140
}
@@ -144,6 +192,55 @@ impl SamplingExt<f32> for AGARCH<f32> {
144192
x
145193
}
146194

195+
#[cfg(feature = "simd")]
196+
fn sample_simd(&self) -> Array1<f32> {
197+
use crate::stats::distr::normal::SimdNormal;
198+
199+
let p = self.alpha.len();
200+
let q = self.beta.len();
201+
202+
// Generate white noise
203+
let z = Array1::random(self.n, SimdNormal::new(0.0, 1.0));
204+
205+
// Arrays for X_t and sigma_t^2
206+
let mut x = Array1::<f32>::zeros(self.n);
207+
let mut sigma2 = Array1::<f32>::zeros(self.n);
208+
209+
// Summation for unconditional variance init
210+
let sum_alpha: f32 = self.alpha.iter().sum();
211+
let sum_delta_half: f32 = self.delta.iter().sum::<f32>() * 0.5;
212+
let sum_beta: f32 = self.beta.iter().sum();
213+
let denom = (1.0 - sum_alpha - sum_delta_half - sum_beta).max(1e-8);
214+
215+
for t in 0..self.n {
216+
if t == 0 {
217+
sigma2[t] = self.omega / denom;
218+
} else {
219+
let mut var_t = self.omega;
220+
// p-lag terms
221+
for i in 1..=p {
222+
if t >= i {
223+
let x_lag = x[t - i];
224+
let indicator = if x_lag < 0.0 { 1.0 } else { 0.0 };
225+
226+
var_t +=
227+
self.alpha[i - 1] * x_lag.powi(2) + self.delta[i - 1] * x_lag.powi(2) * indicator;
228+
}
229+
}
230+
// q-lag terms
231+
for j in 1..=q {
232+
if t >= j {
233+
var_t += self.beta[j - 1] * sigma2[t - j];
234+
}
235+
}
236+
sigma2[t] = var_t;
237+
}
238+
x[t] = sigma2[t].sqrt() * z[t];
239+
}
240+
241+
x
242+
}
243+
147244
fn n(&self) -> usize {
148245
self.n
149246
}

‎src/stochastic/autoregressive/ar.rs‎

Lines changed: 64 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,38 @@ impl SamplingExt<f64> for ARp<f64> {
6464
series
6565
}
6666

67+
#[cfg(feature = "simd")]
68+
fn sample_simd(&self) -> Array1<f64> {
69+
use crate::stats::distr::normal::SimdNormal;
70+
71+
let p = self.phi.len();
72+
let noise = Array1::random(self.n, SimdNormal::new(0.0, self.sigma as f32));
73+
let mut series = Array1::<f64>::zeros(self.n);
74+
75+
// Fill initial conditions if provided
76+
if let Some(init) = &self.x0 {
77+
// Copy up to min(p, n)
78+
for i in 0..p.min(self.n) {
79+
series[i] = init[i];
80+
}
81+
}
82+
83+
// AR recursion
84+
for t in 0..self.n {
85+
let mut val = 0.0;
86+
// Sum over AR lags
87+
for k in 1..=p {
88+
if t >= k {
89+
val += self.phi[k - 1] * series[t - k];
90+
}
91+
}
92+
// Add noise
93+
series[t] += val + noise[t] as f64;
94+
}
95+
96+
series
97+
}
98+
6799
fn n(&self) -> usize {
68100
self.n
69101
}
@@ -105,6 +137,38 @@ impl SamplingExt<f32> for ARp<f32> {
105137
series
106138
}
107139

140+
#[cfg(feature = "simd")]
141+
fn sample_simd(&self) -> Array1<f32> {
142+
use crate::stats::distr::normal::SimdNormal;
143+
144+
let p = self.phi.len();
145+
let noise = Array1::random(self.n, SimdNormal::new(0.0, self.sigma));
146+
let mut series = Array1::<f32>::zeros(self.n);
147+
148+
// Fill initial conditions if provided
149+
if let Some(init) = &self.x0 {
150+
// Copy up to min(p, n)
151+
for i in 0..p.min(self.n) {
152+
series[i] = init[i];
153+
}
154+
}
155+
156+
// AR recursion
157+
for t in 0..self.n {
158+
let mut val = 0.0;
159+
// Sum over AR lags
160+
for k in 1..=p {
161+
if t >= k {
162+
val += self.phi[k - 1] * series[t - k];
163+
}
164+
}
165+
// Add noise
166+
series[t] += val + noise[t];
167+
}
168+
169+
series
170+
}
171+
108172
fn n(&self) -> usize {
109173
self.n
110174
}
@@ -113,21 +177,3 @@ impl SamplingExt<f32> for ARp<f32> {
113177
self.m
114178
}
115179
}
116-
117-
#[cfg(test)]
118-
mod tests {
119-
use ndarray::arr1;
120-
121-
use crate::{
122-
plot_1d,
123-
stochastic::{autoregressive::ar::ARp, SamplingExt},
124-
};
125-
126-
#[test]
127-
fn ar_plot() {
128-
// Suppose p=2 with user-defined coefficients
129-
let phi = arr1(&[0.5, -0.25]);
130-
let ar_model = ARp::new(phi, 1.0, 100, None, None);
131-
plot_1d!(ar_model.sample(), "AR(p) process");
132-
}
133-
}

‎src/stochastic/autoregressive/arch.rs‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,30 @@ impl SamplingExt<f64> for ARCH<f64> {
5252
x
5353
}
5454

55+
#[cfg(feature = "simd")]
56+
fn sample_simd(&self) -> Array1<f64> {
57+
use crate::stats::distr::normal::SimdNormal;
58+
59+
let m = self.alpha.len();
60+
let z = Array1::random(self.n, SimdNormal::new(0.0, 1.0));
61+
let mut x = Array1::<f64>::zeros(self.n);
62+
63+
for t in 0..self.n {
64+
// compute sigma_t^2
65+
let mut var_t = self.omega;
66+
for i in 1..=m {
67+
if t >= i {
68+
let x_lag = x[t - i];
69+
var_t += self.alpha[i - 1] * x_lag.powi(2);
70+
}
71+
}
72+
let sigma_t = var_t.sqrt();
73+
x[t] = sigma_t * z[t] as f64;
74+
}
75+
76+
x
77+
}
78+
5579
fn n(&self) -> usize {
5680
self.n
5781
}
@@ -84,6 +108,30 @@ impl SamplingExt<f32> for ARCH<f32> {
84108
x
85109
}
86110

111+
#[cfg(feature = "simd")]
112+
fn sample_simd(&self) -> Array1<f32> {
113+
use crate::stats::distr::normal::SimdNormal;
114+
115+
let m = self.alpha.len();
116+
let z = Array1::random(self.n, SimdNormal::new(0.0, 1.0));
117+
let mut x = Array1::<f32>::zeros(self.n);
118+
119+
for t in 0..self.n {
120+
// compute sigma_t^2
121+
let mut var_t = self.omega;
122+
for i in 1..=m {
123+
if t >= i {
124+
let x_lag = x[t - i];
125+
var_t += self.alpha[i - 1] * x_lag.powi(2);
126+
}
127+
}
128+
let sigma_t = var_t.sqrt();
129+
x[t] = sigma_t * z[t];
130+
}
131+
132+
x
133+
}
134+
87135
fn n(&self) -> usize {
88136
self.n
89137
}

‎src/stochastic/autoregressive/arima.rs‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,34 @@ impl SamplingExt<f64> for ARIMA<f64> {
5656
result
5757
}
5858

59+
#[cfg(feature = "simd")]
60+
fn sample_simd(&self) -> Array1<f64> {
61+
// 1) Generate an AR(p) series with user-provided coefficients
62+
let ar_model = ARp::new(
63+
self.ar_coefs.clone(),
64+
self.sigma,
65+
self.n,
66+
None, // batch
67+
None, // x0
68+
);
69+
let ar_series = ar_model.sample_simd();
70+
71+
// 2) Generate an MA(q) series with user-provided coefficients
72+
let ma_model = MAq::new(self.ma_coefs.clone(), self.sigma, self.n, None);
73+
let ma_series = ma_model.sample_simd();
74+
75+
// 3) Summation -> ARMA(p,q)
76+
let arma_series = &ar_series + &ma_series;
77+
78+
// 4) Inverse difference d times -> ARIMA(p,d,q)
79+
let mut result = arma_series;
80+
for _ in 0..self.d {
81+
result = inverse_difference(&result);
82+
}
83+
84+
result
85+
}
86+
5987
fn n(&self) -> usize {
6088
self.n
6189
}
@@ -84,6 +112,24 @@ impl SamplingExt<f32> for ARIMA<f32> {
84112
result
85113
}
86114

115+
#[cfg(feature = "simd")]
116+
fn sample_simd(&self) -> Array1<f32> {
117+
let ar_model = ARp::new(self.ar_coefs.clone(), self.sigma, self.n, None, None);
118+
let ar_series = ar_model.sample_simd();
119+
120+
let ma_model = MAq::new(self.ma_coefs.clone(), self.sigma, self.n, None);
121+
let ma_series = ma_model.sample_simd();
122+
123+
let arma_series = &ar_series + &ma_series;
124+
125+
let mut result = arma_series;
126+
for _ in 0..self.d {
127+
result = inverse_difference_f32(&result);
128+
}
129+
130+
result
131+
}
132+
87133
fn n(&self) -> usize {
88134
self.n
89135
}

0 commit comments

Comments
 (0)