Hey everyone. I wanted to share an architectural solution to a mathematical edge case that blew up one of my risk engines a while back, specifically when pricing Knock-Out Barrier options and calculating second-order Greeks.
When we don't have closed-form analytical Greeks (like in most path-dependent exotics), we rely on Central Finite Differences:
Gamma ≈ [V(S_0 + dS) - 2V(S_0) + V(S_0 - dS)] / (dS)^2
To do this efficiently in a Monte Carlo simulation without the variance tearing the derivatives apart, the standard practice is Common Random Numbers (CRN). We apply the exact same stochastic shock Z to the base path, the upper-bumped path, and the lower-bumped path inside the hot loop.
The Discontinuity Problem:
Barrier options possess a step-function discontinuity. Let's say we have a Down-and-Out Put. If the barrier H is 85.0, and our initial spot S_0 drops to 85.0001.
When we calculate the Greeks, the numerical bump -dS forces the lower-bumped path to instantly breach the barrier. The payoff evaluates to strictly 0.0.
Because Gamma divides by (dS)^2 (a microscopically small number), the sudden absolute drop in the V(S_0 - dS) term is interpreted by the algorithm as infinite convexity. Your engine outputs a Gamma of 999,999.0 or -infinity. If you have an automated delta-hedging script hooked to this output, it will violently over-leverage your portfolio trying to hedge a mathematical ghost.
The Algorithmic Solution:
I realized that catching this after the matrix computation was too late and computationally wasteful. The check needs to be embedded directly at the C++ core before the finite difference execution.
If the absolute distance between the Spot and the Barrier is less than or equal to 2 * dS, the boundary is breached by the numerical bump. We must flag the state as unstable and force the engine to yield NaN for Gamma, while preserving the Fair Value and Delta calculations.
Implementation & Testing it out:
I ended up building a dedicated C++ OpenMP pricing engine to handle these massive matrices because Python/NumPy was choking on the GIL when simulating 50M+ paths with barrier logic. I wrapped it behind a Python SDK.
If anyone is backtesting exotic portfolios and wants to see how this discontinuity handling works in practice (or just needs to compute 100 million paths in ~3 seconds), I made a Google Colab notebook demonstrating it.
You can run the stress test directly in the browser here:
https://colab.research.google.com/github/Prometheus-Quant-Engineering/prometheus-quant-examples/blob/main/03_HPC_Asynchronous_Polling_Stress_Test.ipynb
The SDK is open source (pip install prometheus-qengine). Let me know how you guys handle step-function discontinuities in your own proprietary risk engines, always looking to optimize the core loop further.