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..a6d93c7 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 @@ -415,6 +389,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 ────────────────────────────────────────────────── @@ -450,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 @@ -515,20 +570,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/binomial.py b/src/derivatives_pricing/valuation/binomial.py index 70d4014..b77f3cd 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)) @@ -1055,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/src/derivatives_pricing/valuation/contracts.py b/src/derivatives_pricing/valuation/contracts.py index f457aaf..8efc78f 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 — @@ -547,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 13ea9a2..38416c8 100644 --- a/src/derivatives_pricing/valuation/core.py +++ b/src/derivatives_pricing/valuation/core.py @@ -27,11 +27,12 @@ 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 ( AsianAveraging, + BarrierAction, BarrierDirection, BarrierMonitoring, DayCountConvention, @@ -56,7 +57,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 @@ -152,6 +153,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: @@ -518,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): @@ -576,6 +602,19 @@ 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, + ) + + self._reject_barrier_numerical( + method, greek="gamma", engine_constraint=PricingMethod.MONTE_CARLO + ) + if epsilon is None: epsilon = self._underlying.initial_value / 100 if isinstance(self._spec, BarrierSpec): @@ -689,8 +728,17 @@ def theta( if method is not GreekCalculationMethod.NUMERICAL: return float(self._impl.theta()) + self._reject_barrier_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 @@ -924,7 +972,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: @@ -1008,7 +1058,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 @@ -1016,6 +1066,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. @@ -1290,9 +1368,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 @@ -1300,6 +1378,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 @@ -1320,6 +1407,71 @@ def _reject_barrier_binomial_numerical( "bump-and-revalue on the full grid." ) + # 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(). + _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 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_numerical( + self, + method: GreekCalculationMethod, + *, + greek: str, + engine_constraint: PricingMethod | None = None, + monitoring_constraint: BarrierMonitoring | None = None, + ) -> None: + """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 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._BARRIER_NUMERICAL_REASON.get( + greek, + "the bump-and-revalue path is unreliable for this combination", + ) + raise UnsupportedFeatureError( + 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( self, *, @@ -1376,6 +1528,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, diff --git a/src/derivatives_pricing/valuation/monte_carlo.py b/src/derivatives_pricing/valuation/monte_carlo.py index 4fdb856..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 @@ -1127,20 +1142,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) 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/src/derivatives_pricing/valuation/pde.py b/src/derivatives_pricing/valuation/pde.py index e2d956a..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 [] @@ -2612,22 +2621,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 +2800,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 +2817,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 +2833,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 +2911,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): diff --git a/tests/test_barrier.py b/tests/test_barrier.py index 444c0af..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), } @@ -1999,11 +2001,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 +2075,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" @@ -2312,6 +2321,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 # =========================================================================== diff --git a/tests/test_edge_cases.py b/tests/test_edge_cases.py index 73c4dd5..adab554 100644 --- a/tests/test_edge_cases.py +++ b/tests/test_edge_cases.py @@ -10,7 +10,12 @@ import pytest from derivatives_pricing.enums import ( + AsianAveraging, + BarrierAction, + BarrierDirection, + BarrierMonitoring, ExerciseType, + GreekCalculationMethod, OptionType, PDESpaceGrid, PricingMethod, @@ -18,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, @@ -28,11 +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 @@ -141,6 +156,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 +626,216 @@ 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 + + +# ═══════════════════════════════════════════════════════════════════════ +# 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) 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, 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 = { 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."""