我的 Heston 模型程式碼有什麼問題

我正在嘗試編寫 heston 模型定價器。但是，一開始它似乎是正確的，但是當插入極端數據時，我以負機率或負價格檢索自己。

有程式碼：

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

如果您輸入：

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

你應該有 p1 = -3.3306690738754696e-16 和 p2 = -9.103828801926284e-15

哪個不正確

或者如果您嘗試：

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

結果給出的價格為 -4.0850747692973335e-06，這不可能。

然而，對於非極端值（例如貨幣期權），我發現與 black 和 scoles 模型相似的結果，使用平面參數。

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

對於赫斯頓模型的所有計算，我都遵循此處可用的公式：

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

不過，我對 Gil-Palaez 反演積分的數值近似表示懷疑。

謝謝您的幫助。

赫斯頓模型中看漲期權價格的數值近似是出了名的不穩定，很容易導致極端參數的不精確答案。存在幾種不同的公式來計算價格，其中一些公式比其他公式更穩定。您使用的公式可以說是最糟糕的公式之一。

我所知道的最精確的算法是由 Leif Anderson 和 Mark Lake 開發的。你可以在這裡找到一個 python 實現：https ://github.com/tcpedersen/anderson-lake-python

在您的情況下，未驗證 Feller 條件：

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

如果此條件未得到驗證，您可以獲得負變異數，如本維基百科文章中所述：

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

引用自：https://quant.stackexchange.com/questions/63306