寇模型實現Python

嘿，我嘗試在 Python 中實現 Kou 模型。這是我的程式碼：

def Hh(n,x):
   if n<-1: return 0
   elif n==-1:
       return np.exp(-x**2/2)
   elif n==0:
       return math.sqrt(2*np.pi)*scs.norm.cdf(-x)
   else:
       return (Hh(n-2,x)-x*Hh(n-1,x))/n


def P(n,k):
   if k<1 or n<1: return 0
   elif n==k:
       return p**n
   else:
       suma=0
       i=k
       while i<=n-1:
           suma=suma+sc.special.binom(n-k-1, i-k)*sc.special.binom(n, i)*(n_1/(n_1+n_2))**(i-k)*(n_2/(n_1+n_2))**(n-i)*p**i*q**(n-i)
           i+=1
       return suma


def Q(n,k):
if k<1 or n<1: return 0
elif n==k:
   return q**n
else:
   suma=0
   i=k
   while i<=n-1:
       suma=suma+sc.special.binom(n-k-1, i-k)*sc.special.binom(n, i)*(n_1/(n_1+n_2))**(n-i)*(n_2/(n_1+n_2))**(i-k)*p**(n-i)*q**i
       i+=1
   return suma

def Pi(n):
   return (np.exp(-lam*T)*(lam*T)**n)/math.factorial(n)


def I(n,c,a,b,d):
   if b>0 and a!=0:
       suma=0
       i=0
       while i<=n:
           suma=suma+(b/a)**(n-i)*Hh(i,b*c-d)
           i+=1
       return -(np.exp(a*c)/a)*suma+(b/a)**(n+1)*(np.sqrt(2*np.pi)/b)*np.exp((a*d/b)+(a**2/(2*b**2)))*scs.norm.cdf(-b*c+d+a/b)
   elif b<0 and a<0:
       suma=0
       i=0
       while i<=n:
           suma=suma+(b/a)**(n-i)*Hh(i,b*c-d)
           i+=1
       return -(np.exp(a*c)/a)*suma-(b/a)**(n+1)*(np.sqrt(2*np.pi)/b)*np.exp((a*d/b)+(a**2/(2*b**2)))*scs.norm.cdf(b*c-d-a/b)
   else: return 0


def Y(mu,sigma,lam,p, n_1, n_2, a ,T):
   n=1
   suma1=0
   suma2=0
   while n<=10:
       k=1
       suma_1=0
       while k<=n:
           suma_1=suma_1+P(n,k)*(sigma*np.sqrt(T)*n_1)**k*I(k-1,a-mu*T,-n_1, -1/(sigma*np.sqrt(T)), -sigma*n_1*np.sqrt(T))
           k+=1
       suma1=suma1+Pi(n)*suma_1
       n+=1
   n=1
   while n<=10:
       k=1
       suma_2=0
       while k<=n:
           suma_2=suma_2+Q(n,k)*(sigma*np.sqrt(T)*n_2)**k*I(k-1,a-mu*T,n_2, 1/(sigma*np.sqrt(T)), -sigma*n_2*np.sqrt(T))
           k+=1
       suma2=suma2+Pi(n)*suma_2
       n+=1
   return np.exp((sigma*n_1)**2*T/2)/(sigma*np.sqrt(2*np.pi*T))*suma1+np.exp((sigma*n_2)**2*T/2)/(sigma*np.sqrt(2*np.pi*T))*suma2+Pi(0)*scs.norm.cdf(-(a-mu*T)/(sigma*np.sqrt(T)))


def Kou(r,sigma,lam,p,n_1,n_2,S_0,K,T):
   zeta=p*n_1/(n_1-1)+(1-p)*n_2/(n_2+1)-1
   lam2=lam*(zeta+1)
   n_12=n_1-1
   n_22=n_2+1
   p2=p/(1+zeta)*n_1/(n_1-1)
   return S_0*Y(r+1/2*sigma**2-lam*zeta,sigma,lam2,p2,n_12,n_22,math.log(K/S_0),T)-K*np.exp(-r*T)*Y(r-1/2*sigma**2-lam*zeta,sigma,lam,p,n_1,n_2,math.log(K/S_0),T)

並使用以下數據（來自 Kou 2002）

S_0=100
sigma=0.16
r=0.05
lam=1
n_1=10
n_2=5
T=0.5
K=98
p=0.4

不幸的是，我的結果是 6.25，而在 Kou 時應該是 9.14732。有人可以檢查我的程式碼是否正常，或者如果有人有 Kou 模型的程式碼，他是否可以為不同的函式提供多個值，以便我可以檢查哪個函式有錯誤？

這個matlab函式有幫助嗎？

表達式為 $ \Upsilon $ ( Y) 分為三部分。無限和在 10 次迭代後被截斷 ( bound)。我們已經比較了和PFunction的結果。QFunction``IFunction

function Y = Upsilon(x,T,mu,sigma,lambda,etaplus,etaminus,p,q)

   bound = 10;
   
   pi0 = exp(-lambda*T);
   pin = exp(-lambda*T) .*(lambda*T).^(1:bound)./factorial(1:bound);
   
   sump1 = zeros(bound,1);
   sumq1 = zeros(bound,1);
   
   for n=1:bound
       sump2 = zeros(n,1);
       sumq2 = zeros(n,1);
       for k=1:n
           sump2(k) = PFunction(n,k,p,q,etaplus,etaminus) * (sigma*sqrt(T)*etaplus)^k * IFunction(-etaplus,-1/(sigma*sqrt(T)),x-mu*T,-sigma*etaplus*sqrt(T),k-1);
           sumq2(k) = QFunction(n,k,p,q,etaplus,etaminus) * (sigma*sqrt(T)*etaminus)^k * IFunction(etaminus,1/(sigma*sqrt(T)),x-mu*T,-sigma*etaminus*sqrt(T),k-1);
       end
      sump1(n) = pin(n)*sum(sump2);
      sumq1(n) = pin(n)*sum(sumq2);
   end
   
   Y(1) = exp((sigma*etaplus)^2*T/2)/(sigma*sqrt(2*pi*T)) * sum(sump1);
   Y(2) = exp((sigma*etaminus)^2*T/2)/(sigma*sqrt(2*pi*T)) * sum(sumq1);
   Y(3) = pi0*normcdf(-(x-mu*T)/(sigma*sqrt(T)));
   Y = sum(Y);
   
end

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