# What is wrong in my Heston model's code

I am trying to code a heston model pricer.However,it seems correct at the beginning but when inserting extreme data I retrieve myself with negative probabilities or negative prices.

There is the code :

from scipy import *
from math import *
from scipy.stats import norm
from scipy.optimize import fmin_bfgs
from scipy.optimize import brentq
from scipy.optimize import minimize, rosen, rosen_der,least_squares

#public
def call_price(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K,J):
p1 = __p1(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K,J)
p2 = __p2(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K,J)
print(p1,p2)
C = s0 * p1 - K * np.exp(-r*T)*p2
return C

def integrandd(phi):
A = np.exp(-1j *  phi * np.log(K[:]))
C = __fm(phi, kappa, theta, sigma, rho, v0, r, T[:], s0, status)
B = (1j * phi)
return (A * C / B).real

#private
def __p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 , K, status,J):
if J =='No':
integrand = lambda phi: (np.exp(-1j * phi * np.log(K)) * __f(phi, kappa, theta, sigma, rho, v0, r, T, s0, status) /(1j * phi)).real
else :
def integrand(phi):
if len(phi) > 200:
return 0
A = np.array([np.exp(-1j * phi * np.log(i)) for i in K])
C = np.array([__f(phi, kappa, theta, sigma, rho, v0, r,i, s0, status) for i in T])
B = (1j * phi)
return (A * C / B).real

p = 0.50 + 1/pi * (quad(integrand, 0.0, 100))

return p

def __p1(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K,J):
return __p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K, 1,J)
def __p2(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K,J):
return __p(kappa, theta, sigma, rho, v0 ,r ,T ,s0 ,K, 2,J)
def __fm(phi, kappa, theta, sigma, rho, v0, r, T, s0, status):

if status == 1:
u = 0.5
b = kappa - rho * sigma
else:
u = -0.5
b = kappa

a = kappa * theta
x = np.log(s0)
d = np.sqrt((rho * sigma * phi * 1j - b)**2 - sigma**2 * (2 * u * phi * 1j - phi**2))
g = (b - rho * sigma * phi * 1j + d) / (b - rho * sigma * phi * 1j - d)
C = np.array([r * (phi * 1j * i) + (a / sigma**2)*((b - rho * sigma * phi * 1j + d) * i - 2 * np.log((1 - g * np.exp(d * i))/(1 - g))) for i in T])
D = np.array([(b - rho * sigma * phi * 1j + d) / sigma**2 * ((1 - np.exp(d * i)) / (1 - g * np.exp(d * i))) for i in T])
Final = np.exp(C + D * v0 + 1j * phi * x)
return Final

def Optimisor(x,args):
price_,S, K, T, r = args
A = _price_(S, K, T, r, x)
MNM = A - price_
return MNM
#For Matrices

def __f(phi, kappa, theta, sigma, rho, v0, r, T, s0, status):

if status == 1:
u = 0.5
b = kappa - rho * sigma
else:
u = -0.5
b = kappa

a = kappa * theta
x = np.log(s0)
d = np.sqrt((rho * sigma * phi * 1j - b)**2 - sigma**2 * (2 * u * phi * 1j - phi**2))
g = (b - rho * sigma * phi * 1j + d) / (b - rho * sigma * phi * 1j - d)
C = r * ( phi * 1j * T) + (a / sigma**2)*((b - rho * sigma * phi * 1j + d) * T - 2 * np.log((1 - g * np.exp(d * T))/(1 - g)))
D = (b - rho * sigma * phi * 1j + d) / sigma**2 * ((1 - np.exp(d * T)) / (1 - g * np.exp(d * T)))
return np.exp(C + D * v0 + 1j * phi * x)

def implied_volatility(price_ ,S, K, T, r):
args = [price_,S, K, T, r]
try:
res = brentq(Optimisor,0.0001,10,args,maxiter=10000)
except :
res = np.nan
return res

def _price_(S, K, T, r, v):
d1 = (np.log(S/K) + (r + 0.5 * v ** 2) * T)/(v * T**(1/2))
d2 = d1 - v * T **(1/2)
Option_price = S * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
return Option_price


If you input :

S0 = 13780.79296875
Kappa,theta,sigma,rho,v0,r,T = 0.01,0.50,0.01,-1,0.10,0,.1
call_price(Kappa,theta,sigma,rho,v0,r,T,S0,S0*4,'No')


You should have p1 = -3.3306690738754696e-16 and p2 = -9.103828801926284e-15

Which is incorrect

or if you try with :

Kappa,theta,sigma,rho,v0,r,T = 0.01,0.01,0.01,-1,0.10,0,.01


the result gives a price of -4.0850747692973335e-06,which cannot be true.

However for non extreme values(for instance at the money options) I find similar results with the black and scholes model,using flat parameters.

call_price(.01,0.01,0.01,0,0.010,0,.1,S0,S0,'No')# Heston
173.83935148144155

_price_(S0,S0,.1,0,0.1) # Black-Scholes
173.8465909552633


I am following the formulas available here for all the computations of the heston model:

https://frouah.com/finance%20notes/The%20Heston%20model%20short%20version.pdf

I have doubts on the numerical approximation of the Gil-Palaez Inversion integral though.

The numerical approximation of the call option price in the Heston model is notoriously unstable and can easily lead to imprecise answers for extreme parameter. Several different formulas exist for computing the price with some being more stable than others. The formula you are using is arguably one of the worst ones.

The most precise algorithm I know of is developed by Leif Anderson and Mark Lake. You can find a python implementation here: https://github.com/tcpedersen/anderson-lake-python

The Feller condition is not verified in your case:

$$2\kappa\theta>\xi ^ {2}$$

If this condition is not verified, you can get negative variance as explained in this wikipedia article:

https://en.wikipedia.org/wiki/Heston_model

• I used specific parameters that satisfied the Feller conditions in my exemples .
– lays
Apr 9, 2021 at 19:35