# 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+2) :
sigsqX  = L0[i-1]
sigsqY  = L1[i-1]
sigXY   = L2[i-1]

beta[0,0] = sigsqX(X,X)
beta[1,1] = sigsqY(X,X)
beta[0,1] = sigXY(X,X)
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 =  X + alpha*dt + np.dot(bet_sq, W[j,i-1])*np.sqrt(dt)
X =  X + alpha*dt + np.dot(bet_sq, W[j,i-1])*np.sqrt(dt)

if X > x1vec[-1] or X < x1vec:
broke = True
break

if X > x2vec[-1] or X < x2vec:
broke = True
break

if broke == False :
vix_v_t1.append(X)

for i in np.arange(ind+2,ind+2) :

sigsqX  = L0[i-1]
sigsqY  = L1[i-1]
sigXY   = L2[i-1]

beta[0,0] = sigsqX(X,X)
beta[1,1] = sigsqY(X,X)
beta[0,1] = sigXY(X,X)
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 =  X + alpha*dt + np.dot(bet_sq, W[j,i-1])*np.sqrt(dt)
X =  X + alpha*dt + np.dot(bet_sq, W[j,i-1])*np.sqrt(dt)

if X > x1vec[-1] or X < x1vec:
broke = True
break
if X > x2vec[-1] or X < x2vec:
broke = True
break

if broke == False :
spx_t1.append(X)
for i in np.arange(ind+2,ind+2) :

sigsqX  = L0[i-1]
sigsqY  = L1[i-1]
sigXY   = L2[i-1]

beta[0,0] = sigsqX(X,X)
beta[1,1] = sigsqY(X,X)
beta[0,1] = sigXY(X,X)
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 =  X + alpha*dt + np.dot(bet_sq, W[j,i-1])*np.sqrt(dt)
X = X + alpha*dt + np.dot(bet_sq, W[j,i-1])*np.sqrt(dt)
if X > x1vec[-1] or X<x1vec:
broke = True
break
if X > x2vec[-1] or X<x2vec:
broke = True
break

if broke == False :
spx_t2.append(X)

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 !

• To model path dependent options, you will need a lot more than 20k simulations. Industry standard is at least 100k, with many practitioners doing 200k to 500k. 200k simulations should not take more than half a second even on older hardware
– user50421
Oct 16, 2020 at 22:49
• To aid tieing out alternative code suggestions to your implementation, could you provide initial values and a result associated with those inputs. Oct 21, 2020 at 16:00
• Well yes I would love to be able to perform 500k simulations in 0.5 sec, but this code takes 75 secs to run. Oct 26, 2020 at 17:30
• @wgajate I just did thanks :) Oct 26, 2020 at 17:31
• I don't see any changes to the example you provided. The goal would be to provide an executable/runnable example that ensures the correctness of code improvements offered by the community. Oct 28, 2020 at 20:58