Understanding Monte Carlo to solve option price with local volatility

Question

I have read this question pricing using dupire local volatility model which seems to have an answer from here https://www.csie.ntu.edu.tw/~d00922011/python/cases/LocalVol/DUPIRE_FORMULA.PDF

Both of these simply say use monte carlo or finite difference methods once we have the local volatility in order to price options after discussing slightly more in depth how to find the local volatility itself.

I am looking for a more granular description of this last step which is finding the option price itself.

My Attempt

At the point where I am stuck, I have every value needed to solve the dupire equation: $dS_t=\mu(S_t,t) S_tdt+\sigma(S_t;t) S dW$ where $S_t|_{t=0} = S_0$. Assume I have well defined drift and local volatility terms and know the spot price as well. From here, I am able to use a monte carlo solution in order to find $S_t$ at all known $t$. However, $S_t$ is not the option price so I am wondering how the person writing the paper/asking the question linked were able to find these prices.

From my (limited) understanding of monte carlo methods and stochastic differential equations, $S_t$ found using the MC methods is a realization of the random variable that is the underlying price of this equation.

I know that if I find the distribution of this random variable, $p(x)$, at each time step I could price my option at each time step $\big($$C(K,t) = E[(S_{t}-K)^+] = \int_{K}^\infty p(x)(S_{t}-K)dx$$\big)$.

I've read some descriptions on how to solve Ito and Stratonovich integrals using monte carlo methods but these simply tell me how to evaluate these integrals which I already know how to do; I want to find the probability. What is the best way to do this? Are there better methods or python packages I should use instead?

Answer 1

I'm going to go through a coded example to show how you might attempt this, using the python port of the QuantLib library. It will all seem a little mechanical, but hopefully it is instructive. There is quite a bit of setup code required (specifying the interest rates, spot etc., and also a utility function that I use for plotting surfaces), I've pushed this all to the bottom of my answer but you'll need it in the top of your script/notebook.

To build a local vol model, I'm going to need an implied vol surface from market data. I don't know what data you're using so as a very trivial example, I'm instead going to simulate a sequence of SABR smiles at different maturities, and assume that's the market implied vol surface. The code to do that is:

strikes = [70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 130.0]
expirations = [ql.Date(1, 7, 2021), ql.Date(1, 9, 2021), ql.Date(1, 12, 2021), ql.Date(1, 6, 2022)]
vol_matrix = ql.Matrix(len(strikes), len(expirations))

# params are sigma_0, beta, vol_vol, rho
sabr_params = [[0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6]]

for j, expiration in enumerate(expirations):
    for i, strike in enumerate(strikes):
        tte = day_count.yearFraction(today, expiration)
        vol_matrix[i][j] = ql.sabrVolatility(strike, spot, tte, *sabr_params[j])

implied_surface = ql.BlackVarianceSurface(today, calendar, expirations, strikes, vol_matrix, day_count)
implied_surface.setInterpolation('bicubic')
vol_ts = ql.BlackVolTermStructureHandle(implied_surface)

process = ql.BlackScholesMertonProcess(ql.QuoteHandle(ql.SimpleQuote(spot)), dividend_ts, flat_ts, vol_ts)

plot_vol_surface(implied_surface, plot_years=np.arange(0.1, 1, 0.1))
plt.title("Implied Vol")

This also plots the vol surface so we can eyeball it:

Now I'm going to try pricing a vanilla option at 6 month maturity, I'll use a finite difference pricing engine and a monte carlo pricing engine. First, let's query the surface and see what vol we expect:

implied_surface.blackVol(0.5, 90)

returns 0.0787116605540102, or 7.87% IV

In QuantLib, we set up the option and the pricing engines independently. First, the option:

strike = 90.0
option_type = ql.Option.Call

maturity = today + ql.Period(6, ql.Months)
europeanExercise = ql.EuropeanExercise(maturity)

payoff = ql.PlainVanillaPayoff(option_type, strike)

european_option = ql.VanillaOption(payoff, europeanExercise)

and now I'll set up and price using each of two pricing engines (an invaluable resource on these can be found at the ReadTheDocs page maintained by one of the other posters here).

First the FD pricer (note we HAVE TO turn on the final parameter here, which specifies use of Dupire local vol...):

tGrid, xGrid = 3000, 400
fd_engine = ql.FdBlackScholesVanillaEngine(process, tGrid, xGrid, 0, ql.FdmSchemeDesc.Douglas(), True)

european_option.setPricingEngine(fd_engine)
fd_price = european_option.NPV()
fd_price

price is 7.699526511916783

And now the MC pricer (for MC in local vol, a sufficient number of steps is vital, experiment with reducing this and see what happens to the price):

steps = 36
rng = "lowdiscrepancy" # could use "pseudorandom"
numPaths = 2**15

mc_engine = ql.MCEuropeanEngine(process, rng, steps, requiredSamples=numPaths)

european_option.setPricingEngine(mc_engine)
mc_price = european_option.NPV()
mc_price

price is 7.684257656799339

And as a final comparison, we'll calculate the vols that each of these prices imply, to compare with the vol from the surface above:

# Compare calculated volatility to vol from the surface
european_option.impliedVolatility(fd_price, process), european_option.impliedVolatility(mc_price, process)

returns (0.07872046080861067, 0.07707949492498713), both of which are pretty close to the vol we expected.

Hopefully this is helpful! There are hundreds of little assumptions that I've made in the answer however that might not be appropriate to your situation... let me know if you want me to shift anything at all.

Boilerplate code required to run the above snippets:

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
import QuantLib as ql

calc_date = ql.Date(21, ql.December, 2020)

def plot_vol_surface(vol_surface, plot_years=np.arange(0.1, 3, 0.1), plot_strikes=np.arange(70, 130, 1), funct='blackVol'):
    if type(vol_surface) != list:
        surfaces = [vol_surface]
        functs = [funct]
    else:
        surfaces = vol_surface

        if type(funct) != list:
            functs = [funct] * len(surfaces)
        else:
            functs = funct

    fig = plt.figure(figsize=(10, 6))
    ax = fig.gca(projection='3d')
    X, Y = np.meshgrid(plot_strikes, plot_years)
    Z_array, Z_min, Z_max = [], 100, 0

    for surface, funct in zip(surfaces, functs):
        method_to_call = getattr(surface, funct)

        Z = np.array([method_to_call(float(y), float(x)) 
                      for xr, yr in zip(X, Y) 
                          for x, y in zip(xr, yr)]
                     ).reshape(len(X), len(X[0]))

        Z_array.append(Z)
        Z_min, Z_max = min(Z_min, Z.min()), max(Z_max, Z.max())

    # In case of multiple surfaces, need to find universal max and min first for colourmap
    for Z in Z_array:
        N = (Z - Z_min) / (Z_max - Z_min)  # normalize 0 -> 1 for the colormap
        surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0.1, facecolors=cm.coolwarm(N))

    m = cm.ScalarMappable(cmap=cm.coolwarm)
    m.set_array(Z)
    plt.colorbar(m, shrink=0.8, aspect=20)
    ax.view_init(30, 300)

# Simple World State
spot = 100
rate_dom = 0.02
rate_for = 0.05

today = ql.Date(1, 6, 2021)

calendar = ql.NullCalendar()
day_count = ql.Actual365Fixed()

# Set up some risk-free curves
riskFreeCurveDom = ql.FlatForward(today, rate_dom, day_count)
riskFreeCurveFor = ql.FlatForward(today, rate_for, day_count)

flat_ts = ql.YieldTermStructureHandle(riskFreeCurveDom)
dividend_ts = ql.YieldTermStructureHandle(riskFreeCurveFor)

Update - Andraesen-Huge interpolation (nb. I changed the boilerplate graphing function above as well slightly):

strikes = np.linspace(70, 130, 10)
expirations = [ql.Date(1, 7, 2021), ql.Date(1, 8, 2021), ql.Date(1, 9, 2021), ql.Date(1, 12, 2021), ql.Date(1, 6, 2022)]

# params are sigma_0, beta, vol_vol, rho
sabr_params = [[0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6],
               [0.4, 0.6, 0.4, -0.6]]

# Now try Andraeson Huge calibration
calibration_set = ql.CalibrationSet()

for i, strike in enumerate(strikes):
    for j, expiration in enumerate(expirations):
        tte = day_count.yearFraction(today, expiration)
        payoff = ql.PlainVanillaPayoff(ql.Option.Call, strike)
        exercise = ql.EuropeanExercise(expiration)
        vol = ql.sabrVolatility(strike, spot, tte, *sabr_params[j])

        calibration_set.push_back((ql.VanillaOption(payoff, exercise), ql.SimpleQuote(vol)))

ah_interpolation = ql.AndreasenHugeVolatilityInterpl(calibration_set, \
                              ql.QuoteHandle(ql.SimpleQuote(spot)), flat_ts, dividend_ts)
ah_surface = ql.AndreasenHugeVolatilityAdapter(ah_interpolation)

plot_vol_surface([implied_surface, ah_surface], plot_years=np.arange(0.2, 1, 0.05), plot_strikes=np.arange(80, 120, 2))