2
$\begingroup$

I have implemented the Kou option model for pricing vanilla option. I have checked that my program can replicate the price of the option in the original paper of 2002. However, when I use it to price other options, I observe deviation from the prices given by Forsyth (2004) "Robust Numerical Methods for Contingent Claims under Jump Diffusion Processes" or Toivanen "Numerical Valuation of European and American Options under Kou’s Jump-Diffusion Model".

For instance, they price a call at the money at 3.97 whereas with my code, I found 3.69. I attached my code below for your reference.

Just note that I have re-used multiple sources but the principal one was a previous discussion on this forum!

Thanks in advance.

I wonder whether I forgot something in my code.

import numpy as np, scipy as sp, math

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)*sp.stats.norm.cdf(-x)
    else:
        return (Hh(n-2,x)-x*Hh(n-1,x))/n

def P(n, k, p, q, eta_1, eta_2):
    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 + sp.special.binom(n-k-1, i-k) * sp.special.binom(n, i) * (eta_1/(eta_1 + eta_2))**(i-k)*(eta_2/(eta_1 + eta_2))**(n-i)*p**i*q**(n-i)
            i+=1
        return suma

def Q(n, k, p, q, eta_1, eta_2):
    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 + sp.special.binom(n-k-1, i-k) * sp.special.binom(n, i)*(eta_1/(eta_1+eta_2))**(n-i)*(eta_2/(eta_1+eta_2))**(i-k)*p**(n-i)*q**i
            i +=1
        return suma

def I(n,c,a,b,d):
    if n < 0:
        print(n)
    if b>0 and a!=0:
        suma=0
        i=0
        while i<=n:
            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 + 0.5 * a**2 / b**2) * sp.stats.norm.cdf(-b*c+d+a/b)
    elif b<0 and a<0:
        suma=0
        i=0
        while i<=n:
            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 + 0.5 * a**2 / b**2) * sp.stats.norm.cdf(b*c-d-a/b)
    else: 
        return 0

def PSI(mu, Sigma, Lambda, pp, qq, eta_1, eta_2, arg, Time):
    bound=15
    Aa = 0.
    Bb = 0.
    Cc = 0.
    
    pin=np.zeros(bound+1)

    sump1=0.
    sumq1=0.
    
    for i in range(1,bound+1):
        pin[i]=np.exp(-Lambda*Time)*(Lambda*T)**i/math.factorial(i)

    for n in range(1, bound):
        sump2=0.
        sumq2=0.
        
        for k in range(1, n+1):
            sump2 += P(n,k,pp,qq, eta_1, eta_2) * (Sigma * np.sqrt(Time) * eta_1)**k * I(k-1, arg-mu*Time, -eta_1, -1/(Sigma*np.sqrt(Time)), -Sigma * eta_1 * np.sqrt(Time))
            sumq2 += Q(n,k,pp,qq, eta_1, eta_2) * (Sigma * np.sqrt(Time) * eta_2)**k * I(k-1, arg-mu*Time, eta_2,  1/(Sigma*np.sqrt(Time)), -Sigma * eta_2 * np.sqrt(Time))
        
        sump1+=pin[n] * sump2
        sumq1+=pin[n] * sumq2


    Aa = np.exp((Sigma*eta_1)**2*Time/2)/(Sigma*np.sqrt(2*np.pi*Time)) * sump1
    Bb = np.exp((Sigma*eta_2)**2*Time/2)/(Sigma*np.sqrt(2*np.pi*Time)) * sumq1
    Cc = np.exp(-Lambda*Time) * sp.stats.norm.cdf(-(arg-mu*Time)/(Sigma*np.sqrt(Time)))
    
    return Aa + Bb + Cc

def Kou(r,sigma,Lambda,pp,eta_1,eta_2,S_0,K,Time):
    zeta = p * eta_1/(eta_1 - 1.) + (1. - p) * eta_2/(eta_2 + 1.) - 1.
    lam2 = Lambda * (zeta + 1.) 
    eta_12 = eta_1 - 1.
    eta_22 = eta_2 + 1.
    p2   = pp/(1. + zeta) * eta_1/(eta_1 - 1.)   

    return S_0 * PSI(r + 1/2 * sigma**2 - Lambda * zeta, sigma, lam2, p2, 1. - p2, eta_12, eta_22, math.log(K/S_0), Time) - K * np.exp(-r*T) * PSI(r - 1/2 * sigma**2 - Lambda * zeta, sigma, Lambda, pp, 1. - pp, eta_1, eta_2, math.log(K/S_0), Time)
    

S0 = 110.
r = 0.05
K = 100.
Sigma = 0.15
T = 0.25
Lambda_ = .1

Alpha1 = 3.0465
Alpha2 = 3.0775
                                             
p = 0.3445

print(Kou(r,Sigma,Lambda_,p,Alpha1,Alpha2,S0,K,T))

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
q =1.-p

Alpha2

print(Kou(r,sigma,lam,p,n_1,n_2,S_0,K,T))
$\endgroup$
4
  • 1
    $\begingroup$ This is a link to the earlier disscussion (2020) quant.stackexchange.com/questions/57498/… $\endgroup$
    – nbbo2
    Nov 23, 2023 at 11:57
  • $\begingroup$ Thanks, this is one of the discussion I used to write the code. After posting my code, I realized it is actually working. $\endgroup$
    – pierrot
    Nov 23, 2023 at 15:17
  • $\begingroup$ Do you want to close the question ? $\endgroup$
    – nbbo2
    Nov 23, 2023 at 15:20
  • $\begingroup$ Yes please, I don't need any help on this. $\endgroup$
    – pierrot
    Nov 24, 2023 at 1:37

0

Browse other questions tagged or ask your own question.