The Feynman–Kac analysis enables us to define a risk neutral probability in which we can price options. Let \(f_\mathbb{Q}(S_T)\) denote the terminal risk neutral (\(\mathbb{Q}\)-measure) probability at time \(T\) , and let \(F_\mathbb{Q}(S_T)\) denote the cumulative probability. A European call option at time \(t\), expiring at \(T\) , with strike price \(K\), is priced (\(\tau = T – t\)):
\(\large C(K, \tau) = \large E_\mathbb{Q}[e^{-r\tau} (S_T – K)^+]\)
\(\large C(K, \tau) = e^{-r\tau} \int^\infty_K (S_T – K)f_\mathbb{Q}(S_T)dS_T\)
In the case where \(f_\mathbb{Q}\) is the log-normal density and volatility is constant with respect to \(K\), this yields the Black–Scholes formula.
Under basic no-arbitrage restrictions, we can consider more general densities than the log-normal for the underlying. Breeden and Litzenberger (1978) show that the first derivative is a function of the cumulative distribution.
\(\large \frac{\delta C}{\delta K}|{K=S_T} = -e^{-r\tau}(1 – F\mathbb{Q}(S_T))\)
The second derivative then extracts the density
\(\large \frac{\delta^2 C}{\delta K^2}|{K=S_T} = e^{-r\tau}f \mathbb{Q}(S_T)\)
References:
Bruce Mizrach, Chapter 35: Estimating Implied Probabilities from Option Prices and the Underlying
http://econweb.rutgers.edu/mizrach/pubs/%5B42%5D-2010_Handbook.pdf
https://www.repository.utl.pt/bitstream/10400.5/3372/1/Tese_andre_novo_4.0.pdf
import numpy as np import scipy as sc import pandas as pd from scipy import stats import matplotlib.pyplot as plt
European Call Option with Stochastic Volatility
Risk-Neutral Dynamics:
Underlying process:
\(\large dS_{t} = r S_{t} dt + \sqrt{v_{t}} S_{t} dW^\mathbb{Q}_{s,t}\)
Variance process:
\(\large dv_{t} = \kappa (\theta – v_{t})dt + \sigma \sqrt{v_{t}} dW^\mathbb{Q}_{v,t}\)
- \(dW^\mathbb{Q}{S,t}\) Brownian motion of asset price
- \(dW^\mathbb{Q}{v,t}\) Brownian motion of asset’s price variance
- \(\rho^\mathbb{Q}\) correlation between \(dW^\mathbb{Q}{S,t}\) and \(dW^\mathbb{Q}{v,t}\)
# Initialise parameters S0 = 100.0 # initial stock price K = 150.0 # strike price tau = 1.0 # time to maturity in years r = 0.06 # annual risk-free rate # Heston dependent parameters kappa = 3 # rate of mean reversion of variance under risk-neutral dynamics theta = 0.20**2 # long-term mean of variance under risk-neutral dynamics v0 = 0.20**2 # initial variance under risk-neutral dynamics rho = 0.98 # correlation between returns and variances under risk-neutral dynamics sigma = 0.2 # volatility of volatility lambd = 0 # risk premium of variance 2*kappa*theta, ' > ',sigma**2
Breeden-Litzenberger Formula
Let’s determine the risk-neutral probability distribution function (pdf), \(f _\mathbb{Q}\).
\(\large f _\mathbb{Q}(K, \tau) = e^{r \tau} \frac{\delta^2C(K,\tau)}{\delta K^2}\)
Use 2nd order finite difference approximation:
\(\large f_\mathbb{Q}(K, \tau) \approx e^{r\tau} \frac{C(K+\Delta_K,\tau) – 2C(K,\tau) + C(K-\Delta_K,\tau)}{(\Delta_K)^2}\)
Heston Model Characteristic Equation:
\(\large C(S_0, K, v_0, \tau) = \frac{1}{2}(S_0 – Ke^{-r \tau}) + \frac{1}{\pi} \int^\inf_0 \Re [ e^{r\tau} \frac{\varphi(\phi-i)}{i\phi K^{i\phi}} – K\frac{\varphi(\phi)}{i\phi K^{i\phi}} ] d\phi\)
\(\varphi(X_0, K, v_0,\tau; \phi) = e^{r \phi i \tau} S^{i \phi}[\frac{1-ge^{d\tau}}{1-g}]^{\frac{-2a}{\sigma^2}} exp[\frac{a \tau}{\sigma^2} (b_2 -\rho\sigma \phi i + d) + \frac{v_0}{\sigma^2}(b_2 -\rho\sigma \phi i + d)[\frac{1-e^{d\tau}}{1-ge^{d\tau}}]]\)
where:
\(d = \sqrt{(\rho\sigma \phi i – b)^2 + \sigma^2 (\phi i + \phi^2)}\)
\(g = \frac{b -\rho\sigma \phi i + d}{b -\rho\sigma \phi i – d}\)
\(a = \kappa \theta\)
\(b = \kappa + \lambda\)
https://www.maths.univ-evry.fr/pages_perso/crepey/Finance/051111_mikh%20heston.pdf
http://web.math.ku.dk/~rolf/teaching/ctff03/Gatheral.1.pdf
def heston_charfunc(phi, S0, v0, kappa, theta, sigma, rho, lambd, tau, r): # constants a = kappa*theta b = kappa+lambd # common terms w.r.t phi rspi = rho*sigma*phi*1j # define d parameter given phi and b d = np.sqrt( (rspi - b)**2 + (phi*1j+phi**2)*sigma**2 ) # define g parameter given phi, b and d g = (b-rspi+d)/(b-rspi-d) # calculate characteristic function by components exp1 = np.exp(r*phi*1j*tau) term2 = S0**(phi*1j) * ( (1-g*np.exp(d*tau))/(1-g) )**(-2*a/sigma**2) exp2 = np.exp(a*tau*(b-rspi+d)/sigma**2 + v0*(b-rspi+d)*( (1-np.exp(d*tau))/(1-g*np.exp(d*tau)) )/sigma**2) return exp1*term2*exp2 def heston_price_rec(S0, K, v0, kappa, theta, sigma, rho, lambd, tau, r): args = (S0, v0, kappa, theta, sigma, rho, lambd, tau, r) P, umax, N = 0, 100, 650 dphi=umax/N #dphi is width for j in range(1,N): # rectangular integration phi = dphi * (2*j + 1)/2 # midpoint to calculate height numerator = heston_charfunc(phi-1j,*args) - K * heston_charfunc(phi,*args) denominator = 1j*phi*K**(1j*phi) P += dphi * numerator/denominator return np.real((S0 - K*np.exp(-r*tau))/2 + P/np.pi) strikes = np.arange(60, 180, 1.0) option_prices = heston_price_rec(S0, strikes, v0, kappa, theta, sigma, rho, lambd, tau, r) option_prices
Use 2nd order finite difference approximation:
\(\large f_\mathbb{Q}(K, \tau) \approx e^{r\tau} \frac{C(K+\Delta_K,\tau) – 2C(K,\tau) + C(K-\Delta_K,\tau)}{(\Delta_K)^2}\)
prices = pd.DataFrame([strikes, option_prices]).transpose() prices.columns = ['strike', 'price'] prices['curvature'] = (-2 * prices['price'] + prices['price'].shift(1) + prices['price'].shift(-1)) / 1**2 # And plotting... fig = plt.figure() ax = fig.add_subplot(111) plt.ylabel('Call Price ($)') ax2 = ax.twinx() ax.plot(strikes, option_prices, label='Option Prices') ax2.plot(prices['strike'], prices['curvature'], label='$d^2C/dK^2 (\sim pdf)$', color='orange') ax.legend(loc="center right") ax2.legend(loc="upper right") plt.xlabel('Strikes (K)') plt.ylabel('$f_\\tau(K)$') plt.title('Risk-neutral PDF, $f_\mathbb{Q}(K, \\tau)$')
Let’s compare with QuantLib Heston Model
References:
Stack Exchange Question, Computing the Probability Density Function (PDF) for the Heston model
https://quant.stackexchange.com/questions/57827/computing-the-probability-density-function-pdf-for-the-heston-model
import QuantLib as ql today = ql.Date(28, 5, 2022) 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 heston_process = ql.HestonProcess(riskfree_ts, dividend_ts, ql.QuoteHandle(ql.SimpleQuote(S0)), 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 vols = [heston_vol_surface.blackVol(tau, x) for x in strikes] option_prices1 = [] 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_prices1.append(option.NPV()) prices = pd.DataFrame([strikes, option_prices, option_prices1]).transpose() prices.columns = ['strike', 'Rectangular Int','QuantLib'] prices['curvature'] = (-2 * prices['QuantLib'] + prices['QuantLib'].shift(1) + prices['QuantLib'].shift(-1)) / 1**2 # And plotting... fig = plt.figure(figsize=(10,6)) ax = fig.add_subplot(111) plt.ylabel('Call Price ($)') plt.xlabel('Strikes (K)') ax2 = ax.twinx() lns1 = ax.plot(strikes, option_prices1, label='Option Prices') lns2 = ax2.plot(prices['strike'], prices['curvature'], label='$d^2C/dK^2 (\sim pdf)$', color='orange') ax2.fill_between(prices['strike'], prices['curvature'], color='yellow', alpha=0.2) # added these three lines lns = lns1+lns2 labs = [l.get_label() for l in lns] ax.legend(lns, labs, loc=0) # ax.legend(loc="center right") # ax2.legend(loc="upper right") plt.ylabel('$f_\\tau(K)$') plt.title('QuantLib: Risk-neutral PDF, $f_\mathbb{Q}(K, \\tau)$')
What are the differences? – and why
The differences are mainly a result of the method of complex integration over very small increments.
Python has issues with round off errors at small values approximated by binary floating-point numbers.
https://stackoverflow.com/questions/15930381/python-round-off-error
https://docs.python.org/3/tutorial/floatingpoint.html
mse = np.mean( (option_prices - option_prices1)**2 ) print("QuantLib vs. Our Rect Int \n Mean Squared Error: ", mse) prices.dropna() prices.head(40)
Interpolation of the Risk-Neutral density function
inter = prices.dropna() pdf = sc.interpolate.interp1d(inter.strike, np.exp(r*tau)*inter.curvature, kind = 'linear') pdfc = sc.interpolate.interp1d(inter.strike, np.exp(r*tau)*inter.curvature, kind = 'cubic') strikes = np.arange(61, 179, 1.0) plt.plot(strikes, pdfc(strikes), '-+', label='cubic') plt.plot(strikes, pdf(strikes), label='linear') plt.fill_between(strikes, pdf(strikes), color='yellow', alpha=0.2) plt.xlabel('Strikes (K)') plt.ylabel('$f_\\tau(K)$') plt.title('Risk-neutral PDF: $f_\mathbb{Q}(K, \\tau)$') plt.legend() plt.show()
Cumulative distribution function
cdf = sc.interpolate.interp1d(inter.strike, np.cumsum(pdf(strikes)), kind = 'linear') plt.plot(strikes, cdf(strikes)) plt.xlabel('Strikes (K)') plt.ylabel('$F_\\tau(K)$') plt.title('Risk-neutral CDF: $F_\mathbb{Q}(K, \\tau)$')
Using the Risk-neutral PDF to price ‘complex’ derivatives
Remember the RND function \(f^{\tau}_\mathbb{Q}(S_T)\) is defined at a particular time to expiry \(\tau\)
Call
\(\large C(K, \tau) = e^{-r\tau} \int^\infty_K (S_T – K)f^{\tau}_\mathbb{Q}(S_T)dS_T\)
Put
\(\large P(K, \tau) = e^{-r\tau} \int^K_{-\infty} (K – S_T)f^{\tau}_\mathbb{Q}(S_T)dS_T\)
http://econweb.rutgers.edu/mizrach/pubs/%5B42%5D-2010_Handbook.pdf
def integrand_call(x, K): return (x-K)*pdf(x) def integrand_put(x, K): return (K-x)*pdf(x) calls, puts = [], [] for K in strikes: # integral from K to infinity (looking at CDF, 179 last defined value on range) call_int, err = sc.integrate.quad(integrand_call, K, 178, limit=1000, args=K) # integral from -infinity to K (looking at CDF, 61 lowest defined value on range) put_int, err = sc.integrate.quad(integrand_put, 61, K, limit=1000, args=K) call = np.exp(-r*tau)*call_int calls.append( call ) put = np.exp(-r*tau)*put_int puts.append( put ) rnd_prices = pd.DataFrame([strikes, calls, puts]).transpose() rnd_prices.columns = ['strike', 'Calls','Puts'] rnd_prices.tail()