¿Qué es lo que falla en el código de mi modelo Heston?

Question

Estoy intentando codificar un modelo Heston de precios.Sin embargo, parece correcto al principio pero al insertar datos extremos me encuentro con probabilidades negativas o precios negativos.

Ahí está el código :

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

Si introduces :

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

Debería tener p1 = -3,3306690738754696e-16 y p2 = -9,103828801926284e-15

Lo cual es incorrecto

o si lo intentas con :

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

el resultado da un precio de -4,0850747692973335e-06, que no puede ser cierto.

Sin embargo, para los valores no extremos (por ejemplo, las opciones at the money) encuentro resultados similares con el modelo de Black y Scholes, utilizando parámetros planos.

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

Estoy siguiendo las fórmulas disponibles aquí para todos los cálculos del modelo heston:

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

Sin embargo, tengo dudas sobre la aproximación numérica de la integral de inversión de Gil-Palaez.

Gracias por su ayuda.

Answer 1

La aproximación numérica del precio de la opción de compra en el modelo de Heston es notoriamente inestable y puede llevar fácilmente a respuestas imprecisas para un parámetro extremo. Existen varias fórmulas diferentes para calcular el precio, siendo algunas más estables que otras. La fórmula que usted utiliza es posiblemente una de las peores.

El algoritmo más preciso que conozco ha sido desarrollado por Leif Anderson y Mark Lake. Puedes encontrar una implementación en python aquí: https://github.com/tcpedersen/anderson-lake-python

Answer 2

La condición de Feller no se verifica en su caso:

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

Si esta condición no se verifica, se puede obtener una varianza negativa como se explica en este artículo de la wikipedia:

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