Heston Pricing: Distinguishing the Two Characteristic Functions

# Wrong Heston parametrization yields the right result


# Wrong Heston parametrization yields the right result












This is a continuation of the question I asked here, but because there is too much information missing I thought I would rewrite it, and include some code examples. First, we know that the heston model has that $P_1$ and $P_2$ follow the following PDE:

$\begin{equation} -\frac{\partial P_j}{\partial\tau} + \rho \sigma \nu \frac{\partial^2 P_j}{\partial x \partial \nu} +\frac{1}{2}\nu \frac{\partial^2 P_j}{\partial x^2} + \frac{1}{2}\nu \sigma^2 \frac{\partial^2 P_j}{\partial \nu^2}+\Big(r + u_j\nu\Big)\frac{\partial P_j}{\partial x} + \Big(\alpha - b_j\nu \Big)\frac{\partial P_j}{\partial \nu} = 0 \end{equation}$

for j = $1,2$ where $u_1 = \tfrac12,\quad u_2 = -\tfrac12,\quad a = \kappa \theta,\quad b_1 = \kappa + \eta^\nu - \rho\,\sigma,\quad b_2 = \kappa + \eta^\nu. $

This results in the following two ODE's: $\begin{equation} \begin{aligned} &- \frac{\partial C_j}{\partial \tau} + \alpha D_j +r i \phi =0 \\ & -\frac{\partial D_j}{\partial \tau} +\rho \sigma i \phi D_j - \frac{1}{2} \phi^2 + \frac{1}{2} \sigma^2 D_j^2 + u_j i \phi - b_j D_j =0\\ \end{aligned} \end{equation}$

If we then define $A_j = \frac{1}{2}\sigma^2$ and $B_j=b_j - \rho \sigma i \phi $ with $C = -\frac{1}{2}\phi^2 +i u_j\phi $. Then we get $ \begin{equation} \frac{\partial D_j}{\partial \tau} = A_j D_j^2 - B_jD_j + C_j \end{equation}$

One can show that a solution to this Ricatti equation (and the other ODE) becomes: $ \begin{equation} \begin{aligned} & d_j = \sqrt{(-B_j)^2 -4A_jC} = \sqrt{( \rho \sigma i \phi - b_j)^2 - \sigma^2(2i u_j\phi-\phi^2 )}\\ &\lambda_{1,2} = \frac{B_j \pm \sqrt{(-B_j)^2-4A_jC}}{2A_j} = \frac{b_j - \rho \sigma i \phi \pm {d_j}}{\sigma^2}\\ &g_j = \frac{\lambda_1}{\lambda_2} = \frac{b_j -\rho \sigma i \phi + d_j}{b_j-\rho \sigma i \phi-d_j} \\ & D_j = \frac{b_j -\rho \sigma i \phi + d_j}{\sigma^2}\Big(\frac{1-e^{d_j \tau}}{1-g_je^{d_j \tau}} \Big)\\ & C_j(\tau,\phi)= ri\phi \tau + \frac{\alpha}{\sigma^2}\Big[(b_j -\rho \sigma i \phi +d_j)\tau - 2\ln\Big(\ \frac{1- g_j e^{d_j \tau}}{1-g_j} \Big)\Big] \end{aligned} \end{equation} $

Again, following this pdf. However, this parametrization has branch-cut issues. One can show that an equivalent solution is given by:

$\begin{equation} \begin{aligned} % & f_j (\phi, x, \nu) = \exp (C_j(\tau,\phi) +D_j(\tau,\phi)\nu + i\phi x ) \\ & C_j(\tau,\phi)= ri\phi\tau + \frac{\alpha}{\sigma^2}\Big((b_j- \rho\sigma i\phi -d_j)\tau - 2\ln\Big(\ \frac{1- g_j e^{-d_j \tau}}{1-g_j} \Big) \Big)\\ & D_j(\tau,\phi) = \frac{b_j - \rho\sigma i \phi - d_j}{\sigma^2}\Big(\frac{1 - e^{-d_j\tau}}{1-g_je^{-d_j\tau}} \Big) \\ &g_j = \frac{b_j - \rho \sigma i \phi - d_j}{b_j - \rho \sigma i \phi + d_j} ,\quad d_j =\sqrt{( \rho \sigma i \phi - b_j)^2 - \sigma^2(2i u_j\phi-\phi^2 )}\\ \end{aligned} \end{equation}$

So that finally we get: $ \begin{equation} C(t,S_t) = S_t \Big( \frac{1}{2} + \frac{1}{\pi}\int^{\infty}_{0} \mathcal{R} \Big\{ \frac{ \Phi_X{(u-i)}}{iu \Phi_X{(-i)}} e^{-iu \ln K} \Big\} du \Big) - Ke^{-r(T-t)} \Big( \frac{1}{2} + \frac{1}{\pi}\int^{\infty}_{0} \mathcal{R} \Big\{ \frac{ \Phi_X{(u)}}{iu}e^{-iu \ln K} \Big\} du \Big) \end{equation}$

With $\Phi_X= \exp (C_j(\tau,\phi) +D_j(\tau,\phi)\nu + i\phi x ) = E^Q[e^{i\phi X_T}]$

However, when I try to implement this, I can only get it to agree with other implementations 1 2 3 by forcing $P_1$ to follow $j=2$. My understanding was that $P_1$ should have $u_1 = 0.5$ and $b_1 = \kappa + \eta^\nu - \rho\,\sigma,\quad$. See the code below, where if you set version=1 where j=1 (i.e $P_1$) it no longer agrees with e.g the non-analytical Monte Carlo simulation-derived price. What am I misunderstanding? Are you supposed to pick a $j$ for both $P_j$? but in the case I pick $j=1$, the solution also does not agree with the other implementations, only when I pick $j=2$ :/.

```
import numpy as np

from scipy.integrate import quad
import cmath

def heston_d(phi, sigma, rho, u, b):
    i = 1j
    term1 = (i * rho * sigma * phi - b)**2
    term2 = (sigma**2) * (2 * i * u * phi - phi**2)
    return cmath.sqrt(term1 - term2)

def heston_g(b, rho, sigma, phi, d):
    i = 1j
    numerator   = b - i * rho * sigma * phi - d
    denominator = b - i * rho * sigma * phi + d
    return numerator / denominator

def heston_C(phi, tau, r, kappa, theta, sigma, rho, u, b):

    i = 1j
    d = heston_d(phi, sigma, rho, u, b)
    g = heston_g(b, rho, sigma, phi, d)
    prefactor = (kappa * theta) / (sigma**2)
    
    ee      = cmath.e**(-d * tau)
    return i * phi * r * tau + prefactor * ( (b - i * rho * sigma * phi - d) * tau   -2. * cmath.log((1.0 - g * ee) / (1. - g)))

def heston_D(phi, tau, sigma, rho, u, b):
    i = 1j
    d = heston_d(phi, sigma, rho, u, b)
    g = heston_g(b, rho, sigma, phi, d)
    
    numerator = b - i * rho * sigma * phi - d
    ee          = cmath.e**(-d * tau)
    frac_factor = (1. - ee) / (1. - g * ee)
    
    return (numerator / (sigma**2)) * frac_factor

def heston_characteristic(phi, x, v0, tau, r,
                          kappa, theta, sigma, rho,
                          eta_v, version):
    i = 1j

    if version == 1:
        u = 0.5
        b = kappa + eta_v - rho * sigma
    if version==2:
        u = -0.5
        b = (kappa + eta_v)
    Cj = heston_C(phi, tau, r, kappa, theta, sigma, rho, u, b)
    Dj = heston_D(phi, tau, sigma, rho, u, b)
    return cmath.exp(Cj + Dj * v0 + i * phi * cmath.log(x))

def heston_integrand(u, x, v0, tau, r,
                     kappa, theta, sigma, rho,
                     eta_v, K, j):
    i = 1j
    if j == 1:
        Phi_shifted = heston_characteristic(u-i, x, v0, tau, r,kappa, theta, sigma, rho, eta_v, version=2)
        Phi_minus_i = heston_characteristic(-i, x, v0, tau, r,kappa, theta, sigma, rho, eta_v, version=2)
        frac = Phi_shifted / (i * u *Phi_minus_i)
        return (frac * cmath.e**(-i * u * cmath.log(K))).real
    if j == 2:
        Phi_u = heston_characteristic(u, x, v0, tau, r,kappa, theta, sigma, rho, eta_v, version=2)
        frac = Phi_u / (i * u)
        return (frac * cmath.e**(-i * u * cmath.log(K))).real

def heston_Pj(j, r, tau, S0, K,
              kappa, theta, sigma, rho,
              eta_v, v0,):
    x = S0
    integrand = lambda u: heston_integrand(
        u, x, v0, tau, r,
        kappa, theta, sigma, rho,
        eta_v, K, j
    )
    integral_value = quad(integrand, 0.,500)[0]
    return 0.5 + (1.0 / np.pi) * integral_value
def heston_call_price(r, T, S0, K,
                      kappa, theta, sigma, rho,
                      eta_v, v0):
    tau = T
    P1 = heston_Pj(1, r, tau, S0, K,
                   kappa, theta, sigma, rho,
                   eta_v, v0)
    P2 = heston_Pj(2, r, tau, S0, K,
                   kappa, theta, sigma, rho,
                   eta_v, v0)
    return S0 * P1 - K * np.exp(-r * T) * P2

kappa = 4.1
sigma = .3
rho = -0.7
v0 = 0.04
theta = 0.06
r = 0.09
s0 = 1.
T = 1
logMoneyness = 0.1

cp2 = heston_call_price(r,T,s0,s0*np.exp(logMoneyness),kappa,theta,sigma,rho,0,v0)
print('My code, price: ' + str(round(cp2, 5)))

# this is the code of [2]
n_steps = 252   # number of time steps
n_paths = 500  # number of paths
n_blocks = 200 # number of blocks
dt = T/n_steps  # time step
q=0

Vc_list = np.zeros(n_blocks) # call array
Vp_list = np.zeros(n_blocks) # put array

for j in range(n_blocks):
    # Correlated normal random variables
    W1, W2 = np.random.multivariate_normal([0,0], [[1, rho], [rho, 1]], (n_steps, n_paths)).T
    
    # Initialize array for variance
    v = np.zeros((n_steps + 1, n_paths)).T
    v[:, 0] = v0
    
    # Initialize array for stock
    S = np.zeros((n_steps + 1, n_paths)).T
    S[:, 0] = s0
    
    # Compute the paths
    for i in range(1, n_steps + 1):
        S[:, i] = S[:, i-1] * np.exp((r - q - 0.5*v[:, i-1])*dt \
                                     + np.sqrt(v[:, i-1])*np.sqrt(dt)*W2[:, i-1])
        v_prev_pos = np.maximum(v[:, i-1], 0.0)
        d_v        = kappa*(theta - v_prev_pos)*dt \
                    + sigma*np.sqrt(v_prev_pos)*np.sqrt(dt)*W1[:, i-1]
        v[:, i]    = v[:, i-1] + d_v
        v[:, i]    = np.maximum(v[:, i],   0.0)
    
    # Compute the discounted option price for the block    
    Vc_list[j] = np.exp(-r*T)*np.mean(np.maximum(S[:,-1] - s0*np.exp(logMoneyness), 0))
    Vp_list[j] = np.exp(-r*T)*np.mean(np.maximum(s0*np.exp(logMoneyness) - S[:,-1], 0))

# Final option price (mean of the prices from each block)
Vc = np.mean(Vc_list)
Vp = np.mean(Vp_list)

print('Code of [1] price: ' + str(round(Vc, 5)))

# this is the code of [3]
def heston_call_px(S0, K, T, r, kappa, v0, theta, eta, rho_sv):
    def _phi(w, t):
        gamma = eta ** 2 / 2
        beta = kappa - rho_sv * eta * w * 1j
        alpha = -(w ** 2 / 2) - (1j * w / 2)
        h = np.sqrt(beta ** 2 - 4 * alpha * gamma)
        r_plus = (beta + h) / (eta ** 2)
        r_minus = (beta - h) / (eta ** 2)
        g = r_minus / r_plus
        eht = np.exp(-h * t)
        D = r_minus * ((1 - eht) / (1 - g * eht))
        C = kappa * (r_minus * t - (2 / (eta ** 2)) * np.log(
            (1 - g * eht) / (1 - g)))
        return np.exp(
            C * theta + D * v0 + 1j * w * np.log(S0 * np.exp(r * t)))

    def _integrand_1(w):
        f = (np.exp(-1j * w * np.log(K)) * _phi(w - 1j, T)) / (
            1j * w * _phi(-1j, T))
        return f.real

    def _integrand_2(w):
        f = (np.exp(-1j * w * np.log(K)) * _phi(w, T)) / (1j * w)
        return f.real

    p1 = 0.5 + (1 / np.pi) * quad(_integrand_1, 0, 100)[0]
    p2 = 0.5 + (1 / np.pi) * quad(_integrand_2, 0, 100)[0]

    return S0 * p1 - np.exp(-r * T) * K * p2

print("Code of [3] price")
print(heston_call_px(s0, s0*np.exp(logMoneyness),T,r,kappa,v0,theta,sigma,rho))
```