American Option Pricing Engine

GPU-accelerated Bermudan option pricing by Longstaff-Schwartz least-squares Monte Carlo, in CUDA and OpenMP.

Started as a reimplementation of Cvetanoska & Stojanovski, "Using High Performance Computing and Monte Carlo Simulation for Pricing American Options" (arXiv 1205.0106). The paper prices with a perfect-foresight estimator: its algorithm, in its own words, "assume[s] that all future values of the underlying stock price are known" and so "gives the upper bound value of the American option." That upper bound is badly divergent — 56% high at 21 exercise dates, growing without bound — yet the paper presents it as an American-option price. This engine replaces it with Longstaff-Schwartz, the method the paper cites ([12]) but doesn't use, validates the result against two independent references, and runs it 53× faster than a serial CPU baseline.

The headline result

An American call on a non-dividend-paying stock is worth exactly its European counterpart — early exercise is never optimal. That makes Black-Scholes an exact yardstick, and any gap is bias, not noise.

The paper values each path with max(intrinsic, discounted continuation) using that path's own realised future — which, as it acknowledges, gives an upper bound rather than a price. Every path exercises with hindsight, so the estimate is biased high and gets worse the more exercise dates you give it. Longstaff-Schwartz fixes this by estimating continuation from a cross-path regression instead of the realised future.

How it works

flowchart LR
    A["Simulate N paths<br/>GBM, m+1 timesteps"] --> B["Node m<br/>European value of<br/>remaining dt"]
    B --> C{"For i = m-1 … 1"}
    C --> D["Discount every<br/>cashflow one step"]
    D --> E["Regress cashflow on<br/>1, x, x² over<br/>in-the-money paths"]
    E --> F["Exercise where<br/>intrinsic > fitted<br/>continuation"]
    F --> C
    C -->|done| G["Mean × discount<br/>= price ± stderr"]

Two details make it correct rather than merely plausible:

• The regression supplies the exercise decision only. Valuation always uses the realised cashflow, never the fitted value. Substituting the fit back in reintroduces bias.
• Only in-the-money paths enter the regression. Out-of-the-money paths carry no decision, and including them degrades the fit exactly where the exercise boundary lives.

The estimator is low-biased: a fitted policy is suboptimal, and a suboptimal exercise policy under-values. That is the safe direction.

Time grid

With dt = T/(m+1), path value S[i] sits at i·dt. The last grid point is at T − dt, not T: the m exercise opportunities are at dt, 2dt, …, m·dt, and the option still has dt of life left after the final one. Node m therefore needs no regression — holding is worth the European value over the remaining step, which is an exact conditional expectation.

Validation

A call cannot distinguish "LSM works" from "LSM never exercises", because never exercising is the right answer there. So the suite checks both directions:

flowchart TD
    subgraph exact ["Exact references"]
        BS["Black-Scholes<br/>closed form"]
        TREE["CRR binomial tree<br/>same m exercise dates"]
    end
    subgraph mc ["Monte Carlo backends"]
        CALL["Call<br/>never exercises early"]
        PUT["Put<br/>exercise boundary binds"]
    end
    BS -->|"must match exactly"| CALL
    TREE -->|"must match within 0.14%"| PUT
    CALL --> V["17 checks<br/>all passing"]
    PUT --> V
    XB["Cross-backend checks<br/>GPU vs CPU · 1 vs 28 threads"] --> V

The tree is restricted to the same m exercise dates the Monte Carlo grid uses, so both price the identical Bermudan. Comparing against a fully American tree would be apples to oranges — continuous exercise is worth strictly more than m discrete dates.

Performance

Intel i7-14700HX (20C/28T) · RTX 5060 Laptop · CUDA 12.4 · S₀=X=100, T=1, r=5%, σ=20%, m=20

Convergence

QMC reaches three-decimal accuracy at 100,000 paths where the pseudo-random backends need 1,000,000. The LCG curve ticking up between 10⁶ and 4×10⁶ is not a bug — Monte Carlo error is a random walk around zero, not a monotone decline.

Profiling: two wrong guesses, then the real bottleneck

The first correct version ran only 19× serial, which is poor for an embarrassingly parallel problem. Rather than guess, ./build/profile_gpu splits device time into path generation and the LSM backward pass. It killed both obvious hypotheses:

Path generation is transcendental-heavy — one exp per step plus two logs in Moro's tail branch — on a card that runs doubles at 1/64 the FP32 rate. Moving path storage and generation to fp32, while keeping the regression accumulation in double, gives 2.7× for a price change of 8e-6 against a Monte Carlo standard error of 7.1e-3 — a thousandth of the noise it sits inside.

The time-major layout was kept regardless. It is the correct layout and costs nothing; it simply was not the bottleneck.

A bug the fp32 experiment surfaced

At four million paths the single-precision run returned inf. __uint2float_rn(0xFFFFFFFF) rounds up to exactly 2³², so the float LCG can emit exactly 1.0, and Moro then evaluates log(−log(0)). The same hole existed in the double path all along — LCG state 0 gives u = 0 — merely rare enough never to have fired. Both are now clamped inside Φ⁻¹ itself, where the function is genuinely undefined, rather than at the four call sites. Two regression tests cover it.

Build

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j
ctest --test-dir build --output-on-failure

CUDA targets are skipped automatically if no toolkit is present; -DUSE_CUDA=OFF skips them deliberately.

Host compiler: nvcc pins its own supported GCC (CUDA 12.4 tops out at GCC 13). If your default g++ is newer, configure with -DCMAKE_CXX_COMPILER=g++-13, or linking run_cuda_tests fails on an undefined __cxa_call_terminate. CMake warns when it detects the mismatch.

CUDA arch: defaults to sm_80 with PTX emitted alongside, so newer GPUs JIT. Override with -DCMAKE_CUDA_ARCHITECTURES=89 (Ada), 90 (Hopper), etc.

Run

./build/american_serial                 # any of the five backends
./build/american_cuda   call fp32       # [call|put] [fp64|fp32]
./build/american_qmc_cuda put
./build/profile_gpu                     # setup / pathgen / LSM breakdown
./build/validate                        # CPU validation suite

Documentation

Charts regenerate with python3 docs/make_charts.py.

Layout

src/core/     hd_math.hpp     option math compiled for host AND device
              lsm.{hpp,cpp}   Longstaff-Schwartz regression + backward pass
              paths.*         path simulation (shared by serial and OpenMP)
              binomial.*      CRR tree — independent oracle for the tests
              sobol, halton, brownian_bridge, scramble, moro, black_scholes
              backward_induction.hpp   the paper's naive recursion, kept only
                                       so the bias comparison can be made
src/cuda/     path-generation kernels, lsm_gpu.cu (device backward pass),
              cuda_check.h (CUDA_TRY error macros)
src/openmp/   OpenMP backends, standard and QMC
src/benchmark/ benchmark drivers and the GPU profiler
tests/        validate.cpp (no GPU needed) · test_cuda.cpp (GoogleTest)

Known limitations

These are open, not hidden:

• Pseudo-random streams are not independent. Path seeds are (path_id + 1) × 1234567 on a 2³² LCG, so neighbouring paths walk correlated positions on the same lattice. A counter-based generator (cuRAND Philox) is the fix.
• No error bar on the QMC result. sd/√N assumes independent draws and Sobol is deterministic; at 4×10⁶ the QMC backends report a standard error of 0.0071 while the actual error is 0.000017. A real bound needs randomised QMC — average over several independent digital shifts and take the sd of the means.
• LSM is low-biased, so a true confidence interval needs a paired upper bound (the Andersen-Broadie dual).
• Memory bounds the problem. The regression is cross-sectional, so the whole path set must be resident: N × (m+1) values, 336 MB at N=4×10⁶ in fp64.
• The GPU is handicapped here. CUDA 12.4 has no native sm_120 target, so kernels arrive by PTX JIT from sm_80.

References

Glasserman, Monte Carlo Methods in Financial Engineering (Springer, 2004) · Longstaff & Schwartz, "Valuing American Options by Simulation: A Simple Least-Squares Approach", RFS 14(1), 2001 · Joe & Kuo, "Constructing Sobol sequences with better two-dimensional projections", SIAM J. Sci. Comput. 30, 2008 · Cvetanoska & Stojanovski, "Using High Performance Computing and Monte Carlo Simulation for Pricing American Options".