What is the formula for the vanilla option (Call/Put) price in the Heston model?
I only found the bi-variate system of stochastic differential equations of Heston model but no expression for the option prices.
In the Heston Model we have \begin{align} C(t\,,{{S}_{t}},{{v}_{t}},K,T)={{S}_{t}}{{P}_{1}}-K\,{{e}^{-r\tau }}{{P}_{2}} \end{align} where, for $j=1,2$
\begin{align} & {{P}_{j}}({{x}_{t}}\,,\,{{v}_{t}}\,;\,\,{{x}_{T}},\ln K)=\frac{1}{2}+\frac{1}{\pi }\int\limits_{0}^{\infty }{\operatorname{Re}\left( \frac{{{e}^{-i\phi \ln K}}{{f}_{j}}(\phi ;t,x,v)}{i\phi } \right)}\,d\phi \\ & {{f}_{j}}(\phi \,;{{v}_{t}},{{x}_{t}})=\exp [{{C}_{j}}(\tau ,\phi )+{{D}_{j}}(\tau ,\phi ){{v}_{t}}+i\phi {{x}_{t}}] \\ \end{align}
and
\begin{align} & {{C}_{j}}(\tau ,\phi )=(r-q)i\phi \,\tau +\frac{a}{{{\sigma }^{2}}}{{\left( ({{b}_{j}}-\rho \sigma i\phi +{{d}_{j}})\,\tau -2\ln \frac{1-{{g}_{j}}{{e}^{{{d}_{j}}\tau }}}{1-{{g}_{j}}} \right)}_{_{_{_{{}}}}}} \\ & {{D}_{j}}(\tau ,\phi )=\frac{{{b}_{j}}-\rho \sigma i\phi +{{d}_{j}}}{{{\sigma }^{2}}}\left( \frac{1-{{e}^{{{d}_{j}}\tau }}}{1-{{g}_{j}}{{e}^{{{d}_{j}}\tau }}} \right) \\ \end{align}
where \begin{align} & {{g}_{j}}=\frac{{{b}_{j}}-\rho \sigma i\phi +{{d}_{j}}}{{{b}_{j}}-\rho \sigma i\phi -{{d}_{j}}} \\ & {{d}_{j}}=\sqrt{{{({{b}_{j}}-\rho \sigma i\phi )}^{2}}-{{\sigma }^{2}}(2i{{u}_{j}}\phi -{{\phi }^{2}})} \\ & {{u}_{1}}=\frac{1}{2}\,,\,{{u}_{2}}=-\frac{1}{2}\,,\,a=\kappa \theta \,,\,{{b}_{1}}=\kappa +\lambda -\rho \sigma \,,\,{{b}_{2}}=\kappa +\lambda \,,\ {{i}^{2}}=-1 \\ \end{align} Other representations:
The Heston model is given by the following equations in a risk-neutral measure: $$ dS_t = (r - q)S_t dt + \sqrt{V_t} S_t dW_t, \\ dV_t = \kappa(\theta - V_t)dt + \sigma_V \sqrt V_t dZ_t,\\ dS_t dW_t = \rho dt. $$
The call formula from (Lipton, 2002) is as follows:
$$ C(K, T) = S_0 e^{-qT} - \frac{Ke^{-rT}}{\pi} \int\limits_0^\infty \mathfrak{Re}\left[ e^{(iu+\frac 1 2)k} \phi_T(u - \tfrac i 2)\right] \frac{du}{u^2 + \tfrac 1 4}, \\ k = \ln(S_0 / K) + (r-q)T, $$ where $\phi_T$ is the characteristic function. For Heston model: $$ \phi_T(u) = \exp[A(u) + B(u)V_0],\\ A(u) = \frac{\kappa \theta}{\sigma_V^2} \left[(\beta-d)T - 2\ln \left(\frac{ge^{-dT} - 1}{g-1} \right) \right],\\ B(u) = \frac{\beta-d}{\sigma_V^2}\left( \frac{1-e^{-dT}}{1-ge^{-dT}} \right), \\ d = \sqrt{\beta^2 - 4\hat \alpha \gamma},\\ g = \frac{\beta-d}{\beta+d},\\ \hat\alpha = -\frac 1 2 u(u+i), \quad \beta = \kappa - iu\sigma_v \rho, \quad \gamma = \frac 1 2 \sigma_V^2 $$ The code in Python is below:
import numpy as np
import scipy.integrate
r = 0.02
q = 0.01
s = 100.0
v = 0.04
kappa = 1.0
theta = 0.04
sigma_v = 0.2
rho = -0.7
k = 100.0
tau = 1.0
def phi(u, tau):
alpha_hat = -0.5 * u * (u + 1j)
beta = kappa - 1j * u * sigma_v * rho
gamma = 0.5 * sigma_v ** 2
d = np.sqrt(beta**2 - 4 * alpha_hat * gamma)
g = (beta - d) / (beta + d)
h = np.exp(-d*tau)
A_ = (beta-d)*tau - 2*np.log((g*h-1) / (g-1))
A = kappa * theta / (sigma_v**2) * A_
B = (beta - d) / (sigma_v**2) * (1 - h) / (1 - g*h)
return np.exp(A + B * v)
def integral(k, tau):
integrand = (lambda u:
np.real(np.exp((1j*u + 0.5)*k)*phi(u - 0.5j, tau))/(u**2 + 0.25))
i, err = scipy.integrate.quad_vec(integrand, 0, np.inf)
return i
def call(k, tau):
a = np.log(s/k) + (r-q)*tau
i = integral(a, tau)
return s * np.exp(-q*tau) - k * np.exp(-r*tau)/np.pi*i
call(k, tau)