ARCANE: Scalable high-degree cubature formulae for simulating SDEs without Monte Carlo error

Monte Carlo sampling is the standard approach for estimating properties of solutions to stochastic differential equations (SDEs), but accurate estimates require huge sample sizes. Lyons and Victoir (2004) proposed replacing independently sampled Brownian driving paths with “cubature formulae”, deterministic weighted sets of paths that match Brownian “signature moments” up to some degree $D$. They prove that cubature formulae exist for arbitrary $D$, but explicit constructions are difficult and have only reached $D=7$, too small for practical use. We present ARCANE, an algorithm that efficiently and automatically constructs cubature formulae of arbitrary degree. It reproduces the state of the art in seconds and reaches $\boldsymbol{D=19}$ within hours on modest hardware. In simulations across multiple different SDEs and error metrics, our cubature formulae robustly achieve an error orders of magnitude smaller than Monte Carlo with the same number of paths.

💡 Research Summary

The paper tackles a fundamental limitation of Monte Carlo (MC) simulation for stochastic differential equations (SDEs): the 1/√M convergence rate that forces practitioners to use huge sample sizes for accurate statistics. Building on the cubature‑on‑Wiener‑space framework introduced by Lyons and Victoir (2004), the authors propose ARCANE (Algorithm for Recombination of Cubatures from Orthogonal Arrays that match Nested Expected signatures), a scalable method that automatically constructs high‑degree cubature formulae for low‑dimensional SDEs.

A cubature formula is a deterministic weighted set of paths whose signatures (the collection of iterated integrals) match those of Brownian motion up to a prescribed degree D. In theory such formulae exist for any dimension d and any degree D, but explicit constructions have been limited to D≤7 because the naive combinatorial search grows exponentially. ARCANE overcomes this bottleneck through two technical breakthroughs.

First, it uses orthogonal arrays (OAs) to generate an initial pool of candidate piecewise‑linear paths that already have good combinatorial coverage of the multi‑index space required for degree‑D signatures. The OA structure guarantees that each monomial in the signature expansion appears with the correct frequency, turning the moment‑matching problem into a linear system. Second, the algorithm applies a recombination step formulated as a linear programming (LP) problem: it seeks a minimal subset of the OA‑generated paths together with non‑negative weights that exactly satisfy the linear constraints imposed by the signature moments. This LP can be solved in polynomial time in D and d, and the authors further accelerate it with massive GPU parallelism using JAX’s just‑in‑time compilation and automatic differentiation.

The resulting paths are piecewise linear on dyadic time partitions (depth up to 8), which means that the Stratonovich SDE reduces to a sequence of ordinary differential equations on each sub‑interval. Consequently, any standard ODE solver can be used without the need for higher‑order rough‑path machinery, and the computational cost per path is comparable to a standard MC trajectory.

The authors provide rigorous theoretical guarantees: they prove that the recombined weighted set exactly matches all signature tensors up to degree D (Theorems 3.1 and 3.2) and that the algorithm’s runtime scales as a low‑order polynomial in D and d, rather than exponentially. They also discuss numerical stability and the impact of dyadic versus non‑dyadic constructions.

Empirically, ARCANE is evaluated on five widely‑used SDE models: Vasicek (Ornstein‑Uhlenbeck), Cox‑Ingersoll‑Ross (CIR), an Inhomogeneous Geometric Brownian Motion (IGBM), the Wright–Fisher diffusion, and the two‑dimensional Log‑Heston model. Three error metrics are considered: (i) Mean‑Variance Error (MVE) at a fixed terminal time, (ii) Relative error of the initial bond price (closed‑form for Vasicek and CIR), and (iii) Relative error of a European call price (closed‑form for Log‑Heston). Across all experiments, cubature formulae with degrees ranging from 5 to 19 consistently outperform plain MC and two quasi‑Monte‑Carlo baselines (Sobol and Latin hyper‑cube) by several orders of magnitude for the same number of paths (typically 200–8000). For high degrees (D≥15) the error often reaches the numerical tolerance of the ODE solver (≈10⁻⁸), indicating that the cubature error is essentially negligible.

The paper also explores limitations. In the Wright–Fisher case, where the diffusion coefficient becomes non‑smooth near the boundaries 0 and 1, the high‑degree cubature loses some advantage, highlighting the dependence of cubature accuracy on the smoothness of the SDE coefficients. Long‑time horizons (T = 3, 10) degrade accuracy modestly, but the cubature still retains a clear edge over MC/QMC. Dyadic cubatures sometimes exhibit slightly larger degradation than non‑dyadic ones, an effect the authors note for future investigation.

Implementation details are provided: the core algorithm is written in JAX, enabling automatic differentiation, just‑in‑time compilation, and seamless GPU scaling. The authors release both the code repository (github.com/tttc3/ARCANE‑Cubature) and the constructed cubature tables for reproducibility.

In summary, ARCANE delivers a practical pathway from the existence theory of high‑degree Wiener cubature to usable numerical tools. By reducing the combinatorial explosion to a tractable LP, leveraging orthogonal arrays, and exploiting modern GPU hardware, the method constructs deterministic path sets up to degree 19 in a few hours on modest hardware. The resulting cubature formulae achieve dramatically lower estimation error than MC or QMC for a variety of SDEs, while remaining compatible with standard ODE solvers. The work opens the door to high‑precision SDE simulation in finance, physics, biology, and engineering, and suggests several promising extensions: higher dimensions (d > 3), adaptive handling of non‑smooth coefficients, and integration into large‑scale stochastic modeling pipelines.