rustyqlib/equity/
barrier.rs

1//! Analytic pricing of continuously monitored barrier options
2//! (Reiner-Rubinstein 1991), all eight types: up/down x in/out x call/put,
3//! without rebate.
4//!
5//! Implemented as a pure function of the market inputs so Greeks can be
6//! taken by bumping arguments, and so the formulas can be validated
7//! independently of the option object (in-out parity, vanilla limits,
8//! Monte Carlo agreement).
9
10use crate::core::trade::PutOrCall;
11use crate::core::utils::norm_cdf;
12use super::blackscholes::bs_price;
13
14/// Which side the barrier sits on relative to the spot at inception.
15#[derive(Debug, Clone, Copy, PartialEq, Eq)]
16pub enum BarrierDirection {
17    Up,
18    Down,
19}
20
21/// Knock-in options come alive when the barrier is touched; knock-out
22/// options die.
23#[derive(Debug, Clone, Copy, PartialEq, Eq)]
24pub enum KnockType {
25    In,
26    Out,
27}
28
29/// Reiner-Rubinstein price of a European barrier option (no rebate).
30///
31/// If the spot is already at or beyond the barrier the option is treated as
32/// knocked: an `Out` option is worthless, an `In` option is the vanilla.
33#[allow(clippy::too_many_arguments)]
34pub fn barrier_price(
35    s: f64,
36    k: f64,
37    h: f64,
38    r: f64,
39    q: f64,
40    sigma: f64,
41    t: f64,
42    direction: BarrierDirection,
43    knock: KnockType,
44    put_or_call: PutOrCall,
45) -> f64 {
46    assert!(s > 0.0 && k > 0.0 && h > 0.0 && sigma > 0.0 && t > 0.0);
47    let down = direction == BarrierDirection::Down;
48    let knocked_now = if down { s <= h } else { s >= h };
49    if knocked_now {
50        return match knock {
51            KnockType::Out => 0.0,
52            KnockType::In => bs_price(s, k, r, q, sigma, t, put_or_call),
53        };
54    }
55
56    let call = put_or_call == PutOrCall::Call;
57    let phi: f64 = if call { 1.0 } else { -1.0 };
58    let eta: f64 = if down { 1.0 } else { -1.0 };
59    let st = sigma * t.sqrt();
60    let mu = (r - q - 0.5 * sigma * sigma) / (sigma * sigma);
61    let df_q = (-q * t).exp();
62    let df_r = (-r * t).exp();
63    let hs = h / s;
64
65    let x1 = (s / k).ln() / st + (1.0 + mu) * st;
66    let x2 = (s / h).ln() / st + (1.0 + mu) * st;
67    let y1 = (h * h / (s * k)).ln() / st + (1.0 + mu) * st;
68    let y2 = (h / s).ln() / st + (1.0 + mu) * st;
69
70    let a = phi * s * df_q * norm_cdf(phi * x1) - phi * k * df_r * norm_cdf(phi * x1 - phi * st);
71    let b = phi * s * df_q * norm_cdf(phi * x2) - phi * k * df_r * norm_cdf(phi * x2 - phi * st);
72    let c = phi * s * df_q * hs.powf(2.0 * (mu + 1.0)) * norm_cdf(eta * y1)
73        - phi * k * df_r * hs.powf(2.0 * mu) * norm_cdf(eta * y1 - eta * st);
74    let d = phi * s * df_q * hs.powf(2.0 * (mu + 1.0)) * norm_cdf(eta * y2)
75        - phi * k * df_r * hs.powf(2.0 * mu) * norm_cdf(eta * y2 - eta * st);
76
77    let k_above_barrier = k >= h;
78    let knock_in = match (call, down) {
79        (true, true) => if k_above_barrier { c } else { a - b + d },
80        (true, false) => if k_above_barrier { a } else { b - c + d },
81        (false, true) => if k_above_barrier { b - c + d } else { a },
82        (false, false) => if k_above_barrier { a - b + d } else { c },
83    };
84    match knock {
85        KnockType::In => knock_in,
86        // in-out parity (no rebate): out = vanilla - in
87        KnockType::Out => bs_price(s, k, r, q, sigma, t, put_or_call) - knock_in,
88    }
89}
90
91/// When a knock-out rebate is paid.
92#[derive(Debug, Clone, Copy, PartialEq, Eq)]
93pub enum RebateTiming {
94    /// Paid the moment the barrier is touched (Reiner-Rubinstein `F`).
95    AtHit,
96    /// Paid at expiry if the barrier was touched.
97    AtExpiry,
98}
99
100/// Present value of the rebate leg of a single-barrier option.
101///
102/// - **Knock-out**: `rebate` is paid when the barrier is hit — either at
103///   the touch time (`AtHit`, the Reiner-Rubinstein `F` term) or at
104///   expiry (`AtExpiry`, the complement of the survival term);
105/// - **Knock-in**: `rebate` is paid **at expiry** when the option never
106///   knocked in (the `E` term); `timing` is ignored.
107#[allow(clippy::too_many_arguments)]
108pub fn barrier_rebate_value(
109    s: f64,
110    h: f64,
111    rebate: f64,
112    r: f64,
113    q: f64,
114    sigma: f64,
115    t: f64,
116    direction: BarrierDirection,
117    knock: KnockType,
118    timing: RebateTiming,
119) -> f64 {
120    assert!(s > 0.0 && h > 0.0 && sigma > 0.0 && t > 0.0);
121    if rebate == 0.0 {
122        return 0.0;
123    }
124    let down = direction == BarrierDirection::Down;
125    let eta: f64 = if down { 1.0 } else { -1.0 };
126    let knocked_now = if down { s <= h } else { s >= h };
127    let df_r = (-r * t).exp();
128    if knocked_now {
129        return match knock {
130            KnockType::In => 0.0, // knocked in: no rebate
131            KnockType::Out => match timing {
132                RebateTiming::AtHit => rebate,
133                RebateTiming::AtExpiry => rebate * df_r,
134            },
135        };
136    }
137    let st = sigma * t.sqrt();
138    let mu = (r - q - 0.5 * sigma * sigma) / (sigma * sigma);
139    let hs = h / s;
140    // discounted probability of never touching the barrier (the E term)
141    let x2 = (s / h).ln() / st + (1.0 + mu) * st;
142    let y2 = (h / s).ln() / st + (1.0 + mu) * st;
143    let survival_pv =
144        rebate * df_r * (norm_cdf(eta * (x2 - st)) - hs.powf(2.0 * mu) * norm_cdf(eta * (y2 - st)));
145    match knock {
146        KnockType::In => survival_pv,
147        KnockType::Out => match timing {
148            RebateTiming::AtExpiry => rebate * df_r - survival_pv,
149            RebateTiming::AtHit => {
150                // first-touch value (the F term)
151                let lambda = (mu * mu + 2.0 * r / (sigma * sigma)).sqrt();
152                let z = (h / s).ln() / st + lambda * st;
153                rebate
154                    * (hs.powf(mu + lambda) * norm_cdf(eta * z)
155                        + hs.powf(mu - lambda) * norm_cdf(eta * (z - 2.0 * lambda * st)))
156            }
157        },
158    }
159}
160
161/// Single-barrier option with a rebate leg: the Reiner-Rubinstein price
162/// plus [`barrier_rebate_value`].
163#[allow(clippy::too_many_arguments)]
164pub fn barrier_price_with_rebate(
165    s: f64,
166    k: f64,
167    h: f64,
168    rebate: f64,
169    r: f64,
170    q: f64,
171    sigma: f64,
172    t: f64,
173    direction: BarrierDirection,
174    knock: KnockType,
175    timing: RebateTiming,
176    put_or_call: PutOrCall,
177) -> f64 {
178    barrier_price(s, k, h, r, q, sigma, t, direction, knock, put_or_call)
179        + barrier_rebate_value(s, h, rebate, r, q, sigma, t, direction, knock, timing)
180}
181
182/// Double-barrier option (flat lower `l` and upper `u` barriers), by the
183/// Ikeda-Kunitomo image-series expansion (continuous monitoring, no
184/// rebate). Knock-in prices through in-out parity against the vanilla.
185///
186/// The series converges extremely fast; five image pairs are far below
187/// f64 precision for practical inputs. The single-barrier limits
188/// (`l -> 0`, `u -> infinity`) reproduce the Reiner-Rubinstein prices
189/// (tested).
190#[allow(clippy::too_many_arguments)]
191pub fn double_barrier_price(
192    s: f64,
193    k: f64,
194    l: f64,
195    u: f64,
196    r: f64,
197    q: f64,
198    sigma: f64,
199    t: f64,
200    knock: KnockType,
201    put_or_call: PutOrCall,
202) -> f64 {
203    assert!(s > 0.0 && k > 0.0 && sigma > 0.0 && t > 0.0);
204    assert!(l < u, "lower barrier must be below the upper barrier");
205    let knocked_now = s <= l || s >= u;
206    let vanilla = bs_price(s, k, r, q, sigma, t, put_or_call);
207    if knocked_now {
208        return match knock {
209            KnockType::Out => 0.0,
210            KnockType::In => vanilla,
211        };
212    }
213    let b = r - q;
214    let st = sigma * t.sqrt();
215    let mu = 2.0 * b / (sigma * sigma) + 1.0;
216    let df_q = ((b - r) * t).exp();
217    let df_r = (-r * t).exp();
218    let call = put_or_call == PutOrCall::Call;
219    // the effective cap/floor of the payoff region inside the corridor
220    let f = if call { u } else { l };
221    let mut spot_sum = 0.0;
222    let mut strike_sum = 0.0;
223    for n in -5i32..=5 {
224        let un = u.powi(n);
225        let ln = l.powi(n);
226        let ratio1 = (un / ln).powf(mu);
227        let ratio2 = (l.powi(n + 1) / (un * s)).powf(mu);
228        let d1 = ((s * un * un / (k * ln * ln)).ln() + (b + 0.5 * sigma * sigma) * t) / st;
229        let d2 = ((s * un * un / (f * ln * ln)).ln() + (b + 0.5 * sigma * sigma) * t) / st;
230        let d3 = ((l.powi(2 * n + 2) / (k * s * un * un)).ln()
231            + (b + 0.5 * sigma * sigma) * t)
232            / st;
233        let d4 = ((l.powi(2 * n + 2) / (f * s * un * un)).ln()
234            + (b + 0.5 * sigma * sigma) * t)
235            / st;
236        if call {
237            spot_sum += ratio1 * (norm_cdf(d1) - norm_cdf(d2)) - ratio2 * (norm_cdf(d3) - norm_cdf(d4));
238            strike_sum += (un / ln).powf(mu - 2.0) * (norm_cdf(d1 - st) - norm_cdf(d2 - st))
239                - (l.powi(n + 1) / (un * s)).powf(mu - 2.0) * (norm_cdf(d3 - st) - norm_cdf(d4 - st));
240        } else {
241            // put: the payoff region is [l, k], integrated with positive
242            // normal arguments (d2 anchors the floor f = l, d1 the strike)
243            spot_sum += ratio1 * (norm_cdf(d2) - norm_cdf(d1)) - ratio2 * (norm_cdf(d4) - norm_cdf(d3));
244            strike_sum += (un / ln).powf(mu - 2.0) * (norm_cdf(d2 - st) - norm_cdf(d1 - st))
245                - (l.powi(n + 1) / (un * s)).powf(mu - 2.0) * (norm_cdf(d4 - st) - norm_cdf(d3 - st));
246        }
247    }
248    let phi = if call { 1.0 } else { -1.0 };
249    let out = phi * (s * df_q * spot_sum - k * df_r * strike_sum);
250    match knock {
251        KnockType::Out => out.max(0.0),
252        KnockType::In => (vanilla - out.max(0.0)).max(0.0),
253    }
254}
255
256#[cfg(test)]
257mod tests {
258    use super::*;
259
260    const S: f64 = 100.0;
261    const R: f64 = 0.05;
262    const Q: f64 = 0.02;
263    const SIG: f64 = 0.3;
264    const T: f64 = 1.0;
265
266    #[test]
267    fn rebate_identities_hold_exactly() {
268        use BarrierDirection::*;
269        let rebate = 5.0;
270        for (dir, h) in [(Down, 85.0), (Up, 115.0)] {
271            // hit and no-hit rebates at expiry are complementary events
272            let ki = barrier_rebate_value(S, h, rebate, R, Q, SIG, T, dir, KnockType::In,
273                RebateTiming::AtExpiry);
274            let ko_exp = barrier_rebate_value(S, h, rebate, R, Q, SIG, T, dir, KnockType::Out,
275                RebateTiming::AtExpiry);
276            assert!((ki + ko_exp - rebate * (-R * T).exp()).abs() < 1e-12, "{dir:?}");
277            // paying at the touch is worth more than waiting until expiry
278            let ko_hit = barrier_rebate_value(S, h, rebate, R, Q, SIG, T, dir, KnockType::Out,
279                RebateTiming::AtHit);
280            assert!(ko_hit > ko_exp, "{dir:?}: {ko_hit} vs {ko_exp}");
281            // with r = 0 timing is irrelevant
282            let hit0 = barrier_rebate_value(S, h, rebate, 0.0, Q, SIG, T, dir, KnockType::Out,
283                RebateTiming::AtHit);
284            let exp0 = barrier_rebate_value(S, h, rebate, 0.0, Q, SIG, T, dir, KnockType::Out,
285                RebateTiming::AtExpiry);
286            assert!((hit0 - exp0).abs() < 1e-10, "{dir:?}: {hit0} vs {exp0}");
287        }
288        // barrier far away: the knock-out never pays, the knock-in always
289        let far = barrier_rebate_value(S, 1e4, rebate, R, Q, SIG, T, Up, KnockType::Out,
290            RebateTiming::AtHit);
291        assert!(far < 1e-8, "{far}");
292        let sure = barrier_rebate_value(S, 1e4, rebate, R, Q, SIG, T, Up, KnockType::In,
293            RebateTiming::AtExpiry);
294        assert!((sure - rebate * (-R * T).exp()).abs() < 1e-8);
295        // already-touched knock-out pays the rebate immediately
296        let now = barrier_rebate_value(100.0, 100.0, rebate, R, Q, SIG, T, Up, KnockType::Out,
297            RebateTiming::AtHit);
298        assert!((now - rebate).abs() < 1e-12);
299    }
300
301    #[test]
302    fn rebate_at_hit_matches_a_first_touch_monte_carlo() {
303        use crate::core::montecarlo::path_rng;
304        use rand::Rng;
305        let (h, rebate) = (115.0, 10.0);
306        let analytic = barrier_rebate_value(S, h, rebate, R, Q, SIG, T,
307            BarrierDirection::Up, KnockType::Out, RebateTiming::AtHit);
308        // dense-grid first-touch simulation (discrete monitoring misses
309        // some touches, so the MC sits slightly below)
310        let (paths, steps) = (60_000, 2000);
311        let dt = T / steps as f64;
312        let drift = (R - Q - 0.5 * SIG * SIG) * dt;
313        let vol = SIG * dt.sqrt();
314        let mut sum = 0.0;
315        for i in 0..paths {
316            let mut rng = path_rng(97, i as u64);
317            let mut spot = S;
318            for step in 1..=steps {
319                let z: f64 = rng.sample(rand_distr::StandardNormal);
320                spot *= (drift + vol * z).exp();
321                if spot >= h {
322                    sum += rebate * (-R * step as f64 * dt).exp();
323                    break;
324                }
325            }
326        }
327        let mc = sum / paths as f64;
328        assert!(mc < analytic + 0.02, "discrete monitoring should undercount");
329        assert!((mc - analytic).abs() < 0.20, "mc {mc} vs analytic {analytic}");
330    }
331
332    #[test]
333    fn double_barrier_limits_reproduce_single_barriers() {
334        for pc in [PutOrCall::Call, PutOrCall::Put] {
335            for k in [90.0, 100.0, 110.0] {
336                // lower barrier -> 0: pure up-and-out
337                let dko = double_barrier_price(S, k, 10.0, 120.0, R, Q, SIG, T,
338                    KnockType::Out, pc);
339                let uo = barrier_price(S, k, 120.0, R, Q, SIG, T,
340                    BarrierDirection::Up, KnockType::Out, pc);
341                assert!((dko - uo).abs() < 1e-6, "{pc:?} K={k}: {dko} vs UO {uo}");
342                // upper barrier -> infinity: pure down-and-out
343                let dko2 = double_barrier_price(S, k, 80.0, 1000.0, R, Q, SIG, T,
344                    KnockType::Out, pc);
345                let down_out = barrier_price(S, k, 80.0, R, Q, SIG, T,
346                    BarrierDirection::Down, KnockType::Out, pc);
347                assert!((dko2 - down_out).abs() < 1e-6, "{pc:?} K={k}: {dko2} vs DO {down_out}");
348                // both far: the vanilla
349                let wide = double_barrier_price(S, k, 10.0, 1000.0, R, Q, SIG, T,
350                    KnockType::Out, pc);
351                let vanilla = crate::equity::blackscholes::bs_price(S, k, R, Q, SIG, T, pc);
352                assert!((wide - vanilla).abs() < 1e-6, "{pc:?} K={k}");
353            }
354        }
355    }
356
357    #[test]
358    fn double_barrier_parity_and_bounds() {
359        for pc in [PutOrCall::Call, PutOrCall::Put] {
360            let (l, u, k) = (85.0, 120.0, 100.0);
361            let out = double_barrier_price(S, k, l, u, R, Q, SIG, T, KnockType::Out, pc);
362            let inn = double_barrier_price(S, k, l, u, R, Q, SIG, T, KnockType::In, pc);
363            let vanilla = crate::equity::blackscholes::bs_price(S, k, R, Q, SIG, T, pc);
364            assert!((out + inn - vanilla).abs() < 1e-10, "{pc:?} parity");
365            // the corridor is more restrictive than either single barrier
366            let uo = barrier_price(S, k, u, R, Q, SIG, T,
367                BarrierDirection::Up, KnockType::Out, pc);
368            let down_out = barrier_price(S, k, l, R, Q, SIG, T,
369                BarrierDirection::Down, KnockType::Out, pc);
370            assert!(out <= uo + 1e-10 && out <= down_out + 1e-10, "{pc:?} bounds");
371            assert!(out > 0.0);
372        }
373    }
374
375    #[test]
376    fn double_barrier_matches_a_dense_monte_carlo() {
377        use crate::core::montecarlo::path_rng;
378        use rand::Rng;
379        let (l, u, k) = (85.0, 120.0, 100.0);
380        for pc in [PutOrCall::Call, PutOrCall::Put] {
381            let analytic = double_barrier_price(S, k, l, u, R, Q, SIG, T, KnockType::Out, pc);
382            let (paths, steps) = (60_000, 2000);
383            let dt = T / steps as f64;
384            let drift = (R - Q - 0.5 * SIG * SIG) * dt;
385            let vol = SIG * dt.sqrt();
386            let mut sum = 0.0;
387            for i in 0..paths {
388                let mut rng = path_rng(13, i as u64);
389                let mut spot = S;
390                let mut alive = true;
391                for _ in 0..steps {
392                    let z: f64 = rng.sample(rand_distr::StandardNormal);
393                    spot *= (drift + vol * z).exp();
394                    if spot <= l || spot >= u {
395                        alive = false;
396                        break;
397                    }
398                }
399                if alive {
400                    sum += match pc {
401                        PutOrCall::Call => (spot - k).max(0.0),
402                        PutOrCall::Put => (k - spot).max(0.0),
403                    };
404                }
405            }
406            let mc = (-R * T).exp() * sum / paths as f64;
407            // discrete monitoring survives more often -> MC above analytic
408            assert!(mc > analytic - 0.02, "{pc:?}: {mc} vs {analytic}");
409            assert!((mc - analytic).abs() < 0.25, "{pc:?}: mc {mc} vs analytic {analytic}");
410        }
411    }
412
413    #[test]
414    fn matches_independent_oracle_goldens() {
415        use BarrierDirection::*;
416        use KnockType::*;
417        use PutOrCall::*;
418        // values from an independently coded Reiner-Rubinstein implementation
419        let cases = [
420            (Down, In, Call, 90.0, 4.5095197744),
421            (Down, Out, Call, 90.0, 8.5107614943),
422            (Down, In, Put, 90.0, 10.0710164338),
423            (Down, Out, Put, 90.0, 0.0523399543),
424            (Up, In, Call, 120.0, 12.5974705742),
425            (Up, Out, Call, 120.0, 0.4228106946),
426            (Up, In, Put, 120.0, 1.4297711810),
427            (Up, Out, Put, 120.0, 8.6935852071),
428        ];
429        for (direction, knock, pc, h, expected) in cases {
430            let price = barrier_price(S, 100.0, h, R, Q, SIG, T, direction, knock, pc);
431            assert!(
432                (price - expected).abs() < 1e-8,
433                "{direction:?} {knock:?} {pc:?} H={h}: {price} vs {expected}"
434            );
435        }
436    }
437
438    #[test]
439    fn in_plus_out_equals_vanilla() {
440        for pc in [PutOrCall::Call, PutOrCall::Put] {
441            for k in [90.0, 100.0, 110.0] {
442                for (direction, h) in [
443                    (BarrierDirection::Down, 80.0),
444                    (BarrierDirection::Down, 99.0),
445                    (BarrierDirection::Up, 101.0),
446                    (BarrierDirection::Up, 130.0),
447                ] {
448                    let vanilla = bs_price(S, k, R, Q, SIG, T, pc);
449                    let ki = barrier_price(S, k, h, R, Q, SIG, T, direction, KnockType::In, pc);
450                    let ko = barrier_price(S, k, h, R, Q, SIG, T, direction, KnockType::Out, pc);
451                    assert!(
452                        (ki + ko - vanilla).abs() < 1e-10,
453                        "{pc:?} K={k} {direction:?} H={h}: {ki} + {ko} != {vanilla}"
454                    );
455                }
456            }
457        }
458    }
459
460    #[test]
461    fn far_barriers_reduce_to_vanilla_or_zero() {
462        let vanilla_call = bs_price(S, 100.0, R, Q, SIG, T, PutOrCall::Call);
463        // barrier so far away it is never hit: out = vanilla, in = 0
464        let ko = barrier_price(S, 100.0, 1e-4, R, Q, SIG, T, BarrierDirection::Down, KnockType::Out, PutOrCall::Call);
465        let ki = barrier_price(S, 100.0, 1e-4, R, Q, SIG, T, BarrierDirection::Down, KnockType::In, PutOrCall::Call);
466        assert!((ko - vanilla_call).abs() < 1e-9);
467        assert!(ki.abs() < 1e-9);
468        let ko_up = barrier_price(S, 100.0, 1e6, R, Q, SIG, T, BarrierDirection::Up, KnockType::Out, PutOrCall::Call);
469        assert!((ko_up - vanilla_call).abs() < 1e-9);
470    }
471
472    #[test]
473    fn already_knocked_positions() {
474        let vanilla = bs_price(S, 100.0, R, Q, SIG, T, PutOrCall::Call);
475        // spot at the barrier counts as touched
476        let ko = barrier_price(S, 100.0, 100.0, R, Q, SIG, T, BarrierDirection::Down, KnockType::Out, PutOrCall::Call);
477        let ki = barrier_price(S, 100.0, 100.0, R, Q, SIG, T, BarrierDirection::Down, KnockType::In, PutOrCall::Call);
478        assert_eq!(ko, 0.0);
479        assert!((ki - vanilla).abs() < 1e-12);
480    }
481
482    #[test]
483    fn up_out_call_with_strike_above_barrier_is_worthless() {
484        // any payoff requires S_T > K >= H, which forces a knock
485        let price = barrier_price(
486            S, 110.0, 105.0, R, Q, SIG, T,
487            BarrierDirection::Up, KnockType::Out, PutOrCall::Call,
488        );
489        assert!(price.abs() < 1e-12, "{price}");
490    }
491}
rustyqlib/equity/barrier.rs

rustyqlib/equity/
barrier.rs
