如何在Python中重现Kou（2002）双指数模型？

Question

我正在尝试用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

Answer 1

就是I函数，你可以查看论文中的循环和函数。

我实际上正在重复使用您的代码来撰写论文，谢谢。