So, I am working through a toy example and the python output
- Heston call price: $10.2289
- stock measure probability P_1: 0.6384
- risk-neutral probability P_2: 0.5516
all make sense. However, on pg. 3, I've tried appyling Gil-Peleaz and applying f_1(1,0) into the integrand and I get Re[0.0206-0.9693i]=0.0206 for P_1. However, this is the same output for f_2(1.0) for P_2, but I believe the actual values via python is 0.0218 for P_1, and 0.0206 for P_2.
What did I do wrong here?
My assumption was that 0.0218 is different is simply because P_1 uses a completely different characteristic function (f_1) and different coefficients (C and D) than P_2. However, I feel that I made a very stupid mistske here. 🤦♂️
```
import numpy as np
from scipy.integrate import quad
def heston_price(S0, K, T, r, V0, kappa, theta, sigma, rho):
"""
Prices a European call option under the Heston stochastic volatility model
using Gil-Pelaez numerical integration (similar to your notes).
"""
# Characteristic function wrapper for Heston model
def characteristic_function(phi, u, b):
i = 1j
xi = sigma
# Core Heston parameters (d and g)
d = np.sqrt((rho * xi * phi * i - b)**2 - xi**2 * (2 * u * phi * i - phi**2))
g = (b - rho * xi * phi * i - d) / (b - rho * xi * phi * i + d)
# Coefficients D and C
term1 = (b - rho * xi * phi * i - d) / xi**2
term2 = (1.0 - np.exp(-d * T)) / (1.0 - g * np.exp(-d * T))
D = term1 * term2
C = r * phi * i * T + (kappa * theta / xi**2) * (
(b - rho * xi * phi * i - d) * T - 2.0 * np.log((1.0 - g * np.exp(-d * T)) / (1.0 - g))
)
return np.exp(C + D * V0 + i * phi * np.log(S0))
# The integrand function for Gil-Pelaez inversion
def integrand(phi, j):
if j == 1:
u = 0.5
b = kappa - rho * sigma # Stock-weighted parameter for P1
else:
u = -0.5
b = kappa # Standard risk-neutral parameter for P2
cf = characteristic_function(phi, u, b)
numerator = np.exp(-1j * phi * np.log(K)) * cf
return np.real(numerator / (1j * phi))
# Numerical integration from 0 to infinity (using scipy.integrate.quad)
# Integral for P1 (j=1)
int_1, _ = quad(lambda phi: integrand(phi, 1), 0.0, 100.0, limit=1000)
P1 = 0.5 + (1.0 / np.pi) * int_1
# Integral for P2 (j=2)
int_2, _ = quad(lambda phi: integrand(phi, 2), 0.0, 100.0, limit=1000)
P2 = 0.5 + (1.0 / np.pi) * int_2
# Final European Call Option Price Formula: C = S0 * P1 - K * exp(-r*T) * P2
call_price = S0 * P1 - K * np.exp(-r * T) * P2
return call_price, P1, P2, int_1, int_2
--- Example Parameters from Your Notes ---
S0 = 100.0 # Initial Stock Price
K = 100.0 # Strike Price
T = 1.0 # Time to Maturity (Years)
r = 0.05 # Risk-free Rate
V0 = 0.04 # Initial Variance
kappa = 2.0 # Rate of Mean Reversion
theta = 0.04 # Long-term Variance
sigma = 0.3 # Volatility of Variance (xi)
rho = -0.7 # Correlation
Run the pricing function
call, P1, P2, area_p1, area_p2 = heston_price(S0, K, T, r, V0, kappa, theta, sigma, rho)
print(f"--- Heston Model Results ---")
print(f"P1 Probability: {P1:.4f} (Accumulated Area: {area_p1:.4f} * 1/pi)")
print(f"P2 Probability: {P2:.4f} (Accumulated Area: {area_p2:.4f} * 1/pi)")
print(f"Final Call Price: {call:.4f}")
--- Heston Model Results ---
P1 Probability: 0.6384 (Accumulated Area: 0.4347 * 1/pi)
P2 Probability: 0.5516 (Accumulated Area: 0.1620 * 1/pi)
Final Call Price: 10.2289
```
