pathwise_core/process/
markov.rs

1use crate::state::{Diffusion, NoiseIncrement, State};
2use nalgebra::SVector;
3use std::marker::PhantomData;
4
5pub trait Drift<S: State>: Fn(&S, f64) -> S + Send + Sync {}
6impl<S: State, F: Fn(&S, f64) -> S + Send + Sync> Drift<S> for F {}
7
8pub struct SDE<S: State + NoiseIncrement, D: Drift<S>, G: Diffusion<S, S>> {
9    pub drift: D,
10    pub diffusion: G,
11    _s: PhantomData<S>,
12}
13
14impl<S: State + NoiseIncrement, D: Drift<S>, G: Diffusion<S, S>> SDE<S, D, G> {
15    pub fn new(drift: D, diffusion: G) -> Self {
16        Self {
17            drift,
18            diffusion,
19            _s: PhantomData,
20        }
21    }
22
23    pub fn eval_drift(&self, x: &S, t: f64) -> S {
24        (self.drift)(x, t)
25    }
26}
27
28/// Standard Brownian motion: dX = dW
29///
30/// # Example
31///
32/// ```
33/// use pathwise_core::{simulate, bm, euler};
34///
35/// let b = bm();
36/// let scheme = euler();
37/// let paths = simulate(
38///     &b.drift,
39///     &b.diffusion,
40///     &scheme,
41///     0.0,   // x0
42///     0.0,   // t0
43///     1.0,   // t1
44///     5,     // n_paths
45///     100,   // n_steps
46///     0,     // seed
47/// ).expect("simulate failed");
48/// assert_eq!(paths.shape(), &[5, 101]);
49/// ```
50pub fn bm() -> SDE<f64, impl Drift<f64>, impl Diffusion<f64, f64>> {
51    SDE::new(|_x: &f64, _t: f64| 0.0_f64, |_x: f64, _t: f64| 1.0_f64)
52}
53
54/// Geometric Brownian motion: dX = mu*X dt + sigma*X dW
55///
56/// # Example
57///
58/// ```
59/// use pathwise_core::{simulate, gbm, euler};
60///
61/// let g = gbm(0.05, 0.2);
62/// let scheme = euler();
63/// let paths = simulate(&g.drift, &g.diffusion, &scheme, 100.0, 0.0, 1.0, 5, 100, 0).expect("simulate failed");
64/// assert_eq!(paths.shape(), &[5, 101]);
65/// // All starting values equal x0 = 100.0
66/// for i in 0..5 {
67///     assert!((paths[[i, 0]] - 100.0).abs() < 1e-12);
68/// }
69/// ```
70pub fn gbm(mu: f64, sigma: f64) -> SDE<f64, impl Drift<f64>, impl Diffusion<f64, f64>> {
71    SDE::new(
72        move |x: &f64, _t: f64| mu * x,
73        move |x: f64, _t: f64| sigma * x,
74    )
75}
76
77/// Ornstein-Uhlenbeck: dX = theta*(mu - X) dt + sigma dW
78///
79/// # Example
80///
81/// ```
82/// use pathwise_core::{simulate, ou, euler};
83///
84/// // Mean-reverting process: theta=2.0, long-run mean=1.0, vol=0.3
85/// let o = ou(2.0, 1.0, 0.3);
86/// let scheme = euler();
87/// let paths = simulate(&o.drift, &o.diffusion, &scheme, 0.0, 0.0, 1.0, 5, 100, 0).expect("simulate failed");
88/// assert_eq!(paths.shape(), &[5, 101]);
89/// ```
90pub fn ou(theta: f64, mu: f64, sigma: f64) -> SDE<f64, impl Drift<f64>, impl Diffusion<f64, f64>> {
91    SDE::new(
92        move |x: &f64, _t: f64| theta * (mu - x),
93        move |_x: f64, _t: f64| sigma,
94    )
95}
96
97/// Cox-Ingersoll-Ross: dX = kappa*(theta - X) dt + sigma*sqrt(X) dW
98///
99/// Requires `kappa > 0`, `theta > 0`, `sigma > 0`.
100/// Strict Feller condition: `2*kappa*theta > sigma^2`. Returns `Err(FellerViolation)` if not met.
101/// Simulation clips X to 0.0 when discretization produces negative values (full truncation).
102pub fn cir(
103    kappa: f64,
104    theta: f64,
105    sigma: f64,
106) -> Result<
107    SDE<f64, impl Drift<f64>, impl Fn(f64, f64) -> f64 + Send + Sync>,
108    crate::error::PathwiseError,
109> {
110    if kappa <= 0.0 || theta <= 0.0 || sigma <= 0.0 {
111        return Err(crate::error::PathwiseError::InvalidParameters(format!(
112            "CIR requires kappa, theta, sigma > 0; got kappa={}, theta={}, sigma={}",
113            kappa, theta, sigma
114        )));
115    }
116    if 2.0 * kappa * theta <= sigma * sigma {
117        return Err(crate::error::PathwiseError::FellerViolation(
118            format!("2*kappa*theta = {:.4} <= sigma^2 = {:.4}; boundary is reflecting in continuous time but clipping may introduce bias under discretization",
119                2.0 * kappa * theta, sigma * sigma)
120        ));
121    }
122    Ok(SDE::new(
123        move |x: &f64, _t: f64| kappa * (theta - x),
124        move |x: f64, _t: f64| sigma * x.max(0.0).sqrt(),
125    ))
126}
127
128/// NdSDE: N-dimensional SDE with vector state and vector noise.
129pub struct NdSDE<const N: usize, D, G> {
130    pub drift: D,
131    pub diffusion: G,
132}
133
134impl<const N: usize, D, G> NdSDE<N, D, G>
135where
136    D: Fn(&SVector<f64, N>, f64) -> SVector<f64, N> + Send + Sync,
137    G: Diffusion<SVector<f64, N>, SVector<f64, N>>,
138{
139    pub fn new(drift: D, diffusion: G) -> Self {
140        Self { drift, diffusion }
141    }
142}
143
144/// Diffusion term for the Heston model (log-price, variance).
145/// State: [log S, V]. Noise: [dW1, dW2] (independent).
146///
147/// Applies the lower-triangular Cholesky matrix:
148///   d(log S) += sqrt(V) * dW1
149///   dV       += xi * sqrt(V) * (rho * dW1 + sqrt(1-rho^2) * dW2)
150///
151/// Full truncation: V is clipped to 0 in diffusion computation.
152pub struct HestonDiffusion {
153    xi: f64,
154    rho: f64,
155    rho_perp: f64, // sqrt(1 - rho^2)
156}
157
158impl HestonDiffusion {
159    pub fn new(xi: f64, rho: f64) -> Self {
160        Self {
161            xi,
162            rho,
163            rho_perp: (1.0 - rho * rho).sqrt(),
164        }
165    }
166}
167
168impl Diffusion<SVector<f64, 2>, SVector<f64, 2>> for HestonDiffusion {
169    fn apply(&self, x: &SVector<f64, 2>, _t: f64, dw: &SVector<f64, 2>) -> SVector<f64, 2> {
170        let v = x[1].max(0.0); // full truncation
171        let sv = v.sqrt();
172        SVector::from([
173            sv * dw[0],
174            sv * self.xi * (self.rho * dw[0] + self.rho_perp * dw[1]),
175        ])
176    }
177}
178
179/// Heston stochastic volatility model.
180/// State: [log S, V]; use exp(paths[.., .., 0]) to recover S.
181///
182/// d(log S) = (mu - V/2) dt + sqrt(V) dW1
183/// dV       = kappa * (theta - V) dt + xi * sqrt(V) * (rho * dW1 + sqrt(1-rho^2) * dW2)
184///
185/// Parameters:
186/// - `mu`: risk-neutral drift of log-price
187/// - `kappa`: variance mean-reversion speed
188/// - `theta`: long-run variance
189/// - `xi`: volatility of variance (vol of vol)
190/// - `rho`: correlation between price and variance Brownian motions (typically -0.7)
191pub fn heston(
192    mu: f64,
193    kappa: f64,
194    theta: f64,
195    xi: f64,
196    rho: f64,
197) -> NdSDE<2, impl Fn(&SVector<f64, 2>, f64) -> SVector<f64, 2> + Send + Sync, HestonDiffusion> {
198    NdSDE::new(
199        move |x: &SVector<f64, 2>, _t: f64| -> SVector<f64, 2> {
200            let v = x[1].max(0.0);
201            SVector::from([mu - v / 2.0, kappa * (theta - x[1])])
202        },
203        HestonDiffusion::new(xi, rho),
204    )
205}
206
207/// Correlated Ornstein-Uhlenbeck diffusion via Cholesky factor.
208/// apply: L * dW where L = chol(Sigma)
209pub struct CorrOuDiffusion<const N: usize> {
210    l: nalgebra::SMatrix<f64, N, N>,
211}
212
213impl<const N: usize> crate::state::Diffusion<nalgebra::SVector<f64, N>, nalgebra::SVector<f64, N>>
214    for CorrOuDiffusion<N>
215{
216    fn apply(
217        &self,
218        _x: &nalgebra::SVector<f64, N>,
219        _t: f64,
220        dw: &nalgebra::SVector<f64, N>,
221    ) -> nalgebra::SVector<f64, N> {
222        self.l * dw
223    }
224}
225
226/// N-dimensional correlated Ornstein-Uhlenbeck process.
227/// dX = theta*(mu - X) dt + L dW  where L = chol(Sigma)
228///
229/// Returns `Err(DimensionMismatch)` if sigma_mat Cholesky fails (not positive-definite).
230#[allow(clippy::type_complexity)]
231pub fn corr_ou<const N: usize>(
232    theta: f64,
233    mu: nalgebra::SVector<f64, N>,
234    sigma_mat: nalgebra::SMatrix<f64, N, N>,
235) -> Result<
236    NdSDE<
237        N,
238        impl Fn(&nalgebra::SVector<f64, N>, f64) -> nalgebra::SVector<f64, N> + Send + Sync,
239        CorrOuDiffusion<N>,
240    >,
241    crate::error::PathwiseError,
242> {
243    let chol = nalgebra::Cholesky::new(sigma_mat).ok_or_else(|| {
244        crate::error::PathwiseError::DimensionMismatch(
245            "sigma_mat is not positive-definite (Cholesky failed)".into(),
246        )
247    })?;
248    let l = chol.l();
249    Ok(NdSDE::new(
250        move |x: &nalgebra::SVector<f64, N>, _t: f64| -> nalgebra::SVector<f64, N> {
251            (mu - x) * theta
252        },
253        CorrOuDiffusion { l },
254    ))
255}
256
257#[cfg(test)]
258mod tests {
259    use super::*;
260    use crate::state::Increment;
261
262    #[test]
263    fn sde_evaluates_drift_and_diffusion() {
264        let sde = SDE::new(|x: &f64, _t: f64| -0.5 * x, |_x: f64, _t: f64| 1.0_f64);
265        assert_eq!(sde.eval_drift(&2.0, 0.0), -1.0);
266        // diffusion.apply with unit dw=1.0 gives g(x)*1.0 = 1.0
267        assert_eq!(sde.diffusion.apply(&2.0, 0.0, &1.0_f64), 1.0);
268    }
269
270    #[test]
271    fn sde_closures_capture_parameters() {
272        let theta = 0.7_f64;
273        let sde = SDE::new(
274            move |x: &f64, _t: f64| -theta * x,
275            |_x: f64, _t: f64| 1.0_f64,
276        );
277        assert!((sde.eval_drift(&1.0, 0.0) - (-0.7)).abs() < 1e-12);
278    }
279
280    #[test]
281    fn bm_has_zero_drift_unit_diffusion() {
282        let b = bm();
283        assert_eq!(b.eval_drift(&1.5, 0.5), 0.0);
284        // apply with unit dw=1.0 gives g(x,t)*1.0 = 1.0
285        assert_eq!(b.diffusion.apply(&1.5, 0.5, &1.0_f64), 1.0);
286    }
287
288    #[test]
289    fn gbm_drift_and_diffusion() {
290        let g = gbm(0.05, 0.2);
291        assert!((g.eval_drift(&2.0, 0.0) - 0.1).abs() < 1e-12);
292        // sigma*x = 0.2*2.0 = 0.4; apply with dw=1.0
293        assert!((g.diffusion.apply(&2.0, 0.0, &1.0_f64) - 0.4).abs() < 1e-12);
294    }
295
296    #[test]
297    fn ou_drift_and_diffusion() {
298        let o = ou(1.0, 0.0, 0.5);
299        assert!((o.eval_drift(&1.0, 0.0) - (-1.0)).abs() < 1e-12);
300        assert!((o.diffusion.apply(&1.0, 0.0, &1.0_f64) - 0.5).abs() < 1e-12);
301    }
302
303    #[test]
304    fn increment_roundtrip() {
305        let inc = Increment {
306            dw: 0.3_f64,
307            dz: 0.0_f64,
308        };
309        assert!((inc.dw - 0.3).abs() < 1e-12);
310    }
311}
pathwise_core/process/markov.rs

pathwise_core/process/
markov.rs
