From 34e7871005aa6c1277205e8e0f348b4ab5dcb790 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sun, 26 Apr 2026 19:03:58 +0100 Subject: [PATCH 01/12] =?UTF-8?q?Fix=20barrier=20analytical=20zero/tiny-vo?= =?UTF-8?q?l=20NaN=20-=20Relax=20UnderlyingData=20to=20allow=20=CF=83=3D0?= =?UTF-8?q?=20-=20Add=20shared=20coerce=5Fpositive=5Ffloat=20helper=20acro?= =?UTF-8?q?ss=20underlying=20and=20spec=20dataclasses?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../stochastic_processes.py | 42 +++--- src/derivatives_pricing/utils.py | 64 ++++++++- .../valuation/barrier_analytical.py | 107 +++++++++++++-- .../valuation/contracts.py | 94 ++++--------- src/derivatives_pricing/valuation/core.py | 13 +- tests/test_edge_cases.py | 124 ++++++++++++++++++ tests/test_stochastic_processes.py | 19 +++ 7 files changed, 359 insertions(+), 104 deletions(-) diff --git a/src/derivatives_pricing/stochastic_processes.py b/src/derivatives_pricing/stochastic_processes.py index f26f035..4fe1aeb 100644 --- a/src/derivatives_pricing/stochastic_processes.py +++ b/src/derivatives_pricing/stochastic_processes.py @@ -11,7 +11,7 @@ import pandas as pd from .market_environment import MarketData, CorrelationContext from .enums import DayCountConvention -from .utils import calculate_year_fraction, validate_naive_datetime +from .utils import calculate_year_fraction, coerce_positive_float, validate_naive_datetime from .rates import DiscountCurve from .exceptions import ConfigurationError, ValidationError @@ -121,16 +121,16 @@ class GBMParams: dividend_curve: DiscountCurve | None = None def __post_init__(self) -> None: - if self.initial_value is None: - raise ValidationError("GBMParams requires initial_value to be not None") - if self.volatility is None: - raise ValidationError("GBMParams requires volatility to be not None") - if not np.isfinite(float(self.initial_value)): - raise ValidationError("GBMParams requires initial_value to be finite") - if not np.isfinite(float(self.volatility)): - raise ValidationError("GBMParams requires volatility to be finite") - if float(self.volatility) < 0.0: - raise ValidationError("GBMParams requires volatility to be >= 0") + object.__setattr__( + self, + "initial_value", + coerce_positive_float(self.initial_value, name="GBMParams.initial_value"), + ) + object.__setattr__( + self, + "volatility", + coerce_positive_float(self.volatility, name="GBMParams.volatility", strict=False), + ) object.__setattr__( self, "discrete_dividends", @@ -177,16 +177,16 @@ class JDParams: dividend_curve: DiscountCurve | None = None def __post_init__(self) -> None: - if self.initial_value is None: - raise ValidationError("JDParams requires initial_value to be not None") - if self.volatility is None: - raise ValidationError("JDParams requires volatility to be not None") - if not np.isfinite(float(self.initial_value)): - raise ValidationError("JDParams requires initial_value to be finite") - if not np.isfinite(float(self.volatility)): - raise ValidationError("JDParams requires volatility to be finite") - if float(self.volatility) < 0.0: - raise ValidationError("JDParams requires volatility to be >= 0") + object.__setattr__( + self, + "initial_value", + coerce_positive_float(self.initial_value, name="JDParams.initial_value"), + ) + object.__setattr__( + self, + "volatility", + coerce_positive_float(self.volatility, name="JDParams.volatility", strict=False), + ) if self.lambd is None or self.mu is None or self.delta is None: raise ValidationError("JDParams requires lambd, mu, and delta to be not None") diff --git a/src/derivatives_pricing/utils.py b/src/derivatives_pricing/utils.py index b9fb747..5ea9a9c 100644 --- a/src/derivatives_pricing/utils.py +++ b/src/derivatives_pricing/utils.py @@ -12,7 +12,7 @@ import numpy as np from .enums import DayCountConvention, OptionType -from .exceptions import ArbitrageViolationError, ValidationError +from .exceptions import ArbitrageViolationError, ConfigurationError, ValidationError if TYPE_CHECKING: from .rates import DiscountCurve @@ -20,6 +20,7 @@ __all__ = [ "log_timing", "validate_naive_datetime", + "coerce_positive_float", "calculate_year_fraction", "pv_discrete_dividends", "forward_price", @@ -46,6 +47,67 @@ def validate_naive_datetime( raise ValidationError(f"{field_name} must be timezone-naive datetime") +def coerce_positive_float( + value: object, + *, + name: str, + strict: bool = True, + allow_none: bool = False, +) -> float: + """Coerce ``value`` to ``float`` and validate it is finite and positive. + + Centralises the float-coercion / finite / sign-check pattern used by + dataclasses (``UnderlyingData``, ``GBMParams``, ``BarrierSpec``, + etc.) so that all fields raise a consistent ``ConfigurationError`` for + non-numeric input and ``ValidationError`` for non-finite or wrong-sign + input. + + Parameters + ---------- + value + The user-supplied value. May be int, float, or anything ``float()`` + can accept. + name + Human-readable field name used in error messages + (e.g. ``"UnderlyingData.initial_value"``). + strict + If True (default) require ``value > 0``; if False require + ``value >= 0``. Use ``strict=False`` for fields where zero is a + meaningful boundary value (volatility, rebate, strike). + allow_none + If True, ``None`` returns ``None`` (caller handles default). + Otherwise raise ``ValidationError``. + + Returns + ------- + float + The coerced value. + + Raises + ------ + ConfigurationError + If ``value`` cannot be coerced to ``float``. + ValidationError + If ``value`` is None (and ``allow_none`` is False), non-finite, + or violates the sign constraint. + """ + if value is None: + if allow_none: + return None # type: ignore[return-value] + raise ValidationError(f"{name} must be provided") + try: + x = float(value) # type: ignore[arg-type] + except (TypeError, ValueError) as exc: + raise ConfigurationError(f"{name} must be numeric") from exc + if not np.isfinite(x): + raise ValidationError(f"{name} must be finite") + if strict and x <= 0.0: + raise ValidationError(f"{name} must be > 0, got {x}") + if not strict and x < 0.0: + raise ValidationError(f"{name} must be >= 0, got {x}") + return x + + @contextmanager def log_timing(logger, label: str, enabled: bool) -> Iterator[None]: """Context manager that logs elapsed time for a code block. diff --git a/src/derivatives_pricing/valuation/barrier_analytical.py b/src/derivatives_pricing/valuation/barrier_analytical.py index 8e9532c..af7d33b 100644 --- a/src/derivatives_pricing/valuation/barrier_analytical.py +++ b/src/derivatives_pricing/valuation/barrier_analytical.py @@ -415,6 +415,58 @@ def _barrier_price_no_rebate( return B - D if H <= K else A - C +def _deterministic_limit_price( + S: float, + K: float, + H: float, + r: float, + q: float, + T: float, + df_r: float, + option_type: OptionType, + direction: BarrierDirection, + action: BarrierAction, + rebate: float, + rebate_timing: RebateTiming, +) -> float: + """Closed-form barrier price in the sigma -> 0 (deterministic-drift) limit. + + Under GBM with sigma = 0 the spot evolves monotonically from ``S`` to the + forward ``S_T = S * exp((r - q) * T)`` and continuous barrier monitoring + reduces to checking whether the linear path [S, S_T] crosses the barrier. + Caller must have already confirmed the barrier is not triggered at + inception (``_barrier_triggered_at_inception`` returned False). + + Notes + ----- + Discrete monitoring collapses to continuous monitoring in this limit + because the deterministic path is monotone and the Broadie-Glasserman-Kou + shift ``H * exp(beta * sigma * sqrt(dt))`` collapses to ``H``. + """ + S_T = S * np.exp((r - q) * T) + + if direction is BarrierDirection.UP: + hit = S_T >= H + else: + hit = S_T <= H + + if option_type is OptionType.CALL: + intrinsic_T = max(S_T - K, 0.0) + else: + intrinsic_T = max(K - S_T, 0.0) + vanilla_pv = df_r * intrinsic_T + + if action is BarrierAction.OUT: + if hit: + if rebate > 0.0: + return rebate if rebate_timing is RebateTiming.AT_HIT else df_r * rebate + return 0.0 + return vanilla_pv + if hit: + return vanilla_pv + return df_r * rebate if rebate > 0.0 else 0.0 + + # ── Engine class ────────────────────────────────────────────────── @@ -515,20 +567,47 @@ def present_value(self) -> float: H = _broadie_glasserman_adjustment(H, sigma, T, spec.num_observations, spec.direction) # ── No-rebate barrier value ── - value = _barrier_price_no_rebate( - S, - K, - H, - r, - q, - sigma, - T, - df_r, - df_q, - spec.option_type, - spec.direction, - spec.action, - ) + # The Reiner-Rubinstein formulas contain (H/S)**(2*lambda) with + # lambda = (r - q + sigma^2/2)/sigma^2, and d1/x1/y/y1 each divide + # by sigma*sqrt(T). As sigma -> 0 we either divide by zero (when + # computing lambda or any of the d/x/y terms) or blow up at + # (H/S)**(2*lambda). Most operations promote to numpy.float64 (via + # np.log/np.sqrt) so the failure mode is silent inf/nan + + # RuntimeWarning -- np.errstate converts those into + # FloatingPointError so we can fall back to the closed-form + # deterministic-forward price. OverflowError / ZeroDivisionError + # are caught in case any operation stays in Python-float arithmetic. + with np.errstate(over="raise", invalid="raise", divide="raise"): + try: + value = _barrier_price_no_rebate( + S, + K, + H, + r, + q, + sigma, + T, + df_r, + df_q, + spec.option_type, + spec.direction, + spec.action, + ) + except (OverflowError, FloatingPointError, ZeroDivisionError): + return _deterministic_limit_price( + S, + K, + H, + r, + q, + T, + df_r, + spec.option_type, + spec.direction, + spec.action, + spec.rebate, + spec.rebate_timing, + ) # ── Add rebate leg ── if spec.rebate > 0.0: diff --git a/src/derivatives_pricing/valuation/contracts.py b/src/derivatives_pricing/valuation/contracts.py index f457aaf..7be0027 100644 --- a/src/derivatives_pricing/valuation/contracts.py +++ b/src/derivatives_pricing/valuation/contracts.py @@ -18,7 +18,7 @@ RebateTiming, ) from ..exceptions import ConfigurationError, ValidationError -from ..utils import validate_naive_datetime +from ..utils import coerce_positive_float, validate_naive_datetime @dataclass(frozen=True, slots=True) @@ -71,17 +71,11 @@ def __post_init__(self) -> None: f"exercise_type must be ExerciseType enum, got {type(self.exercise_type).__name__}" ) - if self.strike is None: - raise ValidationError("VanillaSpec.strike must be provided") - try: - strike = float(self.strike) - except (TypeError, ValueError) as exc: - raise ConfigurationError("VanillaSpec.strike must be numeric") from exc - if not np.isfinite(strike): - raise ValidationError("VanillaSpec.strike must be finite") - if strike < 0.0: - raise ValidationError("VanillaSpec.strike must be >= 0") - object.__setattr__(self, "strike", strike) + object.__setattr__( + self, + "strike", + coerce_positive_float(self.strike, name="VanillaSpec.strike", strict=False), + ) @dataclass(frozen=True, slots=True) @@ -293,17 +287,11 @@ def __post_init__(self) -> None: if self.option_type not in (OptionType.CALL, OptionType.PUT): raise ValidationError("AsianSpec.option_type must be OptionType.CALL or OptionType.PUT") - if self.strike is None: - raise ValidationError("AsianSpec.strike must be provided") - try: - strike = float(self.strike) - except (TypeError, ValueError) as exc: - raise ConfigurationError("AsianSpec.strike must be numeric") from exc - if not np.isfinite(strike): - raise ValidationError("AsianSpec.strike must be finite") - if strike < 0.0: - raise ValidationError("AsianSpec.strike must be >= 0") - object.__setattr__(self, "strike", strike) + object.__setattr__( + self, + "strike", + coerce_positive_float(self.strike, name="AsianSpec.strike", strict=False), + ) # Exactly one schedule source is required. if (self.fixing_dates is None) == (self.num_observations is None): @@ -347,15 +335,11 @@ def __post_init__(self) -> None: "observed_average and observed_count must both be provided or both omitted." ) if self.observed_average is not None: - try: - obs_avg = float(self.observed_average) - except (TypeError, ValueError) as exc: - raise ConfigurationError("observed_average must be numeric") from exc - if not np.isfinite(obs_avg): - raise ValidationError("observed_average must be finite") - if obs_avg <= 0.0: - raise ValidationError("observed_average must be > 0") - object.__setattr__(self, "observed_average", obs_avg) + object.__setattr__( + self, + "observed_average", + coerce_positive_float(self.observed_average, name="observed_average"), + ) if not isinstance(self.observed_count, int) or self.observed_count < 1: raise ValidationError("observed_count must be a positive integer") @@ -464,41 +448,17 @@ def __post_init__(self) -> None: f"rebate_timing must be RebateTiming enum, got {type(self.rebate_timing).__name__}" ) - # --- strike --- - if self.strike is None: - raise ValidationError("BarrierSpec.strike must be provided") - try: - strike = float(self.strike) - except (TypeError, ValueError) as exc: - raise ConfigurationError("BarrierSpec.strike must be numeric") from exc - if not np.isfinite(strike): - raise ValidationError("BarrierSpec.strike must be finite") - if strike < 0.0: - raise ValidationError("BarrierSpec.strike must be >= 0") - object.__setattr__(self, "strike", strike) - - # --- barrier --- - if self.barrier is None: - raise ValidationError("BarrierSpec.barrier must be provided") - try: - barrier = float(self.barrier) - except (TypeError, ValueError) as exc: - raise ConfigurationError("BarrierSpec.barrier must be numeric") from exc - if not np.isfinite(barrier): - raise ValidationError("BarrierSpec.barrier must be finite") - if barrier <= 0.0: - raise ValidationError("BarrierSpec.barrier must be > 0") - object.__setattr__(self, "barrier", barrier) - - # --- rebate --- - try: - rebate = float(self.rebate) - except (TypeError, ValueError) as exc: - raise ConfigurationError("BarrierSpec.rebate must be numeric") from exc - if not np.isfinite(rebate): - raise ValidationError("BarrierSpec.rebate must be finite") - if rebate < 0.0: - raise ValidationError("BarrierSpec.rebate must be >= 0") + object.__setattr__( + self, + "strike", + coerce_positive_float(self.strike, name="BarrierSpec.strike", strict=False), + ) + object.__setattr__( + self, + "barrier", + coerce_positive_float(self.barrier, name="BarrierSpec.barrier"), + ) + rebate = coerce_positive_float(self.rebate, name="BarrierSpec.rebate", strict=False) object.__setattr__(self, "rebate", rebate) # Knock-in rebate must be paid at expiry (AT_HIT is contradictory — diff --git a/src/derivatives_pricing/valuation/core.py b/src/derivatives_pricing/valuation/core.py index 13ea9a2..a94dd0f 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -27,7 +27,7 @@ import threading import numpy as np import pandas as pd -from ..utils import calculate_year_fraction +from ..utils import calculate_year_fraction, coerce_positive_float from ..stochastic_processes import PathSimulation, GBMProcess from ..exceptions import ConfigurationError, UnsupportedFeatureError, ValidationError from ..enums import ( @@ -152,6 +152,17 @@ class UnderlyingData: dividend_curve: DiscountCurve | None = None def __post_init__(self) -> None: + object.__setattr__( + self, + "initial_value", + coerce_positive_float(self.initial_value, name="UnderlyingData.initial_value"), + ) + object.__setattr__( + self, + "volatility", + coerce_positive_float(self.volatility, name="UnderlyingData.volatility", strict=False), + ) + if self.discrete_dividends is not None: cleaned: list[tuple[dt.datetime, float]] = [] for ex_date, amount in self.discrete_dividends: diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py index 73c4dd5..12d9c3a 100644 --- a/tests/test_edge_cases.py +++ b/tests/test_edge_cases.py @@ -10,6 +10,9 @@ import pytest from derivatives_pricing.enums import ( + BarrierAction, + BarrierDirection, + BarrierMonitoring, ExerciseType, OptionType, PDESpaceGrid, @@ -28,6 +31,7 @@ pv as _pv, ) from derivatives_pricing.valuation import ( + BarrierSpec, BinomialParams, VanillaSpec, PDEParams, @@ -141,6 +145,23 @@ def test_pde_zero_vol_raises(self): with pytest.raises(ValidationError, match="volatility must be positive"): _pv(ud, spec, PricingMethod.PDE_FD) + # --- Construction-time validation: negative vol & non-positive spot rejected --- + + @pytest.mark.parametrize("vol", [-1e-12, -0.2, -1.0]) + def test_underlying_rejects_negative_vol(self, vol): + """UnderlyingData rejects σ < 0; σ = 0 is the deterministic limit and is allowed.""" + with pytest.raises(ValidationError, match=r"volatility must be >= 0"): + _underlying(vol=vol) + + def test_underlying_rejects_nan_vol(self): + with pytest.raises(ValidationError, match=r"volatility must be finite"): + _underlying(vol=float("nan")) + + @pytest.mark.parametrize("spot", [0.0, -1e-12, -100.0]) + def test_underlying_rejects_non_positive_spot(self, spot): + with pytest.raises(ValidationError, match=r"initial_value must be > 0"): + _underlying(spot=spot) + # ═══════════════════════════════════════════════════════════════════════ # Very low volatility (non-zero, but tiny) @@ -594,3 +615,106 @@ def test_american_deep_itm_put_pde(self): ) intrinsic = strike - spot assert pv >= intrinsic - 0.01 + + +# ═══════════════════════════════════════════════════════════════════════ +# Barrier options: tiny / zero σ overflow guard +# ═══════════════════════════════════════════════════════════════════════ + + +class TestBarrierZeroVolatility: + """BSM barrier formula contains ``(H/S)**(2*lambda)`` with + ``lambda = (r-q+sigma^2/2)/sigma^2``. As σ → 0 either λ overflows or + the divide-by-zero kicks in. The engine must fall back to the + deterministic-forward limit instead of returning silent NaN. + """ + + def _barrier_spec( + self, + direction: BarrierDirection, + action: BarrierAction, + barrier: float, + strike: float = 100.0, + option_type: OptionType = OptionType.CALL, + ) -> BarrierSpec: + return BarrierSpec( + option_type=option_type, + exercise_type=ExerciseType.EUROPEAN, + strike=strike, + maturity=MATURITY, + currency="USD", + barrier=barrier, + direction=direction, + action=action, + monitoring=BarrierMonitoring.CONTINUOUS, + ) + + @staticmethod + def _T_and_df(ud) -> tuple[float, float]: + """Year fraction and risk-free DF the engine actually uses.""" + from derivatives_pricing.utils import calculate_year_fraction + + T = calculate_year_fraction(PRICING_DATE, MATURITY) + return T, float(ud.discount_curve.df(T)) + + @pytest.mark.parametrize("vol", [0.0, 1e-12, 1e-10, 1e-8, 1e-4]) + def test_uoc_tiny_vol_returns_deterministic_forward(self, vol): + """UOC with tiny σ: forward = S·exp(rT) < H, option survives → + discounted forward intrinsic.""" + spot, strike, barrier = 100.0, 100.0, 130.0 + ud = _underlying(spot=spot, vol=vol) + spec = self._barrier_spec(BarrierDirection.UP, BarrierAction.OUT, barrier, strike) + pv = _pv(ud, spec, PricingMethod.BSM) + + T, df_r = self._T_and_df(ud) + S_T = spot * np.exp(RATE * T) + assert S_T < barrier # not knocked out by forward + expected = df_r * max(S_T - strike, 0.0) + assert np.isfinite(pv) + assert np.isclose(pv, expected, rtol=1e-10) + + @pytest.mark.parametrize("vol", [0.0, 1e-12, 1e-8]) + def test_uoc_tiny_vol_knocked_out_at_forward(self, vol): + """UOC where forward crosses barrier deterministically → knocked out → 0.""" + spot, strike, barrier = 100.0, 100.0, 102.0 # forward (~102.54) > barrier + ud = _underlying(spot=spot, vol=vol) + spec = self._barrier_spec(BarrierDirection.UP, BarrierAction.OUT, barrier, strike) + pv = _pv(ud, spec, PricingMethod.BSM) + + T, _ = self._T_and_df(ud) + S_T = spot * np.exp(RATE * T) + assert S_T > barrier # knocked out by forward + assert np.isfinite(pv) + assert pv == 0.0 + + @pytest.mark.parametrize("vol", [0.0, 1e-12, 1e-8]) + def test_dop_tiny_vol_returns_deterministic_forward(self, vol): + """Down-and-out put with tiny σ: forward stays above down-barrier → survives.""" + spot, strike, barrier = 100.0, 110.0, 80.0 + ud = _underlying(spot=spot, vol=vol) + spec = self._barrier_spec( + BarrierDirection.DOWN, + BarrierAction.OUT, + barrier, + strike, + option_type=OptionType.PUT, + ) + pv = _pv(ud, spec, PricingMethod.BSM) + + T, df_r = self._T_and_df(ud) + S_T = spot * np.exp(RATE * T) + assert S_T > barrier # not knocked out + expected = df_r * max(strike - S_T, 0.0) + assert np.isfinite(pv) + assert np.isclose(pv, expected, rtol=1e-10) + + @pytest.mark.parametrize("vol", [0.0, 1e-12, 1e-8]) + def test_uic_tiny_vol_not_hit_returns_zero(self, vol): + """Up-and-in call with tiny σ: forward never reaches barrier → KI never + activates, no rebate → 0.""" + spot, strike, barrier = 100.0, 100.0, 130.0 + ud = _underlying(spot=spot, vol=vol) + spec = self._barrier_spec(BarrierDirection.UP, BarrierAction.IN, barrier, strike) + pv = _pv(ud, spec, PricingMethod.BSM) + assert np.isfinite(pv) + assert pv == 0.0 diff --git a/tests/test_stochastic_processes.py b/tests/test_stochastic_processes.py index 49bad20..9b1c8c1 100644 --- a/tests/test_stochastic_processes.py +++ b/tests/test_stochastic_processes.py @@ -179,6 +179,25 @@ def test_gbm_drift_matches_risk_free_rate(self): assert np.abs(mean_return - expected_drift) < 0.01 +class TestGBMParamsValidation: + """GBMParams rejects negative volatility and non-positive initial value.""" + + @pytest.mark.parametrize("vol", [-1e-12, -0.2, -1.0]) + def test_rejects_negative_volatility(self, vol): + with pytest.raises(ValidationError, match=r"volatility must be >= 0"): + GBMParams(initial_value=100.0, volatility=vol) + + @pytest.mark.parametrize("spot", [0.0, -1e-12, -100.0]) + def test_rejects_non_positive_initial_value(self, spot): + with pytest.raises(ValidationError, match=r"initial_value must be > 0"): + GBMParams(initial_value=spot, volatility=0.2) + + def test_zero_volatility_allowed(self): + """σ = 0 is the deterministic-drift limit; GBMParams accepts it.""" + params = GBMParams(initial_value=100.0, volatility=0.0) + assert params.volatility == 0.0 + + class TestSRDProcess: """Tests for the SRDProcess (CIR) class.""" From c1d2b29b9ece8fdd30a1bbdf700751b238c3369d Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sun, 26 Apr 2026 19:13:04 +0100 Subject: [PATCH 02/12] Add BarrierSpec.is_triggered method. - Remove helper from barrier_analytical.py --- .../valuation/barrier_analytical.py | 26 ------------------- src/derivatives_pricing/valuation/binomial.py | 3 +-- .../valuation/contracts.py | 10 +++++++ src/derivatives_pricing/valuation/core.py | 4 +-- 4 files changed, 13 insertions(+), 30 deletions(-) diff --git a/src/derivatives_pricing/valuation/barrier_analytical.py b/src/derivatives_pricing/valuation/barrier_analytical.py index af7d33b..811fe45 100644 --- a/src/derivatives_pricing/valuation/barrier_analytical.py +++ b/src/derivatives_pricing/valuation/barrier_analytical.py @@ -49,32 +49,6 @@ # ── Shared helpers ────────────────────────────────────────────────── -def _is_triggered( - spot: float, - barrier: float, - direction: BarrierDirection, -) -> bool: - """Return ``True`` if the barrier is triggered. - - Parameters - ---------- - spot - Current spot price. - barrier - Barrier level. - direction - UP or DOWN. - - Returns - ------- - bool - ``True`` if the barrier is triggered. - """ - if direction is BarrierDirection.UP: - return spot >= barrier - return spot <= barrier - - _BG_BETA = 0.5826 # Broadie-Glasserman-Kou constant diff --git a/src/derivatives_pricing/valuation/binomial.py b/src/derivatives_pricing/valuation/binomial.py index 70d4014..8e28c6c 100644 --- a/src/derivatives_pricing/valuation/binomial.py +++ b/src/derivatives_pricing/valuation/binomial.py @@ -28,7 +28,6 @@ ValidationError, ) from .params import BinomialParams -from .barrier_analytical import _is_triggered if TYPE_CHECKING: from .core import AsianSpec, BarrierSpec, OptionValuation, UnderlyingData @@ -972,7 +971,7 @@ def _resolve_effective_num_steps(self) -> int: # `_solve_backward` (not the barrier-aware solver), so barrier # alignment is irrelevant. # Skip both the inflation and any warning for those cases. - if _is_triggered(spot, barrier, self.spec.direction): + if self.spec.is_triggered(spot): return base_steps log_distance = np.log(max(spot, barrier) / min(spot, barrier)) diff --git a/src/derivatives_pricing/valuation/contracts.py b/src/derivatives_pricing/valuation/contracts.py index 7be0027..8efc78f 100644 --- a/src/derivatives_pricing/valuation/contracts.py +++ b/src/derivatives_pricing/valuation/contracts.py @@ -507,3 +507,13 @@ def __post_init__(self) -> None: if dates[-1] > self.maturity: raise ValidationError("monitoring_dates must not extend beyond maturity.") object.__setattr__(self, "monitoring_dates", dates) + + def is_triggered(self, spot: float) -> bool: + """Return ``True`` if ``spot`` has crossed the barrier. + + UP barriers trigger when ``spot >= self.barrier``; DOWN barriers + trigger when ``spot <= self.barrier``. + """ + if self.direction is BarrierDirection.UP: + return spot >= self.barrier + return spot <= self.barrier diff --git a/src/derivatives_pricing/valuation/core.py b/src/derivatives_pricing/valuation/core.py index a94dd0f..12c9fcc 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -56,7 +56,7 @@ ) from .bsm import _BSMEuropeanValuation from .asian_analytical import _AnalyticalAsianValuation -from .barrier_analytical import _AnalyticalBarrierValuation, _is_triggered +from .barrier_analytical import _AnalyticalBarrierValuation from .pde import _FDEuropeanValuation, _FDAmericanValuation, _FDBarrierValuation from ..rates import DiscountCurve from ..market_environment import MarketData @@ -1019,7 +1019,7 @@ def _barrier_triggered_at_inception(self) -> bool: "the dispatcher should route BarrierSpecs to barrier engines only." ) spot = float(self._underlying.initial_value) - if not _is_triggered(spot, self._spec.barrier, self._spec.direction): + if not self._spec.is_triggered(spot): return False if self._spec.monitoring is BarrierMonitoring.CONTINUOUS: return True From 34baa86421e7e72e7b7b615ccf640fa15e4dc5dd Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sun, 26 Apr 2026 19:25:27 +0100 Subject: [PATCH 03/12] Handle sigma=0 in _brownian_bridge_hit_prob (MC engine) --- .../valuation/monte_carlo.py | 28 ++++++++++++------- 1 file changed, 18 insertions(+), 10 deletions(-) diff --git a/src/derivatives_pricing/valuation/monte_carlo.py b/src/derivatives_pricing/valuation/monte_carlo.py index 4fdb856..eba3626 100644 --- a/src/derivatives_pricing/valuation/monte_carlo.py +++ b/src/derivatives_pricing/valuation/monte_carlo.py @@ -1127,20 +1127,28 @@ def _brownian_bridge_hit_prob( """ p = np.zeros_like(S_i) + # At sigma == 0 the path is deterministic and monotone, so any in-step + # crossing must show up at an endpoint -- the `crossed` mask catches it + # and the Brownian-bridge correction is both undefined (1/sigma**2) and + # unnecessary (interior hit probability is 0). + sigma_positive = sigma > 0.0 + if direction is BarrierDirection.UP: crossed = (S_i >= barrier) | (S_next >= barrier) - both_below = ~crossed - if np.any(both_below): - log_a = np.log(barrier / S_i[both_below]) - log_b = np.log(barrier / S_next[both_below]) - p[both_below] = np.exp(-2.0 * log_a * log_b / (sigma**2 * dt_step)) + if sigma_positive: + both_below = ~crossed + if np.any(both_below): + log_a = np.log(barrier / S_i[both_below]) + log_b = np.log(barrier / S_next[both_below]) + p[both_below] = np.exp(-2.0 * log_a * log_b / (sigma**2 * dt_step)) else: crossed = (S_i <= barrier) | (S_next <= barrier) - both_above = ~crossed - if np.any(both_above): - log_a = np.log(S_i[both_above] / barrier) - log_b = np.log(S_next[both_above] / barrier) - p[both_above] = np.exp(-2.0 * log_a * log_b / (sigma**2 * dt_step)) + if sigma_positive: + both_above = ~crossed + if np.any(both_above): + log_a = np.log(S_i[both_above] / barrier) + log_b = np.log(S_next[both_above] / barrier) + p[both_above] = np.exp(-2.0 * log_a * log_b / (sigma**2 * dt_step)) p[crossed] = 1.0 return np.clip(p, 0.0, 1.0) From 37ab5d8cae91558a3fac6913f5f9b7c2fab318f7 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sun, 26 Apr 2026 19:56:18 +0100 Subject: [PATCH 04/12] =?UTF-8?q?Add=20=CF=83=3D0=20guard=20to=20MC=20like?= =?UTF-8?q?lihood-ratio=20greeks=20-=20Also=20add=20three=20coverage=20tes?= =?UTF-8?q?ts=20(Asian=20rejection,=20LR-greek=20rejection,=20BSM/MC=20van?= =?UTF-8?q?illa=20agreement)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../valuation/monte_carlo.py | 23 +++- tests/test_edge_cases.py | 121 ++++++++++++++++++ 2 files changed, 140 insertions(+), 4 deletions(-) diff --git a/src/derivatives_pricing/valuation/monte_carlo.py b/src/derivatives_pricing/valuation/monte_carlo.py index eba3626..a7a3e2a 100644 --- a/src/derivatives_pricing/valuation/monte_carlo.py +++ b/src/derivatives_pricing/valuation/monte_carlo.py @@ -435,16 +435,28 @@ def vega_pathwise(self) -> float: # --- likelihood-ratio Greeks ------------------------------------------ + @staticmethod + def _require_positive_sigma_for_lr(sigma: float, greek: str) -> None: + """LR estimators divide by σ; the lognormal density is degenerate at σ=0.""" + if sigma <= 0.0: + raise UnsupportedFeatureError( + f"Likelihood-ratio {greek} is undefined at sigma=0 (the lognormal " + "density degenerates to a Dirac at the forward, so the score " + "function diverges). Use GreekCalculationMethod.NUMERICAL or " + "PATHWISE instead, or call BSM analytical greeks." + ) + def delta_lr(self) -> float: r"""Likelihood-ratio (score-function) delta estimator. :math:`\\Delta = e^{-rT}\\,\\mathbb{E}\\!\\left[\\Phi(S_T)\\, \\frac{Z}{\\sigma\\sqrt{T}\\,S_0}\\right]` """ + sigma = float(self.underlying.volatility) + self._require_positive_sigma_for_lr(sigma, "delta") ST, idx, ttm, df = self._simulate_terminal() Z = self._effective_terminal_z(idx, ttm) S0 = float(self.underlying.initial_value) - sigma = float(self.underlying.volatility) K = self.valuation_ctx.strike payoff = _vanilla_payoff(self.valuation_ctx.option_type, K, ST) return float(np.mean(df * payoff * Z / (sigma * np.sqrt(ttm) * S0))) @@ -462,9 +474,10 @@ def vega_lr(self) -> float: :math:`\\mu = r - q - \\tfrac12\\sigma^2` depends on :math:`\\sigma` via the Itô correction. """ + sigma = float(self.underlying.volatility) + self._require_positive_sigma_for_lr(sigma, "vega") ST, idx, ttm, df = self._simulate_terminal() Z = self._effective_terminal_z(idx, ttm) - sigma = float(self.underlying.volatility) K = self.valuation_ctx.strike payoff = _vanilla_payoff(self.valuation_ctx.option_type, K, ST) score = (Z**2 - 1) / sigma - Z * np.sqrt(ttm) @@ -516,9 +529,10 @@ def theta_lr(self) -> float: where :math:`a = r - q - \tfrac12\sigma^2`. """ + sigma = float(self.underlying.volatility) + self._require_positive_sigma_for_lr(sigma, "theta") ST, idx, ttm, df = self._simulate_terminal() Z = self._effective_terminal_z(idx, ttm) - sigma = float(self.underlying.volatility) K = self.valuation_ctx.strike r, q = self._risk_free_and_div_rates(ttm) a = r - q - 0.5 * sigma**2 @@ -550,9 +564,10 @@ def rho_lr(self) -> float: \rho_{\text{LR}} = \mathrm{df}\,\mathbb{E}\!\left[ \Phi(S_T)\!\left(\frac{\sqrt{T}}{\sigma}\,Z - T\right)\right] """ + sigma = float(self.underlying.volatility) + self._require_positive_sigma_for_lr(sigma, "rho") ST, idx, ttm, df = self._simulate_terminal() Z = self._effective_terminal_z(idx, ttm) - sigma = float(self.underlying.volatility) K = self.valuation_ctx.strike payoff = _vanilla_payoff(self.valuation_ctx.option_type, K, ST) score = np.sqrt(ttm) / sigma * Z - ttm diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py index 12d9c3a..adab554 100644 --- a/tests/test_edge_cases.py +++ b/tests/test_edge_cases.py @@ -10,10 +10,12 @@ import pytest from derivatives_pricing.enums import ( + AsianAveraging, BarrierAction, BarrierDirection, BarrierMonitoring, ExerciseType, + GreekCalculationMethod, OptionType, PDESpaceGrid, PricingMethod, @@ -21,8 +23,14 @@ from derivatives_pricing.exceptions import ( ArbitrageViolationError, NumericalError, + UnsupportedFeatureError, ValidationError, ) +from derivatives_pricing.stochastic_processes import ( + GBMParams, + GBMProcess, + SimulationConfig, +) from helpers import ( flat_curve, market_data, @@ -31,12 +39,15 @@ pv as _pv, ) from derivatives_pricing.valuation import ( + AsianSpec, BarrierSpec, BinomialParams, + MonteCarloParams, VanillaSpec, PDEParams, UnderlyingData, ) +from derivatives_pricing.valuation.core import OptionValuation PRICING_DATE = dt.datetime(2025, 1, 1) MATURITY = dt.datetime(2025, 7, 3) # ~0.5y @@ -718,3 +729,113 @@ def test_uic_tiny_vol_not_hit_returns_zero(self, vol): pv = _pv(ud, spec, PricingMethod.BSM) assert np.isfinite(pv) assert pv == 0.0 + + +# ═══════════════════════════════════════════════════════════════════════ +# Asian options: σ = 0 explicitly rejected by the analytical formula +# ═══════════════════════════════════════════════════════════════════════ + + +class TestAsianZeroVolatility: + """Geometric/arithmetic Asian closed-forms reject σ = 0 explicitly.""" + + @pytest.mark.parametrize( + "averaging", + [AsianAveraging.GEOMETRIC, AsianAveraging.ARITHMETIC], + ) + def test_asian_analytical_rejects_zero_vol(self, averaging): + ud = _underlying(vol=0.0) + spec = AsianSpec( + averaging=averaging, + option_type=OptionType.CALL, + exercise_type=ExerciseType.EUROPEAN, + strike=100.0, + maturity=MATURITY, + currency="USD", + num_observations=12, + ) + with pytest.raises(ValidationError, match=r"volatility must be positive"): + _pv(ud, spec, PricingMethod.BSM) + + +# ═══════════════════════════════════════════════════════════════════════ +# MC likelihood-ratio greeks: σ = 0 rejected (score function diverges) +# ═══════════════════════════════════════════════════════════════════════ + + +class TestMCLikelihoodRatioGreeksZeroVolatility: + """LR estimators divide by σ; the lognormal density is degenerate at σ=0.""" + + @pytest.fixture + def _ov_zero_vol(self): + md = market_data( + pricing_date=PRICING_DATE, + discount_curve=_DEFAULT_RATE_CURVE, + currency="USD", + ) + sim = SimulationConfig(paths=2_000, num_steps=12, end_date=MATURITY) + gbm = GBMProcess(md, GBMParams(initial_value=100.0, volatility=0.0), sim) + spec = make_vanilla_spec( + strike=100.0, + maturity=MATURITY, + option_type=OptionType.CALL, + currency="USD", + ) + return OptionValuation( + gbm, + spec, + PricingMethod.MONTE_CARLO, + params=MonteCarloParams(random_seed=42), + ) + + @pytest.mark.parametrize("greek", ["delta", "vega", "theta", "rho"]) + def test_lr_greek_rejects_zero_vol(self, _ov_zero_vol, greek): + with pytest.raises(UnsupportedFeatureError, match=r"Likelihood-ratio .* sigma=0"): + getattr(_ov_zero_vol, greek)(greek_calc_method=GreekCalculationMethod.LIKELIHOOD_RATIO) + + +# ═══════════════════════════════════════════════════════════════════════ +# Cross-engine vol = 0 agreement (BSM vs MC vanilla) +# ═══════════════════════════════════════════════════════════════════════ + + +class TestVanillaZeroVolCrossEngine: + """At σ = 0 BSM and MC must both return discounted forward intrinsic.""" + + @pytest.mark.parametrize( + "option_type,strike", + [ + (OptionType.CALL, 90.0), # ITM call + (OptionType.PUT, 110.0), # ITM put + (OptionType.CALL, 110.0), # OTM call + (OptionType.PUT, 90.0), # OTM put + ], + ) + def test_bsm_and_mc_agree_at_zero_vol(self, option_type, strike): + md = market_data( + pricing_date=PRICING_DATE, + discount_curve=_DEFAULT_RATE_CURVE, + currency="USD", + ) + ud = UnderlyingData(initial_value=100.0, volatility=0.0, market_data=md) + sim = SimulationConfig(paths=1_000, num_steps=12, end_date=MATURITY) + gbm = GBMProcess(md, GBMParams(initial_value=100.0, volatility=0.0), sim) + spec = make_vanilla_spec( + strike=strike, + maturity=MATURITY, + option_type=option_type, + currency="USD", + ) + + bsm_pv = _pv(ud, spec, PricingMethod.BSM) + mc_pv = OptionValuation( + gbm, + spec, + PricingMethod.MONTE_CARLO, + params=MonteCarloParams(random_seed=42), + ).present_value() + + assert np.isfinite(bsm_pv) and np.isfinite(mc_pv) + # All paths identical at σ=0 → MC variance is exactly zero, so the + # two engines should agree to machine precision (modulo float rounding). + assert np.isclose(bsm_pv, mc_pv, rtol=1e-10, atol=1e-10) From a7665860c4d2611cf91daf9538a0f886e9d44697 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Mon, 27 Apr 2026 16:16:16 +0100 Subject: [PATCH 05/12] Add Binomial barrier stencil guard for greeks (delta, gamma) Co-authored-by: Copilot --- src/derivatives_pricing/valuation/binomial.py | 63 ++++++++++ tests/test_barrier.py | 109 ++++++++++++++++++ 2 files changed, 172 insertions(+) diff --git a/src/derivatives_pricing/valuation/binomial.py b/src/derivatives_pricing/valuation/binomial.py index 8e28c6c..b77f3cd 100644 --- a/src/derivatives_pricing/valuation/binomial.py +++ b/src/derivatives_pricing/valuation/binomial.py @@ -1054,6 +1054,69 @@ def _tree_greeks_data(self) -> tuple[np.ndarray, np.ndarray, float]: dt = T / num_steps return option_lattice, spot_lattice, dt + def _stencil_straddles_barrier(self, step: int) -> bool: + """Return True if any node at ``step`` lies on the wrong side of the barrier. + + The Hull tree-greek formulas assume the option value is locally smooth + in spot. When a node within the extraction stencil sits at or beyond + the barrier, the absorbing-boundary discontinuity (KO → 0 or rebate, + or KI → vanilla) injects a step-function jump that the central-difference + formulas misread as huge curvature. This check is the precondition + for a sensible delta (step=1) or gamma (step=2) extraction. + + Triggered-at-inception barriers are exempt: a triggered KO collapses + to a lattice that is independent of spot (no discontinuity left), and a triggered + KI prices as the vanilla underlying (also smooth). Discretely + monitored barriers are also exempt since the absorbing transition only + happens at monitoring nodes, which don't coincide with intermediate + tree steps in general. + """ + if self.spec.monitoring is not BarrierMonitoring.CONTINUOUS: + return False + if self.valuation_ctx._barrier_triggered_at_inception(): + return False + _, _, spot_lattice = self._setup_binomial_parameters() + nodes = spot_lattice[: step + 1, step] + if self.spec.direction is BarrierDirection.UP: + return bool(np.any(nodes >= self.spec.barrier)) + return bool(np.any(nodes <= self.spec.barrier)) + + def delta(self) -> float: + """Tree delta with a near-barrier safety guard. + + See :meth:`_BinomialValuationBase.delta` for the formula. Raises + :class:`UnsupportedFeatureError` if the step-1 extraction stencil + already crosses the barrier — at that point Hull's central-difference + formula is reading across the absorbing-boundary discontinuity and + the result is unreliable. Use ``PricingMethod.PDE_FD`` instead. + """ + if self._stencil_straddles_barrier(step=1): + raise UnsupportedFeatureError( + "Binomial tree delta is unreliable when the step-1 extraction " + "stencil straddles the barrier (spot is too close to H, so " + "the up- or down-node lies past the absorbing boundary). " + "Use PricingMethod.PDE_FD for greeks in this regime." + ) + return super().delta() + + def gamma(self) -> float: + """Tree gamma with a near-barrier safety guard. + + See :meth:`_BinomialValuationBase.gamma` for the formula. Raises + :class:`UnsupportedFeatureError` if the step-2 extraction stencil + crosses the barrier — Hull's three-node central difference + misreads the absorbing-boundary jump as huge curvature in that + case. Use ``PricingMethod.PDE_FD`` for near-barrier greeks. + """ + if self._stencil_straddles_barrier(step=2): + raise UnsupportedFeatureError( + "Binomial tree gamma is unreliable when the step-2 extraction " + "stencil straddles the barrier (one of the uu/dd nodes lies " + "past the absorbing boundary). Use PricingMethod.PDE_FD for " + "greeks in this regime." + ) + return super().gamma() + def _ko_rebate_values( self, *, diff --git a/tests/test_barrier.py b/tests/test_barrier.py index 444c0af..b5554b5 100644 --- a/tests/test_barrier.py +++ b/tests/test_barrier.py @@ -2312,6 +2312,115 @@ def test_auto_select_rho_uses_numerical(self): assert np.isfinite(rho) +class TestBinomialBarrierStencilGuard: + """Tree delta/gamma must reject when the extraction stencil straddles + the barrier — Hull's central-difference formula is unreliable in that + regime (the absorbing-boundary discontinuity injects false curvature). + """ + + @pytest.fixture(autouse=True) + def _setup(self): + # Spot 72.4 with UP barrier 73 → step-2 up-up node lands ~exactly + # on the barrier under Boyle-Lau alignment. American put deep ITM + # so early-exercise locks pv = K - S = 7.6, true delta = -1, gamma 0. + curve = DiscountCurve.flat(0.05, end_time=2.0) + self.md = MarketData( + PRICING_DATE, + curve, + currency=CURRENCY, + day_count_convention=DayCountConvention.ACT_365F, + ) + self.ud = UnderlyingData( + initial_value=72.4, + volatility=0.205, + market_data=self.md, + ) + self.am_spec = _barrier_spec( + option_type=OptionType.PUT, + exercise_type=ExerciseType.AMERICAN, + strike=80.0, + barrier=73.0, + direction=BarrierDirection.UP, + action=BarrierAction.OUT, + ) + self.ov = OptionValuation( + self.ud, + self.am_spec, + PricingMethod.BINOMIAL, + ) + + def test_pv_and_theta_still_work(self): + """The guard targets greek extraction only — pv/delta/theta are unaffected.""" + assert np.isclose(self.ov.present_value(), 7.6, atol=1e-6) + assert np.isclose(self.ov.delta(), -1.0, atol=1e-6) + # Deep ITM American put with no time value → theta = 0 + assert np.isclose(self.ov.theta(), 0.0, atol=1e-6) + + def test_gamma_rejected_when_stencil_straddles(self): + with pytest.raises(UnsupportedFeatureError, match=r"step-2 .* straddles"): + self.ov.gamma() + + def test_delta_gamma_far_from_barrier_works(self): + """Sanity: when barrier is far from spot, both delta and gamma extract OK.""" + spec = _barrier_spec( + option_type=OptionType.PUT, + exercise_type=ExerciseType.AMERICAN, + strike=80.0, + barrier=120.0, # well above spot 72.4 + direction=BarrierDirection.UP, + action=BarrierAction.OUT, + ) + ov = OptionValuation(self.ud, spec, PricingMethod.BINOMIAL) + # Both should return finite values (no guard trigger) + assert np.isfinite(ov.delta()) + assert np.isfinite(ov.gamma()) + + def test_ko_triggered_at_inception_skips_guard(self): + """KO triggered at inception has no discontinuity around spot → + guard is bypassed, greeks do not raise. + """ + spec = _barrier_spec( + option_type=OptionType.PUT, + exercise_type=ExerciseType.AMERICAN, + strike=80.0, + barrier=72.4, # H == spot → KO triggered at inception + direction=BarrierDirection.UP, + action=BarrierAction.OUT, + ) + ov = OptionValuation(self.ud, spec, PricingMethod.BINOMIAL) + # no rebate so pv and greeks are zero + assert ov.present_value() == 0.0 + assert ov.delta() == 0.0 + assert ov.gamma() == 0.0 + assert ov.theta() == 0.0 + + def test_ki_triggered_at_inception_matches_vanilla(self): + """KI triggered at inception → option becomes vanilla; greeks should + match vanilla, not raise.""" + spec = _barrier_spec( + option_type=OptionType.PUT, + exercise_type=ExerciseType.AMERICAN, + strike=80.0, + barrier=72.4, # H == spot → KI triggered (option activated) + direction=BarrierDirection.UP, + action=BarrierAction.IN, + ) + ov = OptionValuation(self.ud, spec, PricingMethod.BINOMIAL) + van = OptionValuation( + self.ud, + VanillaSpec( + option_type=OptionType.PUT, + exercise_type=ExerciseType.AMERICAN, + strike=80.0, + maturity=MATURITY, + ), + PricingMethod.BINOMIAL, + ) + assert np.isclose(ov.present_value(), van.present_value(), rtol=1e-3) + assert np.isclose(ov.delta(), van.delta(), rtol=1e-3) + assert np.isclose(ov.gamma(), van.gamma(), rtol=1e-3) + + # =========================================================================== # Binomial barrier coverage # =========================================================================== From 5dc115f691e941d0626c47286573f6f3f3e60c25 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Mon, 27 Apr 2026 16:53:29 +0100 Subject: [PATCH 06/12] (test_barrier): replace _BINOMIAL_SKIP_SPOTS with engine-driven skip in BT98 Table6 greeks --- tests/test_barrier.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/tests/test_barrier.py b/tests/test_barrier.py index b5554b5..cefca16 100644 --- a/tests/test_barrier.py +++ b/tests/test_barrier.py @@ -1999,11 +1999,6 @@ class TestBarrierGreeksAgainstBoyleTianTable6: 90.2: (1.2869, -0.0451, -0.1161), } - # Binomial tree greeks are known to degrade very close to the barrier - # (Boyle-Lau retopology + discrete-grid noise). Regression is enforced - # at those spots via BSM numerical and PDE grid greeks only. - _BINOMIAL_SKIP_SPOTS = {90.5, 90.4, 90.3, 90.2} - # Per-engine tolerances for each greek. _TOLS = { "delta": { @@ -2078,10 +2073,22 @@ def test_down_and_out_call_greek_matches_paper(self, spot: float, greek: str): engine_values: dict[PricingMethod, float | None] = {} for method in (PricingMethod.BSM, PricingMethod.BINOMIAL, PricingMethod.PDE_FD): - if method is PricingMethod.BINOMIAL and spot in self._BINOMIAL_SKIP_SPOTS: + try: + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always", RuntimeWarning) + val = self._engine_greek(spot, method, greek) + # Boyle-Lau alignment requires more steps than the cap → engine + # explicitly flagged the result as O(1/√n)-degraded. Skip + # rather than assert paper-truth on a value the engine itself + # said not to trust. + if any("Boyle-Lau step alignment" in str(w.message) for w in caught): + engine_values[method] = None + else: + engine_values[method] = val + except UnsupportedFeatureError: + # The engine itself declared this case unsuitable (e.g. binomial + # tree-greek stencil straddles the barrier near H). engine_values[method] = None - continue - engine_values[method] = self._engine_greek(spot, method, greek) def _fmt(v: float | None) -> str: return f"{v:.6f}" if v is not None else "skipped" From fc7af7099887593ac070b8d8397ad4d7167859dd Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Tue, 28 Apr 2026 15:45:43 +0100 Subject: [PATCH 07/12] =?UTF-8?q?Fix=20barrier=20inception-triggered=20gre?= =?UTF-8?q?eks=20via=20OV-level=20NUMERICAL=20short-circuit=20and=20BSM=20?= =?UTF-8?q?theta=20closed-form=20-=20OV.delta/gamma=20now=20short-circuit?= =?UTF-8?q?=20to=200=20(KO)=20or=20delegate=20the=20bump=20to=20a=20cached?= =?UTF-8?q?=20vanilla=20equivalent=20(KI)=20for=20inception-triggered=20ba?= =?UTF-8?q?rriers.=20-=20BSM=20theta=20short-circuits=20at=20engine=20leve?= =?UTF-8?q?l=20with=20a=20closed-form=20(0=20for=20AT=5FHIT/no-rebate,=20+?= =?UTF-8?q?r=C2=B7V/365=20for=20AT=5FEXPIRY)=20since=20the=20PDE=20identit?= =?UTF-8?q?y=20doesn't=20hold=20for=20at-inception=20triggered=20KOs=20(it?= =?UTF-8?q?=20does=20hold=20for=20KIs=20but=20we=20route=20to=20vanilla.th?= =?UTF-8?q?eta()=20for=20precision)=20-=20Lift=20=5Fvanilla=5Fequivalent?= =?UTF-8?q?=5Fvaluation=20to=20OV=20(memoised)=20so=20PDE's=20existing=20u?= =?UTF-8?q?sages=20share=20the=20cached=20instance.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../valuation/barrier_analytical.py | 37 +++++++++++-- src/derivatives_pricing/valuation/core.py | 52 +++++++++++++++++++ src/derivatives_pricing/valuation/pde.py | 24 ++------- 3 files changed, 89 insertions(+), 24 deletions(-) diff --git a/src/derivatives_pricing/valuation/barrier_analytical.py b/src/derivatives_pricing/valuation/barrier_analytical.py index 811fe45..a6d93c7 100644 --- a/src/derivatives_pricing/valuation/barrier_analytical.py +++ b/src/derivatives_pricing/valuation/barrier_analytical.py @@ -476,16 +476,45 @@ def theta(self) -> float: ``V`` is the closed-form barrier price; ``Δ`` and ``Γ`` come from central-difference bump-and-revalue around the same closed-form evaluator (routed through :attr:`valuation_ctx` so repeated calls - hit the OV-level cache). The identity is exact in the - continuation region (triggered-at-inception cases already - short-circuit in :meth:`present_value`) and delivers better - accuracy than a naive forward-difference time bump. + hit the OV-level cache). The identity holds in the continuation + region; for inception-triggered KOs we short-circuit to a + closed-form θ because the contract is no longer PDE-governed + (it's just paid cash, or a deterministic discounted payment). + For inception-triggered KIs we delegate to the vanilla equivalent. Returned per **calendar day** """ ctx = self.valuation_ctx underlying = self.underlying + # ── Inception-triggered short-circuit ───────────────────────── + # KO triggered: the contract has collapsed to a deterministic + # cashflow and is no longer PDE-governed. The identity then gives + # the wrong answer — e.g. r·V instead of 0 for an AT_HIT rebate + # — so we return the closed-form θ directly. + # KI triggered: the contract is the underlying vanilla, which DOES + # satisfy the BS PDE; the identity would still hold here (with the + # vanilla's δ and γ that OV's NUMERICAL short-circuit already + # provides). We delegate to ``vanilla.theta()`` anyway as a + # precision upgrade — it returns the analytical θ rather than the + # identity evaluated with bumped greeks. + if ctx._barrier_triggered_at_inception(): + spec = self.spec + if spec.action is BarrierAction.IN: + return float(ctx._vanilla_equivalent_valuation().theta()) + # KO triggered. + if spec.rebate <= 0.0 or spec.rebate_timing is RebateTiming.AT_HIT: + # No rebate or rebate paid immediately → pv has no time + # evolution → θ = 0. + return 0.0 + # AT_EXPIRY rebate: pv = R · df_r(T), so dpv/dt = +r · pv; + # per-day θ = r · pv / 365. + T = ctx._maturity_year_fraction() + df_r = float(ctx.discount_curve.df(T)) + pv = float(spec.rebate) * df_r + r = -np.log(df_r) / T + return float(r * pv / 365.0) + S = float(underlying.initial_value) sigma = float(underlying.volatility) sigma2 = sigma**2 diff --git a/src/derivatives_pricing/valuation/core.py b/src/derivatives_pricing/valuation/core.py index 12c9fcc..0dd3711 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -32,6 +32,7 @@ from ..exceptions import ConfigurationError, UnsupportedFeatureError, ValidationError from ..enums import ( AsianAveraging, + BarrierAction, BarrierDirection, BarrierMonitoring, DayCountConvention, @@ -529,6 +530,20 @@ def delta( if method is not GreekCalculationMethod.NUMERICAL: return float(self._impl.delta()) + if isinstance(self._spec, BarrierSpec) and self._barrier_triggered_at_inception(): + # Bump-and-revalue may cross the trigger boundary and price the + # un-triggered contract on the bumped spot — meaningless for an + # already-triggered barrier. Short-circuit: + # - KO triggered → constant-in-spot cashflow → δ = 0. + # - KI triggered → contract IS the underlying vanilla; bump on + # the vanilla equivalent (no state transition). + if self._spec.action is BarrierAction.OUT: + return 0.0 + return self._vanilla_equivalent_valuation().delta( + epsilon=epsilon, + greek_calc_method=GreekCalculationMethod.NUMERICAL, + ) + if epsilon is None: epsilon = self._underlying.initial_value / 100 if isinstance(self._spec, BarrierSpec): @@ -587,6 +602,15 @@ def gamma( if method is not GreekCalculationMethod.NUMERICAL: return float(self._impl.gamma()) + if isinstance(self._spec, BarrierSpec) and self._barrier_triggered_at_inception(): + # See the corresponding short-circuit in :meth:`delta`. + if self._spec.action is BarrierAction.OUT: + return 0.0 + return self._vanilla_equivalent_valuation().gamma( + epsilon=epsilon, + greek_calc_method=GreekCalculationMethod.NUMERICAL, + ) + if epsilon is None: epsilon = self._underlying.initial_value / 100 if isinstance(self._spec, BarrierSpec): @@ -1027,6 +1051,34 @@ def _barrier_triggered_at_inception(self) -> bool: assert mon_dates is not None return any(d == self.pricing_date for d in mon_dates) + @_memoize_result + def _vanilla_equivalent_valuation(self) -> OptionValuation: + """Return the vanilla ``OptionValuation`` that a triggered KI collapses to. + + Cached per-instance so repeated greek calls on a triggered KI reuse + the same vanilla — the vanilla's own solve cache then survives + across delta/gamma/theta and matches the "one solve, three free + greeks" cost profile of a non-triggered native greek. + """ + assert isinstance(self._spec, BarrierSpec), ( + "_vanilla_equivalent_valuation called on non-BarrierSpec valuation." + ) + spec = self._spec + vanilla_spec = VanillaSpec( + option_type=spec.option_type, + exercise_type=spec.exercise_type, + strike=spec.strike, + maturity=spec.maturity, + currency=spec.currency, + contract_size=spec.contract_size, + ) + return OptionValuation( + underlying=self._underlying, + spec=vanilla_spec, + pricing_method=self._pricing_method, + params=self._params, + ) + def _apply_control_variate(self, base_pv: float) -> float: """Apply European control-variate adjustment to American base PV. diff --git a/src/derivatives_pricing/valuation/pde.py b/src/derivatives_pricing/valuation/pde.py index e2d956a..4b4f64b 100644 --- a/src/derivatives_pricing/valuation/pde.py +++ b/src/derivatives_pricing/valuation/pde.py @@ -2612,22 +2612,6 @@ def _resolved_knock_out_value(self) -> float | None: ttm = self.valuation_ctx._maturity_year_fraction() return float(self._spec.rebate) * float(self.valuation_ctx.discount_curve.df(ttm)) - def _vanilla_equivalent_valuation(self) -> OptionValuation: - from .core import OptionValuation - - vanilla_spec = VanillaSpec( - option_type=self._spec.option_type, - exercise_type=self._spec.exercise_type, - strike=self._spec.strike, - maturity=self._spec.maturity, - ) - return OptionValuation( - underlying=self.underlying, - spec=vanilla_spec, - pricing_method=self.valuation_ctx.pricing_method, - params=self.valuation_ctx.params, - ) - def _last_dtau(self) -> float: solve_args = self._base_solve_args() time_to_maturity = float(solve_args["time_to_maturity"]) @@ -2807,7 +2791,7 @@ def delta(self) -> float: if self.valuation_ctx._barrier_triggered_at_inception(): if spec.action is BarrierAction.OUT: return 0.0 - return self._vanilla_equivalent_valuation().delta() + return self.valuation_ctx._vanilla_equivalent_valuation().delta() if self._is_european_ki(): ko_result, van_result = self._solve_european_ki_components() @@ -2824,7 +2808,7 @@ def gamma(self) -> float: if self.valuation_ctx._barrier_triggered_at_inception(): if spec.action is BarrierAction.OUT: return 0.0 - return self._vanilla_equivalent_valuation().gamma() + return self.valuation_ctx._vanilla_equivalent_valuation().gamma() if self._is_european_ki(): ko_result, van_result = self._solve_european_ki_components() @@ -2840,7 +2824,7 @@ def theta(self) -> float: if self.valuation_ctx._barrier_triggered_at_inception(): if spec.action is BarrierAction.OUT: return self._resolved_knock_out_theta() - return self._vanilla_equivalent_valuation().theta() + return self.valuation_ctx._vanilla_equivalent_valuation().theta() if self._is_european_ki(): ko_result, van_result = self._solve_european_ki_components() @@ -2918,7 +2902,7 @@ def present_value(self) -> float: if triggered_value is None: raise ConfigurationError("Resolved knock-out state unexpectedly unavailable") return triggered_value - return self._vanilla_equivalent_valuation().present_value() + return self.valuation_ctx._vanilla_equivalent_valuation().present_value() spec = self._spec label = f"PDE barrier {'American' if spec.exercise_type is ExerciseType.AMERICAN else 'European'}" with log_timing(logger, f"{label} present_value", self.pde_params.log_timings): From d67e8eaebba51e943cd4bc20b9f320dfe271da87 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Wed, 29 Apr 2026 14:35:12 +0100 Subject: [PATCH 08/12] Add mandatory kwarg monitoring to for_barriers cls method - Discrete barriers require finer spatial resolution due to repeated discontinuities at observation dates; default to spot_steps=2400, time_steps=1600 with IMPLICIT (L-stable) scheme to avoid CN ringing. - Continuous monitoring keeps existing CN/1200/800 defaults. - Boyle-Tian Table 8 PDE_FD tolerance loosened to reflect IMPLICIT vs CN paper-truth gap. Co-authored-by: Copilot --- src/derivatives_pricing/valuation/core.py | 4 +- src/derivatives_pricing/valuation/params.py | 49 ++++++++++++++++++--- tests/test_barrier.py | 6 ++- tests/test_params.py | 36 ++++++++++++--- tests/test_pde.py | 1 + tests/test_quantlib_greeks_comparison.py | 3 +- 6 files changed, 83 insertions(+), 16 deletions(-) diff --git a/src/derivatives_pricing/valuation/core.py b/src/derivatives_pricing/valuation/core.py index 0dd3711..2ebdf85 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -959,7 +959,9 @@ def _resolve_params( return BinomialParams(num_steps=num_steps) return BinomialParams() if pricing_method is PricingMethod.PDE_FD: - return PDEParams.for_barriers() if isinstance(spec, BarrierSpec) else PDEParams() + if isinstance(spec, BarrierSpec): + return PDEParams.for_barriers(monitoring=spec.monitoring) + return PDEParams() return None if pricing_method is PricingMethod.MONTE_CARLO: diff --git a/src/derivatives_pricing/valuation/params.py b/src/derivatives_pricing/valuation/params.py index 8714c5d..c327eb2 100644 --- a/src/derivatives_pricing/valuation/params.py +++ b/src/derivatives_pricing/valuation/params.py @@ -10,7 +10,7 @@ from typing import Any import warnings -from ..enums import PDEEarlyExercise, PDEMethod, PDESpaceGrid +from ..enums import BarrierMonitoring, PDEEarlyExercise, PDEMethod, PDESpaceGrid from ..exceptions import ValidationError @@ -237,15 +237,50 @@ class PDEParams: log_timings: bool = False @classmethod - def for_barriers(cls, **overrides: Any) -> PDEParams: - """Create params that mirror the library's internal barrier defaults. + def for_barriers( + cls, + *, + monitoring: BarrierMonitoring, + **overrides: Any, + ) -> PDEParams: + """Create PDE params tuned for barrier pricing. Returns a ``PDEParams`` instance with a finer grid and log-spot - spatial discretization suitable for barrier pricing. Any keyword - argument accepted by the constructor can be passed to override - individual fields. + spatial discretization suitable for barrier pricing. The time- + marching ``method`` is chosen from ``monitoring``: + + - ``BarrierMonitoring.CONTINUOUS`` → ``PDEMethod.CRANK_NICOLSON``. + Continuous monitoring has a single payoff discontinuity at + maturity; the default ``rannacher_steps=2`` startup dampens it + and CN's higher-order time accuracy gives the best PV/greek + quality on the rest of the time march. + - ``BarrierMonitoring.DISCRETE`` → ``PDEMethod.IMPLICIT``. + Discrete monitoring projects ``V(S, t_i) = 0`` past the + barrier at every observation date, introducing a fresh step + discontinuity each time. CN is only A-stable + (Pooley-Forsyth-Vetzal, 2003). ``IMPLICIT`` is L-stable and + empirically performs slightly better than CN at the same grid + resolution for discrete barriers, so it's the default here. + + ``monitoring`` is required (no default) so the dependency is + explicit at the call site. Any other keyword argument accepted + by the constructor can be passed to override individual fields, + including ``method``. """ - defaults = cls(spot_steps=1200, time_steps=800, space_grid=PDESpaceGrid.LOG_SPOT) + if monitoring is BarrierMonitoring.DISCRETE: + method = PDEMethod.IMPLICIT + spot_steps, time_steps = 2400, 1600 + else: + method = PDEMethod.CRANK_NICOLSON + spot_steps, time_steps = 1200, 800 + + defaults = cls( + spot_steps=spot_steps, + time_steps=time_steps, + space_grid=PDESpaceGrid.LOG_SPOT, + method=method, + ) + return dc_replace(defaults, **overrides) if overrides else defaults def __post_init__(self) -> None: diff --git a/tests/test_barrier.py b/tests/test_barrier.py index cefca16..9cc969f 100644 --- a/tests/test_barrier.py +++ b/tests/test_barrier.py @@ -1207,9 +1207,9 @@ def _make_spec(cls, monitoring_kind: str | int) -> BarrierSpec: _TOLS: dict[PricingMethod, dict[str, float]] = { PricingMethod.BSM: dict(rtol=0.022, atol=1.0e-4), PricingMethod.BINOMIAL: dict(rtol=0.013, atol=1.0e-4), - PricingMethod.PDE_FD: dict(rtol=0.001, atol=1.0e-4), + PricingMethod.PDE_FD: dict(rtol=0.002, atol=1.0e-4), } - _PDE_FD_HOURLY_TOL: dict[str, float] = dict(rtol=0.015, atol=1.0e-4) + _PDE_FD_HOURLY_TOL: dict[str, float] = dict(rtol=0.020, atol=1.0e-4) @pytest.mark.parametrize( "frequency,monitoring_kind,paper_pv", @@ -1340,6 +1340,8 @@ class TestBarrierPresentValueAgainstBroadieGlasserman: _TOLS: dict[PricingMethod, dict[str, float]] = { PricingMethod.BSM: dict(rtol=0.0, atol=1.5e-3), PricingMethod.BINOMIAL: dict(rtol=0.035, atol=1.0e-3), + # PDE_FD discrete-monitoring default is IMPLICIT (L-stable); + # slightly less accurate vs paper than CN. PricingMethod.PDE_FD: dict(rtol=0.002, atol=1.0e-3), } diff --git a/tests/test_params.py b/tests/test_params.py index 52759c3..95306ff 100644 --- a/tests/test_params.py +++ b/tests/test_params.py @@ -4,7 +4,7 @@ import pytest -from derivatives_pricing.enums import PDEMethod, PDESpaceGrid +from derivatives_pricing.enums import BarrierMonitoring, PDEMethod, PDESpaceGrid from derivatives_pricing.exceptions import ValidationError from derivatives_pricing.valuation.params import BinomialParams, MonteCarloParams, PDEParams @@ -190,21 +190,47 @@ def test_rejects_bool_for_spot_steps(self): with pytest.raises(ValidationError, match="spot_steps must be an int"): PDEParams(spot_steps=True) - def test_for_barriers_defaults(self): - p = PDEParams.for_barriers() + def test_for_barriers_continuous_defaults(self): + p = PDEParams.for_barriers(monitoring=BarrierMonitoring.CONTINUOUS) assert p.spot_steps == 1200 assert p.time_steps == 800 assert p.space_grid is PDESpaceGrid.LOG_SPOT + # Continuous monitoring → CN (high accuracy on smooth time march + # after rannacher-damped startup at maturity). assert p.method is PDEMethod.CRANK_NICOLSON assert p.control_variate_european is False + def test_for_barriers_discrete_defaults_to_implicit(self): + p = PDEParams.for_barriers(monitoring=BarrierMonitoring.DISCRETE) + assert p.spot_steps == 2400 + assert p.time_steps == 1600 + assert p.space_grid is PDESpaceGrid.LOG_SPOT + # Discrete monitoring → IMPLICIT (L-stable) + assert p.method is PDEMethod.IMPLICIT + def test_for_barriers_with_overrides(self): - p = PDEParams.for_barriers(log_timings=True, control_variate_european=True) + p = PDEParams.for_barriers( + monitoring=BarrierMonitoring.CONTINUOUS, + log_timings=True, + control_variate_european=True, + ) assert p.spot_steps == 1200 assert p.time_steps == 800 assert p.log_timings is True assert p.control_variate_european is True + def test_for_barriers_method_override_wins(self): + """Caller can force CN even on discrete if they really want it.""" + p = PDEParams.for_barriers( + monitoring=BarrierMonitoring.DISCRETE, + method=PDEMethod.CRANK_NICOLSON, + ) + assert p.method is PDEMethod.CRANK_NICOLSON + def test_for_barriers_rejects_invalid_override(self): with pytest.raises(ValidationError, match="spot_steps must be >= 3"): - PDEParams.for_barriers(spot_steps=2) + PDEParams.for_barriers(monitoring=BarrierMonitoring.CONTINUOUS, spot_steps=2) + + def test_for_barriers_requires_monitoring_kwarg(self): + with pytest.raises(TypeError): + PDEParams.for_barriers() # type: ignore[call-arg] diff --git a/tests/test_pde.py b/tests/test_pde.py index ed284a6..35b3bdd 100644 --- a/tests/test_pde.py +++ b/tests/test_pde.py @@ -708,6 +708,7 @@ def test_pde_fd_barrier_european_ki_facade_vs_direct_core_greeks( rebate_timing=RebateTiming.AT_EXPIRY, ) params = PDEParams.for_barriers( + monitoring=BarrierMonitoring.CONTINUOUS, spot_steps=800, time_steps=800, ) diff --git a/tests/test_quantlib_greeks_comparison.py b/tests/test_quantlib_greeks_comparison.py index 15a0cf9..431fac1 100644 --- a/tests/test_quantlib_greeks_comparison.py +++ b/tests/test_quantlib_greeks_comparison.py @@ -21,6 +21,7 @@ AsianAveraging, BarrierAction, BarrierDirection, + BarrierMonitoring, DayCountConvention, ExerciseType, OptionType, @@ -913,7 +914,7 @@ def test_asian_mc_greeks_vs_quantlib( _BARRIER_RATE = 0.05 _BARRIER_DIV = 0.02 _BARRIER_NUMERICAL_SPOT_BUMP_RATIO = 0.025 -_BARRIER_PDE_CFG = PDEParams.for_barriers() +_BARRIER_PDE_CFG = PDEParams.for_barriers(monitoring=BarrierMonitoring.CONTINUOUS) _BARRIER_BINOM_CFG = BinomialParams(num_steps=1000) _QL_BARRIER_TYPE = { From 6d1eb1918fbad63fe67a8b396d5a066247e78390 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Wed, 29 Apr 2026 16:01:54 +0100 Subject: [PATCH 09/12] (pde); Fix American discrete KO handling - For discrete KOs, payoff is the terminal payoff after the KO reset at maturity. The engine previously reused that same array as intrinsic for the American PSOR constraint at every step. - We now assign intrinsic = payoff.copy before reset --- src/derivatives_pricing/valuation/pde.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/derivatives_pricing/valuation/pde.py b/src/derivatives_pricing/valuation/pde.py index 4b4f64b..e576887 100644 --- a/src/derivatives_pricing/valuation/pde.py +++ b/src/derivatives_pricing/valuation/pde.py @@ -1852,6 +1852,16 @@ def _fd_barrier_ko_core( else: payoff = np.maximum(S - strike, 0.0) + # American exercise intrinsic is the holder's exercise value at any + # in-life moment, which is the vanilla payoff — independent of any + # KO-zone modification applied to the maturity payoff below. Under + # discrete monitoring the holder can exercise between observations + # even when sitting in the KO zone, so the PSOR floor must use the + # unmodified vanilla intrinsic. (For continuous KO the grid is + # truncated at the barrier and the KO zone is absent from the grid, + # so the same vanilla array is also correct on the alive side.) + intrinsic = payoff.copy() if early_exercise else None + # For continuous KO: payoff is zero on the barrier side (enforced # by grid truncation since the barrier is at the boundary). # For discrete KO: zero out the payoff on the knocked-out side at maturity @@ -1866,7 +1876,6 @@ def _fd_barrier_ko_core( payoff[S >= barrier] = 0.0 V = payoff.copy() - intrinsic = payoff if early_exercise else None # ── Dividend schedule ───────────────────────────────────────────── schedule = dividend_schedule or [] From c17d64093cceef38b81567532780d2cc7226038d Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sat, 2 May 2026 14:47:49 +0100 Subject: [PATCH 10/12] Block MC barrier gamma and discrete-monitoring theta - Set default barrier theta time bump to 7 days (empirically more accurate than 1 day) --- src/derivatives_pricing/valuation/core.py | 94 ++++++++++++++++++++++- 1 file changed, 90 insertions(+), 4 deletions(-) diff --git a/src/derivatives_pricing/valuation/core.py b/src/derivatives_pricing/valuation/core.py index 2ebdf85..d13ef84 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -611,6 +611,8 @@ def gamma( greek_calc_method=GreekCalculationMethod.NUMERICAL, ) + self._reject_barrier_mc_numerical(method, greek="gamma") + if epsilon is None: epsilon = self._underlying.initial_value / 100 if isinstance(self._spec, BarrierSpec): @@ -724,8 +726,17 @@ def theta( if method is not GreekCalculationMethod.NUMERICAL: return float(self._impl.theta()) + self._reject_barrier_mc_numerical( + method, greek="theta", monitoring_constraint=BarrierMonitoring.DISCRETE + ) + if time_bump_days is None: - time_bump_days = 1.0 + # Barriers use a longer default (7d) — see + # ``_BARRIER_THETA_TIME_BUMP_DAYS`` for rationale. All other + # specs keep the historical 1d default. + time_bump_days = ( + self._BARRIER_THETA_TIME_BUMP_DAYS if isinstance(self._spec, BarrierSpec) else 1.0 + ) bumped_date = self.pricing_date + dt.timedelta(days=time_bump_days) if bumped_date >= self.maturity: return 0.0 @@ -1355,9 +1366,9 @@ def _reject_barrier_binomial_numerical( ) -> None: """Block NUMERICAL bump-and-revalue greeks on binomial barrier specs. - Bumping spot, volatility or time for a barrier option re-invokes - ``_resolve_effective_num_steps`` on each bumped valuation, and the - Boyle-Lau barrier-alignment formula + For continuous barriers, bumping spot, volatility or time for a + barrier option re-invokes ``_resolve_effective_num_steps`` on each + bumped valuation, and the Boyle-Lau barrier-alignment formula ``candidate = i² σ² T / log(H/S)²`` depends on every one of those inputs. The bumped trees therefore end up with *different* step counts from the center tree, so a central difference is comparing @@ -1365,6 +1376,15 @@ def _reject_barrier_binomial_numerical( Rho is exempt because the risk-free rate does not enter the Boyle-Lau formula, so rate bumps reuse the same tree topology and the finite difference is well-defined. + + For discrete barriers, bumping does not amend _effective_num_steps + but empirically, the greeks are noisy. We thus conservatively block + NUMERICAL greeks on barrier-binomial specs (except rho which is + empirically stable). + + For delta, gamma, theta, users should use native TREE Greeks + For rho, bump and revalue is permitted. + For vega, users are advised to switch pricing method to PDE_FD. """ if allow: return @@ -1385,6 +1405,64 @@ def _reject_barrier_binomial_numerical( "bump-and-revalue on the full grid." ) + # Greek-specific rationales for the MC barrier NUMERICAL block. Kept + # alongside the helper so the "what's blocked and why" is in one + # discoverable place rather than spread across gamma()/theta(). + _MC_BARRIER_NUMERICAL_REASON: dict[str, str] = { + "gamma": ( + "with practical path counts the central-difference noise " + "dominates the |Γ| signal and sign flips occur near zero" + ), + "theta": ( + "bump-and-revalue MC theta on discretely-monitored barriers is " + "empirically too noisy and inaccurate to be trusted " + "(continuous-monitoring MC theta is supported)" + ), + } + + def _reject_barrier_mc_numerical( + self, + method: GreekCalculationMethod, + *, + greek: str, + monitoring_constraint: BarrierMonitoring | None = None, + ) -> None: + """Block NUMERICAL bump-and-revalue greeks on MC barrier specs that + are empirically unreliable. + + Two distinct greek failure modes are covered (see + ``_MC_BARRIER_NUMERICAL_REASON``): + + - Gamma (any monitoring): second-derivative MC noise scales as + ~stderr/ε² and at practical path counts the noise floor exceeds + the |Γ| signal across most of parameter space. + - Theta under DISCRETE monitoring: bump-and-revalue MC theta on + discretely-monitored barriers is empirically too noisy and + inaccurate to be trusted. Continuous-monitoring MC theta is + fine and remains supported — callers requesting the + discrete-only block pass + ``monitoring_constraint=BarrierMonitoring.DISCRETE``. + + Other greeks (delta, vega, rho) are not blocked. + """ + if method is not GreekCalculationMethod.NUMERICAL: + return + if not isinstance(self._spec, BarrierSpec): + return + if self._pricing_method is not PricingMethod.MONTE_CARLO: + return + if monitoring_constraint is not None and self._spec.monitoring is not monitoring_constraint: + return + reason = self._MC_BARRIER_NUMERICAL_REASON.get( + greek, + "the bump-and-revalue path is unreliable for this combination", + ) + raise UnsupportedFeatureError( + f"Numerical {greek} is not supported for Monte Carlo barrier " + f"valuations: {reason}. Use PricingMethod.PDE_FD or " + f"PricingMethod.BINOMIAL for accurate barrier {greek}." + ) + def _auto_select_greek_method( self, *, @@ -1441,6 +1519,14 @@ def _validate_mc_greek_method( # at most halfway to the barrier. _BARRIER_BUMP_MAX_FRACTION: float = 0.5 + # Default time-bump (calendar days) for numerical barrier theta. + # 1d (the global default) is too small for MC (CRN noise dominates the + # finite difference) and gives questionable results for FD/binomial too. + # 7d performs well empirically: large enough to dampen MC + # noise, small enough to skip few barrier monitoring dates and stay close + # to instantaneous theta. + _BARRIER_THETA_TIME_BUMP_DAYS: float = 7.0 + @staticmethod def _validate_bump( bump_value: float | None, From 2aa0aa38cf270c96e797d21503b4e65423fbb66c Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sat, 2 May 2026 15:11:25 +0100 Subject: [PATCH 11/12] Refactor _reject_barrier_mc_numerical to _reject_barrier_numerical - Allows for blocking discrete monitoring theta for all pricing methods --- src/derivatives_pricing/valuation/core.py | 69 +++++++++++++---------- 1 file changed, 39 insertions(+), 30 deletions(-) diff --git a/src/derivatives_pricing/valuation/core.py b/src/derivatives_pricing/valuation/core.py index d13ef84..38416c8 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -611,7 +611,9 @@ def gamma( greek_calc_method=GreekCalculationMethod.NUMERICAL, ) - self._reject_barrier_mc_numerical(method, greek="gamma") + self._reject_barrier_numerical( + method, greek="gamma", engine_constraint=PricingMethod.MONTE_CARLO + ) if epsilon is None: epsilon = self._underlying.initial_value / 100 @@ -726,7 +728,7 @@ def theta( if method is not GreekCalculationMethod.NUMERICAL: return float(self._impl.theta()) - self._reject_barrier_mc_numerical( + self._reject_barrier_numerical( method, greek="theta", monitoring_constraint=BarrierMonitoring.DISCRETE ) @@ -1405,62 +1407,69 @@ def _reject_barrier_binomial_numerical( "bump-and-revalue on the full grid." ) - # Greek-specific rationales for the MC barrier NUMERICAL block. Kept + # Greek-specific rationales for the barrier NUMERICAL block. Kept # alongside the helper so the "what's blocked and why" is in one # discoverable place rather than spread across gamma()/theta(). - _MC_BARRIER_NUMERICAL_REASON: dict[str, str] = { + _BARRIER_NUMERICAL_REASON: dict[str, str] = { "gamma": ( "with practical path counts the central-difference noise " "dominates the |Γ| signal and sign flips occur near zero" ), "theta": ( - "bump-and-revalue MC theta on discretely-monitored barriers is " - "empirically too noisy and inaccurate to be trusted " - "(continuous-monitoring MC theta is supported)" + "bump-and-revalue theta on discretely-monitored barriers is " + "unreliable: bumping the pricing date forces re-resolution of " + "the monitoring schedule (a different contract on the bumped " + "side)" ), } - def _reject_barrier_mc_numerical( + def _reject_barrier_numerical( self, method: GreekCalculationMethod, *, greek: str, + engine_constraint: PricingMethod | None = None, monitoring_constraint: BarrierMonitoring | None = None, ) -> None: - """Block NUMERICAL bump-and-revalue greeks on MC barrier specs that - are empirically unreliable. - - Two distinct greek failure modes are covered (see - ``_MC_BARRIER_NUMERICAL_REASON``): - - - Gamma (any monitoring): second-derivative MC noise scales as - ~stderr/ε² and at practical path counts the noise floor exceeds - the |Γ| signal across most of parameter space. - - Theta under DISCRETE monitoring: bump-and-revalue MC theta on - discretely-monitored barriers is empirically too noisy and - inaccurate to be trusted. Continuous-monitoring MC theta is - fine and remains supported — callers requesting the - discrete-only block pass - ``monitoring_constraint=BarrierMonitoring.DISCRETE``. - - Other greeks (delta, vega, rho) are not blocked. + """Block NUMERICAL bump-and-revalue greeks on barrier specs that + are empirically or structurally unreliable. + + Two distinct greek failure modes are currently covered (see + ``_BARRIER_NUMERICAL_REASON``): + + - Gamma on Monte Carlo barriers (any monitoring): second-derivative + MC noise scales as ~stderr/ε² and at practical path counts the + noise floor exceeds the |Γ| signal across most of parameter + space. Caller passes + ``engine_constraint=PricingMethod.MONTE_CARLO``. + - Theta on any discretely-monitored barrier: bumping the pricing + date forces re-resolution of the monitoring schedule, so the + bumped contract is not the same contract — the resulting theta + mixes time decay with a contract-respecification artifact. + Affects every engine equally; caller passes + ``monitoring_constraint=BarrierMonitoring.DISCRETE``. (Binomial + barrier NUMERICAL greeks are already blocked at a finer grain + by ``_reject_barrier_binomial_numerical``.) + + Other greeks (delta, vega, rho) and continuous-monitoring theta + are not blocked. """ if method is not GreekCalculationMethod.NUMERICAL: return if not isinstance(self._spec, BarrierSpec): return - if self._pricing_method is not PricingMethod.MONTE_CARLO: + if engine_constraint is not None and self._pricing_method is not engine_constraint: return if monitoring_constraint is not None and self._spec.monitoring is not monitoring_constraint: return - reason = self._MC_BARRIER_NUMERICAL_REASON.get( + reason = self._BARRIER_NUMERICAL_REASON.get( greek, "the bump-and-revalue path is unreliable for this combination", ) raise UnsupportedFeatureError( - f"Numerical {greek} is not supported for Monte Carlo barrier " - f"valuations: {reason}. Use PricingMethod.PDE_FD or " - f"PricingMethod.BINOMIAL for accurate barrier {greek}." + f"Numerical {greek} is not supported for this barrier " + f"valuation: {reason}. Use PricingMethod.PDE_FD (GRID) or " + f"PricingMethod.BINOMIAL (TREE) for accurate barrier {greek}." ) def _auto_select_greek_method( From faa2f62fa2c3023ec431b7f6ed0f014f13a5cdd1 Mon Sep 17 00:00:00 2001 From: Jogi Sidhu Date: Sun, 3 May 2026 20:59:00 +0100 Subject: [PATCH 12/12] Exclude gamma from barrier MC tests --- tests/test_greeks.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/tests/test_greeks.py b/tests/test_greeks.py index 756634c..49a92d9 100644 --- a/tests/test_greeks.py +++ b/tests/test_greeks.py @@ -1814,6 +1814,11 @@ def test_european_barrier_mc_greeks_vs_bsm(direction, action, option_type, strik which is essentially exact. MC uses NUMERICAL bump-and-revalue on simulated paths with a fixed seed. """ + # Gamma is excluded: NUMERICAL bump-and-revalue gamma on MC barriers is + # blocked at the OV level (second-derivative noise floor exceeds |Γ| + # signal at practical path counts; sign flips occur near zero). See + # ``OptionValuation._reject_barrier_numerical`` for the rationale. + barrier_greeks = ("delta", "vega", "theta", "rho") dp_bsm = _dp_barrier_greeks( pricing_method=PricingMethod.BSM, exercise_type=ExerciseType.EUROPEAN, @@ -1822,6 +1827,7 @@ def test_european_barrier_mc_greeks_vs_bsm(direction, action, option_type, strik barrier=barrier, option_type=option_type, strike=strike, + greeks=barrier_greeks, ) dp_mc = _dp_barrier_mc_greeks( exercise_type=ExerciseType.EUROPEAN, @@ -1830,10 +1836,11 @@ def test_european_barrier_mc_greeks_vs_bsm(direction, action, option_type, strik barrier=barrier, option_type=option_type, strike=strike, + greeks=barrier_greeks, ) # MC bump-and-revalue noise on barrier payoffs is substantial; loose tols. - tols = {"delta": 0.05, "gamma": 0.20, "vega": 0.10, "theta": 0.10, "rho": 0.05} + tols = {"delta": 0.05, "vega": 0.10, "theta": 0.10, "rho": 0.05} assert_greeks_close( lhs=dp_mc, rhs=dp_bsm,