How Broadie–Kaya Exact-Simulation Algorithms Remove Time-Step Bias from Heston Monte Carlo: CIR Transitions, Integrated Variance, Characteristic-Function Inversion and Speed Failure

Reader question: Heston Monte Carlo usually discretises stochastic volatility into time steps. Can we simulate the model without that time-discretisation bias?

Broadie–Kaya exact simulation answers yes for key affine stochastic-volatility models. Instead of approximating the Heston variance process step by step, it samples the terminal variance from its exact CIR transition law, then samples the integrated variance conditional on the variance endpoints, and finally samples the terminal log price conditional on those exact quantities.

This article owns one precise computational job: exact-distribution Monte Carlo transition simulation for Heston-type stochastic volatility. It does not own Heston parameter calibration, Andersen QE, ordinary Euler/full-truncation schemes, Carr–Madan/COS transform pricing, or generic Monte Carlo variance reduction.

This is public mathematical and computational education. It is not financial advice, a trading recommendation, or a claim that an exact transition sampler removes model risk or every numerical error.

1. The Heston dynamics

Under a risk-neutral measure, a standard Heston model is:

dS_t = rS_tdt + √V_tS_tdW_t^S,

dV_t = κ(θ−V_t)dt + σ√V_tdW_t^V,

with:

dW^SdW^V = ρdt.

The variance process is CIR/square-root diffusion. It has a known exact transition distribution, but the stock process depends on the full path of V through the integrated variance:

I = ∫₀^TV_tdt.

That integral is why exact Heston simulation is more difficult than exact simulation of V alone.

2. What discretisation bias means

An Euler-style scheme replaces continuous dynamics by many small steps. Even with infinitely many Monte Carlo paths, a fixed timestep Δt can leave a systematic bias because the simulated process is not exactly the Heston process.

Increasing the number of paths reduces sampling error. It does not remove time-discretisation bias unless the timestep is also refined.

Broadie and Kaya’s central result is that terminal stock price and variance can be sampled from the exact model transition law, allowing an unbiased Monte Carlo estimator for payoffs determined by those sampled dates, apart from ordinary numerical inversion tolerances and random sampling error.

3. Step 1: sample terminal variance exactly

The CIR variance transition has a scaled noncentral chi-square law.

For Δ=T−t, define:

c = σ²(1−e^−κΔ)/(4κ),

d = 4κθ/σ²,

λ = 4κe^−κΔV_t/[σ²(1−e^−κΔ)].

Then:

V_T = cX,

where:

X ~ χ′²_d(λ).

This exact transition avoids the negative-variance repair problem that appears in naive Euler simulation.

4. Step 2: sample integrated variance conditional on endpoints

The difficult random variable is:

I = ∫_t^TV_sds

conditional on V_t and V_T.

Broadie–Kaya derive the conditional characteristic/Laplace-transform machinery for this integrated CIR quantity. The conditional CDF is not available as a simple closed-form inverse. The algorithm therefore evaluates the conditional transform and numerically inverts it to recover the distribution of I.

A uniform random number U is drawn, and the algorithm solves:

F_I|Vt,VT(i) = U

for i.

This numerical inversion is the computational bottleneck.

5. Why exact does not mean “no numerics”

The word exact refers to the stochastic model transition: there is no Euler/Milstein timestep approximation of the SDE over the transition interval.

But a production implementation still has numerical tolerances from:

• transform evaluation;
• Fourier/Laplace inversion;
• quadrature truncation;
• CDF root finding;
• noncentral chi-square random generation;
• floating-point arithmetic.

Those errors should be controlled and benchmarked separately from time-discretisation bias.

6. Step 3: recover the variance Brownian integral

Integrate the variance SDE:

V_T−V_t = κθΔ − κI + σ∫_t^T√V_sdW_s^V.

Therefore:

∫√V dW^V = [V_T−V_t−κθΔ+κI]/σ.

Once V_T and I are known, this stochastic integral is known exactly from the variance equation.

7. Step 4: sample terminal log price

Write:

d log S = (r−½V)dt + √V dW^S.

Decompose the spot Brownian motion:

dW^S = ρdW^V + √(1−ρ²)dW^⊥.

Then conditional on V_t, V_T and I:

log S_T = log S_t + rΔ − ½I + (ρ/σ)[V_T−V_t−κθΔ+κI] + √(1−ρ²)√I Z,

where:

Z ~ N(0,1)

independently of the quantities already conditioned upon.

This is the final exact-transition draw.

8. The complete Broadie–Kaya transition algorithm

1. Input S_t, V_t, model parameters and horizon Δ.
2. Sample V_T from the scaled noncentral chi-square transition.
3. Draw U ~ Uniform(0,1).
4. Evaluate/invert the conditional distribution of I given V_t,V_T until F(I)=U.
5. Compute the correlated variance Brownian integral from the integrated CIR equation.
6. Draw an independent standard normal Z.
7. Construct log S_T from the conditional Gaussian formula.
8. Return S_T,V_T,I and diagnostics.

9. Why the convergence order can improve

With an exact transition sampler, the dominant Monte Carlo error behaves like ordinary sampling error:

O(N^−1/2).

Broadie and Kaya compare total computational-budget convergence with Euler schemes and show that eliminating discretisation bias can improve the rate at which total error falls with computational effort.

The advantage is especially clear when the payoff depends on only a few observation dates, so the expensive conditional integrated-variance inversion is not repeated many times per path.

10. Exact transition versus exact path

Broadie–Kaya samples exact transitions between selected times. It does not magically reveal every point of the continuous path between those times.

For a European payoff observed only at T, one exact transition can be enough.

For an Asian option with 252 monitoring dates, the algorithm may need repeated exact transitions and integrated-variance inversions.

For a continuously monitored barrier, additional path-bridge logic is needed to determine whether a crossing occurred between sampled dates.

This is a crucial scope boundary.

11. Why the integrated-variance inversion is expensive

Each path can require many evaluations of a conditional characteristic function while numerically inverting a CDF.

Compared with a simple QE step, this is computationally heavy.

Later research therefore develops:

• faster inversion;
• precomputed inverse-CDF tables;
• interpolation/caching;
• series or gamma-expansion approximations;
• Poisson-conditioning schemes.

The exact method is therefore an important benchmark even when a faster approximate method wins a production runtime test.

12. Inputs and outputs

Inputs can include:

• S₀,V₀;
• κ,θ,σ,ρ;
• interest/dividend/carry assumptions;
• observation dates;
• conditional transform evaluator;
• CDF inversion tolerance;
• quadrature/integration settings;
• random seeds;
• path count.

Outputs can include:

• exact-transition S and V samples;
• integrated variance I;
• Monte Carlo price;
• sampling standard error;
• inversion iterations;
• transform-evaluation counts;
• numerical inversion residuals;
• runtime per path;
• comparison with QE/Euler benchmarks.

13. Evidence polarity

Evidence for confidence includes:

• sample moments of V_T match CIR analytic moments;
• simulated Heston European prices agree with independent transform prices;
• results are insensitive to tighter inversion tolerance once sufficiently tight;
• Monte Carlo confidence intervals shrink near N^−1/2;
• no systematic timestep bias appears when the transition interval changes for terminal-only payoffs;
• conditional-CDF inversion residuals remain small;
• simulation reproduces known covariance/martingale conditions.

Evidence against confidence includes:

• inversion failures or non-monotone recovered CDFs;
• European price bias relative to characteristic-function benchmarks;
• runtime dominated by root finding without accuracy benefit;
• rare parameter regimes causing numerical overflow/underflow;
• path-dependent results changing strongly with monitoring discretisation;
• exact variance transitions combined with an approximate stock update incorrectly labelled “exact.”

14. Counterexample: a European call benefits strongly

A European call needs only S_T. Broadie–Kaya can sample one exact transition per path.

There is no need for 100 or 1,000 intermediate variance steps merely to control discretisation bias.

This is a favorable use case.

15. Counterexample: a daily Asian can erase the speed advantage

A daily monitored Asian requires many observation dates.

Repeated integrated-variance inversions can make the exact algorithm much slower than a high-quality QE scheme.

Falsifier: benchmark total RMSE per second, not only mathematical exactness. A slower exact estimator can be inferior for a fixed compute budget.

16. Counterexample: continuous barriers still need crossing logic

Exact transition endpoints do not reveal whether S crossed a barrier between dates.

Falsifier: compare against finer monitoring or a validated bridge method. “Exact Heston transition” does not imply exact continuously monitored barrier pricing.

17. Counterexample: numerical inversion error is too loose

If the CDF root tolerance is coarse, sampled I values can be biased by numerical inversion.

Falsifier: tighten the tolerance by orders of magnitude and test whether prices/moments stabilize. The exact probabilistic construction still requires accurate numerical inversion.

18. Counterexample: Feller violation is mishandled

The CIR variance process can approach zero when parameter combinations do not satisfy the usual strict-positivity condition.

The exact transition law remains the correct object, but numerical implementations must handle low degrees of freedom and mass near zero robustly.

Falsifier: benchmark variance samples in boundary-heavy parameter sets against a trusted noncentral-chi-square implementation and analytic moments.

19. Counterexample: exact simulation is used with stale calibration

Perfect transition simulation under badly calibrated parameters gives a precisely simulated wrong model.

Falsifier: keep simulation validation separate from Heston calibration validation.

20. Broadie–Kaya versus Andersen QE

Andersen QE algorithms own a fast approximate discretisation designed to preserve positivity and match conditional moments well.

Broadie–Kaya owns exact model transitions. QE usually wins on speed for dense path simulation; Broadie–Kaya is a gold-standard bias benchmark and can be compelling for sparse observation schedules.

21. Broadie–Kaya versus Heston transform pricing

Heston calibration algorithms and transform pricers use characteristic functions to value vanilla options directly.

Broadie–Kaya is useful when Monte Carlo structure is needed—for path dependence, portfolio simulation, exposures or a common path engine.

22. Alternatives

Euler/full truncation: simple, flexible, biased at finite timestep.

Andersen QE: fast high-quality approximation.

Higher-order weak schemes: trade complexity for lower discretisation bias.

Glasserman–Kim expansions: accelerate conditional integrated-variance simulation.

Transform/PDE pricing: preferable when payoff structure permits direct deterministic pricing.

23. Weak links

• conditional integrated-variance inversion;
• branch/complex transform evaluation;
• CDF bracketing errors;
• slow repeated transitions;
• variance boundary numerics;
• confusing exact endpoint simulation with exact continuous paths;
• inconsistent carry/dividend terms;
• poor random-number generation;
• benchmarking exactness without runtime.

24. What would falsify confidence?

Confidence should be withdrawn if European benchmark prices remain biased; if CIR endpoint moments are wrong; if inversion residuals are unstable; if numerical results change materially with tighter tolerances; if dense path-dependent problems show worse RMSE-per-second than a validated approximation without a compensating accuracy need; or if the implementation labels partially discretised paths as exact.

25. Verification and update triggers

Preserve the model parameters, noncentral-chi-square implementation, conditional transform formulas, inversion algorithm, integration bounds, tolerances, random seeds and benchmark cases.

Revalidate when:

• Heston parameters enter more extreme regimes;
• the payoff monitoring frequency changes;
• the inversion library changes;
• GPU/vectorized implementation changes numerical order;
• jumps or stochastic rates are added;
• benchmark transform prices change implementation.

26. Primary and high-quality references

• Mark Broadie and Özgür Kaya, Exact Simulation of Stochastic Volatility and Other Affine Jump Diffusion Processes, Operations Research, 2006.
• Columbia Business School, faculty publication record for the Broadie–Kaya exact simulation paper.
• Roger Lord, Remmert Koekkoek and Dick van Dijk, A Comparison of Biased Simulation Schemes for Stochastic Volatility Models, Quantitative Finance.
• Simulation Schemes for the Heston Model with Poisson Conditioning, for later discussion of Broadie–Kaya’s exact construction and its computational bottleneck.

Educational boundary: Exact simulation removes one error source—time-discretisation bias for the sampled transitions. It does not remove Monte Carlo error, numerical inversion error, path-monitoring error or model misspecification.