|
| 1 | +use impl_new_derive::ImplNew; |
| 2 | +use ndarray::Array1; |
| 3 | +use ndarray_rand::RandomExt; |
| 4 | +use rand_distr::Normal; |
| 5 | + |
| 6 | +use crate::stochastic::Sampling3D; |
| 7 | + |
| 8 | +#[derive(ImplNew)] |
| 9 | +pub struct HJM { |
| 10 | + pub a: fn(f64) -> f64, |
| 11 | + pub b: fn(f64) -> f64, |
| 12 | + pub p: fn(f64, f64) -> f64, |
| 13 | + pub q: fn(f64, f64) -> f64, |
| 14 | + pub v: fn(f64, f64) -> f64, |
| 15 | + pub alpha: fn(f64, f64) -> f64, |
| 16 | + pub sigma: fn(f64, f64) -> f64, |
| 17 | + pub n: usize, |
| 18 | + pub r0: Option<f64>, |
| 19 | + pub p0: Option<f64>, |
| 20 | + pub f0: Option<f64>, |
| 21 | + pub t: Option<f64>, |
| 22 | + pub m: Option<usize>, |
| 23 | +} |
| 24 | + |
| 25 | +impl Sampling3D<f64> for HJM { |
| 26 | + fn sample(&self) -> [Array1<f64>; 3] { |
| 27 | + let t_max = self.t.unwrap_or(1.0); |
| 28 | + let dt = self.t.unwrap_or(1.0) / (self.n - 1) as f64; |
| 29 | + let mut r = Array1::<f64>::zeros(self.n); |
| 30 | + let mut p = Array1::<f64>::zeros(self.n); |
| 31 | + let mut f = Array1::<f64>::zeros(self.n); |
| 32 | + |
| 33 | + let gn1 = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap()); |
| 34 | + let gn2 = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap()); |
| 35 | + let gn3 = Array1::random(self.n, Normal::new(0.0, dt.sqrt()).unwrap()); |
| 36 | + |
| 37 | + for i in 1..self.n { |
| 38 | + let t = i as f64 * dt; |
| 39 | + |
| 40 | + r[i] = r[i - 1] + (self.a)(t) * dt + (self.b)(t) * gn1[i - 1]; |
| 41 | + p[i] = |
| 42 | + p[i - 1] + (self.p)(t, t_max) * ((self.q)(t, t_max) * dt + (self.v)(t, t_max) * gn2[i - 1]); |
| 43 | + f[i] = f[i - 1] + (self.alpha)(t, t_max) * dt + (self.sigma)(t, t_max) * gn3[i - 1]; |
| 44 | + } |
| 45 | + |
| 46 | + [r, p, f] |
| 47 | + } |
| 48 | + |
| 49 | + fn n(&self) -> usize { |
| 50 | + self.n |
| 51 | + } |
| 52 | + |
| 53 | + fn m(&self) -> Option<usize> { |
| 54 | + self.m |
| 55 | + } |
| 56 | +} |
0 commit comments