Heston Option Pricing: Numerical Instability and the Feller Condition

# What is wrong in my Heston model's code


# 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.integrate import quad
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)[0])
    
    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.

Thank you for your help.

## Answer by tcpedersen (score 4, accepted)

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

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

## Answer by Valometrics.com (score 1)

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

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