Recovering a Risk-Neutral Density from Heston Option Prices

# Computing the Probability Density Function (PDF) for the Heston model


# Computing the Probability Density Function (PDF) for the Heston model












I am trying to compute the PDF for the Heston model using the Breeden Litzenberger formula.

I have calculated the the Heston implied volatilities for a strike range (which i have interpolated using cubic spline interpolation) using python:

In order to get the PDF, I am using the Breeden and Litzenberger forumla:

For the derivative of the implied volatility w.r.t the strike price, I have used numerical differentiation. Doing this together with the Breeden and Litzenberger forumla I get a PDF looking like this:

This does not look right. I have used the same approach but for the SABR model which results in a good looking PDF. I am pretty sure my code is correct so I was wondering if there is something I am missing about the PDF of the Heston model??

For the Heston model I have used the following inputs (S=1 and a strike range K$\in$(0.8,1.3):

I am actually also not sure what parameter to use for the Black Scholes sigma, but it should not make or break the above PDF?

## Answer by StackG (score 3)

https://quant.stackexchange.com/a/57829

I can't quite even re-create your vol smile... when I plug in the parameters you've provided (at $\tau = 0.12$) I get a downward sloping vol smile that doesn't have a minimum at the strikes I looked at

I then backed out the options prices at each of a close-up grid of strikes and calculated the curvature of the prices, which is very close to the rn pdf (just need to correct by a factor of the dcf, which is close to 1 for such short times), and it looks roughly as expected

I've attached my code below, it should be very easy for you to play with the parameters and try to work out what is going wrong in your script (if you share the code, we might be able to help out more)

```
import QuantLib as ql
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

# Your parameters
tau = 0.1219; r = 0.0457; sigma = 0.4433; rho = -0.6175; nu = 0.474; theta = 0.3737; kappa = 1.042

today = ql.Date(1, 9, 2020)
expiry_date = today + ql.Period(int(365*tau), ql.Days)

# Setting up discount curve
risk_free_curve = ql.FlatForward(today, r, ql.Actual365Fixed())
flat_curve = ql.FlatForward(today, 0.0, ql.Actual365Fixed())
riskfree_ts = ql.YieldTermStructureHandle(risk_free_curve)
dividend_ts = ql.YieldTermStructureHandle(flat_curve)

# Setting up a Heston model
spot = 1

# I guess this is the correct mapping?
v0, sigma = nu, sigma

heston_process = ql.HestonProcess(riskfree_ts, dividend_ts, ql.QuoteHandle(ql.SimpleQuote(spot)), v0, kappa, theta, sigma, rho)
heston_model = ql.HestonModel(heston_process)
heston_handle = ql.HestonModelHandle(heston_model)
heston_vol_surface = ql.HestonBlackVolSurface(heston_handle)

# Now doing some pricing and curvature calculations
strikes = np.arange(0.5, 1.6, 0.01)
vols = [heston_vol_surface.blackVol(tau, x) for x in strikes]

option_prices = []

for strike in strikes:
    option = ql.EuropeanOption( ql.PlainVanillaPayoff(ql.Option.Call, strike), ql.EuropeanExercise(expiry_date))

    heston_engine = ql.AnalyticHestonEngine(heston_model)
    option.setPricingEngine(heston_engine)

    option_prices.append(option.NPV())

prices = pd.DataFrame([strikes, option_prices]).transpose()
prices.columns = ['strike', 'price']
prices['curvature'] = (-2 * prices['price'] + prices['price'].shift(1) + prices['price'].shift(-1)) / 0.01**2

# And plotting...
fig = plt.figure()
ax = fig.add_subplot(111)
ax2 = ax.twinx()

ax.plot(strikes, vols, label='Black Vols')
ax2.plot(strikes, option_prices, label='Option Prices', color='orange')
ax2.plot(prices['strike'], prices['curvature'], label='dC/dK (~pdf)', color='purple')

ax.legend(loc="lower left")
ax2.legend(loc="upper right")
ax.grid()
```

## Answer by Nikolai Zaitsev (score 0)

https://quant.stackexchange.com/a/63116

First, couple of corrections (I am not sure, just guessing):

- $X_T$ - is it strike or price of forward underlying?

Let it be strike, $X$ and the underlying is $S_t$ with forward: $Fwd_T=S_t/D_T$.

- Breeden and Litzenberger formula:

No, B&L formula is this: $PDF(S_T)=D_T\cdot\frac{d^2C(X)}{dX^2}$, where $D_T$ is discount factor.

- Finally, my recipe to get RND has two steps:

a. Use a standard B-S function to price European Call with $\sigma(X)$, i.e $Call(S_t,X) = D\cdot (Fwd \cdot N(d_1(\sigma(X))) - K \cdot N(d_2(\sigma(X))))$. Or use Put pricing.

b. Assume smiled volatility and apply B&L formula to B-S Call (or Put) to calculate RND.

Very likely, to get your formula Taylor (2007) twice differentiated B-S by strike and made the same just in one step. Then my recipe can be considered as a cross-check of that derivation (maybe there is typo in formula?).

Is it Heston model or SABR or something else it does not matter much because you start with parameterized volatility shape (smile).