# 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?

• What is the forward at expiry? Usually I would guess it was close to the bottom of the vol smile, but given your large negative rho it could be some way to the left hand side. In that case, you would need to extend your vol smiles down to lower strikes and you might see your pdf peak and then fall. – StackG Sep 6 at 12:14
• @StackG I have used the Spot S=1 such that I can think of the inputs K as the moneyness. I also just tried extrapolating such that the strike range was from 0.2-1.3 and changed the rho to -0.170152. This just made the left tail rise even higher! – Klein Sep 6 at 12:23
• If you have the right parameters (see @JohnDoe's comment), you can evaluate the characteristic function and invert it numerically to backout the density, without any partial derivatives – Alex Sep 6 at 13:12
• @Klein I think JohnDoe meant is that $\theta$ is the long-term of the variance. A reasonable long-term mean of volatility is something like 20%. This however means that $\theta=0.2^2=0.04$, i.e. much smaller than your value. Your number suggests that the long-term volatility is $\sqrt{\theta}\approx60\%$. That's a bit high perhaps? – Alex Sep 6 at 13:33
• @Klein I meant that the Heston closed form solution is $C=SP_1-Ke^{-rT}P_2$. This $P_2$ is the risk-neutral probability of $S_T$ being greater or equal than $K$ and an integral in terms of the characteristic function (Fourier transform of the RN density). You can obtain the RN density this way. – Alex Sep 6 at 13:35

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

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()