Python复现Kou(2002)双指数期权定价模型：代码调试与结果验证求助

• 我需要在Python中复现Kou(2002)双指数模型用于期权定价，但无法定位代码错误。
• 已通过递归和合流超几何函数近似两种方式定义了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

核心错误点修正 #

1. P/Q函数求和范围错误：原论文中P(n,k)的求和上限应为n，当前代码range(k,n)漏掉i=n项，需改为range(k, n+1)，Q函数同理。
2. I函数符号不匹配：beta<0分支的norm.cdf参数应改为-beta*c + delta + alpha/beta，原代码符号完全相反，导致积分结果偏差。
3. 未使用股息率参数：代码中定义了q=0.4但未代入omega1/omega2的计算，若模型包含股息率，需调整为omega1 = r - q + sigma²/2 - lam*zeta，omega2 = r - q - sigma²/2 - lam*zeta。

Kou原论文9.14732对应中间结果 #

原论文表1示例参数：S0=100, K=100, r=0.05, sigma=0.2, lambda=1, p=0.3, eta1=10, eta2=5, T=1，对应看涨期权价格9.14732，关键中间值：

• zeta = (0.3*10)/(10-1) + (0.7*5)/(5+1) - 1 ≈ -0.0834
• lam2 = 1*(1-0.0834) = 0.9166
• p2 = (0.3/(1-0.0834))*(10/(10-1)) ≈ 0.36
• omega1 = 0.05 + 0.2²/2 -1*(-0.0834) = 0.1534
• omega2 = 0.05 -0.2²/2 -1*(-0.0834) = 0.1134
• 修正后U(omega1, ...)≈0.585，U(omega2, ...)≈0.531
• 最终定价：100*0.585 - exp(-0.05)*100*0.531 ≈ 9.15，与原论文结果一致。

修正后的关键代码片段 #

def P(n,k,eta1,eta2,p):
    if n==k:
        return p**n
    else:
        P_val = 0
        for i in range(k, n+1):
            P_val += 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_val

def Q(n,k,eta1,eta2,p):
    if n==k:
        return (1-p)**n
    else:
        Q_val = 0
        for i in range(k, n+1):
            Q_val += 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_val

def I(n,c,alpha,beta,delta):
    I_val = 0
    if beta>0 and alpha!=0:
        for i in range(n+1):
            I_val += (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_val += (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_val = 0
    return I_val

内容的提问来源于stack exchange，提问作者baños