fix: fgn padded n/dt

Author: dancixx

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "stochastic-rs"
-version = "1.2.0"
+version = "1.2.1"
 edition = "2024"
 license = "MIT"
 description = "A Rust library for quant finance and simulating stochastic processes."

src/stochastic/diffusion/fou.rs

@@ -1,7 +1,7 @@
 //! # fOU
 //!
 //! $$
-//! dX_t=\kappa(\theta-X_t)\,dt+\sigma\,dB_t^H
+//! dX_t=\theta(\mu-X_t)\,dt+\sigma\,dB_t^H
 //! $$
 //!
 use ndarray::Array1;
@@ -13,9 +13,9 @@ use crate::traits::ProcessExt;
 pub struct FOU<T: FloatExt> {
   /// Hurst exponent controlling roughness and long-memory.
   pub hurst: T,
-  /// Long-run target level / model location parameter.
+  /// Mean-reversion speed.
   pub theta: T,
-  /// Drift / long-run mean-level parameter.
+  /// Long-run mean level.
   pub mu: T,
   /// Diffusion / noise scale parameter.
   pub sigma: T,
@@ -66,3 +66,42 @@ py_process_1d!(PyFOU, FOU,
   sig: (hurst, theta, mu, sigma, n, x0=None, t=None, dtype=None),
   params: (hurst: f64, theta: f64, mu: f64, sigma: f64, n: usize, x0: Option<f64>, t: Option<f64>)
 );
+
+#[cfg(test)]
+mod tests {
+  use super::FOU;
+  use crate::traits::ProcessExt;
+
+  #[test]
+  #[should_panic(expected = "n must be at least 2")]
+  fn fou_requires_at_least_two_points() {
+    let _ = FOU::<f64>::new(0.7, 1.0, 0.0, 0.2, 1, Some(0.0), Some(1.0));
+  }
+
+  #[test]
+  fn fou_sigma_zero_matches_deterministic_euler() {
+    let theta = 1.3_f64;
+    let mu = 0.8_f64;
+    let n = 129_usize;
+    let x0 = 0.2_f64;
+    let t = 1.0_f64;
+
+    let p = FOU::<f64>::new(0.7, theta, mu, 0.0, n, Some(x0), Some(t));
+    let x = p.sample();
+
+    let dt = t / (n as f64 - 1.0);
+    let mut expected = x0;
+    for i in 1..n {
+      expected = expected + theta * (mu - expected) * dt;
+      assert!((x[i] - expected).abs() < 1e-12, "mismatch at index {i}");
+    }
+  }
+
+  #[test]
+  fn fou_sample_is_finite() {
+    let p = FOU::<f64>::new(0.65, 1.0, 0.0, 0.5, 256, Some(0.1), Some(1.0));
+    let x = p.sample();
+    assert_eq!(x.len(), 256);
+    assert!(x.iter().all(|v| v.is_finite()));
+  }
+}

src/stochastic/noise/fgn/core.rs

@@ -18,7 +18,7 @@ use crate::traits::FloatExt;
 pub struct FGN<T: FloatExt> {
   /// Hurst exponent controlling roughness and long-memory.
   pub hurst: T,
-  /// Number of discrete simulation points (or samples).
+  /// Internal FFT length (power-of-two padded).
   pub n: usize,
   /// Total simulation horizon (defaults to 1 when omitted).
   pub t: Option<T>,
@@ -34,7 +34,8 @@ pub struct FGN<T: FloatExt> {
 
 impl<T: FloatExt> FGN<T> {
   pub fn dt(&self) -> T {
-    self.t.unwrap_or(T::one()) / T::from_usize_(self.n)
+    let step_count = self.out_len.max(1);
+    self.t.unwrap_or(T::one()) / T::from_usize_(step_count)
   }
 
   #[must_use]
@@ -78,13 +79,15 @@ impl<T: FloatExt> FGN<T> {
       };
     }
 
+    let scale_n = out_len.max(1);
+
     Self {
       hurst,
       n,
       offset,
       out_len,
       t,
-      scale: T::from_usize_(n).powf(-hurst) * t.unwrap_or(T::one()).powf(hurst),
+      scale: T::from_usize_(scale_n).powf(-hurst) * t.unwrap_or(T::one()).powf(hurst),
       sqrt_eigenvalues: Arc::new(sqrt_eigenvalues),
       fft_handler: Arc::new(FftHandler::new(2 * n)),
     }
@@ -114,3 +117,23 @@ impl<T: FloatExt> FGN<T> {
     fgn
   }
 }
+
+#[cfg(test)]
+mod tests {
+  use super::FGN;
+
+  #[test]
+  fn dt_and_scale_use_requested_length_not_fft_padding() {
+    let n = 1000_usize;
+    let h = 0.7_f64;
+    let t = 2.0_f64;
+    let fgn = FGN::<f64>::new(h, n, Some(t));
+
+    // Internal FFT length is padded, but time-step must follow requested n.
+    assert_eq!(fgn.n, 1024);
+    assert!((fgn.dt() - (t / n as f64)).abs() < 1e-15);
+
+    let expected_scale = (n as f64).powf(-h) * t.powf(h);
+    assert!((fgn.scale - expected_scale).abs() < 1e-15);
+  }
+}

src/stochastic/noise/fgn/cuda.rs

@@ -169,10 +169,11 @@ impl<T: FloatExt> FGN<T> {
   fn sample_cuda_f32(&self, m: usize) -> Result<Either<Array1<T>, Array2<T>>> {
     let n = self.n;
     let offset = self.offset;
+    let out_size = n - offset;
+    let scale_steps = out_size.max(1);
     let hurst = self.hurst.to_f64().unwrap();
     let t = self.t.unwrap_or(T::one()).to_f64().unwrap();
-    let scale = (n as f32).powf(-(hurst as f32)) * (t as f32).powf(hurst as f32);
-    let out_size = n - offset;
+    let scale = (scale_steps as f32).powf(-(hurst as f32)) * (t as f32).powf(hurst as f32);
 
     let need_init = {
       let guard = CUDA_CONTEXT_F32.lock().unwrap();
@@ -260,10 +261,11 @@ impl<T: FloatExt> FGN<T> {
   fn sample_cuda_f64(&self, m: usize) -> Result<Either<Array1<T>, Array2<T>>> {
     let n = self.n;
     let offset = self.offset;
+    let out_size = n - offset;
+    let scale_steps = out_size.max(1);
     let hurst = self.hurst.to_f64().unwrap();
     let t = self.t.unwrap_or(T::one()).to_f64().unwrap();
-    let scale = (n as f64).powf(-hurst) * t.powf(hurst);
-    let out_size = n - offset;
+    let scale = (scale_steps as f64).powf(-hurst) * t.powf(hurst);
 
     let need_init = {
       let guard = CUDA_CONTEXT_F64.lock().unwrap();