# Heston model with underlying BS dynamics always gives 1/2 of the right value, what am I doing wrong?

Just as an exercise I'm trying to follow this paper: https://arxiv.org/ftp/arxiv/papers/1502/1502.02963.pdf

In the section 2.2 it calculates the value of a Call using the characteristic function of the underlying under Black-Scholes dynamics.

I'm using Python, and strangely my result seems to be exactly 1/2 of the value of the call option, because both the Pi1 and Pi2 I get are also 1/2 of what is supposed to be.

For instance I'm calculating the value of a call option with parameters:

• S = 100
• K = 10
• t = 1
• vol = 0.2
• r = 0

which should give me a Pi1 (delta) equal to 1 (since it's deep ITM) but instead I get a Pi1 = 0.49999999 and a call value of 44.99999

Could you help to to point out where might I have an error? I really hope is not a silly mistake, I've been trying to figure it out for the last 2-3 hours but I cannot see anything wrong. Here is the code I'm using:

import numpy as np
from scipy import stats, special, integrate

def char_func(x,s,vol,t=1,r=0):
mean = np.log(s) + (r - 0.5*vol*vol) * t
var = vol*vol*t
w = np.exp(1j*x*mean - x*x*vol*vol*0.5)
return w.real

def call_value(s,k,vol,t=1,r=0):
def integrand(x,s,k,vol,t=1,r=0):
I = np.exp(-1j*x*np.log(k)) * char_func(x-1j,s,vol,t,r)/(1j*x*char_func(-1j,s,vol,t,r))
return I.real

def integrand2(x,s,k,vol,t=1,r=0):
I = np.exp(-1j*x*np.log(k)) * char_func(x,s,vol,t,r)/(1j*x)
return I.real

int1 = 0.5 + int1[0]/np.pi

int2 = 0.5 + int2[0]/np.pi

return s*int1 - np.exp(-r*t)*k*int2

print(call_value(100,10,0.2))


Thanks

You're implementing the formula $$C=S_0e^{-qT}\Pi_1-Ke^{-rT}\Pi_2$$ where \begin{align*} \Pi_1&=\frac{1}{2}+\frac{1}{\pi}\int_0^\infty \text{Re}\left(\frac{e^{-iu\ln(K)}\varphi_{\ln(S_T)}(u-i)}{iu\varphi_{\ln(S_T)}(-i)}\right)\text{d}u,\\ \Pi_2&=\frac{1}{2}+\frac{1}{\pi}\int_0^\infty \text{Re}\left(\frac{e^{-iu\ln(K)}\varphi_{\ln(S_T)}(u)}{iu}\right)\text{d}u. \end{align*}

A few words about your code. You should not output the real part of the characteristic function, but the complex value. That is, you need return w and not return w.real. Also, you define var but you don't use it.

I get 89.99999999990168 as option value if I run the below code. This makes sense, the option is so deep ITM, its value is essentially $$S_0e^{-qT}-Ke^{-rT}\approx90$$

import numpy as np
from scipy import stats, special, integrate

def char_func(x,s,vol,t=1,r=0):
mean = np.log(s) + (r - 0.5*vol*vol) * t
var = vol*vol*t
w = np.exp(1j*x*mean - x*x*var*0.5)
return w

def call_value(s,k,vol,t=1,r=0):
def integrand(x,s,k,vol,t=1,r=0):
I = np.exp(-1j*x*np.log(k)) * char_func(x-1j,s,vol,t,r)/(1j*x*char_func(-1j,s,vol,t,r))
return I.real

def integrand2(x,s,k,vol,t=1,r=0):
I = np.exp(-1j*x*np.log(k)) * char_func(x,s,vol,t,r)/(1j*x)
return I.real

int1 = 0.5 + int1[0]/np.pi


Note: The option pricing formua you're using is very simple. It was one of the original Fourier methods discovered. Try instead something like \begin{align} c(S_0;K,T) &= c^\mathrm{BS}(S_0;K,T) +\frac{\sqrt{K}e^{-rT}}{\pi}\int_0^\infty \frac{\text{Re}\left( e^{-iku} \left(\varphi_T^\mathrm{BS}(u-0.5i)-\varphi_T(u-0.5i) \right)\right)}{{u^2+0.25}}\mathrm{d}u, \end{align} where $$c^\mathrm{BS}$$ is the Black Scholes option price and $$\varphi_T^\text{BS}$$ the corresponding log-normal characteristic function. That formula uses a control variate and integrates along a different contour.
• Yes! The main problem indeed is that I was returning w.real instead of just w. I'm not very versed in maths since it's not my main field, so I appreciate the help here. Will also have a look at your suggestion to use that formula. Thanks a lot Kevin. Mar 1, 2021 at 12:38