quantark · PDE engine note

Why nodal event masks failed
and how the cell-average projection fixes them

Autocallable PDE prices drifted by up to a few percent of premium, moved when the grid was refined, and flipped sign with a one-ULP nudge of a barrier. The cause was one line of discretization philosophy: applying discrete knock-out, knock-in, and coupon triggers as Boolean masks over grid nodes.

branch worktree-pde-event-projection engines 1D BSM/LV + 2D Heston/SLV 2026-07-24

The operator in question

A snowball or phoenix is priced by backward induction. Between observation dates the solver diffuses the value surface; at each discrete observation it applies an event operator — a hard switch between two branches of the contract, decided by where spot sits relative to a barrier:

V+(x)  =  𝟙breach(x) · Vbreach(x)  +  (1 − 𝟙breach(x)) · Vsurvive(x)

Vbreach is the knocked-out payoff, the knocked-in surface, or the continuation plus a coupon; Vsurvive is the untouched continuation. The operator is exact in the continuum. Everything that went wrong happened when it met a finite grid.

Where nodal breaks

The legacy implementation evaluated the indicator at the nodes: mask = s_vec >= barrier, then overwrote masked nodes with the breach branch. Four distinct failures follow from that single decision.

Half-cell trigger shift

A node speaks for its whole dual cell — the interval reaching halfway to each neighbour. Assigning the node to one branch assigns the entire cell to it, so the effective barrier the discrete system sees is displaced from the contractual one by up to half a cell. That is an O(Δx), one-sided error: on a mesh that always snaps the barrier the same way, refinement shrinks it slowly (O(1/N)) and never averages it out. Measured on a protected 24-observation phoenix: +0.140% / +0.092% / +0.048% of notional at N = 400 / 800 / 1600.

FP roundoff picks the side

The auto-grid deliberately places a node on each barrier — in log space. The node comes back as exp(log(B)), which is B ± 1 ULP with arbitrary sign per barrier. Whether the on-barrier cell counts as breached is then decided by floating-point roundoff, inconsistently across the KO, KI, and coupon barriers of the same product. Switching the comparison between >= and > changed nothing — bitwise. This is the grid-phase lottery: a coarse grid can look accidentally accurate while the refined one converges to the wrong number. The nodal 2D Heston solver plateaus 0.135 below a stride-corrected Monte Carlo reference — about 7% of that product’s premium — while its coarse grids sit near the truth.

Valuation-date blending

An observation falling on the valuation date is not random at all — today’s spot is known. But the event was written onto the grid as a value jump, and price() then linearly interpolated across that jump at an off-node spot. A two-unit coupon due today, spot exactly on the barrier, paid 0.092 / 1.195 / 1.402 / 1.814 / 0.227 units as the uniform grid went N = 100 → 1200. The right answer is 2, at every N.

Undamped jumps

Each event writes a discontinuity onto the surface mid-solve. Crank–Nicolson (and Craig–Sneyd in 2D) are only neutrally stable for the high-frequency modes such a jump excites, so the error rings — visible as gamma oscillating with N. Damping existed in the 1D solver but was accidentally coupled to the auto_grid flag, and the 2D ADI loop ignored the setting entirely.

Watch the cell decide

One knock-out observation, seven grid nodes, the barrier sliding through the dual cell of node 3. The nodal rule flips the whole cell the instant the barrier crosses the node; the projection assigns the node the exact area-weighted average of the two branches over its cell, so the assigned value moves smoothly with the barrier.

Vsurvive (continuation) Vbreach (KO payoff) nodal assignment cell-average projection
Node 3’s dual cell (shaded) spans the midpoints to its neighbours. Away from the barrier the two rules agree; inside the straddled cell only the projection responds to where in the cell the barrier actually sits.
The value assigned to node 3 as the barrier moves through its cell. The nodal rule is a step — a 30-point cliff at the node, where ±1 ULP of exp(log(B)) picks the side. The projection is the smooth ramp the continuum operator actually implies; nudging the barrier by ε moves the result by O(ε).

The fix, in four layers

Dual-cell projection

The event is represented in the finite-volume sense. Every node keeps the pointwise value of its own branch; the one node whose dual cell Ωi straddles the barrier receives the exact cell average of the complete composite function, with both branches read as the grid’s piecewise-linear interpolants:

V+i = (1/|Ωi|) ∫Ωi [ 𝟙breach Vb + (1−𝟙breach) Vs ] dx

Weights live in [0, 1], the result stays inside the min/max envelope of the two branches, and a 1-ULP barrier shift moves it infinitesimally instead of flipping a cell. Smooth convergence is restored, which also makes Richardson extrapolation meaningful again. One subtlety mattered: averaging only the jump and adding it to the pointwise survive value mixes two discretizations — a review counterexample drove a node to −11.875 from entirely non-negative branches. Averaging the complete function is envelope-preserving by construction.

Valuation date is deterministic

At t = 0 the event outcome at the actual spot is known, so nothing is averaged and nothing is blended. The readout interpolates the smooth 0+ branch surface captured before today’s events land on the grid, then applies the transitions pointwise at spot with the product’s own inclusive comparison — the two-unit coupon now pays exactly 2.000 at every N, 1D and 2D. Fixing this exposed a latent 2D bug: the ADI hooks tested is_close(T − τ, 0), but τ accumulates FP step increments, so T − τ lands at ~1e−16 and a relative tolerance against zero can never fire — 2D valuation-date events had been silently cell-averaged all along. All sites now compare τ against T.

Per-event implicit restarts

After every discrete event the scheme restarts with damped (Rannacher-style) implicit steps before returning to second order — the combination the literature shows is required, not optional (Pooley–Vetzal–Forsyth 2003; d’Halluin–Forsyth–Vetzal 2005; Giles–Carter 2006). Damping is decoupled from auto_grid, and the 2D ADI loop now mirrors the 1D schedule exactly: use_rannacher is the master switch, restarts run at event_theta, and maturity-date events are excluded because the terminal Rannacher start-up already owns the payoff discontinuity.

Composite triggers, one pass

A phoenix coupon and a knock-out observed on the same date are one piecewise contractual function with two thresholds, not two events in a row. Projecting them sequentially double-averages any dual cell the thresholds share; project_piecewise_event integrates the K-threshold composite exactly in a single pass (equal thresholds collapse to an empty region and the correct ½ weight). Well-separated barriers reproduce the sequential results.

What stays nodal on purpose. Continuously monitored KI keeps the on-node treatment — for a continuous barrier the value function is continuous there and node alignment is the correct discretization (Boyle–Tian). And event_projection="nodal" survives as an explicit legacy mode, pinned bitwise by the historical characterization oracles.

Evidence

The same protected 24-observation phoenix that exhibited the bias, priced against agreeing QUAD and Monte Carlo references (within +0.01%):

nodal, auto grid cell-average projection (auto and uniform)
PV error vs grid size. Nodal decays O(1/N) from a one-sided offset; the projection sits at +0.010% at every N, on both auto and uniform meshes — and delta converges to 5·10−6 already at N = 400.
grid NPV · cell_averagePV · nodalΓ · grid stencil
40095.537495.6221−0.04335
80095.537795.5768−0.04336
160095.536695.5666−0.04337
dense reference95.5367—−0.04336
Default-certification gate, monthly synthetic snowball. The projected PV is flat in N (spread 8·10−4) while nodal drifts by 0.055; the grid-stencil gamma — the default public greeks path — is stable to 10−5. A separately reviewed “bump gamma” instability turned out to be an estimator artifact (a deprecated 1 bp default bump second-differencing a piecewise-linear readout inside one grid cell) and reproduces identically under nodal — orthogonal to the projection.

In 2D the projected solver was exonerated end-to-end: its Richardson limit sits +0.004 ± 0.003 from a stride-corrected Monte Carlo consensus (the apparent −0.014 gap was the MC’s own observation-stride bias), while nodal plateaus 0.135 below the same reference.