1
$\begingroup$

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 !

$\endgroup$
5
  • $\begingroup$ 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 $\endgroup$
    – user50421
    Oct 16, 2020 at 22:49
  • $\begingroup$ To aid tieing out alternative code suggestions to your implementation, could you provide initial values and a result associated with those inputs. $\endgroup$
    – wgajate
    Oct 21, 2020 at 16:00
  • $\begingroup$ Well yes I would love to be able to perform 500k simulations in 0.5 sec, but this code takes 75 secs to run. $\endgroup$ Oct 26, 2020 at 17:30
  • $\begingroup$ @wgajate I just did thanks :) $\endgroup$ Oct 26, 2020 at 17:31
  • $\begingroup$ 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. $\endgroup$
    – wgajate
    Oct 28, 2020 at 20:58

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy