# Correct Monte Carlo simulation of local volatility models

I am using Monte Carlo simulation to evolve the following SDE over a grid of timepoints $$0,t_1,...,t_N$$.

$$\begin{equation} dS(t)=\sigma(t, S(t))dw(t) \end{equation}$$

Here $$\sigma(t_i,S(t_i)), i=1,...,N$$ has been previously determined from Dupire's formula using European options expiring at $$t_i$$. An Euler discretisation of the SDE from $$t_{i}$$ to $$t_{i+1}$$ gives:

$$\begin{equation} S(t_{i+1}) = S(t_i) + \sigma(t_i, S(t_{i})) \sqrt{\Delta t_i}\zeta(0,1) \end{equation}$$ where $$\zeta(0,1)$$ is a standard Gaussian variate and $$\Delta t_i=t_{i+1}-t_{i}$$. That is, given $$S(t_i)$$, the volatility for the period $$t_i$$ to $$t_{i+1}$$, is interpolated from the local volatility surface at time $$t_i$$, i.e. $$\sigma(t_i,S(t_i))$$ and taken flat over the period $$t_i$$ to $$t_{i+1}$$.

But from a theoretical perspective, this feels wrong. $$\sigma(t_i,S(t_i))$$ comes from Markov projection / Gyongy theorem arguments, i.e. it is the volatility to be taken from time $$0$$ to time $$t_i$$ in order to match the true distribution of $$S(t_i)$$ using an approximate SDE. So, in order to match the distribution of $$S(t_i)$$, I would expect $$\sigma(t_i,S(t_i))$$ to be applied over the entire period from today (time $$0$$) to $$t_i$$, not just from $$t_i$$ to $$t_{i+1}$$. More worryingly, given that in the previous time step, $$t_{i-1}$$ to $$t_i$$, we use volatility from $$\sigma(t_{i-1},S(t_{i-1}))$$ surface, I cannot see how we would ever recover the time $$t_i$$ distribution of $$S(t_i)$$ with this approach.

So I can't help thinking that when simulating from $$t_i$$ to $$t_{i+1}$$, instead of interpolating only from $$\sigma(t_i, S(t_i))$$, we should be interpolating from both $$\sigma(t_i,S(t_i))$$ and $$\sigma(t_{i+1},S(t_{i+1}))$$ to deduce a forward volatility, $$\hat{\sigma}(t_i,S(t_i))$$. Something like (and this is only a guess)

$$\begin{equation} \hat{\sigma}(t_i,S(t_i)) = \sqrt{\frac{\sigma(t_{i+1}, S(t_{i}))^2 t_{i+1}- \sigma(t_i, S(t_i))^2 t_i}{\Delta t_i}} \end{equation}$$

Since all papers I've read only interpolate from $$\sigma(t_i, S(t_i))$$, I am fairly certain I am wrong, but I cannot see why. Any help would be greatly appreciated.

Thanks!

• Two things you might be interested to read up on which are related to this. The first is concerning numerical schemes, there are certain correction based iterative schemes, which may be appropriate in your setting (cf Kleoden and Platten). The second would be different types of stochastic processes. Rather than the usual Ito process construction, see the construction of Stratonovich SDEs/Integrals, where the average of the process is used in the integrand, which seems to be what you're worried about with your numeric scheme. Feb 18 at 13:32

### [I think] the problem is with the SDE, rather than the numerical scheme

At a glance, and as I commented, I think the issue you are coming up against stems more from the underlying SDE rather than the numerical approximation scheme. To explain this a bit more, let's quickly revisit the SDE $$\mathrm{d}S_t = \sigma(t, S_t) \,\mathrm{d}W_t$$ with some given initial condition (e.g. $$S_t = S_0$$ at $$t = 0$$), which is really the solution the integral equation $$S_t = S_0 + \int_0^{t} \sigma(u, S_u) \,\mathrm{d}W_u.$$ The integral as you are using it is interpreted as an Ito integral, defined as the limit $$\int_0^{t} \sigma(u, S_u) \,\mathrm{d}W_u = \lim_{\delta t \to 0} \sum_i \sigma(t_i, S_{t_i}) (W_{t_{i+1}} - W_{t_i}).$$

Recall that the Euler-Maruyama scheme is simply approximating this integral by the same thing appearing in the limit. Thus both the numerical scheme and the underlying stochastic process both suffer from the deficiency your wrestling with. Consequently, I wouldn't blame the numerical scheme for the perceived modelling shortcoming.

### Stratonovich integrals seem a better fit

It seems like what you might be after is a process where the average volatility is seen, especially if you plan on capturing this with a numerical scheme. One way of capturing this is by instead describing the process using a Stratonovich integral/SDE, where instead we have $$\int_0^{t} \sigma(u, S_u) \circ \,\mathrm{d}W_u = \lim_{\delta t \to 0} \sum_i \dfrac{\sigma(t_i, S_{t_i}) + \sigma(t_{i+1}, S_{t_{i+1}})}{2} (W_{t_{i+1}} - W_{t_i}),$$ which would then be written as the SDE $$\mathrm{d}S_t = \sigma(t, S_t) \circ \,\mathrm{d}W_t.$$

You can then see there appears to be a much nicer match between the dynamics you're wishing to model, the stochastic process, and the numerical approximation scheme. See Kloeden and Platen for details about how to approximate solutions to Stratonovich SDEs.

### Some comments on numerical schemes, financial modelling, and stochastic processes

It is worth noting that this is all just a model for financial markets rather than strict mathematical analysis (hence why it's on Quantitative Finance rather than Maths Overflow). In the real world of finance, it is arguably equally (or even more) appropriate to use discrete time models. Furthermore, frequently in mathematical modelling you need to decide what the underlying model is, and whether it is fundamentally discrete or continuous. While the continuous is usually a bit more attractive analytically, if the core model is discrete, then approximating it by a continuous process and then approximating that again by a discrete process is needlessly putting the cart before the horse.

Nonetheless, if you are sure the problem is the shortcoming of the numerical scheme, then there are countless other numeric schemes, many of which you may find more attractive or appropriate for your purposes. Just pick up Kloeden and Platen, and have a look at iterative methods, correction methods, higher order methods, methods for Stratonovich integrals, etc.

• Great answer, +1, but shouldn't the integrals go to $t$ instead of $W_t$ (upper integral limit)? Like the SDE $\text{d}S_t=\sigma(t,S_t)\text{d}W_t$ means $$S_t-S_0=\int_0^\color{red}{t} \sigma(u,S_u)\text{d}W_u$$ Feb 18 at 16:52
• @Kevin correct, will amend it shortly. Feb 18 at 16:56
• Thank you @oliversm, very comprehensive, useful point about the use of Stratonovich integrals, and more sophisticated schemes as found in Kloeden Platen. At the same time, I think there's also a modelling aspect to this question, as I will outline now in an answer. Feb 18 at 17:06

Nothwithstanding @oliversm's excellent response to my question, I think there's also a question on how to interpret $$\sigma(t,S(t))$$ in a local vol model.

Given that $$\sigma(t,S(t))$$ comes from a Markovian projection of a complicated stochastic process on to a simpler one, my interpretation (perhaps wrong) is that it should represent a single volatility (conditional upon knowing $$S(t_i)$$) to be used all the way from time $$0$$ to $$t_i$$ in order to recovery the true distribution of $$S(t_i)$$.

In other words, if I wanted to simulate $$S(t_i)$$ from time $$0$$ to $$t_i$$ and recover the distribution of the complex model, I would, conditional upon knowing $$S(t_i)$$, interpolate the vol from the $$\sigma(t_i,S(t_i))$$ function that I calibrated to European options using the Dupire local volatility approach.

The problem comes with then continuing the simulation from $$t_i$$ to $$t_{i+1}$$. The general consensus in the literature is to interpolate the volatility from $$\sigma(t_{i+1},S(t_{i+1}))$$. But how can this be correct? $$\sigma(t_{i+1},S(t_{i+1}))$$ also results from a Markov projection, (presumably) also meaning that when we interpolate this function (conditional upon knowing $$S(t_{i+1})$$), the resulting volatility should be used from $$0$$ to $$t_{i+1}$$. But we are only using it from $$t_i$$ to $$t_{i+1}$$, i.e. for only part of the timeline. How can we then possibly expect to recover the time $$t_{i+1}$$ distribution of $$S(t_{i+1})$$? This is independent of the numerical scheme we use.

Once we know what volatility is to be used at which time point in the simulation, then we can focus on the available numerical schemes to approximate its evolution in time, taking on board @oliversm's excellent guidance.

I expect that there's a flaw (or perhaps more than one) in my arguments, I'd be grateful to anyone who can shed some light.

Thanks!