Under plain Heston and Hull-White, operator splitting barely outruns a two-point trapezoidal rule on integrated variance. Remark 3 says the Heston trapezoid "performs quite well" and is "only slightly below" the splitting scheme when accuracy is measured against time. Remark 6 reaches much the same view of Hull-White, saying that version "also performs quite well". The Figure 3 discussion gives splitting only a narrow advantage. The framework still earns its place. Lifted Heston and Barndorff-Nielsen-Shephard make the case for it.

The object that disappears

A stochastic volatility simulation normally needs the integrated variance over each step, conditional on the variance endpoint. Much of the existing machinery is built to recover that quantity. Broadie and Kaya invert a characteristic function for every path. Tse and Wan and Abi Jaber instead fit a moment-matched distribution.

This scheme splits the coupled system and avoids constructing the object. A deterministic variance drift sits in one sub-operator, making its integrated variance a closed-form function of the current variance. The other sub-operator contains the variance noise and correlated price term. Under Heston, it reduces to a difference between variance levels. Hull-White and 4/2 add square-root and log terms.

The composition is Strang-style: half step, full step, half step. For Heston, each step requires one driftless CIR draw and two algebraic updates. The appendix reproduces the Alfonsi sampler. The authors extend the same construction to the 4/2 model (3/2 included), Hull-White, an N-factor lifted Heston whose deterministic leg uses a matrix exponential, and four BN-S specifications: Gamma-OU, TS-OU, OU-Gamma, OU-TS.

All results come from simulated paths on published parameter sets. There is no market data or fit to observed option prices. Runs use 10^5 paths and 50 repetitions by default. Accuracy is reported as RMSE = sqrt(StDev^2 + Bias^2), measured against published benchmarks or a Fourier-cosine price. The paper also proves second-order weak convergence. Its global bound is (1/6)C_loc·T·h^2, obtained from a per-step O(h^3) commutator residual.

Do lifted Heston and BN-S carry the claim?

Lifted Heston does. Table 2 reports 0.95e-5 of splitting bias for the out-of-the-money call at K=1.4 and h=1/5. The Bayer-Breneis weak scheme records 0.69e-5 at the same step. Drift-implicit Euler remains at 248.92e-5 with h=1/400 and 63.79e-5 even at h=1/10000.

Its bias also behaves non-monotonically. The 248.92e-5 result at h=1/400 rises to 256.84e-5 at h=1/650, roughly the reciprocal of the largest node 681.37, before declining. Over part of the CPU range, the implicit Euler RMSE curve therefore moves in the wrong direction. Splitting and the weak scheme both deliver biases of order 1e-5. Runtime separates them. For the at-the-money call, Table 3 gives 3.225s versus 6.920s at h=1/800, then 1.828s versus 4.327s at h=1/500.

The BN-S results are the figures I would use. At h=1, the OU-Gamma Asian call has RMSE 5.87e-4, compared with Euler's 2.63e-3 at the same 0.0021s. For the OU-Gamma VIX call at h=1, splitting misses by -0.0294% while Euler misses by 0.7657%. Their times are 0.0031s and 0.0030s.

TS-OU at h=1 is more revealing. Splitting is off by -0.0134%. Cufaro Petroni and Sabino's Algorithm 1, an exact scheme, is off by -0.3805%. Exact sampling of the variance marginal still permits bias in the joint step, and that bias exceeds the splitting error on coarse grids.

The CPU evidence is weaker. Competing-scheme times come from the authors' own runs, and we found no statement of hardware or language. We also could not connect the conclusion's claim that runtime "can be several hundred times faster" to any table entry. The RMSE-versus-CPU figures contain no printed values. In the tables, the largest ratios are about 4.7x for Gamma-OU VIX at h=1/8, where 0.0155s is compared with Sabino's Algorithm 4 at 0.0726s, and up to about 2.4x against Bayer-Breneis, with 1.828s against 4.327s at h=1/500. The introduction uses the same phrase when describing Choi's OU scheme against Li and Wu (2019).

Remark 2 is direct about plain Heston: the Choi-Kwok Poisson-conditioned scheme is more efficient, reaching the benchmark to three significant digits at h=T in under 0.1s.

Choices left to the implementer

Propositions 1 through 4 provide literal update formulas. Most implementation work is transcription. Three choices remain.

The Heston sub-step is CIR(0,0,sigma,h,.), which gives alpha=0 and 4·alpha < sigma^2 always. The applicable Alfonsi branch has a point mass at zero with probability exp(-d_t·V_0/2), followed by a nested draw whose effective alpha is sigma^2/2. A noncentral chi-square sampler produces a different scheme. The 4/2 sub-step, CIR(sigma^2/2,0,sigma,h,.), falls under the other branch.

Lifted Heston uses nodes from 0.63781 to 681.37. At the largest node, 681.37, exp(-Ah/2) with h=1 contains entries of order 1e-148 beside entries of order 1. Meanwhile, I - exp(-Ah/2) approaches zero in the slow direction. The formulas are supplied; the linear algebra is left to the implementer.

The table of background driving Lévy processes, which drive variance jumps, specifies TS(alpha, beta, theta·t) increments for OU-TS. For TS-OU it gives TS(alpha, beta, lambda^alpha·theta·t) plus a compound Gamma term. We did not find a sampler for the tempered stable variate. Both exact benchmarks require that same variate, so the selected sampler affects each side of the comparison.

The VIX benchmark does not reconstruct

The printed information was insufficient for us to reproduce their VIX benchmark. The stated parameters are lambda=0.5, rho=-10, r=0.03, K=0.15, V0=0.04, T=1. We found no value for the VIX horizon tau_vix, although it enters both A and B in the affine representation.

There is also a sign issue. We derive the OU-Gamma integral in Eq. (36b) as alpha·[rho/beta + ln((beta-rho)/beta)], the negative of the entry tabulated in Corollary 1. With (alpha,beta)=(2,2) and rho=-10, our expression gives -6.42 while the printed expression gives +6.42. The difference moves B by about 25.7 in VIX^2 units. An A·V_T term is of order V_0 = 0.04, while the strike is K=0.15 (VIX^2 = 0.0225). The benchmark price level depends entirely on that sign. Reported relative errors of order 0.03% therefore cannot be checked externally until the discrepancy is resolved.

The convergence theory has another boundary. Proposition 8 assumes sup-norm bounds on double commutators applied to the smoothed test function. The conclusion restricts the result to payoffs "whose payoff functions satisfy these assumptions". We did not find a demonstration that a call or a knocked-out barrier payoff meets them.

Empirically, the rate plots follow the h^2 reference. Yet the Asian convergence benchmark is a fine-grid splitting run using 240 substeps per monitoring interval and shared random increments. It measures the rate rather than the distance from the true price. For BN-S Asians, the RMSE benchmark is itself a biased scheme: Euler at h=1/100 with 10^8 paths.

On a desk

For discretely monitored exotics under BN-S or a lifted Heston approximation, splitting delivers the best accuracy per second among the schemes described here. At the same h, a drift-implicit Euler step costs about a third as much, yet its bias is 248.92e-5 at h=1/400. Coarse-grid behaviour makes the splitting method usable because steps can coincide with monitoring dates. Their BN-S Asian has M=4 monitoring dates. Since error accumulates per step, the monitoring calendar determines the grid.

Four restrictions matter. Every experiment we found uses T=1. The 4/2 scheme is stated to require the Feller condition 2·k·theta >= sigma^2. For lifted Heston, non-negativity in the deterministic leg is proved only when weights, nodes and initial factors are non-negative. The Normal Tempered Stable specification is excluded because it permits negative jumps.

We did not test these results on our own data. The paper contains no trading rule. Converting its pricing engine into one would require option quotes to price against, specifically bid and ask, which we do not hold for the listed option universe. The claim here remains confined to the pricing engine.

Before trusting the convergence result at trading horizons, I would want a parameter set away from T=1 and away from at-the-money.