Skip to content

Commit 365c996

Browse files
dancixxdancixx
committed
feat: improve fractional process par generation
1 parent 88f13ce commit 365c996

File tree

6 files changed

+378
-129
lines changed

6 files changed

+378
-129
lines changed

‎src/ai/fou/fou_lstm_datasets.rs‎

Lines changed: 3 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ use ndarray::{s, Array1};
88
use ndarray_rand::RandomExt;
99
use rand_distr::Uniform;
1010

11-
use crate::stochastic::{diffusion::fou::FOU, noise::fgn::FGN, SamplingExt};
11+
use crate::stochastic::{diffusion::fou::FOU, SamplingExt};
1212

1313
pub fn test_vasicek_1_d(
1414
epoch_size: usize,
@@ -34,16 +34,7 @@ pub fn test_vasicek_1_d(
3434
for idx in 0..epoch_size {
3535
let hurst = hursts[idx];
3636
let theta = thetas[idx];
37-
let fou = FOU::new(
38-
theta,
39-
mu,
40-
sigma,
41-
n,
42-
Some(0.0),
43-
Some(16.0),
44-
None,
45-
FGN::<f64>::new(hurst, n - 1, Some(1.0), None),
46-
);
37+
let fou = FOU::<f64>::new(hurst, theta, mu, sigma, n, Some(0.0), Some(16.0), None);
4738
let mut path = fou.sample();
4839
let mean = path.mean().unwrap();
4940
let std = path.std(0.0);
@@ -88,16 +79,7 @@ pub fn test_vasicek_2_d(
8879
for idx in 0..epoch_size {
8980
let hurst = hursts[idx];
9081
let theta = thetas[idx];
91-
let fou = FOU::new(
92-
theta,
93-
mu,
94-
sigma,
95-
n,
96-
Some(0.0),
97-
Some(16.0),
98-
None,
99-
FGN::<f64>::new(hurst, n - 1, Some(1.0), None),
100-
);
82+
let fou = FOU::<f64>::new(hurst, theta, mu, sigma, n, Some(0.0), Some(16.0), None);
10183
let mut path = fou.sample();
10284
let mean = path.mean().unwrap();
10385
let std = path.std(0.0);

‎src/stats/fou_estimator.rs‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -428,15 +428,14 @@ impl FOUParameterEstimationV3 {
428428
#[cfg(test)]
429429
mod tests {
430430
use super::*;
431-
use crate::stochastic::{diffusion::fou::FOU, noise::fgn::FGN, SamplingExt};
431+
use crate::stochastic::{diffusion::fou::FOU, SamplingExt};
432432

433433
#[test]
434434
fn test_fou_parameter_estimation_v1() {
435435
const N: usize = 10000;
436436
const X0: f64 = 0.0;
437437

438-
let fgn = FGN::<f64>::new(0.70, 4095, Some(1.0), None);
439-
let fou = FOU::new(5.0, 2.8, 1.0, 4096, Some(X0), Some(16.0), None, fgn);
438+
let fou = FOU::<f64>::new(0.70, 5.0, 2.8, 1.0, 4096, Some(X0), Some(16.0), None);
440439
let path = fou.sample();
441440
let mut estimator = FOUParameterEstimationV1::new(path, FilterType::Daubechies, None);
442441

@@ -457,8 +456,7 @@ mod tests {
457456
const X0: f64 = 0.0;
458457
let delta = 1.0 / 256.0;
459458

460-
let fgn = FGN::<f64>::new(0.70, N - 1, Some(1.0), None);
461-
let fou = FOU::new(5.0, 2.8, 2.0, N, Some(X0), Some(16.0), None, fgn);
459+
let fou = FOU::<f64>::new(0.70, 5.0, 2.8, 2.0, N, Some(X0), Some(16.0), None);
462460
let path = fou.sample();
463461
let mut estimator = FOUParameterEstimationV2::new(path, Some(delta), N);
464462

‎src/stochastic/diffusion/fcir.rs‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,8 @@ impl FCIR<f64> {
4949
use_sym,
5050
m,
5151
fgn,
52+
#[cfg(feature = "cuda")]
53+
cuda: false,
5254
}
5355
}
5456
}
@@ -176,6 +178,8 @@ impl FCIR<f32> {
176178
use_sym,
177179
m,
178180
fgn,
181+
#[cfg(feature = "cuda")]
182+
cuda: false,
179183
}
180184
}
181185
}

‎src/stochastic/diffusion/fgbm.rs‎

Lines changed: 113 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,10 @@
1-
use ndarray::Array1;
1+
use ndarray::{Array1, Array2, Axis, Zip};
2+
use rayon::iter::{IntoParallelIterator, ParallelIterator};
23

34
use crate::stochastic::{noise::fgn::FGN, SamplingExt};
45

56
pub struct FGBM<T> {
7+
pub hurst: T,
68
pub mu: T,
79
pub sigma: T,
810
pub n: usize,
@@ -17,6 +19,32 @@ pub struct FGBM<T> {
1719

1820
#[cfg(feature = "f64")]
1921
impl FGBM<f64> {
22+
#[must_use]
23+
pub fn new(
24+
hurst: f64,
25+
mu: f64,
26+
sigma: f64,
27+
n: usize,
28+
x0: Option<f64>,
29+
t: Option<f64>,
30+
m: Option<usize>,
31+
) -> Self {
32+
let fgn = FGN::<f64>::new(hurst, n - 1, t, m);
33+
34+
Self {
35+
hurst,
36+
mu,
37+
sigma,
38+
n,
39+
x0,
40+
t,
41+
m,
42+
fgn,
43+
#[cfg(feature = "cuda")]
44+
cuda: false,
45+
}
46+
}
47+
2048
fn fgn(&self) -> Array1<f64> {
2149
#[cfg(feature = "cuda")]
2250
if self.cuda {
@@ -42,12 +70,38 @@ impl SamplingExt<f64> for FGBM<f64> {
4270
fgbm[0] = self.x0.unwrap_or(0.0);
4371

4472
for i in 1..self.n {
45-
fgbm[i] = fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1]
73+
fgbm[i] = fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1];
4674
}
4775

4876
fgbm
4977
}
5078

79+
/// Sample the Fractional Geometric Brownian Motion (FGBM) process in parallel
80+
fn sample_par(&self) -> ndarray::Array2<f64> {
81+
let n = self.n();
82+
let m = self.m().unwrap();
83+
let dt = self.t.unwrap_or(1.0) / (n - 1) as f64;
84+
let mut xs = Array2::zeros((m, n));
85+
let fgn = self.fgn.sample_par();
86+
87+
debug_assert_eq!(fgn.nrows(), m);
88+
debug_assert_eq!(fgn.ncols(), n - 1);
89+
90+
Zip::from(xs.axis_iter_mut(Axis(0)))
91+
.and(fgn.axis_iter(Axis(0)))
92+
.into_par_iter()
93+
.for_each(|(mut fgbm, fgn)| {
94+
fgbm[0] = self.x0.unwrap_or(0.0);
95+
96+
for i in 1..n {
97+
fgbm[i] =
98+
fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1];
99+
}
100+
});
101+
102+
xs
103+
}
104+
51105
/// Number of time steps
52106
fn n(&self) -> usize {
53107
self.n
@@ -67,6 +121,32 @@ impl SamplingExt<f64> for FGBM<f64> {
67121

68122
#[cfg(feature = "f32")]
69123
impl FGBM<f32> {
124+
#[must_use]
125+
pub fn new(
126+
hurst: f32,
127+
mu: f32,
128+
sigma: f32,
129+
n: usize,
130+
x0: Option<f32>,
131+
t: Option<f32>,
132+
m: Option<usize>,
133+
) -> Self {
134+
let fgn = FGN::<f32>::new(hurst, n - 1, t, m);
135+
136+
Self {
137+
hurst,
138+
mu,
139+
sigma,
140+
n,
141+
x0,
142+
t,
143+
m,
144+
fgn,
145+
#[cfg(feature = "cuda")]
146+
cuda: false,
147+
}
148+
}
149+
70150
fn fgn(&self) -> Array1<f32> {
71151
self.fgn.sample()
72152
}
@@ -83,12 +163,38 @@ impl SamplingExt<f32> for FGBM<f32> {
83163
fgbm[0] = self.x0.unwrap_or(0.0);
84164

85165
for i in 1..self.n {
86-
fgbm[i] = fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1]
166+
fgbm[i] = fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1];
87167
}
88168

89169
fgbm
90170
}
91171

172+
/// Sample the Fractional Geometric Brownian Motion (FGBM) process in parallel
173+
fn sample_par(&self) -> ndarray::Array2<f32> {
174+
let n = self.n();
175+
let m = self.m().unwrap();
176+
let dt = self.t.unwrap_or(1.0) / (n - 1) as f32;
177+
let mut xs = Array2::zeros((m, n));
178+
let fgn = self.fgn.sample_par();
179+
180+
debug_assert_eq!(fgn.nrows(), m);
181+
debug_assert_eq!(fgn.ncols(), n - 1);
182+
183+
Zip::from(xs.axis_iter_mut(Axis(0)))
184+
.and(fgn.axis_iter(Axis(0)))
185+
.into_par_iter()
186+
.for_each(|(mut fgbm, fgn)| {
187+
fgbm[0] = self.x0.unwrap_or(0.0);
188+
189+
for i in 1..n {
190+
fgbm[i] =
191+
fgbm[i - 1] + self.mu * fgbm[i - 1] * dt + self.sigma * fgbm[i - 1] * fgn[i - 1];
192+
}
193+
});
194+
195+
xs
196+
}
197+
92198
/// Number of time steps
93199
fn n(&self) -> usize {
94200
self.n
@@ -105,52 +211,28 @@ impl SamplingExt<f32> for FGBM<f32> {
105211
mod tests {
106212
use crate::{
107213
plot_1d,
108-
stochastic::{noise::fgn::FGN, SamplingExt, N, X0},
214+
stochastic::{SamplingExt, N, X0},
109215
};
110216

111217
use super::*;
112218

113219
#[test]
114220
fn fgbm_length_equals_n() {
115-
let fgbm = FGBM::new(
116-
1.0,
117-
0.8,
118-
N,
119-
Some(X0),
120-
Some(1.0),
121-
None,
122-
FGN::<f64>::new(0.7, N - 1, Some(1.0), None),
123-
);
221+
let fgbm = FGBM::<f64>::new(0.7, 1.0, 0.8, N, Some(X0), Some(1.0), None);
124222

125223
assert_eq!(fgbm.sample().len(), N);
126224
}
127225

128226
#[test]
129227
fn fgbm_starts_with_x0() {
130-
let fgbm = FGBM::new(
131-
1.0,
132-
0.8,
133-
N,
134-
Some(X0),
135-
Some(1.0),
136-
None,
137-
FGN::<f64>::new(0.7, N - 1, Some(1.0), None),
138-
);
228+
let fgbm = FGBM::<f64>::new(0.7, 1.0, 0.8, N, Some(X0), Some(1.0), None);
139229

140230
assert_eq!(fgbm.sample()[0], X0);
141231
}
142232

143233
#[test]
144234
fn fgbm_plot() {
145-
let fgbm = FGBM::new(
146-
1.0,
147-
0.8,
148-
N,
149-
Some(X0),
150-
Some(1.0),
151-
None,
152-
FGN::<f64>::new(0.7, N - 1, Some(1.0), None),
153-
);
235+
let fgbm = FGBM::<f64>::new(0.7, 1.0, 0.8, N, Some(X0), Some(1.0), None);
154236

155237
plot_1d!(
156238
fgbm.sample(),

0 commit comments

Comments
 (0)