Feature/barriers

README.md

@@ -41,7 +41,6 @@ knock-in and knock-out structures, and rebates.
 - **Stochastic processes** — Geometric Brownian Motion, Jump Diffusion (Merton), Square-Root Diffusion (CIR)
 - **Discount curves** — log-linear interpolation on arbitrary term structures; deterministic time-varying forward rate and dividend curves
 - **Discrete dividends** — supported across all pricing methods
-- **Barrier options** — continuous and discrete monitoring, knock-in/knock-out, rebates (at-hit and at-expiry)
 - **Control variates** — European analytical control variates for American pricing variance reduction
 - **Custom payoffs** — user-defined payoff functions via `PayoffSpec`
 
@@ -93,7 +92,7 @@ import derivatives_pricing as dp
 pricing_date = dt.datetime(2025, 1, 1)
 maturity = dt.datetime(2025, 7, 1)
 
-dc = dp.DiscountCurve.flat(rate=0.05, end_time=1.0)
+dc = dp.DiscountCurve.flat(rate=0.05)
 md = dp.MarketData(
     pricing_date=pricing_date,
     discount_curve=dc,
@@ -164,10 +163,6 @@ examples/                     # API usage notebooks
 tutorials/                    # Theory deep-dive notebooks
 ```
 
-## Roadmap
-
-Planned: stochastic volatility models.
-
 Found a bug or have a feature request? [Open an issue](https://github.com/jsidhu06/derivatives-pricing/issues).
 
 ## Disclaimer

src/derivatives_pricing/rates.py

@@ -114,22 +114,17 @@ def from_zero_rates(
         return cls(times=times, dfs=dfs)
 
     @classmethod
-    def flat(
-        cls,
-        rate: float,
-        end_time: float,
-        steps: int = 1,
-    ) -> DiscountCurve:
+    def flat(cls, rate: float, end_time: float = 100.0) -> DiscountCurve:
         """Build a flat continuously-compounded discount curve.
 
         Parameters
         ----------
         rate
             Flat continuously-compounded annual rate.
         end_time
-            Final maturity in years.
-        steps
-            Number of intervals used to discretize ``[0, end_time]``.
+            Final maturity in years.  Defaults to 100, large enough to
+            cover any realistic option maturity.  Override only if you
+            have a specific reason to truncate the curve domain.
 
         Returns
         -------
@@ -138,9 +133,7 @@ def flat(
         """
         if end_time <= 0.0:
             raise ValidationError("end_time must be positive")
-        if steps < 1:
-            raise ValidationError("steps must be >= 1")
-        times = np.linspace(0.0, float(end_time), int(steps) + 1)
+        times = np.array([0.0, float(end_time)])
         dfs = np.exp(-float(rate) * times)
         return cls(times=times, dfs=dfs)
 

src/derivatives_pricing/valuation/__init__.py

@@ -1,15 +1,18 @@
 """Option valuation and pricing engines.
 
-This module provides a unified interface for pricing vanilla and custom options
-using various methods: Monte Carlo simulation, Binomial trees, Black-Scholes-Merton
-analytical formulas, and PDE finite difference methods.
+This module provides a unified interface for pricing vanilla, custom-payoff,
+Asian, and barrier options using various methods: Monte Carlo simulation,
+Binomial trees, Black-Scholes-Merton analytical formulas, and PDE finite
+difference methods.
 
 Public API
 ----------
 Core classes:
     OptionValuation: Main dispatcher for option pricing
     VanillaSpec: Contract specification for vanilla options
     PayoffSpec: Contract specification for custom payoffs
+    AsianSpec: Contract specification for Asian options
+    BarrierSpec: Contract specification for barrier options
     UnderlyingData: Minimal underlying data container
 
 Parameter classes:

src/derivatives_pricing/valuation/barrier_analytical.py

@@ -543,16 +543,12 @@ def present_value(self) -> float:
             H = _broadie_glasserman_adjustment(H, sigma, T, spec.num_observations, spec.direction)
 
         # ── No-rebate barrier value ──
-        # 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.
+        # Reiner-Rubinstein contains ``(H/S)**(2*lambda)`` and divides by
+        # ``sigma*sqrt(T)``, so ``sigma -> 0`` produces inf/nan via numpy
+        # promotion or OverflowError/ZeroDivisionError in pure-Python
+        # paths.  Catch both (``np.errstate`` turns the silent
+        # floating-point case into ``FloatingPointError``) and fall back
+        # to the deterministic-forward closed form below.
         with np.errstate(over="raise", invalid="raise", divide="raise"):
             try:
                 value = _barrier_price_no_rebate(

src/derivatives_pricing/valuation/binomial.py

@@ -1,7 +1,14 @@
 """Binomial-tree valuation engines (Cox-Ross-Rubinstein).
 
-Implements European and American vanilla option pricing, plus Asian-option
-extensions used by the core dispatcher.
+Implements European and American pricing for:
+
+- vanilla call/put
+- Asian options (Hull-style tree averages or MC sampling on the tree)
+- barrier options (continuous and discrete monitoring, KO and KI, with
+  Boyle-Lau on-node alignment for continuous and a half-step CRR-layer
+  analog for discrete)
+
+Plugged into the registry dispatcher in ``core.py``.
 """
 
 from __future__ import annotations
@@ -958,37 +965,26 @@ def _resolve_effective_num_steps(self) -> int:
         if spot <= 0.0 or barrier <= 0.0 or sigma <= 0.0 or ttm <= 0.0:
             return base_steps
 
-        # If the option is already triggered at inception, `present_value()`
-        # will short-circuit through `_inception_short_circuit_value` without
-        # running the barrier tree. Neither Boyle-Lau inflation nor the
-        # cap-bind warning is meaningful in that case:
-        #   - KO: the option is dead at t=0 and PV is a closed-form rebate
-        #     discount; the barrier-aware solver never runs.
-        #   - KI: the option becomes vanilla immediately and is priced via
-        #     `_solve_backward` (not the barrier-aware solver), so barrier
-        #     alignment is irrelevant.
+        # Inception-triggered specs short-circuit upstream in
+        # `_inception_short_circuit_value`; the barrier solver never runs,
+        # alignment is moot.
         if self.valuation_ctx._barrier_triggered_at_inception():
             return base_steps
 
         log_distance = np.log(max(spot, barrier) / min(spot, barrier))
         divisor = log_distance * log_distance
-        # Only bail if log_distance is exactly zero (barrier coincides with
-        # spot — a triggered-at-inception case handled upstream by
-        # `_inception_short_circuit_value`). We must NOT use `np.isclose` here
-        # with its default atol=1e-8 because that would silently disable the
-        # Boyle-Lau step adjustment for all near-spot barriers (any barrier
-        # closer to spot than ~1e-4 in log-space would fall inside the
-        # tolerance), exactly the regime where BL alignment matters most.
+        # Strict `<= 0` (not `np.isclose`): default atol=1e-8 would
+        # silently disable alignment for any barrier within ~1e-4 of spot
+        # in log-space — the regime where alignment matters most.
         if divisor <= 0.0:
             return base_steps
 
-        # Continuous monitoring: place a CRR layer exactly ON H (Boyle-Lau 1994).
-        # Discrete monitoring: place H midway between two CRR layers — the
-        # binomial probability mass spans half a layer on each side of each
-        # layer's nominal position, so on-H placement biases the effective
-        # kill threshold by -Δ/2.  Half-step placement (analog of Cheuk-Vorst 1996 /
-        # Boyle-Tian 1998 / in the FD setting) gives an unbiased
-        # discrete kill probability.
+        # Continuous: place a CRR layer on H (Boyle-Lau 1994, factor = i).
+        # Discrete:   place H midway between two CRR layers (factor = i - 0.5).
+        # Under discrete monitoring, placing H exactly on a CRR layer introduces
+        # a first-order grid-placement bias in the effective kill threshold;
+        # half-step placement is the binomial analog of Cheuk-Vorst 1996 /
+        # Boyle-Tian 1998.
         shift = 0.0 if self.spec.monitoring is BarrierMonitoring.CONTINUOUS else 0.5
         max_steps = max(1000, base_steps * 5)
         optimum_steps = base_steps
@@ -997,14 +993,9 @@ def _resolve_effective_num_steps(self) -> int:
             candidate = int((factor * factor * sigma * sigma * ttm) / divisor)
             if candidate >= base_steps:
                 if candidate > max_steps:
-                    # Boyle-Lau alignment requires a tree with `candidate` time
-                    # steps to place a layer of CRR nodes exactly on the barrier
-                    # (typically because the barrier sits very close to spot).
-                    # The `max_steps` cap prevents runaway memory, so the final
-                    # tree will run at `max_steps` without barrier alignment and
-                    # will converge extremely slowly (classic Boyle-Lau bias,
-                    # O(1/√n) with a large constant). Warn the user so they
-                    # know to switch engines.
+                    # Alignment unattainable within memory budget — warn
+                    # and fall back to unaligned at `max_steps` (the warning
+                    # message itself spells out the implications).
                     log_distance_pct = 100.0 * log_distance / max(np.log(spot), 1.0e-12)
                     warnings.warn(
                         (

src/derivatives_pricing/valuation/core.py

@@ -2,7 +2,7 @@
 
 This module is the central orchestration layer for pricing:
 
-- Spec dataclasses (`VanillaSpec`, `PayoffSpec`, `AsianSpec`)
+- Spec dataclasses (`VanillaSpec`, `PayoffSpec`, `AsianSpec`, `BarrierSpec`)
 - Underlying data container (`UnderlyingData`)
 - Registry-based dispatcher (`OptionValuation`) that maps
     `(PricingMethod, ExerciseType)` to a private implementation engine
@@ -531,15 +531,10 @@ def delta(
             return float(self._impl.delta())
 
         if isinstance(self._spec, BarrierSpec) and self._barrier_triggered_at_inception():
-            # Inception-triggered short-circuit (NUMERICAL only):
-            # bump-and-revalue would cross the trigger boundary and price
-            # the un-triggered contract on the bumped spot, which is
-            # meaningless for an already-triggered barrier.  Engines handle
-            # TREE/GRID triggered greeks natively, so this branch is only
-            # reached for NUMERICAL.
-            #   • KO triggered → cashflow constant in spot → δ = 0.
-            #   • KI triggered → contract IS the vanilla equivalent; bump
-            #     there (no state transition).
+            # Inception-triggered NUMERICAL short-circuit: bumping spot
+            # would cross the trigger boundary and price the un-triggered
+            # contract.  KO → constant cashflow → δ = 0;  KI → delegate
+            # to the vanilla equivalent.
             if self._spec.action is BarrierAction.OUT:
                 return 0.0
             return self._vanilla_equivalent_valuation().delta(
@@ -1041,14 +1036,11 @@ def _resolve_params(
                 return MonteCarloParams()
             if pricing_method is PricingMethod.BINOMIAL:
                 if isinstance(spec, BarrierSpec):
-                    # Continuous barriers get Boyle-Lau on-node step inflation
-                    # (place a CRR layer ON H); discrete barriers get half-step
-                    # inflation (place H midway between two CRR layers — the
-                    # binomial-tree analog of Cheuk-Vorst 1996 / Boyle-Tian 1998
-                    # for the FD/trinomial setting).  Discrete defaults to 3000
-                    # base steps because CRR's intrinsic O(1/N) finite-step error
-                    # remains visible on small-price reverse-barrier cases (UOC
-                    # tighter, DOP tighter, UIP wider).
+                    # Continuous: on-node alignment (Boyle-Lau 1994).  Discrete:
+                    # half-step alignment (binomial analog of Cheuk-Vorst 1996 /
+                    # Boyle-Tian 1998).  Discrete defaults to 3000 base steps to
+                    # keep CRR's O(1/N) finite-step error tight on small-price
+                    # reverse-barrier cases.
                     num_steps = 1000 if spec.monitoring is BarrierMonitoring.CONTINUOUS else 3000
                     return BinomialParams(num_steps=num_steps)
                 return BinomialParams()
@@ -1449,25 +1441,20 @@ def _reject_barrier_binomial_numerical(
     ) -> None:
         """Block NUMERICAL bump-and-revalue greeks on binomial barrier specs.
 
-        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
-        two unrelated tree topologies rather than approximating ``∂V/∂x``.
-        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.
+        Both continuous (Boyle-Lau 1994 on-node) and discrete (half-step
+        Cheuk-Vorst / Boyle-Tian analog) alignment pick the tree step
+        count from ``candidate ≈ factor² σ² T / log(H/S)²``, which depends
+        on spot, volatility, and time-to-maturity.  Bumping any of those
+        re-resolves alignment, so a central difference compares trees of
+        *different topology* rather than approximating ``∂V/∂x``.
+
+        Rho is exempt — the risk-free rate doesn't enter the alignment
+        formula, so rate bumps reuse the same tree topology and the
+        finite difference is well-defined.
+
+        Guidance: use ``GreekCalculationMethod.TREE`` for Δ/Γ/Θ; switch
+        to ``PricingMethod.PDE_FD`` for vega (and for any NUMERICAL
+        bump-and-revalue, where the grid topology is bump-invariant).
         """
         if allow:
             return
@@ -1516,28 +1503,23 @@ def _reject_barrier_numerical(
         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.
+        """Block NUMERICAL bump-and-revalue on barrier greek/engine
+        combinations known to be unreliable.
 
-        Two distinct greek failure modes are currently covered (see
-        ``_BARRIER_NUMERICAL_REASON``):
+        Currently:
 
-        - 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
+        - **Γ on MC barriers** (any monitoring): second-derivative MC
+          noise scales ~ stderr/ε²; at practical path counts the noise
+          floor exceeds the |Γ| signal.  Pass
           ``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.
+        - **Θ on discretely-monitored barriers** (any engine): bumping
+          the pricing date forces re-resolution of the monitoring
+          schedule, so the bumped contract is a different contract.
+          Pass ``monitoring_constraint=BarrierMonitoring.DISCRETE``.
+          (Binomial barrier NUMERICAL greeks are blocked at a finer
+          grain in ``_reject_barrier_binomial_numerical``.)
+
+        Per-greek reasons are kept in ``_BARRIER_NUMERICAL_REASON``.
         """
         if method is not GreekCalculationMethod.NUMERICAL:
             return
@@ -1630,25 +1612,19 @@ def _validate_bump(
     ) -> None:
         """Validate a numerical-greek bump argument against the resolved method.
 
-        Two checks are bundled here so every greek entry point can run a
-        single line:
-
-        1. The bump must be strictly positive when supplied. Negative
-           bumps are almost certainly user mistakes — central differences
-           are sign-symmetric so a negative spot/vol/rate bump silently
-           gives the same magnitude with confusing semantics, while a
-           negative ``time_bump_days`` flips the forward-difference theta.
-        2. The bump must be compatible with the *resolved* greek method.
-           Passing ``epsilon=...`` together with a method that doesn't use
-           a finite-difference bump (ANALYTICAL / GRID / TREE / PATHWISE
-           where applicable / LR) is a contradiction — the bump would be
-           silently ignored, letting users believe they are controlling
-           it when they aren't.  We therefore require callers to resolve
-           the method first (via ``_resolve_greek_method``) and pass the
-           resolved method here; that way the same rule applies whether
-           the user supplied an explicit method or relied on auto-select.
-           ``extra_allowed_methods`` lets e.g. gamma also accept PATHWISE
-           (which uses a finite-difference epsilon under the hood).
+        Two checks, bundled so every greek entry point is one line:
+
+        1. **Sign**: bump must be strictly positive — negative bumps are
+           silent footguns (sign-symmetric central diff for spot/vol/rate
+           gives the same magnitude; negative ``time_bump_days`` flips
+           the forward-difference theta).
+        2. **Method compatibility**: passing a bump with a non-bumping
+           method (ANALYTICAL / GRID / TREE / LR / PATHWISE where it
+           doesn't use ε) is a contradiction — the bump would be
+           silently ignored.  Callers must pre-resolve via
+           ``_resolve_greek_method`` so this check is invariant to
+           explicit-method vs auto-select.  ``extra_allowed_methods``
+           lets e.g. gamma also accept PATHWISE (which uses ε internally).
         """
         if bump_value is None:
             return