qualia_core_db/domains/financial/economics/
stochastic.rs

1//! Stochastic finance/economics kernels.
2//!
3//! The legacy public helpers are convenience APIs. The seeded `*_into` variants
4//! are the deterministic, caller-buffered kernels to prefer for tests, WASM, and
5//! any future hot-path integration.
6
7use rand_distr::{Distribution, StandardNormal};
8
9#[cfg(all(not(target_os = "android"), not(target_arch = "wasm32")))]
10use rayon::prelude::*;
11
12/// Default seed for deterministic Monte Carlo callers that do not supply their
13/// own scenario seed.
14pub const DEFAULT_MONTE_CARLO_SEED: u64 = 0x5144_4245_434f_4e31;
15
16#[derive(Debug, Clone, Copy, PartialEq, Eq)]
17pub enum StochasticError {
18    InvalidSteps,
19    InvalidPaths,
20    OutputBufferTooSmall,
21    NonFiniteInput,
22}
23
24#[derive(Debug, Clone, Copy)]
25struct SplitMix64 {
26    state: u64,
27}
28
29impl SplitMix64 {
30    fn new(seed: u64) -> Self {
31        Self { state: seed }
32    }
33
34    fn next_u64(&mut self) -> u64 {
35        self.state = self.state.wrapping_add(0x9E37_79B9_7F4A_7C15);
36        let mut z = self.state;
37        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
38        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
39        z ^ (z >> 31)
40    }
41
42    fn unit_open(&mut self) -> f64 {
43        // Map to (0, 1) with 53 random mantissa bits.
44        let bits = self.next_u64() >> 11;
45        ((bits as f64) + 0.5) * (1.0 / ((1u64 << 53) as f64))
46    }
47
48    fn gaussian(&mut self) -> f64 {
49        let u1 = self.unit_open();
50        let u2 = self.unit_open();
51        (-2.0 * u1.ln()).sqrt() * (core::f64::consts::TAU * u2).cos()
52    }
53}
54
55fn valid_gbm_inputs(
56    initial_price: f64,
57    drift: f64,
58    volatility: f64,
59    time_horizon: f64,
60    steps: usize,
61) -> bool {
62    steps > 0
63        && initial_price.is_finite()
64        && drift.is_finite()
65        && volatility.is_finite()
66        && time_horizon.is_finite()
67        && initial_price >= 0.0
68        && volatility >= 0.0
69        && time_horizon >= 0.0
70}
71
72fn gbm_step(current_price: f64, drift: f64, volatility: f64, dt: f64, z: f64) -> f64 {
73    current_price
74        * f64::exp((drift - 0.5 * volatility * volatility) * dt + volatility * dt.sqrt() * z)
75}
76
77/// Simulates a single path of Geometric Brownian Motion (GBM) using the ambient
78/// random source. Kept for compatibility with existing CLI/WASM surfaces.
79pub fn simulate_gbm_path(
80    initial_price: f64,
81    drift: f64,
82    volatility: f64,
83    time_horizon: f64,
84    steps: usize,
85) -> f64 {
86    let dt = time_horizon / steps as f64;
87    let mut current_price = initial_price;
88    let mut rng = rand::rng();
89
90    for _ in 0..steps {
91        let z: f64 = StandardNormal.sample(&mut rng);
92        current_price = gbm_step(current_price, drift, volatility, dt, z);
93    }
94
95    current_price
96}
97
98/// Deterministic single-path GBM simulation.
99pub fn simulate_gbm_path_seeded(
100    initial_price: f64,
101    drift: f64,
102    volatility: f64,
103    time_horizon: f64,
104    steps: usize,
105    seed: u64,
106) -> Result<f64, StochasticError> {
107    if !valid_gbm_inputs(initial_price, drift, volatility, time_horizon, steps) {
108        return if steps == 0 {
109            Err(StochasticError::InvalidSteps)
110        } else {
111            Err(StochasticError::NonFiniteInput)
112        };
113    }
114    let dt = time_horizon / steps as f64;
115    let mut current_price = initial_price;
116    let mut rng = SplitMix64::new(seed);
117
118    for _ in 0..steps {
119        current_price = gbm_step(current_price, drift, volatility, dt, rng.gaussian());
120    }
121    Ok(current_price)
122}
123
124/// Writes each GBM step price into `out`. Returns the number of prices written.
125pub fn simulate_gbm_steps_into(
126    initial_price: f64,
127    drift: f64,
128    volatility: f64,
129    time_horizon: f64,
130    steps: usize,
131    seed: u64,
132    out: &mut [f64],
133) -> Result<usize, StochasticError> {
134    if steps == 0 {
135        return Err(StochasticError::InvalidSteps);
136    }
137    if out.len() < steps {
138        return Err(StochasticError::OutputBufferTooSmall);
139    }
140    if !valid_gbm_inputs(initial_price, drift, volatility, time_horizon, steps) {
141        return Err(StochasticError::NonFiniteInput);
142    }
143
144    let dt = time_horizon / steps as f64;
145    let mut current_price = initial_price;
146    let mut rng = SplitMix64::new(seed);
147    for slot in out.iter_mut().take(steps) {
148        current_price = gbm_step(current_price, drift, volatility, dt, rng.gaussian());
149        *slot = current_price;
150    }
151    Ok(steps)
152}
153
154/// Runs a deterministic Monte Carlo VaR calculation into a caller-owned buffer.
155///
156/// `final_prices_out` is sorted ascending before returning, so callers may reuse
157/// it for quantile diagnostics. Returns `(paths, mean_final_price, var_95_loss)`.
158pub fn run_monte_carlo_var_seeded_into(
159    initial_price: f64,
160    drift: f64,
161    volatility: f64,
162    time_horizon: f64,
163    steps: usize,
164    paths: usize,
165    seed: u64,
166    final_prices_out: &mut [f64],
167) -> Result<(usize, f64, f64), StochasticError> {
168    if paths == 0 {
169        return Err(StochasticError::InvalidPaths);
170    }
171    if final_prices_out.len() < paths {
172        return Err(StochasticError::OutputBufferTooSmall);
173    }
174    if !valid_gbm_inputs(initial_price, drift, volatility, time_horizon, steps) {
175        return if steps == 0 {
176            Err(StochasticError::InvalidSteps)
177        } else {
178            Err(StochasticError::NonFiniteInput)
179        };
180    }
181
182    let mut sum = 0.0;
183    for (i, slot) in final_prices_out.iter_mut().take(paths).enumerate() {
184        let path_seed = seed ^ ((i as u64).wrapping_mul(0xD1B5_4A32_D192_ED03));
185        let final_price = simulate_gbm_path_seeded(
186            initial_price,
187            drift,
188            volatility,
189            time_horizon,
190            steps,
191            path_seed,
192        )?;
193        *slot = final_price;
194        sum += final_price;
195    }
196
197    let prices = &mut final_prices_out[..paths];
198    prices.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
199    let mean = sum / paths as f64;
200    let var_index = ((paths as f64 * 0.05).floor() as usize).min(paths - 1);
201    let var_95 = initial_price - prices[var_index];
202    Ok((paths, mean, var_95))
203}
204
205/// Runs a Monte Carlo simulation to calculate the expected end value and the
206/// Value at Risk (VaR) at a 95% confidence interval.
207///
208/// This is the compatibility facade. Prefer `run_monte_carlo_var_seeded_into`
209/// for deterministic, caller-buffered execution.
210pub fn run_monte_carlo_var(
211    initial_price: f64,
212    drift: f64,
213    volatility: f64,
214    time_horizon: f64,
215    steps: usize,
216    paths: usize,
217) -> (f64, f64) {
218    if paths == 0 || steps == 0 {
219        return (f64::NAN, f64::NAN);
220    }
221
222    #[cfg(all(not(target_os = "android"), not(target_arch = "wasm32")))]
223    let mut final_prices: Vec<f64> = (0..paths)
224        .into_par_iter()
225        .map(|_| simulate_gbm_path(initial_price, drift, volatility, time_horizon, steps))
226        .collect();
227
228    #[cfg(any(target_os = "android", target_arch = "wasm32"))]
229    let mut final_prices: Vec<f64> = (0..paths)
230        .into_iter()
231        .map(|_| simulate_gbm_path(initial_price, drift, volatility, time_horizon, steps))
232        .collect();
233
234    final_prices.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
235    let mean: f64 = final_prices.iter().sum::<f64>() / paths as f64;
236    let var_index = ((paths as f64 * 0.05).floor() as usize).min(paths.saturating_sub(1));
237    let var_95 = initial_price - final_prices[var_index];
238
239    (mean, var_95)
240}
241
242#[cfg(test)]
243mod tests {
244    use super::*;
245
246    #[test]
247    fn seeded_gbm_is_reproducible() {
248        let a = simulate_gbm_path_seeded(100.0, 0.05, 0.2, 1.0, 32, 7).unwrap();
249        let b = simulate_gbm_path_seeded(100.0, 0.05, 0.2, 1.0, 32, 7).unwrap();
250        assert_eq!(a, b);
251    }
252
253    #[test]
254    fn caller_buffered_path_writes_steps() {
255        let mut out = [0.0f64; 4];
256        let n = simulate_gbm_steps_into(100.0, 0.0, 0.0, 1.0, 4, 1, &mut out).unwrap();
257        assert_eq!(n, 4);
258        for price in out {
259            assert!(price > 0.0);
260        }
261    }
262
263    #[test]
264    fn seeded_var_uses_caller_buffer() {
265        let mut prices = [0.0f64; 16];
266        let (n, mean, var) =
267            run_monte_carlo_var_seeded_into(100.0, 0.02, 0.1, 1.0, 12, 16, 42, &mut prices)
268                .unwrap();
269        assert_eq!(n, 16);
270        assert!(mean.is_finite());
271        assert!(var.is_finite());
272        assert!(prices.windows(2).all(|w| w[0] <= w[1]));
273    }
274}
qualia_core_db/domains/financial/economics/stochastic.rs

qualia_core_db/domains/financial/economics/
stochastic.rs
