From-scratch numerical methods for stochastic differential equations — built for correctness, speed, and reproducibility.
A growing, reproducible collection of numerical methods for stochastic differential equations (SDEs), written from scratch in Julia and C++. Each topic pairs a clean implementation with a benchmarked, written-up analysis — the goal is not just code that works, but code whose allocation profile, parallel scaling, and numerical behaviour are measured and understood.
Built as an applied companion to my PhD research, with an eye toward quantitative finance and research applications.
| Status | Topic | Focus | Write-up |
|---|---|---|---|
| ✅ | Brownian motion | path generation, allocation profiling, SIMD vs. scalar loops, safe multithreading | Read → |
| ✅ | BAOAB Langevin integrator | C++ ⇄ Julia, exact Gibbs validation, a hand-written RNG engine, and a benchmark that is honest about what it does not prove | Read → |
| 🔜 | Monte Carlo option pricing | GBM paths, European options, variance reduction | planned |
| 🔜 | Euler–Maruyama / Milstein | strong & weak convergence, time-stepping schemes | planned |
Adding a topic: drop a note in
docs/<Topic>/, add ascripts/<Topic>/experiment, and add one row to this table.
The same underdamped Langevin integrator written in C++ and in Julia, validated against the exact Gibbs distribution before either implementation is timed.
The sharp correctness test is not var(q) = kT/k. It is the momentum deficit: the
end-of-step momentum satisfies var(p) = m kT (1 - (ω dt)²/4) = 2.8125, not the textbook
m kT = 3.0. That is a falsifiable prediction with a coefficient — a mistyped port can
land on 0.6 by luck, but not also on 2.8125. Both implementations hit it to within 0.1%.
8×10⁸ steps, Apple M2 Pro, one sitting (2026-08-15):
| implementation | RNG engine | normal transform | M steps/s | ns/draw |
|---|---|---|---|---|
| C++ | mt19937_64 |
Marsaglia polar | 46.39 | 13.104 |
| C++ | Xoshiro256pp (written for this) |
Marsaglia polar | 85.42 | 6.627 |
| Julia | Xoshiro |
ziggurat | 48.53 | 3.480 |
Which language wins depends entirely on what is held constant. On library defaults Julia leads by 4.6%. Give both sides the same RNG algorithm — a from-scratch xoshiro256++ in C++, ~10 lines — and C++ leads by 76%, while still using the slower of the two normal transforms.
Key findings:
- Per-step time does not decompose into RNG + arithmetic: the same B/A/O arithmetic leaves a 5.1 ns remainder in C++ and a 17.1 ns remainder in Julia. Both remainders are reported as invalid rather than as measurements.
- How much of a step the engine owns does not transfer between setups. In Julia a 2×
engine difference moved the loop by 8 ps; in C++, swapping
mt19937_64for xoshiro256++ moved it by 9.85 ns — more than it moved the isolated draw. - Two predictions were written down before that run. Both were wrong, and the write-up says so, says why, and marks the replacement explanation as an unverified hypothesis.
- ~33% of remaining time sits in two divisions by mass, left unapplied so that both implementations carry the identical inefficiency.
📖 Full write-up: validation, benchmark, and what it does not prove →
Generating ensembles of Brownian paths is the foundation of every Monte Carlo SDE method. This first topic profiles two kernels — a scalar loop vs. an in-place vectorized version — across three execution strategies, measuring where the time and memory actually go.
Ensemble of 10,000 paths (N = 500, 6 threads):
| Strategy | Median |
|---|---|
| serial + loop | 29.32 ms |
| serial + vectorized | 23.21 ms |
| threaded + vectorized | 7.77 ms |
≈ 3.8× end-to-end (serial loop → threaded vectorized).
Key findings:
- Both kernels do exactly one heap allocation — loop vs. vectorized is a compute win (~1.4× from SIMD), not a memory win.
- Multithreading is safe through buffer ownership, not the
!convention; data races come from sharing a passed-in buffer, never from in-place mutation itself. - Threading caps near ~3× on 6 threads because every task allocates, making the garbage collector a synchronization point. To scale further, kill the per-path allocation (pre-allocated matrix + disjoint columns).
📖 Full write-up: allocation & threading → · C++ vs Julia performance →
git clone https://github.com/Louhokseson/NumericalSDEs.git
cd NumericalSDEs
julia setup.jlsetup.jl adds the required registries (General + HokseonRegistry),
resolves, and instantiates the exact environment pinned in Manifest.toml — so it runs
out of the box on a fresh clone.
Scripts auto-activate the project via DrWatson:
using DrWatson
@quickactivate "NumericalSDEs"Then run a topic, for example:
julia scripts/Brownianmotion/Brownianpath_benchmark.jlThe C++ topics need no Julia setup — CMake ≥ 3.16 and a C++17 compiler:
cmake -S scripts/BAOAB -B scripts/BAOAB/build-release -DCMAKE_BUILD_TYPE=Release
cmake --build scripts/BAOAB/build-release
./scripts/BAOAB/build-release/baoab_serial # correctness: ALL CHECKS PASSED
./scripts/BAOAB/build-release/bench_baoab # ns/stepscripts/ # runnable experiments & benchmarks, organised by topic
# (C++ topics carry their own CMakeLists.txt and README)
docs/ # written analyses, one folder per topic
src/ # shared library code
plots/ # generated figures
test/ # test suite (run in CI)
This project is licensed under the MIT License. See LICENSE for details.
Copyright © 2026 Xuexun Lu.
Xuexun Lu (Hokseon)
PhD Candidate
The Maurer Computational Surface Science Group
University of Warwick, UK