STOCHLABSTOCHASTIC SYSTEMS WORKBENCH

REFERENCE / NUMERICAL METHODS

Inspect the assumptions.

Every experiment runs the same C++20 implementation in the native executable and in WebAssembly. Empirical quantities come from seeded simulations; the theory column is computed independently from analytic formulas.

01 / Symmetric random walk

Independent signs take values −1 and +1 with equal probability. With n steps over a horizon T, each increment is scaled by √(T/n). The terminal mean is zero and variance is T for every n. Increasing the number of steps demonstrates diffusive scaling toward Brownian motion.

X(kT/n) = √(T/n) ∑ⱼ₌₁ᵏ εⱼ

All Bernoulli increments are simulated, including undisplayed realizations. There is no Gaussian replacement for the finite random walk.

02 / Brownian motion

B(t + Δ) − B(t) ∼ N(0, Δ) · B(0) = 0

Plotted paths use independent, exact Gaussian increments on a uniform grid. Lines between observations are a rendering convention. Undisplayed endpoints are sampled from N(0,T), avoiding unnecessary trajectory allocations. The endpoint mean and variance are compared with 0 and T.

03 / First passage

τₐ = inf {t ≥ 0 : B(t) = a}
P(τₐ ≤ T) = erfc(a / √(2T)), a ≥ 0

A discrete grid misses some crossings. StochLab samples the continuous maximum conditional on the endpoints of a Brownian bridge. For endpoints x and y separated by Δ, and U uniform on (0,1):

M = ½[x + y + √((y − x)² − 2Δ log U)]

Rendered paths use a bridge maximum on each interval; undisplayed realizations use one exact terminal/maximum pair over the whole horizon. The Monte Carlo crossing probability includes excursions between grid observations. The displayed crossing time is a grid diagnostic, not an exact continuous hitting-time sample.

04 / Reflection principle

After the first hit of a positive barrier a, reflect the remaining trajectory about that barrier: B̃(t) = 2a − B(t). Brownian symmetry leaves the reflected path law unchanged. This gives the factor of two in the first-passage formula.

Use Reflect tail and Next crossing to inspect the construction. The browser reflects a sampled, piecewise-linear path; its visible crossing is an approximation to the continuous hitting time. Aggregate probabilities still use exact bridge maxima.

05 / Quadratic variation

QVₙ = ∑(ΔB)² ∼ (T/n) χ²ₙ
E[QVₙ] = T · Var[QVₙ] = 2T²/n

The lower plot refines nested partitions of the same rendered path. Individual realizations need not converge monotonically. The table also compares aggregate QV moments with theory. For undisplayed realizations, a Gaussian projection and an independent chi-square residual sample the exact joint endpoint/QV law without generating every increment. The chi-square sampler uses the Marsaglia–Tsang gamma algorithm.

06 / Ornstein–Uhlenbeck

dX = θ(μ − X)dt + σdB
X(t+Δ) = μ + (X(t)−μ)e⁻ᶿΔ + σ√((1−e⁻²ᶿΔ)/(2θ)) Z

The engine uses the exact Gaussian transition, with expm1 for small-interval accuracy. It does not use Euler discretization. A deterministic initial state x₀ gives terminal mean μ + (x₀ − μ)e⁻ᶿᵀ and variance σ²(1 − e⁻²ᶿᵀ)/(2θ).

07 / OU stationarity

X∞ ∼ N(μ, σ²/(2θ)), θ > 0

The stationarity experiment starts away from the mean. Increase T to observe relaxation. The table distinguishes finite-horizon theory from the stationary target. Horizons shorter than five relaxation times receive a warning; a finite-time discrepancy from the stationary law can be transient bias, not a simulation error.

08 / Poisson arrivals

Δtᵢ ∼ Exp(λ) · E[N(T)] = Var[N(T)] = λT

The engine advances through exact event times. There are no fixed time buckets. Counting paths contain pre- and post-event points; the raster uses the original arrival times. A zero rate gives an empty event process.

09 / Compound Poisson

Each Poisson event carries an independent Gaussian mark J with configurable mean m and deviation s. Their sum is a jump process with terminal mean λTm and variance λT(s² + m²). Marks can be negative; only the underlying event count is monotone.

10 / Finite-state CTMC

Enter a square generator matrix Q with nonnegative off-diagonal rates and rows summing to zero. In state i, draw an exponential holding time at rate −Qᵢᵢ, then choose j ≠ i in proportion to Qᵢⱼ. Zero exit rates are absorbing. Up to 16 states are supported.

Occupation uses actual holding times, including the final interval. A pivoted linear solver finds πQ = 0 subject to ∑πᵢ = 1. If multiple closed classes or ill-conditioning prevent a reliable unique solution, the comparison is omitted with a warning. Finite-horizon occupation includes the initial transient.

11 / Self-exciting Hawkes arrivals

λ(t) = μ + ∑ₜᵢ<ₜ α exp(−β(t − tᵢ))

Ogata adaptive thinning samples exact event times. The current intensity bounds the decaying intensity until the next candidate; a candidate is accepted by the ratio of the decayed intensity to that bound. Accepted events raise intensity by α.

Here α is the jump in intensity. The integrated kernel mass, or branching ratio, is α/β. The workbench requires α < β and warns near criticality. Intensity begins at μ with an empty event history.

d = β − α
E[N(T)] = (μβ/d)T − (μα/d²)(1 − e⁻ᵈᵀ)
E[λ(T)] = μ + (μα/d)(1 − e⁻ᵈᵀ)

The comparison uses these finite-horizon values. The stationary intensity μ/(1 − α/β) is a separate long-run diagnostic. Intensity curves retain pre-/post-event jumps and sample the intervening exponential decay; the event raster shows exact accepted times.

12 / Synthetic exchange

The exchange stores individual FIFO orders at integer price levels. Limit orders match eligible opposing levels, market orders consume available depth, and cancellations remove resting quantities. Prices emerge from book events. They are not drawn from a stock-price SDE.

External eventProbabilityBehavior
Limit order52%Join, improve, or rest near the book; 1–24 units
Market order30%Match immediately; discard unfilled remainder
Cancellation18%Remove a uniformly selected external order

Arrivals use Poisson or univariate Hawkes clocks. Buy and sell sides have equal probability. The book is initialized near 100.00. A synthetic reservoir replenishes an empty side, and capacity maintenance bounds the resting book at 4,096 orders. These interventions are labeled in the tape and materially shape the price dynamics.

Inventory-skew market maker

A single agent posts passive, ordinary FIFO quotes around an inventory-adjusted reservation price. Every eight external events, or after a quote disappears, it replaces its quotes. The reservation price is the external midpoint minus skew × inventory, measured in ticks. Quotes have at most eight units and inventory is limited to ±100.

Fills arise from actual resting orders. Cash is accounted in integer tick units. Marked-to-market P&L is cash + inventory × midpoint. Spread capture measures signed fill value relative to the pre-event midpoint; it is not realized profit. There are no fees, latency, or market calibration. This sandbox demonstrates stochastic interaction and makes no claim about real-market profitability.

Validation & reproducibility

The native suite checks deterministic invariants and 18 seeded empirical-versus-theoretical comparisons. Most statistical tolerances use five or six sampling standard errors. Exchange tests exercise FIFO, partial fills, cancellations, accounting, book bounds, flow clocks, and independence from batch size.

All realizations contribute to moments and histograms. At most 24 continuous paths or 8 event paths are transferred to the browser. The RNG is an explicit SplitMix64 stream with Box–Muller normals and inverse-transform exponential waits. Repeating a seed and parameter set reproduces the experiment in the same build, excluding timing. Transcendental rounding may differ slightly across platforms.

Variance uses n − 1 in the denominator. SE denotes one Monte Carlo standard error. A single realization cannot estimate variance reliably and produces a warning. Work budgets reject oversized or pathological requests explicitly; partial aggregates are never silently shown as complete.

make test
make validate
npm test
npm run test:browser

Further reading