rustyqlib/equity/
rainbow.rs

1//! Rainbow (multi-asset) options: best-of, worst-of, spread, basket and
2//! exchange payoffs on n correlated lognormal assets.
3//!
4//! Engines:
5//! - **Analytic**: Margrabe (exchange, exact), Kirk's approximation
6//!   (spread), moment-matched lognormal (basket). Best-of / worst-of have
7//!   no analytic pricer yet (Stulz for n = 2 is future work) and price on
8//!   Monte Carlo.
9//! - **Monte Carlo**: correlated terminal GBM (Cholesky), low-discrepancy
10//!   or antithetic pseudo-random sampling, deterministic parallel
11//!   reduction, standard errors.
12//!
13//! Greeks: per-asset `deltas` and `vegas` by common-random-number bumps;
14//! scalar theta and rho. Each asset carries a flat vol; per-asset smiles
15//! for multi-asset payoffs are future work.
16
17use chrono::NaiveDate;
18use libm::exp;
19use rayon::prelude::*;
20use serde::{Deserialize, Serialize};
21
22use crate::core::curves::{Compounding, YieldCurve};
23use crate::core::daycount::DayCountConvention;
24use crate::core::linalg::{cholesky, nearest_correlation};
25use crate::core::trade::PutOrCall;
26use crate::core::utils::norm_cdf;
27use crate::equity::montecarlo::{McStats, MonteCarloConfig, Sampler};
28use crate::equity::utils::PricingEngine;
29use crate::core::montecarlo::{path_normals, QmcSequence};
30use crate::core::errors::RustyQLibError;
31use crate::core::results::{Greeks, PricingResult};
32use crate::core::traits::Instrument;
33
34const PATH_CHUNK: usize = 4096;
35
36// ── Contract data (JSON) ────────────────────────────────────────────────
37
38#[derive(Clone, Debug, Deserialize, Serialize)]
39pub struct RainbowAssetData {
40    pub symbol: String,
41    pub spot: f64,
42    pub volatility: f64,
43    pub dividend: Option<f64>,
44}
45
46#[derive(Clone, Debug, Deserialize, Serialize)]
47pub struct RainbowOptionData {
48    pub symbol: String,
49    /// "best_of" | "worst_of" | "spread" | "basket" | "exchange"
50    pub rainbow_type: String,
51    pub put_or_call: Option<String>,
52    pub assets: Vec<RainbowAssetData>,
53    /// Full correlation matrix, n x n.
54    pub correlations: Vec<Vec<f64>>,
55    pub strike_price: Option<f64>,
56    /// Basket weights (defaults to equal weights).
57    pub weights: Option<Vec<f64>>,
58    pub maturity: String,
59    pub risk_free_rate: Option<f64>,
60    pub discount_curve: Option<crate::core::curves::CurveInput>,
61    pub pricer: Option<String>,
62    pub simulation: Option<u64>,
63    pub mc_sampler: Option<String>,
64    pub mc_seed: Option<u64>,
65    /// Pricing as-of date (`YYYY-MM-DD`); defaults to today.
66    pub valuation_date: Option<String>,
67}
68
69#[derive(Debug, Clone, Copy, PartialEq, Eq)]
70pub enum RainbowType {
71    BestOf,
72    WorstOf,
73    Spread,
74    Basket,
75    Exchange,
76}
77
78// ── Instrument ──────────────────────────────────────────────────────────
79
80#[derive(Debug)]
81pub struct RainbowOption {
82    pub symbol: String,
83    pub rainbow_type: RainbowType,
84    pub put_or_call: PutOrCall,
85    pub spots: Vec<f64>,
86    pub vols: Vec<f64>,
87    pub dividends: Vec<f64>,
88    pub correlations: Vec<Vec<f64>>,
89    pub strike_price: f64,
90    pub weights: Vec<f64>,
91    pub maturity_date: NaiveDate,
92    pub valuation_date: NaiveDate,
93    pub discount_curve: YieldCurve,
94    /// The numerical method with its settings. Rainbow payoffs price on
95    /// the analytic engine (Margrabe / Kirk / moment matching) or Monte
96    /// Carlo (terminal correlated GBM: `paths`, `sampler` and `seed` from
97    /// the config; `time_steps`/`scheme` do not apply).
98    pub engine: PricingEngine,
99    /// Cholesky factor of the correlation matrix (lower triangular).
100    chol: Vec<Vec<f64>>,
101}
102
103impl Instrument for RainbowOption {
104    fn try_npv(&self) -> Result<f64, RustyQLibError> {
105        self.check_engine_support()?;
106        Ok(match self.engine {
107            PricingEngine::BlackScholes => self.analytic_npv_with(&self.params()),
108            _ => self.mc_stats_with(&self.params()).pv,
109        })
110    }
111
112    /// Value, scalar theta/rho and (under Monte Carlo) the standard
113    /// error. Spot Greeks are per-asset for rainbows — see
114    /// [`RainbowOption::deltas`] and [`RainbowOption::vegas`] — so the
115    /// scalar delta/gamma/vega slots stay zero.
116    fn price(&self) -> Result<PricingResult, RustyQLibError> {
117        self.check_engine_support()?;
118        let (pv, std_err) = match self.engine {
119            PricingEngine::MonteCarlo(_) => {
120                let stats = self.mc_stats_with(&self.params());
121                (stats.pv, Some(stats.std_err))
122            }
123            _ => (self.try_npv()?, None),
124        };
125        Ok(PricingResult {
126            pv,
127            greeks: Greeks { theta: self.theta(), rho: self.rho(), ..Default::default() },
128            std_err,
129        })
130    }
131}
132
133/// Market snapshot bumped by the Greeks (common random numbers).
134#[derive(Clone)]
135struct Params {
136    spots: Vec<f64>,
137    vols: Vec<f64>,
138    r: f64,
139    t: f64,
140}
141
142impl RainbowOption {
143    /// Monte Carlo settings. Invariant: only called on the Monte Carlo
144    /// code paths (the engine dispatch guarantees it).
145    fn mc_cfg(&self) -> &MonteCarloConfig {
146        match &self.engine {
147            PricingEngine::MonteCarlo(cfg) => cfg,
148            _ => unreachable!("Monte Carlo code path reached on a non-MC engine"),
149        }
150    }
151
152    /// Build from contract data, panicking on any invalid field. Fallible
153    /// callers should use [`RainbowOption::try_from_json`].
154    pub fn from_json(data: &RainbowOptionData) -> Box<RainbowOption> {
155        Self::try_from_json(data).unwrap_or_else(|e| panic!("{e}"))
156    }
157
158    pub fn try_from_json(data: &RainbowOptionData) -> Result<Box<RainbowOption>, RustyQLibError> {
159        let valuation_date =
160            crate::core::data_models::parse_valuation_date(data.valuation_date.as_deref())?;
161        let n = data.assets.len();
162        if n < 2 {
163            return Err(RustyQLibError::invalid_input("assets", "rainbow options need at least two assets"));
164        }
165        let rainbow_type = match data.rainbow_type.trim().to_lowercase().as_str() {
166            "best_of" | "bestof" | "max" => RainbowType::BestOf,
167            "worst_of" | "worstof" | "min" => RainbowType::WorstOf,
168            "spread" => RainbowType::Spread,
169            "basket" => RainbowType::Basket,
170            "exchange" | "margrabe" => RainbowType::Exchange,
171            other => return Err(RustyQLibError::invalid_input(
172                "rainbow_type",
173                format!("invalid rainbow_type '{other}'"),
174            )),
175        };
176        if matches!(rainbow_type, RainbowType::Spread | RainbowType::Exchange) && n != 2 {
177            return Err(RustyQLibError::invalid_input(
178                "assets",
179                "spread and exchange options take exactly two assets",
180            ));
181        }
182        let put_or_call = match data.put_or_call.as_deref().unwrap_or("C").trim() {
183            "C" | "c" | "Call" | "call" => PutOrCall::Call,
184            "P" | "p" | "Put" | "put" => PutOrCall::Put,
185            other => return Err(RustyQLibError::invalid_input(
186                "put_or_call",
187                format!("invalid put_or_call '{other}' (use 'C' or 'P')"),
188            )),
189        };
190        let strike_price = data.strike_price.unwrap_or(0.0);
191        if rainbow_type != RainbowType::Exchange && data.strike_price.is_none() {
192            return Err(RustyQLibError::invalid_input("strike_price", "strike_price is required"));
193        }
194        let weights = match &data.weights {
195            Some(w) => {
196                if w.len() != n {
197                    return Err(RustyQLibError::invalid_input(
198                        "weights",
199                        "weights must match the number of assets",
200                    ));
201                }
202                w.clone()
203            }
204            None => vec![1.0 / n as f64; n],
205        };
206        if data.correlations.len() != n || data.correlations.iter().any(|row| row.len() != n) {
207            return Err(RustyQLibError::invalid_input(
208                "correlations",
209                "correlations must be an n x n matrix",
210            ));
211        }
212        // an empirical / hand-stressed matrix that fails PSD is repaired
213        // with Higham's nearest-correlation projection; asymmetry or a
214        // non-unit diagonal is a data error and still rejected
215        let chol = match cholesky(&data.correlations) {
216            Ok(l) => l,
217            Err(RustyQLibError::NumericalError(ref msg))
218                if msg.contains("positive semi-definite") =>
219            {
220                log::warn!(
221                    "correlation matrix is not PSD; \
222                     projecting to the nearest correlation matrix (Higham)"
223                );
224                let repaired = nearest_correlation(&data.correlations, 1e-12, 200)?;
225                cholesky(&repaired)?
226            }
227            Err(e) => return Err(e),
228        };
229        let discount_curve = match &data.discount_curve {
230            Some(input) => YieldCurve::from_input(input, valuation_date)?,
231            None => YieldCurve::flat(
232                data.risk_free_rate.unwrap_or(0.0),
233                valuation_date,
234                DayCountConvention::Act365,
235                Compounding::Continuous,
236            )?,
237        };
238        let maturity_date = NaiveDate::parse_from_str(&data.maturity, "%Y-%m-%d")
239            .map_err(|_| RustyQLibError::invalid_input(
240                "maturity",
241                format!("invalid date '{}' (expected YYYY-MM-DD)", data.maturity),
242            ))?;
243        Ok(Box::new(RainbowOption {
244            symbol: data.symbol.clone(),
245            rainbow_type,
246            put_or_call,
247            spots: data.assets.iter().map(|a| a.spot).collect(),
248            vols: data.assets.iter().map(|a| a.volatility).collect(),
249            dividends: data.assets.iter().map(|a| a.dividend.unwrap_or(0.0)).collect(),
250            correlations: data.correlations.clone(),
251            strike_price,
252            weights,
253            maturity_date,
254            valuation_date,
255            discount_curve,
256            engine: match data.pricer.as_deref().map_or("MC", |v| v).trim() {
257                "Analytical" | "analytical" => PricingEngine::BlackScholes,
258                "MonteCarlo" | "montecarlo" | "MC" | "mc" => {
259                    PricingEngine::MonteCarlo(MonteCarloConfig {
260                        paths: data.simulation.unwrap_or(100_000) as usize,
261                        sampler: data
262                            .mc_sampler
263                            .as_deref()
264                            .map(|s| {
265                                s.parse::<Sampler>().map_err(|_| {
266                                    RustyQLibError::invalid_input(
267                                        "mc_sampler",
268                                        format!("invalid mc_sampler '{s}'"),
269                                    )
270                                })
271                            })
272                            .transpose()?
273                            .unwrap_or(Sampler::Sobol),
274                        seed: data.mc_seed.unwrap_or(42),
275                        ..Default::default()
276                    })
277                }
278                other => return Err(RustyQLibError::invalid_input(
279                    "pricer",
280                    format!("invalid pricer '{other}' for rainbow (Analytical or MC)"),
281                )),
282            },
283            chol,
284        }))
285    }
286
287    pub fn time_to_maturity(&self) -> f64 {
288        (self.maturity_date - self.valuation_date).num_days() as f64 / 365.0
289    }
290
291    fn params(&self) -> Params {
292        let t = self.time_to_maturity();
293        Params {
294            spots: self.spots.clone(),
295            vols: self.vols.clone(),
296            r: self.discount_curve.zero_rate_with(t, Compounding::Continuous),
297            t,
298        }
299    }
300
301    /// Terminal payoff on realized asset levels.
302    fn payoff(&self, terminal: &[f64]) -> f64 {
303        let phi = match self.put_or_call {
304            PutOrCall::Call => 1.0,
305            PutOrCall::Put => -1.0,
306        };
307        let k = self.strike_price;
308        match self.rainbow_type {
309            RainbowType::BestOf => {
310                let best = terminal.iter().cloned().fold(f64::MIN, f64::max);
311                (phi * (best - k)).max(0.0)
312            }
313            RainbowType::WorstOf => {
314                let worst = terminal.iter().cloned().fold(f64::MAX, f64::min);
315                (phi * (worst - k)).max(0.0)
316            }
317            RainbowType::Spread => (phi * (terminal[0] - terminal[1] - k)).max(0.0),
318            RainbowType::Basket => {
319                let basket: f64 =
320                    self.weights.iter().zip(terminal).map(|(w, s)| w * s).sum();
321                (phi * (basket - k)).max(0.0)
322            }
323            RainbowType::Exchange => match self.put_or_call {
324                PutOrCall::Call => (terminal[0] - terminal[1]).max(0.0),
325                PutOrCall::Put => (terminal[1] - terminal[0]).max(0.0),
326            },
327        }
328    }
329
330    // ── Pricing ─────────────────────────────────────────────────────────
331
332    /// Reject engine/payoff combinations the library cannot price,
333    /// with an error naming the engine that can.
334    pub(crate) fn check_engine_support(&self) -> Result<(), RustyQLibError> {
335        match &self.engine {
336            PricingEngine::BlackScholes => {
337                if matches!(self.rainbow_type, RainbowType::BestOf | RainbowType::WorstOf) {
338                    return Err(RustyQLibError::UnsupportedEngine(
339                        "best-of / worst-of rainbows have no analytic pricer yet \
340                         (Stulz for two assets is future work); use MonteCarlo"
341                            .to_string(),
342                    ));
343                }
344                Ok(())
345            }
346            PricingEngine::MonteCarlo(_) => Ok(()),
347            other => Err(RustyQLibError::UnsupportedEngine(format!(
348                "rainbow options price on the Analytical or MonteCarlo engines, not {:?}",
349                other.kind()
350            ))),
351        }
352    }
353
354    pub fn npv_with_stats(&self) -> Option<McStats> {
355        match self.engine {
356            PricingEngine::MonteCarlo(_) => Some(self.mc_stats_with(&self.params())),
357            _ => None,
358        }
359    }
360
361    fn price_with(&self, p: &Params) -> f64 {
362        match self.engine {
363            PricingEngine::BlackScholes => self.analytic_npv_with(p),
364            _ => self.mc_stats_with(p).pv,
365        }
366    }
367
368    /// Per-asset spot deltas (central bumps, common random numbers).
369    pub fn deltas(&self) -> Vec<f64> {
370        let base = self.params();
371        (0..self.spots.len())
372            .map(|i| {
373                let h = base.spots[i] * 0.01;
374                let mut up = base.clone();
375                up.spots[i] += h;
376                let mut dn = base.clone();
377                dn.spots[i] -= h;
378                (self.price_with(&up) - self.price_with(&dn)) / (2.0 * h)
379            })
380            .collect()
381    }
382
383    /// Per-asset vegas (central bumps of each asset's vol).
384    pub fn vegas(&self) -> Vec<f64> {
385        let base = self.params();
386        (0..self.vols.len())
387            .map(|i| {
388                let h = 0.01;
389                let mut up = base.clone();
390                up.vols[i] += h;
391                let mut dn = base.clone();
392                dn.vols[i] = (dn.vols[i] - h).max(1e-6);
393                (self.price_with(&up) - self.price_with(&dn)) / (2.0 * h)
394            })
395            .collect()
396    }
397
398    pub fn theta(&self) -> f64 {
399        let base = self.params();
400        let h = (1.0 / 365.0_f64).min(0.5 * base.t);
401        let mut up = base.clone();
402        up.t += h;
403        let mut dn = base.clone();
404        dn.t -= h;
405        -(self.price_with(&up) - self.price_with(&dn)) / (2.0 * h)
406    }
407
408    pub fn rho(&self) -> f64 {
409        let base = self.params();
410        let h = 1e-4;
411        let mut up = base.clone();
412        up.r += h;
413        let mut dn = base.clone();
414        dn.r -= h;
415        (self.price_with(&up) - self.price_with(&dn)) / (2.0 * h)
416    }
417
418    // ── Analytic pricers ────────────────────────────────────────────────
419
420    fn analytic_npv_with(&self, p: &Params) -> f64 {
421        match self.rainbow_type {
422            RainbowType::Exchange => self.margrabe(p),
423            RainbowType::Spread => self.kirk(p),
424            RainbowType::Basket => self.basket_moment_match(p),
425            // invariant: check_engine_support refuses these before pricing
426            RainbowType::BestOf | RainbowType::WorstOf => unreachable!(
427                "best-of/worst-of on the analytic engine is rejected before pricing"
428            ),
429        }
430    }
431
432    /// Margrabe (1978), exact: exchange option pays (S1 - S2)^+.
433    fn margrabe(&self, p: &Params) -> f64 {
434        let (i, j) = match self.put_or_call {
435            PutOrCall::Call => (0, 1),
436            PutOrCall::Put => (1, 0),
437        };
438        let rho = self.correlations[0][1];
439        let sigma = (p.vols[i] * p.vols[i] + p.vols[j] * p.vols[j]
440            - 2.0 * rho * p.vols[i] * p.vols[j])
441            .sqrt();
442        let (q_i, q_j) = (self.dividends[i], self.dividends[j]);
443        let st = sigma * p.t.sqrt();
444        if st < 1e-12 {
445            // perfectly correlated identical dynamics: the exchange is
446            // deterministic — discounted positive forward difference
447            return (p.spots[i] * exp(-q_i * p.t) - p.spots[j] * exp(-q_j * p.t)).max(0.0);
448        }
449        let d1 = ((p.spots[i] / p.spots[j]).ln() + (q_j - q_i + 0.5 * sigma * sigma) * p.t) / st;
450        let d2 = d1 - st;
451        p.spots[i] * exp(-q_i * p.t) * norm_cdf(d1) - p.spots[j] * exp(-q_j * p.t) * norm_cdf(d2)
452    }
453
454    /// Kirk's (1995) approximation for spread options (S1 - S2 - K)^+.
455    fn kirk(&self, p: &Params) -> f64 {
456        let f1 = p.spots[0] * exp((p.r - self.dividends[0]) * p.t);
457        let f2 = p.spots[1] * exp((p.r - self.dividends[1]) * p.t);
458        let k = self.strike_price;
459        let rho = self.correlations[0][1];
460        let w = f2 / (f2 + k);
461        let sigma = (p.vols[0] * p.vols[0] - 2.0 * rho * p.vols[0] * p.vols[1] * w
462            + p.vols[1] * p.vols[1] * w * w)
463            .sqrt();
464        let st = sigma * p.t.sqrt();
465        let d1 = ((f1 / (f2 + k)).ln() + 0.5 * sigma * sigma * p.t) / st;
466        let d2 = d1 - st;
467        let df = exp(-p.r * p.t);
468        match self.put_or_call {
469            PutOrCall::Call => df * (f1 * norm_cdf(d1) - (f2 + k) * norm_cdf(d2)),
470            PutOrCall::Put => df * ((f2 + k) * norm_cdf(-d2) - f1 * norm_cdf(-d1)),
471        }
472    }
473
474    /// Lognormal moment matching for basket options (Levy / Turnbull-Wakeman
475    /// style): match the basket forward's first two moments, price with
476    /// Black's formula.
477    fn basket_moment_match(&self, p: &Params) -> f64 {
478        let n = p.spots.len();
479        let fwds: Vec<f64> = (0..n)
480            .map(|i| self.weights[i] * p.spots[i] * exp((p.r - self.dividends[i]) * p.t))
481            .collect();
482        let m1: f64 = fwds.iter().sum();
483        let mut m2 = 0.0;
484        for i in 0..n {
485            for j in 0..n {
486                m2 += fwds[i]
487                    * fwds[j]
488                    * exp(self.correlations[i][j] * p.vols[i] * p.vols[j] * p.t);
489            }
490        }
491        let log_var = (m2 / (m1 * m1)).ln().max(1e-12);
492        let sqrt_v = log_var.sqrt();
493        let k = self.strike_price;
494        let d1 = ((m1 / k).ln() + 0.5 * log_var) / sqrt_v;
495        let d2 = d1 - sqrt_v;
496        let df = exp(-p.r * p.t);
497        match self.put_or_call {
498            PutOrCall::Call => df * (m1 * norm_cdf(d1) - k * norm_cdf(d2)),
499            PutOrCall::Put => df * (k * norm_cdf(-d2) - m1 * norm_cdf(-d1)),
500        }
501    }
502
503    // ── Monte Carlo (correlated terminal GBM) ───────────────────────────
504
505    fn mc_stats_with(&self, p: &Params) -> McStats {
506        let n = self.spots.len();
507        let t = p.t;
508        let df = exp(-p.r * t);
509        let sqrt_t = t.sqrt();
510        let drifts: Vec<f64> = (0..n)
511            .map(|i| (p.r - self.dividends[i] - 0.5 * p.vols[i] * p.vols[i]) * t)
512            .collect();
513        let cfg = self.mc_cfg();
514        let qmc = match cfg.sampler {
515            Sampler::Sobol => Some(QmcSequence::new(n, cfg.seed)),
516            Sampler::PseudoRandom => None,
517        };
518        let chunks = cfg.paths.div_ceil(PATH_CHUNK);
519        let partials: Vec<(f64, f64)> = (0..chunks)
520            .into_par_iter()
521            .map(|chunk| {
522                let mut eps = vec![0.0; n];
523                let mut terminal = vec![0.0; n];
524                let (mut sum, mut sum_sq) = (0.0, 0.0);
525                for path in chunk * PATH_CHUNK..((chunk + 1) * PATH_CHUNK).min(cfg.paths) {
526                    match &qmc {
527                        Some(seq) => seq.normals(path as u64 + 1, &mut eps),
528                        None => {
529                            // antithetic pairs from per-pair streams
530                            path_normals(cfg.seed, (path / 2) as u64, &mut eps);
531                            if path % 2 == 1 {
532                                for e in eps.iter_mut() {
533                                    *e = -*e;
534                                }
535                            }
536                        }
537                    }
538                    for i in 0..n {
539                        // z_i = sum_j L[i][j] eps_j (Cholesky-correlated)
540                        let z: f64 =
541                            (0..=i).map(|j| self.chol[i][j] * eps[j]).sum();
542                        terminal[i] = p.spots[i] * exp(drifts[i] + p.vols[i] * sqrt_t * z);
543                    }
544                    let v = df * self.payoff(&terminal);
545                    sum += v;
546                    sum_sq += v * v;
547                }
548                (sum, sum_sq)
549            })
550            .collect();
551        let (sum, sum_sq) =
552            partials.into_iter().fold((0.0, 0.0), |a, b| (a.0 + b.0, a.1 + b.1));
553        let nf = cfg.paths as f64;
554        let mean = sum / nf;
555        let var = (sum_sq / nf - mean * mean).max(0.0);
556        McStats { pv: mean, std_err: (var / nf).sqrt(), paths: cfg.paths, steps: 1 }
557    }
558}
559
560#[cfg(test)]
561mod tests {
562    use super::*;
563    use chrono::Local;
564    use crate::equity::blackscholes::bs_price;
565    use crate::equity::utils::Engine;
566
567    fn two_asset(rainbow_type: &str, pc: &str, strike: Option<f64>, rho: f64) -> Box<RainbowOption> {
568        RainbowOption::from_json(&RainbowOptionData {
569            symbol: "RB".to_string(),
570            rainbow_type: rainbow_type.to_string(),
571            put_or_call: Some(pc.to_string()),
572            assets: vec![
573                RainbowAssetData {
574                    symbol: "A".into(),
575                    spot: 100.0,
576                    volatility: 0.3,
577                    dividend: Some(0.02),
578                },
579                RainbowAssetData {
580                    symbol: "B".into(),
581                    spot: 95.0,
582                    volatility: 0.25,
583                    dividend: Some(0.01),
584                },
585            ],
586            correlations: vec![vec![1.0, rho], vec![rho, 1.0]],
587            strike_price: strike,
588            weights: None,
589            maturity: maturity_1y(),
590            risk_free_rate: Some(0.05),
591            discount_curve: None,
592            pricer: Some("MC".to_string()),
593            simulation: Some(100_000),
594            mc_sampler: None,
595            mc_seed: None,
596            valuation_date: None,
597        })
598    }
599
600    fn maturity_1y() -> String {
601        let d = Local::now().date_naive() + chrono::Duration::days(365);
602        d.format("%Y-%m-%d").to_string()
603    }
604
605    #[test]
606    fn unsupported_engines_error_instead_of_panicking() {
607        // best-of has no analytic pricer: typed refusal, not a panic
608        let mut option = two_asset("best_of", "C", Some(100.0), 0.6);
609        option.engine = PricingEngine::BlackScholes;
610        match option.try_npv() {
611            Err(RustyQLibError::UnsupportedEngine(msg)) => {
612                assert!(msg.contains("MonteCarlo"), "should name the right engine: {msg}")
613            }
614            other => panic!("expected UnsupportedEngine, got {other:?}"),
615        }
616        // engines that never apply to rainbows are refused too
617        option.engine = PricingEngine::from_kind(Engine::Binomial);
618        assert!(matches!(
619            option.try_npv(),
620            Err(RustyQLibError::UnsupportedEngine(_))
621        ));
622        // and price() carries the same guarantee
623        assert!(option.price().is_err());
624    }
625
626    #[test]
627    fn price_reports_theta_rho_and_mc_std_err() {
628        let option = two_asset("exchange", "C", None, 0.6);
629        let result = option.price().unwrap();
630        let se = result.std_err.expect("MC rainbow must report a standard error");
631        assert!(se > 0.0 && se.is_finite());
632        assert!((result.pv - option.npv()).abs() < 1e-12);
633        assert_eq!(result.greeks.theta, option.theta());
634        assert_eq!(result.greeks.rho, option.rho());
635        assert_eq!(result.greeks.delta, 0.0, "spot Greeks are per-asset");
636    }
637
638    #[test]
639    fn margrabe_matches_monte_carlo() {
640        let mut option = two_asset("exchange", "C", None, 0.6);
641        option.engine = PricingEngine::BlackScholes;
642        let analytic = option.npv();
643        option.engine = PricingEngine::from_kind(Engine::MonteCarlo);
644        let mc = option.npv();
645        assert!((mc - analytic).abs() < 0.05, "mc={mc} margrabe={analytic}");
646        assert!(analytic > 0.0);
647    }
648
649    #[test]
650    fn margrabe_vanishes_for_identical_assets() {
651        let mut option = two_asset("exchange", "C", None, 1.0);
652        option.spots = vec![100.0, 100.0];
653        option.vols = vec![0.3, 0.3];
654        option.dividends = vec![0.02, 0.02];
655        option.engine = PricingEngine::BlackScholes;
656        assert!(option.npv().abs() < 1e-10);
657    }
658
659    #[test]
660    fn kirk_close_to_monte_carlo() {
661        let mut option = two_asset("spread", "C", Some(5.0), 0.6);
662        option.engine = PricingEngine::BlackScholes;
663        let kirk = option.npv();
664        option.engine = PricingEngine::from_kind(Engine::MonteCarlo);
665        let mc = option.npv();
666        // Kirk is an approximation: agreement at the few-cents level
667        assert!((mc - kirk).abs() < 0.10, "mc={mc} kirk={kirk}");
668    }
669
670    #[test]
671    fn spread_with_zero_strike_equals_margrabe() {
672        let mut spread = two_asset("spread", "C", Some(0.0), 0.6);
673        spread.engine = PricingEngine::BlackScholes;
674        let mut exchange = two_asset("exchange", "C", None, 0.6);
675        exchange.engine = PricingEngine::BlackScholes;
676        assert!((spread.npv() - exchange.npv()).abs() < 1e-10);
677    }
678
679    #[test]
680    fn basket_moment_match_close_to_monte_carlo() {
681        let data = RainbowOptionData {
682            symbol: "BK".into(),
683            rainbow_type: "basket".into(),
684            put_or_call: Some("C".into()),
685            assets: vec![
686                RainbowAssetData { symbol: "A".into(), spot: 100.0, volatility: 0.3, dividend: None },
687                RainbowAssetData { symbol: "B".into(), spot: 90.0, volatility: 0.25, dividend: None },
688                RainbowAssetData { symbol: "C".into(), spot: 110.0, volatility: 0.35, dividend: None },
689            ],
690            correlations: vec![
691                vec![1.0, 0.5, 0.3],
692                vec![0.5, 1.0, 0.4],
693                vec![0.3, 0.4, 1.0],
694            ],
695            strike_price: Some(100.0),
696            weights: None,
697            maturity: maturity_1y(),
698            risk_free_rate: Some(0.05),
699            discount_curve: None,
700            pricer: Some("Analytical".into()),
701            simulation: Some(100_000),
702            mc_sampler: None,
703            mc_seed: None,
704            valuation_date: None,
705        };
706        let mut option = RainbowOption::from_json(&data);
707        let analytic = option.npv();
708        option.engine = PricingEngine::from_kind(Engine::MonteCarlo);
709        let mc = option.npv();
710        assert!((mc - analytic).abs() < 0.15, "mc={mc} moment-match={analytic}");
711    }
712
713    #[test]
714    fn best_of_plus_worst_of_equals_sum_of_vanillas() {
715        // max + min = S1 + S2 pathwise, so (max-K)+ + (min-K)+ = (S1-K)+ + (S2-K)+
716        let k = 100.0;
717        let best = two_asset("best_of", "C", Some(k), 0.6).npv();
718        let worst = two_asset("worst_of", "C", Some(k), 0.6).npv();
719        let t = two_asset("best_of", "C", Some(k), 0.6).time_to_maturity();
720        let vanillas = bs_price(100.0, k, 0.05, 0.02, 0.3, t, PutOrCall::Call)
721            + bs_price(95.0, k, 0.05, 0.01, 0.25, t, PutOrCall::Call);
722        assert!(
723            (best + worst - vanillas).abs() < 0.1,
724            "best {best} + worst {worst} vs vanillas {vanillas}"
725        );
726    }
727
728    #[test]
729    fn worst_of_call_at_zero_strike_is_forward_minus_margrabe() {
730        // min(S1, S2) = S2 - (S2 - S1)^+
731        let worst = two_asset("worst_of", "C", Some(1e-9), 0.6);
732        let t = worst.time_to_maturity();
733        let worst_pv = worst.npv();
734        let mut exchange_21 = two_asset("exchange", "P", None, 0.6); // pays (S2 - S1)^+
735        exchange_21.engine = PricingEngine::BlackScholes;
736        let expected = 95.0 * (-0.01 * t as f64).exp() - exchange_21.npv();
737        assert!((worst_pv - expected).abs() < 0.05, "{worst_pv} vs {expected}");
738    }
739
740    #[test]
741    fn correlation_orders_worst_of_prices() {
742        // higher correlation raises the worst-of call (the min rises)
743        let low = two_asset("worst_of", "C", Some(100.0), 0.0).npv();
744        let high = two_asset("worst_of", "C", Some(100.0), 0.9).npv();
745        assert!(high > low, "high-corr {high} must exceed low-corr {low}");
746    }
747
748    #[test]
749    fn monte_carlo_is_reproducible_and_reports_stats() {
750        let option = two_asset("worst_of", "C", Some(100.0), 0.6);
751        assert_eq!(option.npv(), option.npv());
752        let stats = option.npv_with_stats().unwrap();
753        assert!(stats.std_err > 0.0 && stats.std_err < 0.5);
754    }
755
756    #[test]
757    fn deltas_and_vegas_have_sensible_signs() {
758        let option = two_asset("spread", "C", Some(5.0), 0.6);
759        let deltas = option.deltas();
760        assert!(deltas[0] > 0.0, "long asset 1: {deltas:?}");
761        assert!(deltas[1] < 0.0, "short asset 2: {deltas:?}");
762        let vegas = option.vegas();
763        assert!(vegas[0] > 0.0);
764    }
765
766    #[test]
767    fn cholesky_rejects_invalid_correlations() {
768        assert!(cholesky(&[vec![1.0, 0.5], vec![0.4, 1.0]]).is_err()); // asymmetric
769        assert!(cholesky(&[vec![2.0, 0.0], vec![0.0, 1.0]]).is_err()); // diagonal != 1
770        // correlation > 1 in disguise: not positive definite
771        assert!(cholesky(&[
772            vec![1.0, 0.9, -0.9],
773            vec![0.9, 1.0, 0.9],
774            vec![-0.9, 0.9, 1.0]
775        ])
776        .is_err());
777        assert!(cholesky(&[vec![1.0, 0.5], vec![0.5, 1.0]]).is_ok());
778    }
779}
rustyqlib/equity/rainbow.rs

rustyqlib/equity/
rainbow.rs
