Simulation of square-root processes made simple: applications to the Heston model

We introduce a simple, efficient and accurate nonnegative preserving numerical scheme for simulating the square-root process. The novel idea is to simulate the integrated square-root process first instead of the square-root process itself. Numerical experiments on realistic parameter sets, applied for the integrated process and the Heston model, display high precision with a very low number of time steps. As a bonus, our scheme yields the exact limiting Inverse Gaussian distributions of the integrated square-root process with only one single time-step in two scenarios: (i) for high mean-reversion and volatility-of-volatility regimes, regardless of maturity; and (ii) for long maturities, independent of the other parameters.