[Python] Complex integration of option pricing

I'm finishing my economics degree and need to present monograph. The theme I've chosen is option pricing in jum-diffusion models. To find the solution I'm using the Mellin transform and found the integral,



where c is a complex variable, and



I tried to use python to calculate it, but my results doesn't make sense. The option price is higher than the underlying asset. And comparing with the approach of finite series of Merton (1976), for my inputs the output should be around 20, but it's 177!

The inputs I'm using are:

# Models parameters
# BSM parameters
S = 100 # Preço do ativo no momento de compra da put
K = 115 # Strike
T = 1 # Tempo até o vencimento
r = 0.125 # Taxa livre de risco
sigma = 0.15 # Volatilidade

# Merton JD parameters
lambda_ = 0.5 # Intensidade dos saltos
mu_j = 0.2 # Média dos saltos
sigma_j = 0.45 # Desvio padrão dos saltos

# Kou JD parameters
p = 0.65
q = (1-p)
eta1 = 40
eta2 = 50


And the function for the integration is (ignore the Greek letters, it's the easy part):

def JD_put_value(S, K, T, t, r, sigma, jump_dist):
'''
Inputs:
S : float - Asset price
K : float - Strike price
T : float - Time
r : float - Free-risk rate
sigma : float - Asset volatility
jump_dist : str - Jumps-distribution

Outputs:
put_price : float
delta : float
gamma : float
theta : float
vega : float
rho : float'''

E_y_LN = np.exp(mu_j + 0.5 * sigma_j**2)
def E_yw_LN(omega):
return np.exp(-omega * mu_j + 0.5 * (omega**2) * (sigma_j**2)) # Esperança matemática de y^{-omega}

E_y_DE = p/eta1 - q/eta2
def E_yw_DE(omega): return p*((eta1**omega)/(eta1 - omega)) + omega*((eta2**omega)/(eta2-omega))

def Omega(omega):
if jump_dist == 'log-normal':
return (1/2)*sigma**2*(omega+1)*omega - (r - lambda_)*(omega+1) + lambda_*(omega*E_y_LN + E_yw_LN(omega))
elif jump_dist == 'exponencial-dupla':
return (1/2)*sigma**2*(omega+1)*omega - (r - lambda_)*(omega+1) + lambda_*(omega*E_y_DE + E_yw_DE(omega))

def putValueIntegral(omega):
c_omega = 1 + omega*1j
return K**(c_omega+1)/(c_omega*(c_omega+1))*np.exp(Omega(c_omega)*(T-t))*S**(-c_omega)
c = 1
result, error = quad(putValueIntegral, -np.inf, np.inf, complex_func=True)
return (result/(2*np.pi), error)

JD_put_value(S, K, T, 0, r, sigma, 'log-normal')


The method used in Merton (1976) is coded as:

def merton_put_price(S, K, T, r, sigma, lambda_j, mu_j, sigma_j, max_jumps=20):
tau = T
put_price = 0.0

for n in range(max_jumps + 1):
# Ajuste da volatilidade e taxa de juros
sigma_n = sqrt(sigma**2 + n * sigma_j**2 / tau)
kappa = exp(mu_j + 0.5 * sigma_j**2) - 1
r_n = r - lambda_j * kappa + n * mu_j / tau

# Preço da put de Black-Scholes para n saltos
bs_put = BSMPutValue(S, K, T, r_n, sigma_n)

# Probabilidade de n saltos (Poisson)
prob = exp(-lambda_j * tau) * (lambda_j * tau)**n / factorial(n)

# Soma ponderada
put_price += prob * bs_put[0]

return put_price

# Calculando o preço da put
put_price = merton_put_price(S, K, T, r, sigma, lambda_, mu_j, sigma_j)
print(f"Preço da Put no Modelo de Merton: {put_price:.4f}")

Continue reading...