Skip to main content

optionstratlib_simulation/simulation/
walk_driver.rs

1/******************************************************************************
2   Author: Joaquín Béjar García
3   Email: jb@taunais.com
4   Date: 12/7/26
5******************************************************************************/
6//! Generic walk driver shared by every step generator.
7//!
8//! The three public generators (`generator_positive` here, plus
9//! `optionstratlib_market::chains::generator_optionchain` and `optionstratlib_market::series::generator_optionseries`,
10//! both behind the `synthetic` feature) used to duplicate the same
11//! algorithm: dispatch on [`WalkType`], extract the walk volatility, skip the
12//! first walker value, advance the x-step until expiration, build the next
13//! y-value and cap the result at `WalkParams::size`. [`walk_steps`] is the
14//! single implementation of that algorithm; each generator only supplies the
15//! closure that builds its `Y` value from the new price.
16
17use crate::error::SimulationError;
18use crate::simulation::steps::{Step, Xstep, Ystep};
19use crate::simulation::{WalkParams, WalkType};
20use optionstratlib_core::model::Positive;
21use optionstratlib_core::model::decimal::d_sqrt;
22use optionstratlib_core::utils::TimeFrame;
23use optionstratlib_core::utils::numeric::calculate_log_returns;
24use optionstratlib_pricing::volatility::{adjust_volatility, constant_volatility};
25use rust_decimal::Decimal;
26use std::convert::TryInto;
27use std::fmt::Display;
28use tracing::debug;
29
30/// Annualized volatility driving the walk.
31///
32/// For the nine synthetic walk types this is the variant's `volatility`
33/// parameter (see [`WalkType::volatility`]). For [`WalkType::Historical`] it
34/// is the constant volatility estimated from the log returns of the full
35/// price history, annualized from the history's timeframe.
36///
37/// # Errors
38///
39/// Propagates errors from the historical helpers (`calculate_log_returns`,
40/// `constant_volatility`, `adjust_volatility`) as [`SimulationError`], so
41/// the driver stays free of chain-layer error types.
42fn walk_volatility<Y>(
43    walk_params: &WalkParams<Positive, Y>,
44) -> Result<Option<Positive>, SimulationError>
45where
46    Y: TryInto<Positive> + Display + Clone,
47{
48    match &walk_params.walk_type {
49        WalkType::Historical {
50            timeframe, prices, ..
51        } => {
52            let log_returns: Vec<Decimal> = calculate_log_returns(prices)?;
53            let constant_volatility = constant_volatility(&log_returns)?;
54            let implied_volatility =
55                adjust_volatility(constant_volatility, *timeframe, TimeFrame::Year)?;
56            Ok(Some(implied_volatility))
57        }
58        walk_type => Ok(walk_type.volatility()),
59    }
60}
61
62/// Per-step expanding-window volatility estimates for a price path, free of
63/// look-ahead bias: the estimate at index `i` uses only the log returns of
64/// `prices[..=i]`.
65///
66/// This is the estimator [`walk_steps`] and [`walk_steps_par`] use for
67/// [`WalkType::Historical`], whose volatility is not a model parameter and has
68/// to be estimated from the series ([`WalkType::volatility`] returns `None` for
69/// it). It is public so that a consumer pricing a historical walk outside the
70/// driver — generating the path with
71/// [`crate::simulation::WalkTypeAble::generate_with_vol`] and building its own
72/// chains, say — uses the same numbers rather than a second copy of the
73/// mathematics that can drift release to release.
74///
75/// Uses the same sample-variance convention as
76/// [`optionstratlib_pricing::volatility::constant_volatility`], computed incrementally with
77/// checked arithmetic, and annualizes each estimate from `timeframe`.
78///
79/// # Contract
80///
81/// * The result has exactly `prices.len()` elements, so it indexes in step with
82///   the price path.
83/// * Element `i` is a function of `prices[..=i]` only. Appending prices never
84///   changes an earlier element.
85/// * Indices with fewer than two returns (0 and 1) cannot carry a sample
86///   variance and are backfilled with the first computable estimate, which
87///   keeps the series usable from step 0 without inventing a number from the
88///   future.
89/// * `Ok(None)` means no index had enough data — fewer than three prices.
90///
91/// # Errors
92///
93/// Propagates errors from [`optionstratlib_core::utils::numeric::calculate_log_returns`] /
94/// [`optionstratlib_pricing::volatility::adjust_volatility`] and surfaces arithmetic overflow as
95/// [`SimulationError`], keeping the driver free of chain-layer error types.
96///
97/// # Examples
98///
99/// ```
100/// use optionstratlib_simulation::simulation::expanding_window_vols;
101/// use optionstratlib_core::utils::TimeFrame;
102/// use optionstratlib_core::model::Positive;
103/// use optionstratlib_core::pos_or_panic;
104///
105/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
106/// // A flat prefix followed by a jump: the early estimates cannot see it.
107/// let mut prices = vec![Positive::HUNDRED; 6];
108/// prices.push(pos_or_panic!(200.0));
109///
110/// let vols = expanding_window_vols(&prices, TimeFrame::Day)?
111///     .expect("seven prices are enough for a sample variance");
112///
113/// assert_eq!(vols.len(), prices.len());
114/// assert_eq!(vols[4], Positive::ZERO); // flat window so far
115/// assert!(vols[6] > Positive::ZERO); // the jump is now in the window
116/// # Ok(())
117/// # }
118/// ```
119pub fn expanding_window_vols(
120    prices: &[Positive],
121    timeframe: TimeFrame,
122) -> Result<Option<Vec<Positive>>, SimulationError> {
123    let n = prices.len();
124    if n < 3 {
125        // Fewer than two returns anywhere: no sample variance is computable.
126        return Ok(None);
127    }
128    let log_returns: Vec<Decimal> = calculate_log_returns(prices)?;
129
130    let overflow =
131        |what: &str| SimulationError::walk_error(&format!("expanding-window {what} overflowed"));
132
133    let mut sum = Decimal::ZERO;
134    let mut sq_sum = Decimal::ZERO;
135    let mut raw: Vec<Option<Positive>> = Vec::with_capacity(n);
136    raw.push(None); // index 0 has no returns yet
137
138    for (i, r) in log_returns.iter().copied().enumerate() {
139        sum = sum
140            .checked_add(r)
141            .ok_or_else(|| overflow("running sum of log returns"))?;
142        let r_sq = r
143            .checked_mul(r)
144            .ok_or_else(|| overflow("squared log return"))?;
145        sq_sum = sq_sum
146            .checked_add(r_sq)
147            .ok_or_else(|| overflow("running sum of squared log returns"))?;
148
149        // `r` is the return between prices[i] and prices[i + 1]; after
150        // consuming it, `count` returns are available at price index i + 1.
151        let count = i.checked_add(1).ok_or_else(|| overflow("return count"))?;
152        if count < 2 {
153            raw.push(None);
154            continue;
155        }
156        let count_dec = Decimal::from(count as u64);
157        let denom = count_dec
158            .checked_sub(Decimal::ONE)
159            .ok_or_else(|| overflow("sample size minus one"))?;
160        // Sample variance via prefix sums: (Σr² − (Σr)²/n) / (n − 1),
161        // algebraically identical to the two-pass form in
162        // `constant_volatility`.
163        let sum_sq = sum
164            .checked_mul(sum)
165            .ok_or_else(|| overflow("squared running sum"))?;
166        let mean_sq_total = sum_sq
167            .checked_div(count_dec)
168            .ok_or_else(|| overflow("mean projection"))?;
169        let variance = sq_sum
170            .checked_sub(mean_sq_total)
171            .ok_or_else(|| overflow("variance numerator"))?
172            .checked_div(denom)
173            .ok_or_else(|| overflow("variance"))?
174            .max(Decimal::ZERO);
175        let std_dev = d_sqrt(variance, "simulation::walk_driver::volatility_sqrt")
176            .map_err(|_| overflow("volatility sqrt"))?;
177        let std_dev = Positive::new_decimal(std_dev).unwrap_or(Positive::ZERO);
178        let annualized = adjust_volatility(std_dev, timeframe, TimeFrame::Year)?;
179        raw.push(Some(annualized));
180    }
181
182    let first_computable = raw.iter().flatten().next().copied();
183    match first_computable {
184        None => Ok(None),
185        Some(fill) => Ok(Some(
186            raw.into_iter().map(|vol| vol.unwrap_or(fill)).collect(),
187        )),
188    }
189}
190
191/// Runs a walk and materializes it as a vector of [`Step`]s, delegating the
192/// construction of each y-value to `next_y`.
193///
194/// This is the shared driver behind [`generator_positive`] and, behind the
195/// `synthetic` feature, `optionstratlib_market::chains::generator_optionchain` and
196/// `optionstratlib_market::series::generator_optionseries`.
197///
198/// # Contract
199///
200/// * The returned vector always starts with `walk_params.init_step`.
201/// * If the walker yields no values beyond the initial one (e.g. a size-1
202///   walk), only the initial step is returned.
203/// * The walk is truncated when the x-step reaches expiration
204///   (`SimulationError::ExpirationReached`) or when `next_y` returns
205///   `Ok(None)`; any other step-advance error is propagated.
206/// * The result is truncated to at most `walk_params.size` steps.
207/// * With `walk_params.seed` set, the path, and therefore the result for a
208///   pure `next_y`, is identical on every run (see
209///   [`WalkParams`]); [`walk_steps_par`] returns the same steps.
210///
211/// # Parameters
212///
213/// * `walk_params` - The walk configuration; `walk_params.walk_type` selects
214///   the stochastic process via
215///   [`crate::simulation::WalkTypeAble::generate_with_vol`].
216/// * `next_y` - Builds the y-value for one step from `(new_price,
217///   step_volatility, advanced_x_step)`. Returning `Ok(None)` ends the walk
218///   gracefully at the previous step.
219///
220/// # Volatility
221///
222/// The volatility passed to `next_y` is per-step wherever the model provides
223/// one: the simulated vol path for stochastic-volatility walk types
224/// (`Garch`, `Heston`, `Custom`, `Telegraph`), an expanding-window estimate
225/// free of look-ahead bias for `Historical` walks, and the walk type's
226/// constant `volatility` parameter otherwise.
227///
228/// # Returns
229///
230/// * `Ok(Vec<Step<Positive, Y>>)` - The simulated walk.
231/// * `Err(E)` - If the walker, the volatility estimation, the x-step
232///   advance, or `next_y` fail.
233///
234/// # Errors
235///
236/// The driver is generic over its error type: every internal failure is a
237/// [`SimulationError`] lifted into `E` via `From`, so the simulation layer
238/// carries no chain-specific error types — the chain and series generators
239/// simply instantiate `E = ChainError`. Walker errors include
240/// `SimulationError::InsufficientHistoricalData` when a `Historical` walk
241/// has fewer prices than `walk_params.size`; `next_y` errors are returned
242/// verbatim.
243pub fn walk_steps<Y, E, F>(
244    walk_params: &WalkParams<Positive, Y>,
245    mut next_y: F,
246) -> Result<Vec<Step<Positive, Y>>, E>
247where
248    Y: TryInto<Positive> + Display + Clone,
249    E: From<SimulationError>,
250    F: FnMut(&Positive, Option<Positive>, &Xstep<Positive>) -> Result<Option<Y>, E>,
251{
252    debug!("{}", walk_params);
253    let path = walk_params
254        .walker
255        .generate_with_vol(walk_params)
256        .map_err(E::from)?;
257    let y_steps = path.prices;
258    if y_steps.len() <= 1 {
259        // Preserve the init-step invariant when the walker produces no
260        // values beyond the initial one; downstream consumers expect at
261        // least the initial step to be present.
262        return Ok(vec![walk_params.init_step.clone()]);
263    }
264
265    // Per-step volatilities: from the walker when the model simulates a vol
266    // path (Garch/Heston/Custom/Telegraph); for Historical walks, an
267    // expanding-window estimate over the walked prices that uses no future
268    // data. Otherwise a single constant from the walk type.
269    let step_vols: Option<Vec<Positive>> = match path.vols {
270        Some(vols) => Some(vols),
271        None => match &walk_params.walk_type {
272            WalkType::Historical { timeframe, .. } => {
273                expanding_window_vols(&y_steps, *timeframe).map_err(E::from)?
274            }
275            _ => None,
276        },
277    };
278    let constant_vol = if step_vols.is_some() {
279        None
280    } else {
281        walk_volatility(walk_params).map_err(E::from)?
282    };
283
284    let mut steps: Vec<Step<Positive, Y>> = vec![walk_params.init_step.clone()];
285    let mut previous_x_step = walk_params.init_step.x;
286    let mut y_index = *walk_params.ystep_ref().index();
287
288    // The first walker value duplicates the init step, so skip it.
289    for (i, y_step) in y_steps.iter().enumerate().skip(1) {
290        previous_x_step = match previous_x_step.next() {
291            Ok(x_step) => x_step,
292            // Reaching expiration is the normal end of a walk: truncate.
293            Err(SimulationError::ExpirationReached) => break,
294            // Any other step-advance failure is a real error: propagate.
295            Err(e) => return Err(E::from(e)),
296        };
297        let volatility = step_vols
298            .as_ref()
299            .and_then(|vols| vols.get(i).copied())
300            .or(constant_vol);
301        let Some(y_value) = next_y(y_step, volatility, &previous_x_step)? else {
302            break;
303        };
304        // The initial index is caller-supplied, so a walk seeded at
305        // `i32::MAX` overflows on its first step.
306        y_index = y_index.checked_add(1).ok_or_else(|| {
307            E::from(SimulationError::step_error(
308                "step index overflowed advancing the walk",
309            ))
310        })?;
311        steps.push(Step {
312            x: previous_x_step,
313            y: Ystep::new(y_index, y_value),
314        });
315    }
316
317    if steps.len() > walk_params.size {
318        debug!(
319            "walk produced {} steps, truncating to configured size {}",
320            steps.len(),
321            walk_params.size
322        );
323        steps.truncate(walk_params.size);
324    }
325    Ok(steps)
326}
327
328/// Parallel variant of [`walk_steps`]: identical contract and output, with
329/// the per-step `next_y` calls executed on the rayon thread pool.
330///
331/// The walk itself (walker output, per-step volatilities, x-step sequence)
332/// is computed serially and is identical to [`walk_steps`]; only the y-value
333/// construction — the dominant cost when `next_y` rebuilds an option chain
334/// or series — is fanned out. Output is deterministic and equal to the
335/// serial driver's for the same inputs because every per-step input is
336/// precomputed and order is preserved on collection.
337///
338/// # Sequential/parallel equivalence
339///
340/// * The price path is drawn once, serially, from a single generator before
341///   anything is fanned out: `deterministic_rng(seed)` when
342///   [`WalkParams::seed`] is `Some(seed)`, the thread RNG otherwise. No
343///   random number is drawn on a rayon worker, so the thread count and the
344///   scheduling cannot reach the random stream.
345/// * Nothing is reduced across steps in parallel; the y-values are assembled
346///   in step order, so no `Decimal` or `f64` result depends on the order in
347///   which workers finish.
348/// * Hence, for a seeded walk and a `next_y` that is a pure function of its
349///   arguments, `walk_steps_par(p, f)` equals `walk_steps(p, f)` bit for bit,
350///   on every run.
351/// * An unseeded walk draws a new path on every call, serial or parallel, so
352///   two calls agree only in distribution.
353///
354/// Semantics notes versus the serial driver:
355///
356/// * `next_y` must be `Fn + Sync` (stateless or internally synchronized)
357///   instead of `FnMut`.
358/// * `Ok(None)` still ends the walk at the previous step; steps past the
359///   first `None` are computed speculatively and their results — values AND
360///   errors — are discarded, exactly as the serial driver never evaluates
361///   them.
362/// * On failure, the first error in step order that precedes any stop
363///   signal is returned.
364///
365/// # Errors
366///
367/// Same as [`walk_steps`].
368pub fn walk_steps_par<Y, E, F>(
369    walk_params: &WalkParams<Positive, Y>,
370    next_y: F,
371) -> Result<Vec<Step<Positive, Y>>, E>
372where
373    Y: TryInto<Positive> + Display + Clone + Send,
374    E: From<SimulationError> + Send,
375    F: Fn(&Positive, Option<Positive>, &Xstep<Positive>) -> Result<Option<Y>, E> + Sync,
376{
377    use rayon::prelude::*;
378
379    debug!("{}", walk_params);
380    let path = walk_params
381        .walker
382        .generate_with_vol(walk_params)
383        .map_err(E::from)?;
384    let y_steps = path.prices;
385    if y_steps.len() <= 1 {
386        // Preserve the init-step invariant when the walker produces no
387        // values beyond the initial one; downstream consumers expect at
388        // least the initial step to be present.
389        return Ok(vec![walk_params.init_step.clone()]);
390    }
391
392    // Per-step volatilities: same policy as the serial driver.
393    let step_vols: Option<Vec<Positive>> = match path.vols {
394        Some(vols) => Some(vols),
395        None => match &walk_params.walk_type {
396            WalkType::Historical { timeframe, .. } => {
397                expanding_window_vols(&y_steps, *timeframe).map_err(E::from)?
398            }
399            _ => None,
400        },
401    };
402    let constant_vol = if step_vols.is_some() {
403        None
404    } else {
405        walk_volatility(walk_params).map_err(E::from)?
406    };
407
408    // Precompute the x-step sequence serially, honoring the
409    // expiration-truncation contract; Xstep::next is deterministic.
410    let mut x_steps: Vec<Xstep<Positive>> = Vec::with_capacity(y_steps.iter().skip(1).len());
411    let mut current_x = walk_params.init_step.x;
412    for _ in 1..y_steps.len() {
413        current_x = match current_x.next() {
414            Ok(x_step) => x_step,
415            // Reaching expiration is the normal end of a walk: truncate.
416            Err(SimulationError::ExpirationReached) => break,
417            // Any other step-advance failure is a real error: propagate.
418            Err(e) => return Err(E::from(e)),
419        };
420        x_steps.push(current_x);
421    }
422
423    // Every step's inputs are known up front: build the y-values in
424    // parallel, preserving order. Collect WITHOUT short-circuiting so a
425    // speculative error past a stop signal cannot mask the stop: the
426    // ordered assembly below decides what the serial driver would have
427    // observed.
428    let built: Vec<Result<Option<Y>, E>> = x_steps
429        .par_iter()
430        .enumerate()
431        .map(|(offset, x_step)| {
432            let price_index = offset
433                .checked_add(1)
434                .ok_or_else(|| E::from(SimulationError::walk_error("x-step offset overflowed")))?;
435            let y_step = y_steps.get(price_index).copied().ok_or_else(|| {
436                E::from(SimulationError::walk_error(
437                    "walker path shorter than x-steps",
438                ))
439            })?;
440            let volatility = step_vols
441                .as_ref()
442                .and_then(|vols| vols.get(price_index).copied())
443                .or(constant_vol);
444            next_y(&y_step, volatility, x_step)
445        })
446        .collect();
447
448    // Assemble serially in step order, mirroring the serial driver exactly:
449    // an `Err` is propagated only if it occurs before any `Ok(None)` stop
450    // signal; results past the stop are speculative and discarded.
451    let mut steps: Vec<Step<Positive, Y>> = vec![walk_params.init_step.clone()];
452    let mut y_index = *walk_params.ystep_ref().index();
453    for (x_step, y_result) in x_steps.into_iter().zip(built) {
454        let y_value = match y_result {
455            Ok(Some(y_value)) => y_value,
456            // The serial driver would have stopped here and never evaluated
457            // the remaining steps.
458            Ok(None) => break,
459            Err(e) => return Err(e),
460        };
461        // The initial index is caller-supplied, so a walk seeded at
462        // `i32::MAX` overflows on its first step.
463        y_index = y_index.checked_add(1).ok_or_else(|| {
464            E::from(SimulationError::step_error(
465                "step index overflowed advancing the walk",
466            ))
467        })?;
468        steps.push(Step {
469            x: x_step,
470            y: Ystep::new(y_index, y_value),
471        });
472    }
473
474    if steps.len() > walk_params.size {
475        debug!(
476            "walk produced {} steps, truncating to configured size {}",
477            steps.len(),
478            walk_params.size
479        );
480        steps.truncate(walk_params.size);
481    }
482    Ok(steps)
483}
484
485/// Generates a vector of `Step`s containing `Positive` x-values and `Positive` y-values.
486///
487/// Simulates the stochastic process selected by `walk_params.walk_type` for a plain
488/// positive value (typically an underlying price) without rebuilding option chains.
489///
490/// # Contract (shared with `optionstratlib_market::chains::generator_optionchain` and
491/// `optionstratlib_market::series::generator_optionseries`, both behind the `synthetic` feature)
492///
493/// * The returned vector always starts with `walk_params.init_step`.
494/// * If the walker yields no values beyond the initial one (e.g. a size-1 walk),
495///   only the initial step is returned.
496/// * The walk is truncated when the x-step reaches expiration
497///   (`SimulationError::ExpirationReached`); any other step-advance error is propagated.
498/// * The result is truncated to at most `walk_params.size` steps.
499///
500/// # Arguments
501///
502/// * `walk_params` - A reference to the `WalkParams` struct containing the walk parameters.
503///
504/// # Returns
505///
506/// * `Ok(Vec<Step<Positive, Positive>>)` - A vector of `Step`s representing the simulated walk.
507/// * `Err(ChainError)` - If the underlying simulator fails for any reason.
508///
509/// # Errors
510///
511/// Returns [`SimulationError`] if the random-walk generator fails — including
512/// `SimulationError::InsufficientHistoricalData` when a `Historical` walk has fewer
513/// prices than `walk_params.size`.
514pub fn generator_positive(
515    walk_params: &WalkParams<Positive, Positive>,
516) -> Result<Vec<Step<Positive, Positive>>, SimulationError> {
517    // The generator reports the error of its own layer (#511); a market
518    // market caller maps it into its own error where it uses it.
519    walk_steps(walk_params, |new_price, _volatility, _x_step| {
520        Ok(Some(*new_price))
521    })
522}
523
524#[cfg(test)]
525mod tests {
526    use super::*;
527    use crate::simulation::WalkTypeAble;
528    use crate::simulation::randomwalk::RandomWalk;
529    use optionstratlib_core::model::ExpirationDate;
530    use optionstratlib_core::pos_or_panic;
531    use optionstratlib_core::utils::time::convert_time_frame;
532    use optionstratlib_core::utils::{Len, TimeFrame};
533    use rust_decimal_macros::dec;
534
535    /// Seed of the simulated walks, so every run draws the same stream.
536    const SEED: u64 = 685;
537
538    #[derive(Clone)]
539    struct Walker {}
540    impl Walker {
541        fn new() -> Self {
542            Walker {}
543        }
544    }
545    impl WalkTypeAble<Positive, Positive> for Walker {}
546
547    #[test]
548    fn test_generator_positive() -> Result<(), Box<dyn std::error::Error>> {
549        let n_steps = 100;
550        let initial_price = Positive::HUNDRED;
551        let std_dev = pos_or_panic!(20.0);
552        let walker = Box::new(Walker::new());
553        let days = pos_or_panic!(30.0);
554
555        let walk_params = WalkParams {
556            size: n_steps,
557            init_step: Step {
558                x: Xstep::new(Positive::ONE, TimeFrame::Minute, ExpirationDate::Days(days)),
559                y: Ystep::new(0, initial_price),
560            },
561            walk_type: WalkType::GeometricBrownian {
562                dt: convert_time_frame(Positive::ONE / days, &TimeFrame::Minute, &TimeFrame::Day)?,
563                drift: dec!(0.0),
564                volatility: std_dev,
565            },
566            walker,
567            seed: Some(SEED),
568        };
569        let random_walk =
570            RandomWalk::new("Random Walk".to_string(), &walk_params, generator_positive)
571                .expect("random walk construction");
572        assert_eq!(random_walk.len(), n_steps);
573        Ok(())
574    }
575
576    #[test]
577    fn test_generator_positive_early_return() {
578        let initial_price = Positive::HUNDRED;
579        let walker = Box::new(Walker::new());
580
581        let walk_params = WalkParams {
582            size: 1, // Just one step to trigger early return
583            init_step: Step {
584                x: Xstep::new(
585                    Positive::ONE,
586                    TimeFrame::Minute,
587                    ExpirationDate::Days(pos_or_panic!(30.0)),
588                ),
589                y: Ystep::new(0, initial_price),
590            },
591            walk_type: WalkType::GeometricBrownian {
592                dt: pos_or_panic!(0.01),
593                drift: dec!(0.0),
594                volatility: pos_or_panic!(0.2),
595            },
596            walker,
597            seed: Some(SEED),
598        };
599
600        let steps = match generator_positive(&walk_params) {
601            Ok(steps) => steps,
602            Err(e) => panic!("generator_positive failed: {e}"),
603        };
604
605        // We should just get the initial step back
606        assert_eq!(steps.len(), 1);
607    }
608
609    /// A walker that returns an empty path; the driver must fall back to the
610    /// init-only walk instead of panicking or erroring.
611    #[test]
612    fn test_walk_steps_empty_walker_output() {
613        #[derive(Clone)]
614        struct EmptyWalker {}
615        impl WalkTypeAble<Positive, Positive> for EmptyWalker {
616            fn brownian(
617                &self,
618                _params: &WalkParams<Positive, Positive>,
619            ) -> Result<Vec<Positive>, SimulationError> {
620                Ok(Vec::new())
621            }
622        }
623
624        let walk_params = WalkParams {
625            size: 5,
626            init_step: Step {
627                x: Xstep::new(
628                    Positive::ONE,
629                    TimeFrame::Day,
630                    ExpirationDate::Days(pos_or_panic!(30.0)),
631                ),
632                y: Ystep::new(0, Positive::HUNDRED),
633            },
634            walk_type: WalkType::Brownian {
635                dt: pos_or_panic!(0.01),
636                drift: dec!(0.0),
637                volatility: pos_or_panic!(0.2),
638            },
639            walker: Box::new(EmptyWalker {}),
640            seed: None,
641        };
642
643        let steps = match generator_positive(&walk_params) {
644            Ok(steps) => steps,
645            Err(e) => panic!("empty walker output must not error: {e}"),
646        };
647        assert_eq!(steps.len(), 1, "empty walker output yields init-only walk");
648    }
649
650    /// Stochastic-vol walkers must expose a vol path aligned with prices
651    /// and actually varying when vol-of-vol is high.
652    #[test]
653    fn test_heston_with_vol_exposes_varying_vols() {
654        use rust_decimal_macros::dec;
655        let walker = Walker::new();
656        let params: WalkParams<Positive, Positive> = WalkParams {
657            size: 50,
658            init_step: Step {
659                x: Xstep::new(
660                    Positive::ONE,
661                    TimeFrame::Day,
662                    ExpirationDate::Days(pos_or_panic!(365.0)),
663                ),
664                y: Ystep::new(0, Positive::HUNDRED),
665            },
666            walk_type: WalkType::Heston {
667                dt: pos_or_panic!(0.004),
668                drift: dec!(0.0),
669                volatility: pos_or_panic!(0.2),
670                kappa: Positive::TWO,
671                theta: pos_or_panic!(0.04),
672                xi: Positive::TWO, // high vol-of-vol so the path must move
673                rho: dec!(-0.5),
674            },
675            walker: Box::new(Walker::new()),
676            seed: Some(SEED),
677        };
678        let path = match walker.heston_with_vol(&params) {
679            Ok(path) => path,
680            Err(e) => panic!("heston_with_vol failed: {e}"),
681        };
682        let vols = match path.vols {
683            Some(vols) => vols,
684            None => panic!("heston must expose a vol path"),
685        };
686        assert_eq!(vols.len(), path.prices.len());
687        assert_eq!(vols.first(), Some(&pos_or_panic!(0.2)));
688        let varying = vols.windows(2).any(|w| w[0] != w[1]);
689        assert!(varying, "high vol-of-vol Heston vol path cannot be flat");
690    }
691
692    /// Constant-vol walk types expose no vol path.
693    #[test]
694    fn test_generate_with_vol_none_for_gbm() {
695        use rust_decimal_macros::dec;
696        let walker = Walker::new();
697        let params: WalkParams<Positive, Positive> = WalkParams {
698            size: 10,
699            init_step: Step {
700                x: Xstep::new(
701                    Positive::ONE,
702                    TimeFrame::Day,
703                    ExpirationDate::Days(pos_or_panic!(30.0)),
704                ),
705                y: Ystep::new(0, Positive::HUNDRED),
706            },
707            walk_type: WalkType::GeometricBrownian {
708                dt: pos_or_panic!(0.01),
709                drift: dec!(0.0),
710                volatility: pos_or_panic!(0.2),
711            },
712            walker: Box::new(Walker::new()),
713            seed: Some(SEED),
714        };
715        let path = match walker.generate_with_vol(&params) {
716            Ok(path) => path,
717            Err(e) => panic!("generate_with_vol failed: {e}"),
718        };
719        assert!(path.vols.is_none());
720        assert_eq!(path.prices.len(), 10);
721    }
722
723    /// The expanding-window Historical estimator must not see future data:
724    /// with a flat prefix followed by a jump, the early estimates are zero
725    /// while the post-jump estimates are strictly positive.
726    #[test]
727    fn test_expanding_window_vols_no_look_ahead() {
728        let mut prices = vec![Positive::HUNDRED; 8];
729        prices.push(pos_or_panic!(200.0));
730        prices.push(pos_or_panic!(150.0));
731
732        let vols = match expanding_window_vols(&prices, TimeFrame::Day) {
733            Ok(Some(vols)) => vols,
734            Ok(None) => panic!("estimator returned no vols"),
735            Err(e) => panic!("estimator failed: {e}"),
736        };
737        assert_eq!(vols.len(), prices.len());
738        // Index 5 sees only the flat prefix: zero variance, zero vol.
739        assert_eq!(vols.get(5), Some(&Positive::ZERO));
740        // The last index has seen the jump: strictly positive vol.
741        match vols.last() {
742            Some(last) => assert!(*last > Positive::ZERO),
743            None => panic!("empty vols"),
744        }
745    }
746
747    /// The public contract of the estimator: element `i` depends on
748    /// `prices[..=i]` only, so extending the series must leave every earlier
749    /// estimate byte-identical. This is the property a consumer pricing a
750    /// historical walk outside the driver relies on (#423).
751    #[test]
752    fn test_expanding_window_vols_prefix_is_stable() {
753        let prices = vec![
754            Positive::HUNDRED,
755            pos_or_panic!(101.5),
756            pos_or_panic!(99.25),
757            pos_or_panic!(104.0),
758            pos_or_panic!(97.5),
759            pos_or_panic!(110.0),
760            pos_or_panic!(108.75),
761        ];
762
763        let full = match expanding_window_vols(&prices, TimeFrame::Day) {
764            Ok(Some(vols)) => vols,
765            Ok(None) => panic!("estimator returned no vols"),
766            Err(e) => panic!("estimator failed: {e}"),
767        };
768
769        // Every prefix of length >= 3 must reproduce the head of `full`.
770        for len in 3..=prices.len() {
771            let prefix = match prices.get(..len) {
772                Some(prefix) => prefix,
773                None => panic!("prefix of len {len} missing"),
774            };
775            let vols = match expanding_window_vols(prefix, TimeFrame::Day) {
776                Ok(Some(vols)) => vols,
777                Ok(None) => panic!("prefix of len {len} returned no vols"),
778                Err(e) => panic!("estimator failed on prefix {len}: {e}"),
779            };
780            assert_eq!(vols.len(), len, "one estimate per price");
781            match full.get(..len) {
782                Some(head) => assert_eq!(vols, head, "prefix of len {len} drifted"),
783                None => panic!("full series shorter than {len}"),
784            }
785        }
786    }
787
788    /// Fewer than three prices cannot carry a sample variance, and the
789    /// estimator says so instead of inventing one.
790    #[test]
791    fn test_expanding_window_vols_needs_three_prices() {
792        for len in 0..3 {
793            let prices = vec![Positive::HUNDRED; len];
794            match expanding_window_vols(&prices, TimeFrame::Day) {
795                Ok(None) => {}
796                Ok(Some(vols)) => panic!("len {len} unexpectedly produced {vols:?}"),
797                Err(e) => panic!("estimator failed on len {len}: {e}"),
798            }
799        }
800    }
801
802    /// The final estimate must equal the whole-series constant volatility,
803    /// annualized the same way: by the last index the expanding window is the
804    /// full window, so the incremental prefix-sum form and the two-pass form in
805    /// `constant_volatility` have to agree.
806    #[test]
807    fn test_expanding_window_vols_last_matches_constant_volatility() {
808        let prices = vec![
809            Positive::HUNDRED,
810            pos_or_panic!(102.0),
811            pos_or_panic!(98.5),
812            pos_or_panic!(103.25),
813            pos_or_panic!(101.0),
814        ];
815
816        let vols = match expanding_window_vols(&prices, TimeFrame::Day) {
817            Ok(Some(vols)) => vols,
818            Ok(None) => panic!("estimator returned no vols"),
819            Err(e) => panic!("estimator failed: {e}"),
820        };
821        let last = match vols.last() {
822            Some(last) => *last,
823            None => panic!("empty vols"),
824        };
825
826        let log_returns = match calculate_log_returns(&prices) {
827            Ok(returns) => returns,
828            Err(e) => panic!("log returns failed: {e}"),
829        };
830        let expected = match constant_volatility(&log_returns)
831            .and_then(|vol| adjust_volatility(vol, TimeFrame::Day, TimeFrame::Year))
832        {
833            Ok(expected) => expected,
834            Err(e) => panic!("reference volatility failed: {e}"),
835        };
836
837        let diff = (last.to_dec() - expected.to_dec()).abs();
838        assert!(
839            diff < dec!(0.000_000_1),
840            "expanding window {last} != constant volatility {expected}"
841        );
842    }
843
844    /// `next_y` returning `Ok(None)` must end the walk gracefully.
845    #[test]
846    fn test_walk_steps_next_y_none_stops_walk() {
847        let walker = Box::new(Walker::new());
848        let walk_params = WalkParams {
849            size: 10,
850            init_step: Step {
851                x: Xstep::new(
852                    Positive::ONE,
853                    TimeFrame::Day,
854                    ExpirationDate::Days(pos_or_panic!(30.0)),
855                ),
856                y: Ystep::new(0, Positive::HUNDRED),
857            },
858            walk_type: WalkType::GeometricBrownian {
859                dt: pos_or_panic!(0.01),
860                drift: dec!(0.0),
861                volatility: pos_or_panic!(0.2),
862            },
863            walker,
864            seed: Some(SEED),
865        };
866
867        let mut calls = 0;
868        let steps = match walk_steps(
869            &walk_params,
870            |price, _vol, _x| -> Result<Option<Positive>, SimulationError> {
871                calls += 1;
872                if calls > 3 {
873                    Ok(None)
874                } else {
875                    Ok(Some(*price))
876                }
877            },
878        ) {
879            Ok(steps) => steps,
880            Err(e) => panic!("walk_steps failed: {e}"),
881        };
882        // init + 3 accepted steps
883        assert_eq!(steps.len(), 4);
884    }
885
886    /// Multi-step behavior: prices, indices and time-to-expiry must all be
887    /// exact under a deterministic ramp walker.
888    #[test]
889    fn test_generator_positive_multi_step_behavior() {
890        use crate::simulation::walk_test_support::RampWalker;
891        use rust_decimal_macros::dec;
892
893        let size = 5;
894        let walk_params = WalkParams {
895            size,
896            init_step: Step {
897                x: Xstep::new(
898                    Positive::ONE,
899                    TimeFrame::Day,
900                    ExpirationDate::Days(pos_or_panic!(30.0)),
901                ),
902                y: Ystep::new(0, Positive::HUNDRED),
903            },
904            walk_type: WalkType::GeometricBrownian {
905                dt: pos_or_panic!(0.01),
906                drift: dec!(0.0),
907                volatility: pos_or_panic!(0.2),
908            },
909            walker: Box::new(RampWalker {
910                delta: Positive::TWO,
911            }),
912            seed: None,
913        };
914
915        let steps = match generator_positive(&walk_params) {
916            Ok(steps) => steps,
917            Err(e) => panic!("generator_positive failed: {e}"),
918        };
919        assert_eq!(steps.len(), size);
920        for (i, step) in steps.iter().enumerate() {
921            // Price follows the ramp exactly.
922            let expected_price = Positive::HUNDRED + Positive::TWO * i as f64;
923            assert_eq!(*step.y.value(), expected_price, "price at step {i}");
924            // Y index increments from the init step's index.
925            assert_eq!(*step.y.index(), i as i32, "y index at step {i}");
926            // Time-to-expiry decreases by exactly one day per step.
927            let days = match step.x.days_left() {
928                Ok(days) => days,
929                Err(e) => panic!("days_left failed at step {i}: {e}"),
930            };
931            assert_eq!(days, pos_or_panic!(30.0 - i as f64), "days at step {i}");
932        }
933    }
934
935    /// The walk must truncate exactly at expiration: with a 3-day expiry and
936    /// daily steps, a size-10 request yields init + 3 steps.
937    #[test]
938    fn test_generator_positive_truncates_at_expiration() {
939        use crate::simulation::walk_test_support::RampWalker;
940        use rust_decimal_macros::dec;
941
942        let walk_params = WalkParams {
943            size: 10,
944            init_step: Step {
945                x: Xstep::new(
946                    Positive::ONE,
947                    TimeFrame::Day,
948                    ExpirationDate::Days(pos_or_panic!(3.0)),
949                ),
950                y: Ystep::new(0, Positive::HUNDRED),
951            },
952            walk_type: WalkType::GeometricBrownian {
953                dt: pos_or_panic!(0.01),
954                drift: dec!(0.0),
955                volatility: pos_or_panic!(0.2),
956            },
957            walker: Box::new(RampWalker {
958                delta: Positive::ONE,
959            }),
960            seed: None,
961        };
962
963        let steps = match generator_positive(&walk_params) {
964            Ok(steps) => steps,
965            Err(e) => panic!("generator_positive failed: {e}"),
966        };
967        // init(3d) -> 2d -> 1d -> 0d, then ExpirationReached breaks.
968        assert_eq!(steps.len(), 4);
969        match steps.last() {
970            Some(last) => match last.x.days_left() {
971                Ok(days) => assert_eq!(days, Positive::ZERO),
972                Err(e) => panic!("days_left failed: {e}"),
973            },
974            None => panic!("empty steps"),
975        }
976    }
977
978    /// size = 0 yields the init-only walk from the shared driver.
979    #[test]
980    fn test_generator_positive_size_zero() {
981        use crate::simulation::walk_test_support::RampWalker;
982        use rust_decimal_macros::dec;
983
984        let walk_params = WalkParams {
985            size: 0,
986            init_step: Step {
987                x: Xstep::new(
988                    Positive::ONE,
989                    TimeFrame::Day,
990                    ExpirationDate::Days(pos_or_panic!(30.0)),
991                ),
992                y: Ystep::new(0, Positive::HUNDRED),
993            },
994            walk_type: WalkType::GeometricBrownian {
995                dt: pos_or_panic!(0.01),
996                drift: dec!(0.0),
997                volatility: pos_or_panic!(0.2),
998            },
999            walker: Box::new(RampWalker {
1000                delta: Positive::ONE,
1001            }),
1002            seed: None,
1003        };
1004
1005        let steps = match generator_positive(&walk_params) {
1006            Ok(steps) => steps,
1007            Err(e) => panic!("size-0 walk must not error: {e}"),
1008        };
1009        assert_eq!(steps.len(), 1);
1010    }
1011
1012    /// The parallel driver must produce exactly the serial driver's output
1013    /// for identical inputs.
1014    #[test]
1015    fn test_walk_steps_par_matches_serial() {
1016        use crate::simulation::walk_test_support::RampWalker;
1017        use rust_decimal_macros::dec;
1018
1019        let walk_params = WalkParams {
1020            size: 8,
1021            init_step: Step {
1022                x: Xstep::new(
1023                    Positive::ONE,
1024                    TimeFrame::Day,
1025                    ExpirationDate::Days(pos_or_panic!(30.0)),
1026                ),
1027                y: Ystep::new(0, Positive::HUNDRED),
1028            },
1029            walk_type: WalkType::GeometricBrownian {
1030                dt: pos_or_panic!(0.01),
1031                drift: dec!(0.0),
1032                volatility: pos_or_panic!(0.2),
1033            },
1034            walker: Box::new(RampWalker {
1035                delta: Positive::TWO,
1036            }),
1037            seed: None,
1038        };
1039
1040        let double =
1041            |price: &Positive,
1042             _vol: Option<Positive>,
1043             _x: &Xstep<Positive>|
1044             -> Result<Option<Positive>, SimulationError> { Ok(Some(*price * 2.0)) };
1045
1046        let serial = match walk_steps(&walk_params, double) {
1047            Ok(steps) => steps,
1048            Err(e) => panic!("serial driver failed: {e}"),
1049        };
1050        let parallel = match walk_steps_par(&walk_params, double) {
1051            Ok(steps) => steps,
1052            Err(e) => panic!("parallel driver failed: {e}"),
1053        };
1054
1055        assert_eq!(serial.len(), parallel.len());
1056        for (i, (s, p)) in serial.iter().zip(parallel.iter()).enumerate() {
1057            assert_eq!(s.y.value(), p.y.value(), "y value differs at step {i}");
1058            assert_eq!(s.y.index(), p.y.index(), "y index differs at step {i}");
1059            let (sd, pd) = match (s.x.days_left(), p.x.days_left()) {
1060                (Ok(sd), Ok(pd)) => (sd, pd),
1061                _ => panic!("days_left failed at step {i}"),
1062            };
1063            assert_eq!(sd, pd, "x days differ at step {i}");
1064        }
1065    }
1066
1067    /// Review regression (#419): an error from a speculative step past an
1068    /// `Ok(None)` stop must be discarded (the serial driver never evaluates
1069    /// those steps), while an error before the stop still propagates.
1070    #[test]
1071    fn test_walk_steps_par_error_after_stop_is_discarded() {
1072        use crate::simulation::walk_test_support::RampWalker;
1073        use rust_decimal_macros::dec;
1074
1075        let make_params = || WalkParams {
1076            size: 8,
1077            init_step: Step {
1078                x: Xstep::new(
1079                    Positive::ONE,
1080                    TimeFrame::Day,
1081                    ExpirationDate::Days(pos_or_panic!(30.0)),
1082                ),
1083                y: Ystep::new(0, Positive::HUNDRED),
1084            },
1085            walk_type: WalkType::GeometricBrownian {
1086                dt: pos_or_panic!(0.01),
1087                drift: dec!(0.0),
1088                volatility: pos_or_panic!(0.2),
1089            },
1090            walker: Box::new(RampWalker {
1091                delta: Positive::ONE, // prices 100, 101, 102, ...
1092            }),
1093            seed: None,
1094        };
1095
1096        // Stop at price 102 (step 2), error at price 103 (step 3): the
1097        // speculative error must be discarded and the walk succeed.
1098        let steps = match walk_steps_par(
1099            &make_params(),
1100            |price: &Positive, _vol, _x| -> Result<Option<Positive>, SimulationError> {
1101                if *price == pos_or_panic!(103.0) {
1102                    Err(SimulationError::invalid_parameters("speculative error"))
1103                } else if *price == pos_or_panic!(102.0) {
1104                    Ok(None)
1105                } else {
1106                    Ok(Some(*price))
1107                }
1108            },
1109        ) {
1110            Ok(steps) => steps,
1111            Err(e) => panic!("error past the stop signal must be discarded: {e}"),
1112        };
1113        // init + step 1 (price 101); price 102 stopped the walk.
1114        assert_eq!(steps.len(), 2);
1115
1116        // Error at price 101 (step 1), stop at price 102 (step 2): the
1117        // error precedes the stop and must propagate, as in the serial
1118        // driver.
1119        let result = walk_steps_par(
1120            &make_params(),
1121            |price: &Positive, _vol, _x| -> Result<Option<Positive>, SimulationError> {
1122                if *price == pos_or_panic!(101.0) {
1123                    Err(SimulationError::invalid_parameters("early error"))
1124                } else if *price == pos_or_panic!(102.0) {
1125                    Ok(None)
1126                } else {
1127                    Ok(Some(*price))
1128                }
1129            },
1130        );
1131        assert!(result.is_err(), "error before the stop must propagate");
1132    }
1133}