Monte Carlo Simulation of a Two-Factor Local Volatility Model
Article Quant Q&A · Author: Jordi Lecoch
Summary
The document describes a two-factor Markov diffusion used to price options by Monte Carlo. Its drift and covariance are obtained from a finite-difference solution to a Hamilton–Jacobi–Bellman equation, then interpolated at each simulated state. The simulation advances paths over several maturity intervals and records the relevant state variables for pricing.
The author reports that 20,000 paths take about 75 seconds, while interpolation accounts for only a small fraction of runtime. The code checks whether the covariance matrix is positive definite, repairs it when necessary, and stops paths that leave the interpolation grid. The question asks how to speed up the simulation, but provides no solution or comparison of acceleration methods. Its performance observations describe one implementation and setup, so they do not establish how other implementations or hardware would perform.
Key ideas
- The model simulates two state variables with state-dependent drift and covariance.
- Covariance values are interpolated from grids produced by a finite-difference HJB solution.
- The implementation factors each path sequentially through time and handles multiple maturities in separate intervals.
- Positive-definiteness repairs and stopping paths that leave the grid are part of the simulation logic.
- The reported runtime suggests interpolation is not the main cost in this implementation.
Tags
Full text
# Fast Monte Carlo of Local Volatility Model
# Fast Monte Carlo of Local Volatility Model
I want to compute option prices via a Monte Carlo simulation. The model implemented is a Markov process, following the SDE :
```
d X_t = alpha * dt + beta^(1/2) * d W_t
```
where alpha lives in R^2, beta in R^(2*2) and is Positive semi definite. (^(1/2) -> cholesky) Therefore, X_t is in R^2.
I obtain value for beta and alpha after a resolution of an HJB equation using finite differences on a mesh grid (yielding arrays `sigmasqX`,`sigmasqY` and `sigmaXY`, which I then interpolate to obtain value of beta and alpha at `(t,X_t)`).
My problem is, I then run a Monte Carlo simulation to obtain option prices, to derive implied volatility from, but the simulation takes forever on a very limited 20 000 path simulation.
Here are the broad strokes of my strategy :
- I start with the initial `X_0`
- I generate a normal sample of size `t_size*sample`, where `tsize` is the length of my time vector, and `sample` the number of paths generated
- I interpolate using `spicy.RectBivariateSpline` over the `sigma` arrays. Here are the nested `for` loops :
```
for j in np.arange(sample):
X = np.array([np.log(var.S0), var.x2_0])
broke = False
for i in np.arange(1,ind[0]+2) :
sigsqX = L0[i-1]
sigsqY = L1[i-1]
sigXY = L2[i-1]
beta[0,0] = sigsqX(X[0],X[1])
beta[1,1] = sigsqY(X[0],X[1])
beta[0,1] = sigXY(X[0],X[1])
beta[1,0] = beta[0,1]
if not isPD(beta) :
beta = nearestPD(beta)
alpha = r_rate-q_rate-0.5*beta[0,0]*np.ones(2)
bet_sq = np.linalg.cholesky(beta)
X[0] = X[0] + alpha[0]*dt + np.dot(bet_sq, W[j,i-1])[0]*np.sqrt(dt)
X[1] = X[1] + alpha[1]*dt + np.dot(bet_sq, W[j,i-1])[1]*np.sqrt(dt)
if X[0] > x1vec[-1] or X[0] < x1vec[0]:
broke = True
break
if X[1] > x2vec[-1] or X[1] < x2vec[0]:
broke = True
break
if broke == False :
vix_v_t1.append(X[1])
for i in np.arange(ind[0]+2,ind[1]+2) :
sigsqX = L0[i-1]
sigsqY = L1[i-1]
sigXY = L2[i-1]
beta[0,0] = sigsqX(X[0],X[1])
beta[1,1] = sigsqY(X[0],X[1])
beta[0,1] = sigXY(X[0],X[1])
beta[1,0] = beta[0,1]
if not isPD(beta) :
beta = nearestPD(beta)
alpha = r_rate-q_rate-0.5*beta[0,0]*np.ones(2)
bet_sq = np.linalg.cholesky(beta)
X[0] = X[0] + alpha[0]*dt + np.dot(bet_sq, W[j,i-1])[0]*np.sqrt(dt)
X[1] = X[1] + alpha[1]*dt + np.dot(bet_sq, W[j,i-1])[1]*np.sqrt(dt)
if X[0] > x1vec[-1] or X[0] < x1vec[0]:
broke = True
break
if X[1] > x2vec[-1] or X[1] < x2vec[0]:
broke = True
break
if broke == False :
spx_t1.append(X[0])
for i in np.arange(ind[1]+2,ind[2]+2) :
sigsqX = L0[i-1]
sigsqY = L1[i-1]
sigXY = L2[i-1]
beta[0,0] = sigsqX(X[0],X[1])
beta[1,1] = sigsqY(X[0],X[1])
beta[0,1] = sigXY(X[0],X[1])
beta[1,0] = beta[0,1]
if not isPD(beta) :
beta = nearestPD(beta)
alpha = r_rate-q_rate-0.5*beta[0,0]*np.ones(2)
bet_sq = np.linalg.cholesky(beta)
X[0] = X[0] + alpha[0]*dt + np.dot(bet_sq, W[j,i-1])[0]*np.sqrt(dt)
X[1] = X[1] + alpha[1]*dt + np.dot(bet_sq, W[j,i-1])[1]*np.sqrt(dt)
if X[0] > x1vec[-1] or X[0]<x1vec[0]:
broke = True
break
if X[1] > x2vec[-1] or X[1]<x2vec[0]:
broke = True
break
if broke == False :
spx_t2.append(X[0])
spot = [spx_t1, spx_t2, vix_v_t1]
return spot
```
A few things to note :
- I check whether `beta`, and adjust it if need be, because of numerical approximation it sometimes has a tiny negative eigenvalue and corrupts the cholesky decomposition
- There are 3 linear `for` loops to walk through the entire time set because there a three maturities for options, it is easier to manage that way (not impactful here)
- Sometimes the value for `X_t` wanders outside of the grid, and breaks the simulation (because of extrapolation errors), so I implemented all these verifications (`if broke : break`)
- The `L0`, `L1` and `L2` variables are python `list` type, that hold the interpolation of the grid, at each date of the time vector.
I know why this code is slow :
- many `if` statements slow it down a lot
- finding nearest PD for beta can be costly
- using local volatility means, I can't parallelize according to NumPy's rules, the Monte Carlo process
I checked, and the interpolation is negligible (0.5 sec on average). But this path simulation takes 75 sec for 20 000 paths...
My question : is there a way to speed up the whole process ? I can't imagine that I am stuck with such poor execution times... (comparable to solving backward pricing PDEs which are notorious for being expansive computation wise)
Thanks !Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.