1
$\begingroup$

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.

$\endgroup$
4
$\begingroup$

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

$\endgroup$
1
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ I used specific parameters that satisfied the Feller conditions in my exemples . $\endgroup$
    – lays
    Apr 9 at 19:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.