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

Question

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 = integrate.quad(integrand,0,np.inf,args=(s,k,vol,t,r))
    int1 = 0.5 + int1[0]/np.pi

    int2 = integrate.quad(integrand2,0,np.inf,args=(s,k,vol,t,r)))
    int2 = 0.5 + int2[0]/np.pi

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

print(call_value(100,10,0.2))

Thanks

Answer 1

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 = integrate.quad(integrand,0,np.inf,args=(s,k,vol,t,r))
    int1 = 0.5 + int1[0]/np.pi

    int2 = integrate.quad(integrand2,0,np.inf,args=(s,k,vol,t,r))
    int2 = 0.5 + int2[0]/np.pi
    
    return s*int1 - np.exp(-r*t)*k*int2

print(call_value(100,10,0.2))

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.