pathwise_core/process/
markov.rs

1use std::marker::PhantomData;
2use nalgebra::SVector;
3use crate::state::{Diffusion, NoiseIncrement, State};
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<SDE<f64, impl Drift<f64>, impl Fn(f64, f64) -> f64 + Send + Sync>, crate::error::PathwiseError> {
107    if kappa <= 0.0 || theta <= 0.0 || sigma <= 0.0 {
108        return Err(crate::error::PathwiseError::InvalidParameters(
109            format!("CIR requires kappa, theta, sigma > 0; got kappa={}, theta={}, sigma={}", kappa, theta, sigma)
110        ));
111    }
112    if 2.0 * kappa * theta <= sigma * sigma {
113        return Err(crate::error::PathwiseError::FellerViolation(
114            format!("2*kappa*theta = {:.4} <= sigma^2 = {:.4}; boundary is reflecting in continuous time but clipping may introduce bias under discretization",
115                2.0 * kappa * theta, sigma * sigma)
116        ));
117    }
118    Ok(SDE::new(
119        move |x: &f64, _t: f64| kappa * (theta - x),
120        move |x: f64, _t: f64| sigma * x.max(0.0).sqrt(),
121    ))
122}
123
124/// NdSDE: N-dimensional SDE with vector state and vector noise.
125pub struct NdSDE<const N: usize, D, G> {
126    pub drift: D,
127    pub diffusion: G,
128}
129
130impl<const N: usize, D, G> NdSDE<N, D, G>
131where
132    D: Fn(&SVector<f64, N>, f64) -> SVector<f64, N> + Send + Sync,
133    G: Diffusion<SVector<f64, N>, SVector<f64, N>>,
134{
135    pub fn new(drift: D, diffusion: G) -> Self {
136        Self { drift, diffusion }
137    }
138}
139
140/// Diffusion term for the Heston model (log-price, variance).
141/// State: [log S, V]. Noise: [dW1, dW2] (independent).
142///
143/// Applies the lower-triangular Cholesky matrix:
144///   d(log S) += sqrt(V) * dW1
145///   dV       += xi * sqrt(V) * (rho * dW1 + sqrt(1-rho^2) * dW2)
146///
147/// Full truncation: V is clipped to 0 in diffusion computation.
148pub struct HestonDiffusion {
149    xi: f64,
150    rho: f64,
151    rho_perp: f64,  // sqrt(1 - rho^2)
152}
153
154impl HestonDiffusion {
155    pub fn new(xi: f64, rho: f64) -> Self {
156        Self { xi, rho, rho_perp: (1.0 - rho * rho).sqrt() }
157    }
158}
159
160impl Diffusion<SVector<f64, 2>, SVector<f64, 2>> for HestonDiffusion {
161    fn apply(&self, x: &SVector<f64, 2>, _t: f64, dw: &SVector<f64, 2>) -> SVector<f64, 2> {
162        let v = x[1].max(0.0);  // full truncation
163        let sv = v.sqrt();
164        SVector::from([
165            sv * dw[0],
166            sv * self.xi * (self.rho * dw[0] + self.rho_perp * dw[1]),
167        ])
168    }
169}
170
171/// Heston stochastic volatility model.
172/// State: [log S, V]; use exp(paths[.., .., 0]) to recover S.
173///
174/// d(log S) = (mu - V/2) dt + sqrt(V) dW1
175/// dV       = kappa * (theta - V) dt + xi * sqrt(V) * (rho * dW1 + sqrt(1-rho^2) * dW2)
176///
177/// Parameters:
178/// - `mu`: risk-neutral drift of log-price
179/// - `kappa`: variance mean-reversion speed
180/// - `theta`: long-run variance
181/// - `xi`: volatility of variance (vol of vol)
182/// - `rho`: correlation between price and variance Brownian motions (typically -0.7)
183pub fn heston(
184    mu: f64,
185    kappa: f64,
186    theta: f64,
187    xi: f64,
188    rho: f64,
189) -> NdSDE<2, impl Fn(&SVector<f64, 2>, f64) -> SVector<f64, 2> + Send + Sync, HestonDiffusion> {
190    NdSDE::new(
191        move |x: &SVector<f64, 2>, _t: f64| -> SVector<f64, 2> {
192            let v = x[1].max(0.0);
193            SVector::from([mu - v / 2.0, kappa * (theta - x[1])])
194        },
195        HestonDiffusion::new(xi, rho),
196    )
197}
198
199/// Correlated Ornstein-Uhlenbeck diffusion via Cholesky factor.
200/// apply: L * dW where L = chol(Sigma)
201pub struct CorrOuDiffusion<const N: usize> {
202    l: nalgebra::SMatrix<f64, N, N>,
203}
204
205impl<const N: usize> crate::state::Diffusion<nalgebra::SVector<f64, N>, nalgebra::SVector<f64, N>> for CorrOuDiffusion<N> {
206    fn apply(&self, _x: &nalgebra::SVector<f64, N>, _t: f64, dw: &nalgebra::SVector<f64, N>) -> nalgebra::SVector<f64, N> {
207        self.l * dw
208    }
209}
210
211/// N-dimensional correlated Ornstein-Uhlenbeck process.
212/// dX = theta*(mu - X) dt + L dW  where L = chol(Sigma)
213///
214/// Returns `Err(DimensionMismatch)` if sigma_mat Cholesky fails (not positive-definite).
215pub fn corr_ou<const N: usize>(
216    theta: f64,
217    mu: nalgebra::SVector<f64, N>,
218    sigma_mat: nalgebra::SMatrix<f64, N, N>,
219) -> Result<
220    NdSDE<N, impl Fn(&nalgebra::SVector<f64, N>, f64) -> nalgebra::SVector<f64, N> + Send + Sync, CorrOuDiffusion<N>>,
221    crate::error::PathwiseError,
222> {
223    let chol = nalgebra::Cholesky::new(sigma_mat).ok_or_else(|| {
224        crate::error::PathwiseError::DimensionMismatch(
225            "sigma_mat is not positive-definite (Cholesky failed)".into()
226        )
227    })?;
228    let l = chol.l();
229    Ok(NdSDE::new(
230        move |x: &nalgebra::SVector<f64, N>, _t: f64| -> nalgebra::SVector<f64, N> {
231            (mu - x) * theta
232        },
233        CorrOuDiffusion { l },
234    ))
235}
236
237#[cfg(test)]
238mod tests {
239    use super::*;
240    use crate::state::Increment;
241
242    #[test]
243    fn sde_evaluates_drift_and_diffusion() {
244        let sde = SDE::new(|x: &f64, _t: f64| -0.5 * x, |_x: f64, _t: f64| 1.0_f64);
245        assert_eq!(sde.eval_drift(&2.0, 0.0), -1.0);
246        // diffusion.apply with unit dw=1.0 gives g(x)*1.0 = 1.0
247        assert_eq!(sde.diffusion.apply(&2.0, 0.0, &1.0_f64), 1.0);
248    }
249
250    #[test]
251    fn sde_closures_capture_parameters() {
252        let theta = 0.7_f64;
253        let sde = SDE::new(
254            move |x: &f64, _t: f64| -theta * x,
255            |_x: f64, _t: f64| 1.0_f64,
256        );
257        assert!((sde.eval_drift(&1.0, 0.0) - (-0.7)).abs() < 1e-12);
258    }
259
260    #[test]
261    fn bm_has_zero_drift_unit_diffusion() {
262        let b = bm();
263        assert_eq!(b.eval_drift(&1.5, 0.5), 0.0);
264        // apply with unit dw=1.0 gives g(x,t)*1.0 = 1.0
265        assert_eq!(b.diffusion.apply(&1.5, 0.5, &1.0_f64), 1.0);
266    }
267
268    #[test]
269    fn gbm_drift_and_diffusion() {
270        let g = gbm(0.05, 0.2);
271        assert!((g.eval_drift(&2.0, 0.0) - 0.1).abs() < 1e-12);
272        // sigma*x = 0.2*2.0 = 0.4; apply with dw=1.0
273        assert!((g.diffusion.apply(&2.0, 0.0, &1.0_f64) - 0.4).abs() < 1e-12);
274    }
275
276    #[test]
277    fn ou_drift_and_diffusion() {
278        let o = ou(1.0, 0.0, 0.5);
279        assert!((o.eval_drift(&1.0, 0.0) - (-1.0)).abs() < 1e-12);
280        assert!((o.diffusion.apply(&1.0, 0.0, &1.0_f64) - 0.5).abs() < 1e-12);
281    }
282
283    #[test]
284    fn increment_roundtrip() {
285        let inc = Increment { dw: 0.3_f64, dz: 0.0_f64 };
286        assert!((inc.dw - 0.3).abs() < 1e-12);
287    }
288}
pathwise_core/process/markov.rs

pathwise_core/process/
markov.rs
