optionstratlib_simulation/simulation/
quadratic_rough_heston.rs

1/******************************************************************************
2   Author: Joaquín Béjar García
3   Email: jb@taunais.com
4   Date: 6/11/26
5******************************************************************************/
6
7//! The multi-factor quadratic rough Heston model: its kernel, its simulation by the
8//! explicit-implicit Euler scheme (A.1) and the Monte Carlo price of a vanilla on it.
9//! [`QuadraticRoughHestonParams`] documents the model for consumers.
10
11use crate::error::QuadraticRoughHestonError;
12use optionstratlib_core::model::Positive;
13use optionstratlib_core::model::decimal::{
14    d_div, d_sub, decimal_to_f64, f64_to_decimal, finite_decimal,
15};
16use optionstratlib_core::model::types::OptionStyle;
17use optionstratlib_core::utils::deterministic_rng;
18use rand::RngExt;
19use rand::rngs::StdRng;
20use rand_distr::{Distribution, StandardNormal};
21use rayon::prelude::*;
22use rust_decimal::Decimal;
23use rust_decimal_macros::dec;
24use std::f64::consts::PI;
25use std::num::NonZeroUsize;
26
27/// The excluded lower bound of `alpha`.
28const ALPHA_LOWER_BOUND: Decimal = dec!(0.5);
29
30/// The fewest paths with a standard error.
31const MIN_PATHS: usize = 2;
32
33/// An empty buffer able to hold `len` entries without growing.
34fn reserve<T>(len: usize, buffer: &'static str) -> Result<Vec<T>, QuadraticRoughHestonError> {
35    let mut reserved = Vec::new();
36    reserved
37        .try_reserve_exact(len)
38        .map_err(|_| QuadraticRoughHestonError::BufferTooLarge { buffer, len })?;
39    Ok(reserved)
40}
41
42#[cold]
43#[inline(never)]
44fn invalid_parameter(
45    parameter: &'static str,
46    value: Decimal,
47    reason: &'static str,
48) -> QuadraticRoughHestonError {
49    QuadraticRoughHestonError::InvalidParameter {
50        parameter,
51        value,
52        reason,
53    }
54}
55
56/// `raw` as a positive `Decimal`, or the typed error of a value that has none.
57fn kernel_decimal(
58    raw: f64,
59    quantity: &'static str,
60    index: usize,
61) -> Result<Decimal, QuadraticRoughHestonError> {
62    finite_decimal(raw)
63        .filter(|value| *value > Decimal::ZERO)
64        .ok_or(QuadraticRoughHestonError::KernelOutOfRange { quantity, index })
65}
66
67/// The factor weights and discounts of the geometric partition of the fractional kernel.
68///
69/// With `eta_i = x^(i - n/2)` for `i = 1..=n`, the measure
70/// `mu(dg) = g^(-alpha) dg / (Gamma(alpha) Gamma(1 - alpha))` is lumped on each
71/// `[eta_(i-1), eta_i]` into one atom, whose mass is the weight and whose mean is the
72/// discount:
73///
74/// ```text
75/// c_i     = x^((1-alpha)(i-n/2)) (1 - x^(alpha-1)) / ((1-alpha) Gamma(alpha) Gamma(1-alpha))
76/// gamma_i = (1-alpha) x^(i-1-n/2) (x^(2-alpha) - 1) / ((2-alpha) (x^(1-alpha) - 1))
77/// ```
78///
79/// [`value_at`](Self::value_at) is `K^n(t) = sum c_i exp(-gamma_i t)`, which approximates
80/// `K(t) = t^(alpha-1) / Gamma(alpha)` on a range of `t` set by `n` and `x`. The ratio that
81/// suits a given `n` is the caller's choice: there is no default.
82#[derive(Debug, Clone, PartialEq)]
83pub struct QuadraticRoughHestonKernel {
84    alpha: Decimal,
85    ratio: Decimal,
86    weights: Vec<Decimal>,
87    discounts: Vec<Decimal>,
88    raw_weights: Vec<f64>,
89    raw_discounts: Vec<f64>,
90}
91
92impl QuadraticRoughHestonKernel {
93    /// Builds the kernel of `factors` factors from `alpha` and the geometric `ratio`.
94    ///
95    /// # Errors
96    ///
97    /// - [`QuadraticRoughHestonError::AlphaOutOfRange`] when `alpha` is outside `(1/2, 1)`.
98    /// - [`QuadraticRoughHestonError::NoFactors`] when `factors` is zero.
99    /// - [`QuadraticRoughHestonError::RatioNotAboveOne`] when `ratio` is at most one.
100    /// - [`QuadraticRoughHestonError::KernelOutOfRange`] when a weight or a discount has no
101    ///   positive value in a `Decimal`.
102    /// - [`QuadraticRoughHestonError::BufferTooLarge`] when `factors` cannot be allocated.
103    /// - [`QuadraticRoughHestonError::Decimal`] when an input has no `f64` value.
104    pub fn try_new(
105        alpha: Decimal,
106        factors: usize,
107        ratio: Decimal,
108    ) -> Result<Self, QuadraticRoughHestonError> {
109        if alpha <= ALPHA_LOWER_BOUND || alpha >= Decimal::ONE {
110            return Err(QuadraticRoughHestonError::AlphaOutOfRange { alpha });
111        }
112        if factors == 0 {
113            return Err(QuadraticRoughHestonError::NoFactors);
114        }
115        if ratio <= Decimal::ONE {
116            return Err(QuadraticRoughHestonError::RatioNotAboveOne { ratio });
117        }
118
119        // `1 - alpha` and `ratio - 1` come from the `Decimal`, so an `alpha` just below one
120        // and a ratio just above it keep their distance from one, which an `f64` of the
121        // whole value would round away. `1 - alpha` is in `(0, 1/2)`,
122        // `Gamma(alpha) Gamma(1 - alpha)` is `pi / sin(pi (1 - alpha))`, and `exp_m1` keeps
123        // the differences of powers accurate for a ratio close to one.
124        let beta = decimal_to_f64(d_sub(
125            Decimal::ONE,
126            alpha,
127            "quadratic_rough_heston::kernel",
128        )?)?;
129        let log_ratio = decimal_to_f64(d_sub(
130            ratio,
131            Decimal::ONE,
132            "quadratic_rough_heston::kernel",
133        )?)?
134        .ln_1p(); // scan-banned: allow -- f64 `ln_1p` of a positive value: positive and finite, and it does not abort
135        let sin_pi_beta = (PI * beta).sin(); // scan-banned: allow -- f64 `sin` of a finite value: finite, and it does not abort
136        let weight_scale = -(-beta * log_ratio).exp_m1() * sin_pi_beta / (beta * PI);
137        let discount_scale = beta * ((1.0 + beta) * log_ratio).exp_m1()
138            / ((1.0 + beta) * (beta * log_ratio).exp_m1());
139        let half_factors = factors as f64 / 2.0;
140
141        let mut weights = reserve(factors, "kernel weights")?;
142        let mut discounts = reserve(factors, "kernel discounts")?;
143        let mut raw_weights = reserve(factors, "kernel weights")?;
144        let mut raw_discounts = reserve(factors, "kernel discounts")?;
145        for index in 1..=factors {
146            let exponent = index as f64 - half_factors;
147            let raw_weight = (beta * exponent * log_ratio).exp() * weight_scale; // scan-banned: allow -- f64 `exp`: returns inf on overflow, it does not abort; `kernel_decimal` rejects it
148            let raw_discount = ((exponent - 1.0) * log_ratio).exp() * discount_scale; // scan-banned: allow -- f64 `exp`: returns inf on overflow, it does not abort; `kernel_decimal` rejects it
149            weights.push(kernel_decimal(raw_weight, "weight", index)?);
150            discounts.push(kernel_decimal(raw_discount, "discount", index)?);
151            raw_weights.push(raw_weight);
152            raw_discounts.push(raw_discount);
153        }
154
155        Ok(Self {
156            alpha,
157            ratio,
158            weights,
159            discounts,
160            raw_weights,
161            raw_discounts,
162        })
163    }
164
165    /// The roughness exponent `alpha`, in `(1/2, 1)`.
166    #[must_use]
167    #[inline]
168    pub fn alpha(&self) -> Decimal {
169        self.alpha
170    }
171
172    /// The geometric ratio `x_n` of the partition.
173    #[must_use]
174    #[inline]
175    pub fn ratio(&self) -> Decimal {
176        self.ratio
177    }
178
179    /// The number of factors `n`.
180    #[must_use]
181    #[inline]
182    pub fn factors(&self) -> usize {
183        self.weights.len()
184    }
185
186    /// The weights `c_i`, from `i = 1` to `n`, rounded to a `Decimal`.
187    #[must_use]
188    #[inline]
189    pub fn weights(&self) -> &[Decimal] {
190        &self.weights
191    }
192
193    /// The discounts `gamma_i`, per year, from `i = 1` to `n`, rounded to a `Decimal`.
194    #[must_use]
195    #[inline]
196    pub fn discounts(&self) -> &[Decimal] {
197        &self.discounts
198    }
199
200    /// The approximated kernel `K^n(t) = sum c_i exp(-gamma_i t)` at `t` years.
201    ///
202    /// # Errors
203    ///
204    /// Returns [`QuadraticRoughHestonError::Decimal`] when `t` has no `f64` value or the
205    /// result does not fit a `Decimal`.
206    pub fn value_at(&self, t: Positive) -> Result<Decimal, QuadraticRoughHestonError> {
207        let t = decimal_to_f64(t.to_dec())?;
208        let value: f64 = self
209            .raw_weights
210            .iter()
211            .zip(&self.raw_discounts)
212            .map(|(weight, discount)| weight * (-discount * t).exp()) // scan-banned: allow -- f64 `exp` of a non-positive value: finite, and it does not abort
213            .sum();
214        Ok(f64_to_decimal(value)?)
215    }
216}
217
218/// The parameters of the model and its initial state.
219///
220/// # Model
221///
222/// The multi-factor approximation of the quadratic rough Heston model, driven by one
223/// Brownian motion `W` under the risk-neutral measure:
224///
225/// ```text
226/// dS   = S sqrt(V) dW,   V = a (Z - b)^2 + c,   Z = sum c_i Z^i
227/// dZ^i = (-gamma_i Z^i - lambda Z) dt + eta sqrt(V) dW,   Z^i_0 given
228/// ```
229///
230/// `(c_i, gamma_i)` come from the [`QuadraticRoughHestonKernel`]. The initial factors are
231/// free parameters, which replace the function `theta_0` of the rough model and set the
232/// initial term structure of `Z`.
233///
234/// # Scheme
235///
236/// [`simulate_quadratic_rough_heston_path`] and [`price_vanilla_quadratic_rough_heston`]
237/// use the explicit-implicit Euler scheme (A.1) of the calibration paper, with `dt = T / N`
238/// and `dW ~ N(0, dt)`:
239///
240/// ```text
241/// V_k       = a (Z_k - b)^2 + c,   Z_k = sum c_i Z^i_k
242/// S_(k+1)   = S_k + S_k sqrt(V_k) dW_k
243/// Z^i_(k+1) = (Z^i_k - lambda Z_k dt + eta sqrt(V_k) dW_k) / (1 + gamma_i dt)
244/// ```
245///
246/// The implicit denominator keeps the scheme stable for a large `gamma_i`.
247///
248/// # What this implementation does and does not do
249///
250/// - **The spot can go negative.** Scheme (A.1) is an arithmetic Euler step, so `S_(k+1)` is
251///   not above zero when `sqrt(V_k) dW_k <= -1`. It is not clamped, since a clamp would break
252///   the martingale property. [`QuadraticRoughHestonEstimate::non_positive_paths`] counts the
253///   affected paths, and on one of them a put can pay more than its strike. A low count does
254///   not show that the step is fine enough, as the last point below shows.
255/// - **The variance is not truncated.** The calibration paper's footnote 1 mentions a
256///   truncated square that secures the martingale property. This is the plain square, so a
257///   path whose state stops being finite returns [`QuadraticRoughHestonError::NonFinite`]
258///   with its step, and no value is replaced.
259/// - **There are no rates or dividends.** `S` is driftless, and a price is the undiscounted
260///   mean payoff, as with `r = q = 0`.
261/// - **No parameter set is a default.** The papers' values are for SPX and VIX, and nothing
262///   here is calibrated to any underlying.
263/// - **One run does not show that a price has converged.** A strong feedback gives heavy
264///   tails, and neither a small standard error, a low count of non positive paths nor a mean
265///   of `S_T` close to `S_0` rules out a wrong price. Take `alpha = 0.51`, `n = 10`, `x = 3.92`,
266///   `a = 1`, `b = 0.1`, `c = 0.01`, `lambda = 1.5`, `eta = 1`, `Z^i_0 = -0.05`,
267///   `S_0 = K = 100`, `T = 0.25` and 100,000 paths, as an illustration. At `N = 100`, about
268///   0.15% of the paths are non positive, and an ATM call is `150589 +/- 150584` at seed 1 and
269///   `5.32 +/- 0.02` at seed 3, while `N = 800` gives about `4.89`. So the small standard error
270///   at seed 3 hides a bias of about 20 of them, and the blown-up run has a mean of `S_T` about
271///   one standard error from `S_0`, because the standard error grows with the mean. **Compare
272///   several seeds and several step counts before trusting a price.**
273///
274/// # References
275///
276/// - Rosenbaum and Zhang, *Deep calibration of the quadratic rough Heston model*,
277///   [arXiv:2107.01611v2](https://arxiv.org/abs/2107.01611v2), Section 2 and Appendix A.
278/// - *Multi-asset market making under the quadratic rough Heston*,
279///   [arXiv:2212.10164v1](https://arxiv.org/abs/2212.10164v1), Section 2.1.
280///
281/// # Examples
282///
283/// ```rust
284/// use optionstratlib_core::model::types::OptionStyle;
285/// use optionstratlib_core::pos_or_panic;
286/// use optionstratlib_simulation::simulation::{
287///     QuadraticRoughHestonGrid, QuadraticRoughHestonKernel, QuadraticRoughHestonParams,
288///     price_vanilla_quadratic_rough_heston,
289/// };
290/// use rust_decimal_macros::dec;
291/// use std::num::NonZeroUsize;
292///
293/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
294/// // Illustrative values, calibrated to nothing.
295/// let kernel = QuadraticRoughHestonKernel::try_new(dec!(0.6), 5, dec!(5))?;
296/// let params = QuadraticRoughHestonParams::try_new(
297///     dec!(0.05),
298///     dec!(0.1),
299///     dec!(0.04),
300///     dec!(1),
301///     dec!(0.1),
302///     kernel,
303///     vec![dec!(0); 5],
304///     pos_or_panic!(100.0),
305/// )?;
306/// let grid = QuadraticRoughHestonGrid::try_new(dec!(0.25), 50)?;
307/// let paths = NonZeroUsize::try_from(200_usize)?;
308/// let estimate = price_vanilla_quadratic_rough_heston(
309///     &params,
310///     &grid,
311///     OptionStyle::Call,
312///     pos_or_panic!(100.0),
313///     paths,
314///     42,
315/// )?;
316/// assert!(estimate.price() > dec!(0));
317/// # Ok(())
318/// # }
319/// ```
320#[derive(Debug, Clone, PartialEq)]
321pub struct QuadraticRoughHestonParams {
322    a: Decimal,
323    b: Decimal,
324    c: Decimal,
325    lambda: Decimal,
326    eta: Decimal,
327    kernel: QuadraticRoughHestonKernel,
328    initial_factors: Vec<Decimal>,
329    spot: Positive,
330    initial_z: Decimal,
331    initial_variance: Decimal,
332    raw: RawParams,
333    raw_initial: State,
334}
335
336/// The parameters in the `f64` of the kernel.
337#[derive(Debug, Clone, Copy, PartialEq)]
338struct RawParams {
339    a: f64,
340    b: f64,
341    c: f64,
342    lambda: f64,
343    eta: f64,
344}
345
346impl QuadraticRoughHestonParams {
347    /// Validates the parameters and reads the initial state.
348    ///
349    /// `a >= 0`, `c > 0`, `lambda > 0`, `eta > 0` and a positive `spot`; `b` is free. The
350    /// `initial_factors` are the `Z^i_0`, one per factor of the `kernel`.
351    ///
352    /// # Errors
353    ///
354    /// - [`QuadraticRoughHestonError::InvalidParameter`] when a parameter breaks its condition.
355    /// - [`QuadraticRoughHestonError::InitialFactorCount`] when `initial_factors` does not hold
356    ///   one value per factor.
357    /// - [`QuadraticRoughHestonError::Decimal`] when the initial `Z` or `V` does not fit a
358    ///   `Decimal`, or a parameter has no `f64` value.
359    #[allow(clippy::too_many_arguments)]
360    pub fn try_new(
361        a: Decimal,
362        b: Decimal,
363        c: Decimal,
364        lambda: Decimal,
365        eta: Decimal,
366        kernel: QuadraticRoughHestonKernel,
367        initial_factors: Vec<Decimal>,
368        spot: Positive,
369    ) -> Result<Self, QuadraticRoughHestonError> {
370        if a < Decimal::ZERO {
371            return Err(invalid_parameter("a", a, "must not be negative"));
372        }
373        if c <= Decimal::ZERO {
374            return Err(invalid_parameter("c", c, "must be positive"));
375        }
376        if lambda <= Decimal::ZERO {
377            return Err(invalid_parameter("lambda", lambda, "must be positive"));
378        }
379        if eta <= Decimal::ZERO {
380            return Err(invalid_parameter("eta", eta, "must be positive"));
381        }
382        if spot.to_dec() <= Decimal::ZERO {
383            return Err(invalid_parameter("spot", spot.to_dec(), "must be positive"));
384        }
385        if initial_factors.len() != kernel.factors() {
386            return Err(QuadraticRoughHestonError::InitialFactorCount {
387                expected: kernel.factors(),
388                actual: initial_factors.len(),
389            });
390        }
391
392        let raw = RawParams {
393            a: decimal_to_f64(a)?,
394            b: decimal_to_f64(b)?,
395            c: decimal_to_f64(c)?,
396            lambda: decimal_to_f64(lambda)?,
397            eta: decimal_to_f64(eta)?,
398        };
399        let raw_factors = initial_factors
400            .iter()
401            .map(|factor| decimal_to_f64(*factor))
402            .collect::<Result<Vec<f64>, _>>()?;
403        let (z, variance) = read(&kernel.raw_weights, &raw, &raw_factors);
404        let raw_initial = State {
405            spot: decimal_to_f64(spot.to_dec())?,
406            z,
407            variance,
408            factors: raw_factors,
409        };
410        let initial_z = f64_to_decimal(z)?;
411        let initial_variance = f64_to_decimal(variance)?;
412
413        Ok(Self {
414            a,
415            b,
416            c,
417            lambda,
418            eta,
419            kernel,
420            initial_factors,
421            spot,
422            initial_z,
423            initial_variance,
424            raw,
425            raw_initial,
426        })
427    }
428
429    /// The strength `a` of the feedback on the variance.
430    #[must_use]
431    #[inline]
432    pub fn a(&self) -> Decimal {
433        self.a
434    }
435
436    /// The asymmetry `b` of the feedback.
437    #[must_use]
438    #[inline]
439    pub fn b(&self) -> Decimal {
440        self.b
441    }
442
443    /// The base level `c` of the variance, per year.
444    #[must_use]
445    #[inline]
446    pub fn c(&self) -> Decimal {
447        self.c
448    }
449
450    /// The mean reversion `lambda` of `Z`.
451    #[must_use]
452    #[inline]
453    pub fn lambda(&self) -> Decimal {
454        self.lambda
455    }
456
457    /// The volatility `eta` of the factors.
458    #[must_use]
459    #[inline]
460    pub fn eta(&self) -> Decimal {
461        self.eta
462    }
463
464    /// The kernel.
465    #[must_use]
466    #[inline]
467    pub fn kernel(&self) -> &QuadraticRoughHestonKernel {
468        &self.kernel
469    }
470
471    /// The initial factors `Z^i_0`.
472    #[must_use]
473    #[inline]
474    pub fn initial_factors(&self) -> &[Decimal] {
475        &self.initial_factors
476    }
477
478    /// The initial spot `S_0`.
479    #[must_use]
480    #[inline]
481    pub fn spot(&self) -> Positive {
482        self.spot
483    }
484
485    /// The initial `Z_0 = sum c_i Z^i_0`.
486    #[must_use]
487    #[inline]
488    pub fn initial_z(&self) -> Decimal {
489        self.initial_z
490    }
491
492    /// The initial variance `V_0 = a (Z_0 - b)^2 + c`, per year.
493    #[must_use]
494    #[inline]
495    pub fn initial_variance(&self) -> Decimal {
496        self.initial_variance
497    }
498}
499
500/// A uniform time grid: `steps` steps over `horizon` years.
501#[derive(Debug, Clone, Copy, PartialEq, Eq)]
502pub struct QuadraticRoughHestonGrid {
503    horizon: Decimal,
504    steps: usize,
505    dt: Decimal,
506}
507
508impl QuadraticRoughHestonGrid {
509    /// Builds the grid of `steps` steps over `horizon` years.
510    ///
511    /// # Errors
512    ///
513    /// - [`QuadraticRoughHestonError::NonPositiveStep`] when `horizon` is not above zero,
514    ///   `steps` is zero, or the step rounds to zero in a `Decimal`.
515    /// - [`QuadraticRoughHestonError::Decimal`] when the division fails.
516    pub fn try_new(horizon: Decimal, steps: usize) -> Result<Self, QuadraticRoughHestonError> {
517        let non_positive = || QuadraticRoughHestonError::NonPositiveStep { horizon, steps };
518        if horizon <= Decimal::ZERO || steps == 0 {
519            return Err(non_positive());
520        }
521        let dt = d_div(
522            horizon,
523            Decimal::from(steps),
524            "quadratic_rough_heston::grid",
525        )?;
526        if dt.is_zero() {
527            return Err(non_positive());
528        }
529        Ok(Self { horizon, steps, dt })
530    }
531
532    /// The horizon `T`, in years.
533    #[must_use]
534    #[inline]
535    pub fn horizon(&self) -> Decimal {
536        self.horizon
537    }
538
539    /// The number of steps `N`.
540    #[must_use]
541    #[inline]
542    pub fn steps(&self) -> usize {
543        self.steps
544    }
545
546    /// The step `dt = T / N`, in years, rounded to a `Decimal`.
547    #[must_use]
548    #[inline]
549    pub fn dt(&self) -> Decimal {
550        self.dt
551    }
552}
553
554/// One simulated path, by step of the grid from 0 to `N`.
555///
556/// `spot` and `z` can be negative. See [`QuadraticRoughHestonParams`].
557#[derive(Debug, Clone, PartialEq, Eq)]
558pub struct QuadraticRoughHestonPath {
559    spot: Vec<Decimal>,
560    variance: Vec<Decimal>,
561    z: Vec<Decimal>,
562    factors: Vec<Decimal>,
563    factor_count: usize,
564}
565
566impl QuadraticRoughHestonPath {
567    /// The spot `S_k` at each step.
568    #[must_use]
569    #[inline]
570    pub fn spot(&self) -> &[Decimal] {
571        &self.spot
572    }
573
574    /// The variance `V_k`, per year, at each step.
575    #[must_use]
576    #[inline]
577    pub fn variance(&self) -> &[Decimal] {
578        &self.variance
579    }
580
581    /// `Z_k` at each step.
582    #[must_use]
583    #[inline]
584    pub fn z(&self) -> &[Decimal] {
585        &self.z
586    }
587
588    /// The factors `Z^i_k` at `step`, or `None` past the end of the path.
589    #[must_use]
590    pub fn factors_at(&self, step: usize) -> Option<&[Decimal]> {
591        let start = step.checked_mul(self.factor_count)?;
592        let end = start.checked_add(self.factor_count)?;
593        self.factors.get(start..end)
594    }
595}
596
597/// The Monte Carlo price of a vanilla, with its diagnostics.
598#[derive(Debug, Clone, Copy, PartialEq, Eq)]
599pub struct QuadraticRoughHestonEstimate {
600    price: Decimal,
601    standard_error: Decimal,
602    mean_terminal_spot: Decimal,
603    terminal_spot_standard_error: Decimal,
604    paths: usize,
605    non_positive_paths: usize,
606}
607
608impl QuadraticRoughHestonEstimate {
609    /// The mean payoff per unit of underlying, undiscounted.
610    #[must_use]
611    #[inline]
612    pub fn price(&self) -> Decimal {
613        self.price
614    }
615
616    /// The standard error of [`price`](Self::price).
617    #[must_use]
618    #[inline]
619    pub fn standard_error(&self) -> Decimal {
620        self.standard_error
621    }
622
623    /// The mean of `S_T`, which estimates `S_0`.
624    #[must_use]
625    #[inline]
626    pub fn mean_terminal_spot(&self) -> Decimal {
627        self.mean_terminal_spot
628    }
629
630    /// The standard error of [`mean_terminal_spot`](Self::mean_terminal_spot).
631    #[must_use]
632    #[inline]
633    pub fn terminal_spot_standard_error(&self) -> Decimal {
634        self.terminal_spot_standard_error
635    }
636
637    /// The number of simulated paths.
638    #[must_use]
639    #[inline]
640    pub fn paths(&self) -> usize {
641        self.paths
642    }
643
644    /// The paths with `S_k <= 0` at some step. See [`QuadraticRoughHestonParams`].
645    #[must_use]
646    #[inline]
647    pub fn non_positive_paths(&self) -> usize {
648        self.non_positive_paths
649    }
650}
651
652/// `Z = sum c_i Z^i` and `V = a (Z - b)^2 + c`.
653#[inline]
654fn read(weights: &[f64], raw: &RawParams, factors: &[f64]) -> (f64, f64) {
655    let z: f64 = weights
656        .iter()
657        .zip(factors)
658        .map(|(weight, factor)| weight * factor)
659        .sum();
660    let shifted = z - raw.b;
661    (z, raw.a * shifted * shifted + raw.c)
662}
663
664/// The state of the scheme at one step, with its `Z` and `V`.
665#[derive(Debug, Clone, PartialEq)]
666struct State {
667    spot: f64,
668    z: f64,
669    variance: f64,
670    factors: Vec<f64>,
671}
672
673/// A value of the state that stopped being finite: the kernel's small form of
674/// [`QuadraticRoughHestonError::NonFinite`].
675#[derive(Debug, Clone, Copy, PartialEq, Eq)]
676struct NonFiniteValue {
677    quantity: &'static str,
678    step: usize,
679}
680
681impl NonFiniteValue {
682    #[cold]
683    #[inline(never)]
684    fn into_error(self) -> QuadraticRoughHestonError {
685        QuadraticRoughHestonError::NonFinite {
686            quantity: self.quantity,
687            step: self.step,
688        }
689    }
690}
691
692impl State {
693    /// Every weight is positive, so a factor that is not finite leaves `z` not finite.
694    #[inline]
695    fn ensure_finite(&self, step: usize) -> Result<(), NonFiniteValue> {
696        let quantity = if !self.spot.is_finite() {
697            "spot"
698        } else if !self.z.is_finite() {
699            "z"
700        } else if !self.variance.is_finite() {
701            "variance"
702        } else {
703            return Ok(());
704        };
705        Err(NonFiniteValue { quantity, step })
706    }
707}
708
709/// One step of the scheme (A.1) from `state` into `next`, for the increment `dw` over `dt`.
710///
711/// `next` is scratch of the same shape: nothing of it is read, and nothing allocates.
712#[inline(always)]
713fn step(params: &QuadraticRoughHestonParams, state: &State, dw: f64, dt: f64, next: &mut State) {
714    let raw = &params.raw;
715    let diffusion = state.variance.sqrt() * dw; // scan-banned: allow -- f64 `sqrt`: returns NaN for negative input, it does not abort; `ensure_finite` rejects it
716    let forcing = -raw.lambda * state.z * dt + raw.eta * diffusion;
717    next.spot = state.spot + state.spot * diffusion;
718    for ((next_factor, factor), discount) in next
719        .factors
720        .iter_mut()
721        .zip(&state.factors)
722        .zip(&params.kernel.raw_discounts)
723    {
724        *next_factor = (factor + forcing) / (1.0 + discount * dt);
725    }
726    (next.z, next.variance) = read(&params.kernel.raw_weights, raw, &next.factors);
727}
728
729/// The increments of one path, and the `next` scratch of its steps.
730struct Walk<'a> {
731    params: &'a QuadraticRoughHestonParams,
732    dt: f64,
733    sqrt_dt: f64,
734    rng: StdRng,
735    state: State,
736    scratch: State,
737}
738
739impl<'a> Walk<'a> {
740    fn new(params: &'a QuadraticRoughHestonParams, (dt, sqrt_dt): (f64, f64), seed: u64) -> Self {
741        Self {
742            params,
743            dt,
744            sqrt_dt,
745            rng: deterministic_rng(seed),
746            state: params.raw_initial.clone(),
747            scratch: params.raw_initial.clone(),
748        }
749    }
750
751    /// Moves to grid step `to`, drawing its normal increment.
752    #[inline]
753    fn advance(&mut self, to: usize) -> Result<(), NonFiniteValue> {
754        let normal: f64 = StandardNormal.sample(&mut self.rng);
755        self.apply(self.sqrt_dt * normal, to)
756    }
757
758    /// Moves to grid step `to` on the increment `dw`.
759    #[inline]
760    fn apply(&mut self, dw: f64, to: usize) -> Result<(), NonFiniteValue> {
761        step(self.params, &self.state, dw, self.dt, &mut self.scratch);
762        std::mem::swap(&mut self.state, &mut self.scratch);
763        self.state.ensure_finite(to)
764    }
765}
766
767/// `dt` and its square root, in `f64`.
768fn raw_step(grid: &QuadraticRoughHestonGrid) -> Result<(f64, f64), QuadraticRoughHestonError> {
769    let dt = decimal_to_f64(grid.dt)?;
770    Ok((dt, dt.sqrt())) // scan-banned: allow -- f64 `sqrt` of a positive finite value: finite, and it does not abort
771}
772
773/// Simulates one path of the model by the scheme (A.1), from `seed`.
774///
775/// The same arguments give the same path on the same build and platform. Across platforms the
776/// `f64` `exp`, `ln` and `sin` may round differently. See [`QuadraticRoughHestonParams`] for
777/// the model and the scheme.
778///
779/// # Errors
780///
781/// - [`QuadraticRoughHestonError::NonFinite`] when the state stops being finite, with its step.
782/// - [`QuadraticRoughHestonError::BufferTooLarge`] when the path cannot be allocated.
783/// - [`QuadraticRoughHestonError::Decimal`] when a value does not fit a `Decimal`.
784#[tracing::instrument(level = "debug", skip(params))]
785pub fn simulate_quadratic_rough_heston_path(
786    params: &QuadraticRoughHestonParams,
787    grid: &QuadraticRoughHestonGrid,
788    seed: u64,
789) -> Result<QuadraticRoughHestonPath, QuadraticRoughHestonError> {
790    let factor_count = params.kernel.factors();
791    let points = grid
792        .steps
793        .checked_add(1)
794        .ok_or(QuadraticRoughHestonError::BufferTooLarge {
795            buffer: "path",
796            len: grid.steps,
797        })?;
798    let factor_points =
799        points
800            .checked_mul(factor_count)
801            .ok_or(QuadraticRoughHestonError::BufferTooLarge {
802                buffer: "path factors",
803                len: points,
804            })?;
805    let mut path = QuadraticRoughHestonPath {
806        spot: reserve(points, "path spots")?,
807        variance: reserve(points, "path variances")?,
808        z: reserve(points, "path z")?,
809        factors: reserve(factor_points, "path factors")?,
810        factor_count,
811    };
812
813    let mut walk = Walk::new(params, raw_step(grid)?, seed);
814    record(&mut path, &walk.state)?;
815    for to in 1..=grid.steps {
816        walk.advance(to).map_err(NonFiniteValue::into_error)?;
817        record(&mut path, &walk.state)?;
818    }
819    Ok(path)
820}
821
822/// Appends `state` to `path`.
823fn record(
824    path: &mut QuadraticRoughHestonPath,
825    state: &State,
826) -> Result<(), QuadraticRoughHestonError> {
827    path.spot.push(f64_to_decimal(state.spot)?);
828    path.variance.push(f64_to_decimal(state.variance)?);
829    path.z.push(f64_to_decimal(state.z)?);
830    for factor in &state.factors {
831        path.factors.push(f64_to_decimal(*factor)?);
832    }
833    Ok(())
834}
835
836/// What the price keeps of one path.
837#[derive(Debug, Clone, Copy, PartialEq)]
838struct PathOutcome {
839    terminal_spot: f64,
840    non_positive: bool,
841}
842
843/// Runs one path without recording it.
844fn run_path(
845    params: &QuadraticRoughHestonParams,
846    grid: &QuadraticRoughHestonGrid,
847    step_sizes: (f64, f64),
848    seed: u64,
849) -> Result<PathOutcome, NonFiniteValue> {
850    let mut walk = Walk::new(params, step_sizes, seed);
851    let mut non_positive = false;
852    for to in 1..=grid.steps {
853        walk.advance(to)?;
854        non_positive |= walk.state.spot <= 0.0;
855    }
856    Ok(PathOutcome {
857        terminal_spot: walk.state.spot,
858        non_positive,
859    })
860}
861
862/// The running mean and squared deviations of Welford's algorithm.
863#[derive(Debug, Clone, Copy, Default, PartialEq)]
864struct Moments {
865    count: f64,
866    mean: f64,
867    squared_deviations: f64,
868}
869
870impl Moments {
871    #[inline]
872    fn add(&mut self, value: f64) {
873        self.count += 1.0;
874        let deviation = value - self.mean;
875        self.mean += deviation / self.count;
876        self.squared_deviations += deviation * (value - self.mean);
877    }
878
879    /// The mean and its standard error, as `Decimal`s. Needs two values or more.
880    fn estimate(&self) -> Result<(Decimal, Decimal), QuadraticRoughHestonError> {
881        let variance_of_mean = self.squared_deviations / (self.count - 1.0) / self.count;
882        Ok((
883            f64_to_decimal(self.mean)?,
884            f64_to_decimal(variance_of_mean.sqrt())?, // scan-banned: allow -- f64 `sqrt`: returns NaN for negative input, it does not abort; `f64_to_decimal` rejects it
885        ))
886    }
887}
888
889/// Prices a European call or put on `S_T` by Monte Carlo, undiscounted, as with `r = q = 0`.
890///
891/// One path is drawn per seed of a master generator seeded with `seed`. The seeds are drawn
892/// serially and the paths summed in their order, so the result does not depend on the number
893/// of threads of the rayon pool it runs on. A path that fails makes the price fail, and the
894/// failure reported is the one of the first such path. See [`QuadraticRoughHestonParams`] for
895/// the model, the scheme and the treatment of a spot that goes negative.
896///
897/// # Errors
898///
899/// - [`QuadraticRoughHestonError::InvalidParameter`] for a single path, which has no standard
900///   error.
901/// - [`QuadraticRoughHestonError::NonFinite`] when a path stops being finite, with its step.
902/// - [`QuadraticRoughHestonError::BufferTooLarge`] when `paths` cannot be allocated.
903/// - [`QuadraticRoughHestonError::Decimal`] when `strike` has no `f64` value or a result does
904///   not fit a `Decimal`.
905#[tracing::instrument(level = "debug", skip(params))]
906pub fn price_vanilla_quadratic_rough_heston(
907    params: &QuadraticRoughHestonParams,
908    grid: &QuadraticRoughHestonGrid,
909    style: OptionStyle,
910    strike: Positive,
911    paths: NonZeroUsize,
912    seed: u64,
913) -> Result<QuadraticRoughHestonEstimate, QuadraticRoughHestonError> {
914    let path_count = paths.get();
915    if path_count < MIN_PATHS {
916        return Err(invalid_parameter(
917            "paths",
918            Decimal::from(path_count),
919            "a standard error needs at least two paths",
920        ));
921    }
922    let raw_strike = decimal_to_f64(strike.to_dec())?;
923    let step_sizes = raw_step(grid)?;
924
925    let mut master = deterministic_rng(seed);
926    let mut seeds = reserve::<u64>(path_count, "path seeds")?;
927    seeds.extend((0..path_count).map(|_| master.random::<u64>()));
928
929    let mut outcomes = reserve(path_count, "path outcomes")?;
930    seeds
931        .par_iter()
932        .map(|path_seed| run_path(params, grid, step_sizes, *path_seed))
933        .collect_into_vec(&mut outcomes);
934
935    let non_positive_paths = outcomes
936        .iter()
937        .filter(|outcome| matches!(outcome, Ok(path) if path.non_positive))
938        .count();
939    let mut payoff = Moments::default();
940    let mut terminal = Moments::default();
941    for outcome in outcomes {
942        let outcome = outcome.map_err(NonFiniteValue::into_error)?;
943        payoff.add(match style {
944            OptionStyle::Call => (outcome.terminal_spot - raw_strike).max(0.0),
945            OptionStyle::Put => (raw_strike - outcome.terminal_spot).max(0.0),
946        });
947        terminal.add(outcome.terminal_spot);
948    }
949
950    let (price, standard_error) = payoff.estimate()?;
951    let (mean_terminal_spot, terminal_spot_standard_error) = terminal.estimate()?;
952    Ok(QuadraticRoughHestonEstimate {
953        price,
954        standard_error,
955        mean_terminal_spot,
956        terminal_spot_standard_error,
957        paths: path_count,
958        non_positive_paths,
959    })
960}
961
962#[cfg(test)]
963mod tests {
964    use super::*;
965    use optionstratlib_core::pos_or_panic;
966
967    /// Largest difference accepted against a value computed by hand with `f64`.
968    const HAND_TOLERANCE: f64 = 1e-12;
969
970    /// `n = 2`, `alpha = 3/4`, `x = 4`: `eta = (1/4, 1, 4)`. By hand, with
971    /// `G = Gamma(3/4) Gamma(1/4) = pi sqrt(2)` and `mu(dg) = g^(-3/4) dg / G`:
972    /// `c_i = 4 (eta_i^(1/4) - eta_(i-1)^(1/4)) / G`, and
973    /// `gamma_i = (4/5) (eta_i^(5/4) - eta_(i-1)^(5/4)) / (G c_i)`.
974    const SMALL_WEIGHTS: [f64; 2] = [0.26369654378952456, 0.37292322857805654];
975    const SMALL_DISCOUNTS: [f64; 2] = [0.5621320343559644, 2.2485281374238566];
976
977    fn small_kernel() -> QuadraticRoughHestonKernel {
978        QuadraticRoughHestonKernel::try_new(dec!(0.75), 2, dec!(4)).unwrap()
979    }
980
981    /// `a = 0.5`, `b = 0.1`, `c = 0.04`, `lambda = 2`, `eta = 0.5`, `S_0 = 100`, `Z^i_0 = (0.1, 0.2)`.
982    fn small_params() -> QuadraticRoughHestonParams {
983        QuadraticRoughHestonParams::try_new(
984            dec!(0.5),
985            dec!(0.1),
986            dec!(0.04),
987            dec!(2),
988            dec!(0.5),
989            small_kernel(),
990            vec![dec!(0.1), dec!(0.2)],
991            Positive::HUNDRED,
992        )
993        .unwrap()
994    }
995
996    fn close(actual: f64, expected: f64) -> bool {
997        (actual - expected).abs() < HAND_TOLERANCE
998    }
999
1000    #[test]
1001    fn test_kernel_small_case_matches_hand_values() {
1002        let kernel = small_kernel();
1003        assert_eq!(kernel.alpha(), dec!(0.75));
1004        assert_eq!(kernel.ratio(), dec!(4));
1005        assert_eq!(kernel.factors(), 2);
1006        for i in 0..2 {
1007            assert!(close(kernel.raw_weights[i], SMALL_WEIGHTS[i]), "c_{i}");
1008            assert!(close(kernel.raw_discounts[i], SMALL_DISCOUNTS[i]), "g_{i}");
1009            assert!(close(
1010                decimal_to_f64(kernel.weights()[i]).unwrap(),
1011                SMALL_WEIGHTS[i]
1012            ));
1013            assert!(close(
1014                decimal_to_f64(kernel.discounts()[i]).unwrap(),
1015                SMALL_DISCOUNTS[i]
1016            ));
1017        }
1018    }
1019
1020    /// `n = 3` has a half-integer `n / 2`: `eta = (4^(-3/2), 4^(-1/2), 4^(1/2), 4^(3/2))`.
1021    #[test]
1022    fn test_kernel_odd_factors_match_hand_values() {
1023        let kernel = QuadraticRoughHestonKernel::try_new(dec!(0.75), 3, dec!(4)).unwrap();
1024        let weights = [0.2217414783874054, 0.3135898060761293, 0.4434829567748108];
1025        let discounts = [0.2810660171779822, 1.1242640687119285, 4.497056274847715];
1026        for i in 0..3 {
1027            assert!(close(kernel.raw_weights[i], weights[i]), "c_{i}");
1028            assert!(close(kernel.raw_discounts[i], discounts[i]), "g_{i}");
1029        }
1030    }
1031
1032    /// The one figure of the calibration paper that reproduces: `gamma_10 = 542.32`, printed
1033    /// to two decimals. It is a check on the formula, and not a suggested ratio.
1034    #[test]
1035    fn test_kernel_discount_matches_the_printed_figure() {
1036        let kernel = QuadraticRoughHestonKernel::try_new(dec!(0.51), 10, dec!(3.92)).unwrap();
1037        let discount = kernel.discounts().last().copied().unwrap();
1038        assert!((discount - dec!(542.32)).abs() < dec!(0.005), "{discount}");
1039    }
1040
1041    /// `K^n(t)` from the hand values, at zero and away from it.
1042    #[test]
1043    fn test_kernel_value_at_is_the_sum_of_exponentials() {
1044        let kernel = small_kernel();
1045        let at_zero = kernel.value_at(Positive::ZERO).unwrap();
1046        assert!(close(
1047            decimal_to_f64(at_zero).unwrap(),
1048            SMALL_WEIGHTS[0] + SMALL_WEIGHTS[1]
1049        ));
1050        let t = 0.5_f64;
1051        let value = kernel.value_at(pos_or_panic!(0.5)).unwrap();
1052        let expected = SMALL_WEIGHTS[0] * (-SMALL_DISCOUNTS[0] * t).exp()
1053            + SMALL_WEIGHTS[1] * (-SMALL_DISCOUNTS[1] * t).exp();
1054        assert!(close(decimal_to_f64(value).unwrap(), expected));
1055        let far = kernel.value_at(Positive::MAX).unwrap();
1056        assert_eq!(far, Decimal::ZERO);
1057    }
1058
1059    #[test]
1060    fn test_kernel_invalid_alpha_returns_typed_error() {
1061        for alpha in [dec!(-1), dec!(0), dec!(0.5), dec!(1), dec!(2)] {
1062            assert!(
1063                matches!(
1064                    QuadraticRoughHestonKernel::try_new(alpha, 2, dec!(4)),
1065                    Err(QuadraticRoughHestonError::AlphaOutOfRange { alpha: rejected })
1066                        if rejected == alpha
1067                ),
1068                "{alpha}"
1069            );
1070        }
1071    }
1072
1073    #[test]
1074    fn test_kernel_alpha_next_to_its_bounds_is_accepted() {
1075        let tiny = Decimal::new(1, 20);
1076        for alpha in [dec!(0.5) + tiny, Decimal::ONE - tiny] {
1077            let kernel = QuadraticRoughHestonKernel::try_new(alpha, 2, dec!(4)).unwrap();
1078            assert!(kernel.raw_weights.iter().all(|weight| *weight > 0.0));
1079        }
1080    }
1081
1082    #[test]
1083    fn test_kernel_ratio_just_above_one_keeps_its_distance_from_one() {
1084        // An `f64` of the whole ratio rounds 1 + 1e-16 to one, and 1 + 3e-16 to 1 + 2.2e-16.
1085        let lowest = QuadraticRoughHestonKernel::try_new(dec!(0.75), 2, dec!(1.0000000000000001));
1086        assert!(lowest.is_ok());
1087
1088        // With `alpha = 0.75`, `n = 2`, `i = 1`, the exponent is zero and the first weight is
1089        // `(1 - x^(-1/4)) sin(pi / 4) / (pi / 4)`, about `(x - 1) sin(pi / 4) / pi`.
1090        let kernel =
1091            QuadraticRoughHestonKernel::try_new(dec!(0.75), 2, dec!(1.0000000000000003)).unwrap();
1092        let weight = decimal_to_f64(kernel.weights()[0]).unwrap();
1093        let expected = 3e-16 * (PI / 4.0).sin() / PI;
1094        assert!(
1095            (weight / expected - 1.0).abs() < 1e-6,
1096            "weight {weight}, expected {expected}"
1097        );
1098    }
1099
1100    #[test]
1101    fn test_kernel_no_factors_returns_typed_error() {
1102        assert!(matches!(
1103            QuadraticRoughHestonKernel::try_new(dec!(0.75), 0, dec!(4)),
1104            Err(QuadraticRoughHestonError::NoFactors)
1105        ));
1106    }
1107
1108    #[test]
1109    fn test_kernel_ratio_not_above_one_returns_typed_error() {
1110        for ratio in [dec!(-2), dec!(0), dec!(1)] {
1111            assert!(
1112                matches!(
1113                    QuadraticRoughHestonKernel::try_new(dec!(0.75), 2, ratio),
1114                    Err(QuadraticRoughHestonError::RatioNotAboveOne { ratio: rejected })
1115                        if rejected == ratio
1116                ),
1117                "{ratio}"
1118            );
1119        }
1120    }
1121
1122    /// The first factor of a very wide partition has no positive value in a `Decimal`, the
1123    /// last has none that fits, and a ratio that reads as one in `f64` has no partition.
1124    #[test]
1125    fn test_kernel_out_of_decimal_range_returns_typed_error() {
1126        assert!(matches!(
1127            QuadraticRoughHestonKernel::try_new(dec!(0.75), 1000, dec!(4)),
1128            Err(QuadraticRoughHestonError::KernelOutOfRange {
1129                quantity: "weight",
1130                index: 1
1131            })
1132        ));
1133        assert!(matches!(
1134            QuadraticRoughHestonKernel::try_new(dec!(0.75), 100, dec!(4)),
1135            Err(QuadraticRoughHestonError::KernelOutOfRange {
1136                quantity: "discount",
1137                index: 1
1138            })
1139        ));
1140        assert!(matches!(
1141            QuadraticRoughHestonKernel::try_new(dec!(0.75), 2, Decimal::ONE + Decimal::new(1, 28)),
1142            Err(QuadraticRoughHestonError::KernelOutOfRange { index: 1, .. })
1143        ));
1144    }
1145
1146    #[test]
1147    fn test_kernel_unallocatable_factors_return_typed_error() {
1148        assert!(matches!(
1149            QuadraticRoughHestonKernel::try_new(dec!(0.75), usize::MAX, dec!(4)),
1150            Err(QuadraticRoughHestonError::BufferTooLarge { .. })
1151        ));
1152    }
1153
1154    #[test]
1155    fn test_params_read_the_initial_state_and_keep_the_inputs() {
1156        let params = small_params();
1157        assert_eq!(params.a(), dec!(0.5));
1158        assert_eq!(params.b(), dec!(0.1));
1159        assert_eq!(params.c(), dec!(0.04));
1160        assert_eq!(params.lambda(), dec!(2));
1161        assert_eq!(params.eta(), dec!(0.5));
1162        assert_eq!(params.kernel(), &small_kernel());
1163        assert_eq!(params.initial_factors(), &[dec!(0.1), dec!(0.2)]);
1164        assert_eq!(params.spot(), Positive::HUNDRED);
1165        // By hand: `Z_0 = c_1 0.1 + c_2 0.2` and `V_0 = 0.5 (Z_0 - 0.1)^2 + 0.04`.
1166        assert!(close(
1167            decimal_to_f64(params.initial_z()).unwrap(),
1168            0.10095430009456377
1169        ));
1170        assert!(close(
1171            decimal_to_f64(params.initial_variance()).unwrap(),
1172            0.040000455344335245
1173        ));
1174    }
1175
1176    #[test]
1177    fn test_params_invalid_values_return_typed_errors() {
1178        let build = |a, c, lambda, eta, spot| {
1179            QuadraticRoughHestonParams::try_new(
1180                a,
1181                dec!(0.1),
1182                c,
1183                lambda,
1184                eta,
1185                small_kernel(),
1186                vec![dec!(0.1), dec!(0.2)],
1187                spot,
1188            )
1189        };
1190        let one = Decimal::ONE;
1191        let hundred = Positive::HUNDRED;
1192        let cases = [
1193            ("a", build(dec!(-0.1), one, one, one, hundred)),
1194            ("c", build(one, Decimal::ZERO, one, one, hundred)),
1195            ("c", build(one, dec!(-1), one, one, hundred)),
1196            ("lambda", build(one, one, Decimal::ZERO, one, hundred)),
1197            ("eta", build(one, one, one, Decimal::ZERO, hundred)),
1198            ("spot", build(one, one, one, one, Positive::ZERO)),
1199        ];
1200        for (expected, result) in cases {
1201            assert!(
1202                matches!(
1203                    result,
1204                    Err(QuadraticRoughHestonError::InvalidParameter { parameter, .. })
1205                        if parameter == expected
1206                ),
1207                "{expected}"
1208            );
1209        }
1210    }
1211
1212    #[test]
1213    fn test_params_wrong_factor_count_returns_typed_error() {
1214        let result = QuadraticRoughHestonParams::try_new(
1215            dec!(0.5),
1216            dec!(0.1),
1217            dec!(0.04),
1218            dec!(2),
1219            dec!(0.5),
1220            small_kernel(),
1221            vec![dec!(0.1)],
1222            Positive::HUNDRED,
1223        );
1224        assert!(matches!(
1225            result,
1226            Err(QuadraticRoughHestonError::InitialFactorCount {
1227                expected: 2,
1228                actual: 1
1229            })
1230        ));
1231    }
1232
1233    /// `V_0` is finite in `f64` and overflows a `Decimal`.
1234    #[test]
1235    fn test_params_initial_variance_out_of_decimal_range_returns_typed_error() {
1236        let result = QuadraticRoughHestonParams::try_new(
1237            dec!(1000000000000000000000),
1238            dec!(1000000000000000000000),
1239            dec!(0.04),
1240            dec!(2),
1241            dec!(0.5),
1242            small_kernel(),
1243            vec![dec!(0.1), dec!(0.2)],
1244            Positive::HUNDRED,
1245        );
1246        assert!(matches!(result, Err(QuadraticRoughHestonError::Decimal(_))));
1247    }
1248
1249    #[test]
1250    fn test_grid_keeps_the_step_and_rejects_a_non_positive_one() {
1251        let grid = QuadraticRoughHestonGrid::try_new(dec!(0.5), 4).unwrap();
1252        assert_eq!(grid.horizon(), dec!(0.5));
1253        assert_eq!(grid.steps(), 4);
1254        assert_eq!(grid.dt(), dec!(0.125));
1255        for (horizon, steps) in [
1256            (Decimal::ZERO, 4),
1257            (dec!(-1), 4),
1258            (dec!(1), 0),
1259            (Decimal::new(1, 28), 1000),
1260        ] {
1261            assert!(
1262                matches!(
1263                    QuadraticRoughHestonGrid::try_new(horizon, steps),
1264                    Err(QuadraticRoughHestonError::NonPositiveStep { .. })
1265                ),
1266                "{horizon} over {steps}"
1267            );
1268        }
1269    }
1270
1271    /// One step against `dW = 0.01`, `dt = 0.25`, computed by hand from `SMALL_WEIGHTS`,
1272    /// `SMALL_DISCOUNTS` and `small_params`:
1273    /// `S_1 = 100 + 100 sqrt(V_0) dW`, and
1274    /// `Z^i_1 = (Z^i_0 - 2 Z_0 0.25 + 0.5 sqrt(V_0) dW) / (1 + gamma_i 0.25)`.
1275    #[test]
1276    fn test_step_matches_a_hand_computation_with_the_implicit_denominator() {
1277        let params = small_params();
1278        let state = params.raw_initial.clone();
1279        let mut next = state.clone();
1280        step(&params, &state, 0.01, 0.25, &mut next);
1281        assert!(close(next.spot, 100.20000113835759));
1282        assert!(close(next.factors[0], 0.044297583028316205));
1283        assert!(close(next.factors[1], 0.09635731956970187));
1284        assert!(close(next.z, 0.04761500225385726));
1285        assert!(close(next.variance, 0.04137209399443169));
1286    }
1287
1288    /// A zero increment leaves the spot, and each factor keeps only its drift.
1289    #[test]
1290    fn test_step_zero_increment_keeps_the_spot() {
1291        let params = small_params();
1292        let state = params.raw_initial.clone();
1293        let mut next = state.clone();
1294        step(&params, &state, 0.0, 0.25, &mut next);
1295        assert_eq!(next.spot, state.spot);
1296        let expected = (0.1 - 2.0 * state.z * 0.25) / (1.0 + SMALL_DISCOUNTS[0] * 0.25);
1297        assert!(close(next.factors[0], expected));
1298    }
1299
1300    /// Halving the step on one Brownian path, without the noise of two independent estimates:
1301    /// each coarse `dW` is the sum of two fine ones.
1302    ///
1303    /// `alpha = 0.6`, `n = 5`, `x = 5`, `a = 0.5`, `b = 0.1`, `c = 0.01`, `lambda = 1.5`,
1304    /// `eta = 1`, `Z^i_0 = -0.05`, `S_0 = K = 100`, `T = 0.25`, an ATM call, `N = 100` against
1305    /// `2N = 200`, on 20,000 paths. The change measures `0.0087`, with a standard error of
1306    /// `0.0013`, against a tolerance of `0.02`.
1307    ///
1308    /// The coarse side goes through [`Walk`], and the fine side drives [`step`] with its own
1309    /// `dt`. If both went through `Walk`, a wrong `dt` there would pass: it is a consistent
1310    /// scheme for another model, and halving its step changes the price just as little.
1311    #[test]
1312    fn test_walk_halving_the_step_on_coupled_increments_stays_within_its_tolerance() {
1313        const COARSE_STEPS: usize = 100;
1314        const PATHS: usize = 20_000;
1315        const TOLERANCE: f64 = 0.02;
1316        let kernel = QuadraticRoughHestonKernel::try_new(dec!(0.6), 5, dec!(5)).unwrap();
1317        let params = QuadraticRoughHestonParams::try_new(
1318            dec!(0.5),
1319            dec!(0.1),
1320            dec!(0.01),
1321            dec!(1.5),
1322            dec!(1),
1323            kernel,
1324            vec![dec!(-0.05); 5],
1325            Positive::HUNDRED,
1326        )
1327        .unwrap();
1328        let coarse_grid = QuadraticRoughHestonGrid::try_new(dec!(0.25), COARSE_STEPS).unwrap();
1329        let coarse_step = raw_step(&coarse_grid).unwrap();
1330        let fine_dt = 0.25 / (2.0 * COARSE_STEPS as f64);
1331        let fine_sqrt_dt = fine_dt.sqrt();
1332
1333        let mut rng = deterministic_rng(958);
1334        let mut change = Moments::default();
1335        for _ in 0..PATHS {
1336            let mut coarse = Walk::new(&params, coarse_step, 0);
1337            let mut fine = params.raw_initial.clone();
1338            let mut scratch = fine.clone();
1339            for to in 1..=COARSE_STEPS {
1340                let first: f64 = StandardNormal.sample(&mut rng);
1341                let second: f64 = StandardNormal.sample(&mut rng);
1342                let (first, second) = (fine_sqrt_dt * first, fine_sqrt_dt * second);
1343                step(&params, &fine, first, fine_dt, &mut scratch);
1344                step(&params, &scratch, second, fine_dt, &mut fine);
1345                coarse.apply(first + second, to).unwrap();
1346            }
1347            let payoff = |spot: f64| (spot - 100.0).max(0.0);
1348            change.add(payoff(coarse.state.spot) - payoff(fine.spot));
1349        }
1350        assert!(change.mean.abs() < TOLERANCE, "{change:?}");
1351    }
1352
1353    #[test]
1354    fn test_state_not_finite_names_its_quantity_and_step() {
1355        let finite = small_params().raw_initial;
1356        assert_eq!(finite.ensure_finite(3), Ok(()));
1357        for (quantity, state) in [
1358            (
1359                "spot",
1360                State {
1361                    spot: f64::NAN,
1362                    ..finite.clone()
1363                },
1364            ),
1365            (
1366                "z",
1367                State {
1368                    z: f64::INFINITY,
1369                    ..finite.clone()
1370                },
1371            ),
1372            (
1373                "variance",
1374                State {
1375                    variance: f64::INFINITY,
1376                    ..finite.clone()
1377                },
1378            ),
1379        ] {
1380            assert_eq!(
1381                state.ensure_finite(3),
1382                Err(NonFiniteValue { quantity, step: 3 })
1383            );
1384        }
1385    }
1386
1387    #[test]
1388    fn test_moments_match_the_sample_mean_and_standard_error() {
1389        let mut moments = Moments::default();
1390        for value in [1.0, 2.0, 3.0, 4.0] {
1391            moments.add(value);
1392        }
1393        let (mean, standard_error) = moments.estimate().unwrap();
1394        assert_eq!(mean, dec!(2.5));
1395        // Sample variance 5/3, over 4 values.
1396        let expected = (5.0_f64 / 3.0 / 4.0).sqrt();
1397        assert!(close(decimal_to_f64(standard_error).unwrap(), expected));
1398    }
1399}
optionstratlib_simulation/simulation/quadratic_rough_heston.rs

optionstratlib_simulation/simulation/
quadratic_rough_heston.rs
