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?

  • 3
    $\begingroup$ What options are you trying to price? QuantLib-Python has many methods for pricing options under local volatility. If you can tell us a bit more about your problem I'll provide some sample code. From there you can dig in to the implementation and understand what is going on 'under the hood'. $\endgroup$
    – StackG
    Commented Jul 19, 2021 at 23:45
  • $\begingroup$ Hi @StackG. I am trying to price CDX options which are vanilla options. I have followed an algorithm that helps interpolate prices across strike and tenor to help find local volatility. This interpolation was done so that the derivatives that define local vol according to the dupire equation make sense. Let me know if you need to learn more $\endgroup$ Commented Jul 20, 2021 at 12:55

1 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)
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: Not a very interesting surface, but never mind!!

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)

fd_price = european_option.NPV()

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)

mc_price = european_option.NPV()

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]
        surfaces = vol_surface

        if type(funct) != list:
            functs = [funct] * len(surfaces)
            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_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)
    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))

Note the slightly funky behaviour near to the edges of the interpolation region

  • $\begingroup$ Thanks so much for your answer! QuantLib is very helpful. I looked through QuantLib and it looks like you even have a calculator to find the volatility surface using Andreeson and Huge's method which is what I had been doing $\endgroup$ Commented Jul 21, 2021 at 13:20
  • 1
    $\begingroup$ That is correct - let me know if you would like and I can add a snippet above to do the QuantLib interpolation using the AH method $\endgroup$
    – StackG
    Commented Jul 21, 2021 at 13:25
  • $\begingroup$ I would appreciate that if you wouldn't mind $\endgroup$ Commented Jul 21, 2021 at 13:27
  • 2
    $\begingroup$ Added, you might need to play with the interpolation knots a little bit to get a smoother function near the edge of the interpolation region! $\endgroup$
    – StackG
    Commented Jul 21, 2021 at 13:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.