我正在尝试用Python实现Kou双指数模型来进行期权定价,但是我找不到哪里出错了。我已经递归地定义了
Hh
函数,并通过汇合超几何函数进行逼近,两者都产生相似的结果。有人可以帮我找到错误,并提供 Kou 原始论文中 9.14732 的中间结果吗?请参阅下面的代码
import numpy as np
import scipy as sp
from scipy.stats import norm
from scipy.special import hyp1f1, gamma, comb, factorial
def Hh(n,x):
if n<-1: return 0
elif n==-1:
return np.exp(-x**2/2)
elif n==0:
return np.sqrt(2*np.pi)*norm.cdf(-x)
else:
return (Hh(n-2,x)-x*Hh(n-1,x))/n
def Hh1(n,x):
p1 = 2**(-n/2)*np.sqrt(np.pi)*np.exp(-x**2/2)
p2 = hyp1f1(n/2+1/2,1/2,x**2/2)/(np.sqrt(2)*gamma(1+n/2))
p3 = x*hyp1f1(n/2+1,3/2,x**2/2)/gamma(1/2+n/2)
return p1*(p2-p3)
def P(n,k,eta1,eta2,p):
if n==k:
return p**n
else:
P = 0
for i in range(k,n):
P += comb(n-k-1,i-k)*comb(n,i)*(eta1/(eta1+eta2))**(i-k)*(eta2/(eta1+eta2))**(n-i)*p**i*(1-p)**(n-i)
return P
def Q(n,k,eta1,eta2,p):
if n==k:
return (1-p)**n
else:
Q = 0
for i in range(k,n):
Q += comb(n-k-1,i-k)*comb(n,i)*(eta1/(eta1+eta2))**(n-i)*(eta2/(eta1+eta2))**(i-k)*p**(n-i)*(1-p)**i
return Q
def I(n,c,alpha,beta,delta):
I = 0
if beta>0 and alpha!=0:
for i in range(n+1):
I += (beta/alpha)**(n-i) * Hh1(i,beta*c-delta) + (beta/alpha)**(n+1) * (np.sqrt(2*np.pi)/beta) * np.exp(alpha*delta/beta+alpha**2/(2*beta**2)) * norm.cdf(-beta*c+delta+alpha/beta)
elif beta<0 and alpha<0:
for i in range(n+1):
I += (beta/alpha)**(n-i) * Hh1(i,beta*c-delta) - (beta/alpha)**(n+1) * (np.sqrt(2*np.pi)/beta) * np.exp(alpha*delta/beta+alpha**2/(2*beta**2)) * norm.cdf(beta*c-delta-alpha/beta)
else:
I = 0
return I
def U(mu,sigma,lambd,p,eta1,eta2,a,T, bound=13):
def Pi(n):
return np.exp(-lambd*T)*(lambd*T)**n/factorial(n)
exp1 = np.exp((sigma*eta1)**2*T/2) / (sigma*np.sqrt(2*np.pi*T))
exp2 = np.exp((sigma*eta2)**2*T/2) / (sigma*np.sqrt(2*np.pi*T))
sum1 = 0
sum2 = 0
for n in range(1,bound):
sumP = 0
sumQ = 0
for k in range(1,n+1):
sumP += P(n,k,eta1,eta2,p) * (sigma*np.sqrt(T)*eta1)**k * I(k-1,a-mu*T,-eta1,-1/(sigma*np.sqrt(T)),-sigma*eta1*np.sqrt(T))
sumQ += Q(n,k,eta1,eta2,p) * (sigma*np.sqrt(T)*eta2)**k * I(k-1,a-mu*T,eta2,1/(sigma*np.sqrt(T)),-sigma*eta2*np.sqrt(T))
sum1 += Pi(n)*sumP
sum2 += Pi(n)*sumQ
return exp1*sum1 + exp2*sum2 + Pi(0)*norm.cdf(-(a-mu*T)/(sigma*np.sqrt(T)))
def Kou(r,sigma,lam,p,eta1,eta2,S0,K,T):
zeta = (p*eta1)/(eta1-1)+((1-p)*eta2)/(eta2+1)-1
lam2 = lam * (zeta+1)
eta12 = eta1 - 1
eta22 = eta2 + 1
p2 = (p/(1+zeta))*(eta1/(eta1-1))
omega1 = r+sigma**2/2-lam*zeta
omega2 = r-sigma**2/2-lam*zeta
return S0*U(omega1,sigma,lam2,p2,eta12,eta22,np.log(K/S0),T) - np.exp(-r*T)*K*U(omega2,sigma,lam,p,eta1,eta2,np.log(K/S0),T)
S0=100
sigma=0.16
r=0.05
lam=1
n1=10
n2=5
tau=0.5
K=98
q=0.4
就是I函数,你可以查看论文中的循环和函数。
我实际上正在重复使用您的代码来撰写论文,谢谢。