Kou model implementation Python

Hey I try to implement Kou model in Python. This is my code:

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)


and use following data (from 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


Unfortunately, my result is 6.25, when in Kou should be 9.14732. Can someone check if my code is OK, or if someone has code of Kou model, would he be able to give several values ​​for different functions so I can check in which function I have mistake?

• Kou (2002, Footnote 9) says that’s enough to truncate the infinite sum at $n=10$ or $n=15$. I can confirm this from experimental testing. So you can reduce the number of iterations :) Commented Aug 20, 2020 at 14:53
• I change from $100$ to $10$ but I get result $0$. Something is wrong :( Commented Aug 20, 2020 at 17:14
• @KeSchn do you use Kou model? Can you tell me what values for different functions (for example Y) you get from your code? It would help me to find mistake in my code, because now I have no Idea what is wrong. Commented Aug 20, 2020 at 21:25
• Sorry, yes of course. I get $\Upsilon(0)=0.5932$. Bloody hell, must've been yet another typo in my comment. I'm so sorry. The code below yields $\Upsilon(0)=0.5932$. So, you're right. Well done! Commented Aug 21, 2020 at 19:34
• Please have a look at this question Commented Aug 22, 2020 at 11:12

Does this matlab function help?

The expression for $$\Upsilon$$ (Y) is split up in three parts. The infinite sum is truncated after 10 iterations (bound). We already compared our results for PFunction, QFunction and 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

• Can I ask you how do you calibrate Kou model? It's 5 parameters to estimate and it takes a long time. Commented Aug 23, 2020 at 18:48
• Simply ask Python to find the parameters which minimise the sum of squared errors. This shouldn’t take too long if you follow best practices and write the code efficiently. You should use reasonable initial guesses (look at Kou’s papers). Commented Aug 23, 2020 at 18:54
• Ok thanks, I wil try and probably improve my code becuase now it takes a looong time :D Commented Aug 23, 2020 at 21:36
• is it normal that from calibration I get $\lambda=50$ and $p=1$ and $\eta_1=50$, $\eta_2=155$? Are these values ​​not too high? You write that I should use reasonable initial guesses but in which paper they are given? Commented Sep 12, 2020 at 19:34
• @Mr.Price Sorry, fully forgot to come back to you. The numbers are completely off. Something is going wrong. Have a look at how to implement calibration''. Your implementation my be wrong? You should simply minimise the sum of squared (relative) errors between your model and the market. The market quotes should be free of arbitrage of course. Commented Sep 15, 2020 at 11:08