Stochastic processes

The stochastic-rs-stochastic crate ships 120+ processes organised into six families. Every process implements ProcessExt<T>: sample() returns one path, sample_map(m, f) folds over m paths in parallel, and sample_par(m) returns Vec<Output> (one ndarray::Array1<T> per path for single-factor models).

Families

Examples

Geometric Brownian Motion

The textbook log-normal diffusion $dS_t = \mu S_t\,dt + \sigma S_t\,dW_t$.

use stochastic_rs::prelude::*;
use stochastic_rs::stochastic::diffusion::gbm::Gbm;

let p = Gbm::<f64>::new(0.05, 0.2, 1_000, Some(100.0), Some(1.0));
let path = p.sample();              // Array1<f64>, length 1000
let paths = p.sample_par(10_000);   // Vec<Array1<f64>>, 10_000 paths of length 1000

import stochastic_rs as srs

p = srs.Gbm(mu=0.05, sigma=0.2, n=1000, x0=100.0, t=1.0)
path = p.sample()                   # numpy.ndarray, shape (1000,)
paths = p.sample_par(10_000)        # shape (10_000, 1000)

Heston stochastic volatility

Two-factor model $dS_t = \mu S_t\,dt + \sqrt{V_t}\,S_t\,dW^1_t$, $dV_t = \kappa(\theta - V_t)\,dt + \sigma\sqrt{V_t}\,dW^2_t$ with $d\langle W^1, W^2\rangle_t = \rho\,dt$.

use stochastic_rs::simd_rng::Unseeded;
use stochastic_rs::stochastic::volatility::HestonPow;
use stochastic_rs::stochastic::volatility::heston::Heston;

let p = Heston::<f64, _>::new(
    Some(100.0),       // s0
    Some(0.04),        // v0
    2.0,               // kappa
    0.04,              // theta
    0.3,               // sigma (vol-of-vol)
    -0.7,              // rho
    0.03,              // mu
    1_000,             // n
    Some(1.0),         // t
    HestonPow::Sqrt,   // 0.5 = classic Heston; ThreeHalves = 3/2 model
    Some(true),        // use_sym: reflect variance to stay non-negative
    Unseeded,
);
let [s_path, v_path] = p.sample();   // [Array1<f64>; 2]

import stochastic_rs as srs

p = srs.Heston(kappa=2.0, theta=0.04, sigma=0.3, rho=-0.7, mu=0.03,
               n=1000, s0=100.0, v0=0.04, t=1.0)
s, v = p.sample()                   # both numpy arrays

Discretisation schemes

sample() integrates the variance with the Euler full-truncation scheme by default. For large κ, high |ρ|, or long maturities Euler carries a noticeable discretisation bias; switch to the Andersen (2008) Quadratic-Exponential (QE) scheme at compile time with .qe():

use stochastic_rs::simd_rng::Unseeded;
use stochastic_rs::stochastic::volatility::HestonPow;
use stochastic_rs::stochastic::volatility::heston::Heston;

let qe = Heston::<f64, _>::new(
    Some(100.0), Some(0.04), 2.0, 0.04, 0.3, -0.7, 0.03,
    1_000, Some(1.0), HestonPow::Sqrt, Some(true), Unseeded,
)
.qe();                               // -> Heston<f64, _, AndersenQe>
let [s_path, v_path] = qe.sample();

The scheme is a zero-sized type parameter (Heston<T, S, Sch>), so the choice is monomorphised with no per-step branch. QE is defined for the square-root (CIR) variance only, so keep HestonPow::Sqrt. Reference: Andersen, L. (2008), Efficient simulation of the Heston stochastic volatility model, Journal of Computational Finance 11(3). (.qe() is Rust-only; the Python Heston uses the Euler scheme.)

Fractional Brownian motion

Roughness controlled by the Hurst parameter $H \in (0, 1)$. $H = 0.5$ recovers standard Brownian motion; $H < 0.5$ is rough, $H > 0.5$ is persistent.

use stochastic_rs::simd_rng::Deterministic;
use stochastic_rs::stochastic::noise::fgn::Fgn;

let fgn = Fgn::<f64, _>::new(/* hurst */ 0.3, /* n */ 4096,
                             /* t */ Some(1.0), Deterministic::new(42));
let increments = fgn.sample();

import stochastic_rs as srs
import numpy as np

fgn = srs.Fgn(hurst=0.3, n=4096, t=1.0, seed=42)
increments = fgn.sample()
fbm = np.cumsum(increments)         # fBM = cumsum(fGN)

Common patterns

Construction

Every process follows the same new(args, seed) pattern (mirrors the distribution constructors — see the seeding concept page). The seed argument is a value implementing SeedExt:

use stochastic_rs::simd_rng::{Deterministic, Unseeded};

let p = Foo::<T, _>::new(/* params */, n, x0, t, Unseeded);                  // auto-seeded
let p = Foo::<T, _>::new(/* params */, n, x0, t, Deterministic::new(42));    // reproducible

Use Deterministic in tests so they replay bit-exactly; pass a shared Deterministic to several processes when chaining correlated factors — each call to seed.rng() / seed.rng_ext() atomically advances the internal counter, so the streams diverge despite the shared root seed. SeedExt::reseed(u64) swaps a Deterministic source in place to sweep seeds without rebuilding the process.

Sampling

let path  = p.sample();          // Array1<T>, length n
let paths = p.sample_par(m);     // Array2<T>, shape (m, n)

sample_par(m) uses Rayon. Each path gets a deterministic seed derived from the master seed plus the path index, so the result is thread-count-independent.

Acceleration

CPU SIMD (f64x4 / f32x8) is on by default. The fractional / fGN family (Fgn, Fbm, Fou, Fcir, Fgbm, FJacobi, Cfou, Cfgns, JumpFou) additionally picks a sampling backend at compile time with .on::<B>():

use stochastic_rs::simd_rng::Unseeded;
use stochastic_rs::stochastic::device::CudaNative;
use stochastic_rs::stochastic::noise::fgn::Fgn;

let fgn = Fgn::<f32, _>::new(0.7, 65_536, None, Unseeded);
let path = fgn.on::<CudaNative>().sample();   // FFT on the GPU, zero runtime branch

Backends: Cpu (default), CudaNative (cuda-native), CubeCl (gpu / gpu-cuda / gpu-wgpu), MetalNative (metal), Accelerate (accelerate) — each marker exists only when its feature is compiled, so an unavailable backend is a compile error, not a runtime fallback. See the Backends concept page and the Feature flags matrix.

Adding a new process

If you contribute, the four relevant SKILLs are:

• add-diffusion-process — GBM / OU / CIR / Heston-style
• add-jump-process — Merton-jump / Kou / Bates / compound-Poisson
• add-fractional-process — fBM / rough Bergomi / fractional CIR
• add-gpu-sampler — CUDA / Metal port

Each contains the file-by-file recipe.