from hestonpy.models.heston import Heston

import matplotlib.pyplot as plt
import time
import numpy as np

Pricing with Heston

Initialisation of the model

S0 = 100
V0 = 0.06
r = 0.05

params = {
    "kappa": 1,
    "theta": 0.06,
    "sigma": 0.3,
    "rho": -0.5,
    "drift_emm": 0.00,
}

heston = Heston(spot=S0, vol_initial=V0, r=r, **params)

heston.price_surface()

There are three different ways to price with an Heston object: Monte Carlo, Heston original formula and Carr-Madan

Price via Monte Carlo

Parameters

nbr_points = 100
nbr_simulations = 10**3

Via Euler-Maruyama scheme

start_time = time.time()
result = heston.monte_carlo_price(nbr_points=nbr_points, nbr_simulations=nbr_simulations, strike=100, time_to_maturity=1, scheme="euler")
time_delta = round(time.time() - start_time,4)
price_euler = round(result.price, 2)
std_euler = round(result.std, 2)
infinum_euler = round(result.infinum, 2)
supremum_euler = round(result.supremum, 2)
print(f"Monte Carlo Euler scheme in {time_delta}s : price ${price_euler}, std {std_euler}, and Confidence interval [{infinum_euler},{supremum_euler}]\n")

Variance has been null 0 times over the 100000 iterations (0.0%)
Monte Carlo Euler scheme in 0.0445s : price $12.06, std 0.49, and Confidence interval [12.02,12.11]

Via Milstein scheme

start_time = time.time()
result = heston.monte_carlo_price(nbr_points=nbr_points, nbr_simulations=nbr_simulations, strike=100, time_to_maturity=1, scheme="milstein")
time_delta = round(time.time() - start_time,4)
price_milstein = round(result.price, 2)
std_milstein = round(result.std, 2)
infinum_milstein = round(result.infinum, 2)
supremum_milstein = round(result.supremum, 2)
print(f"Monte Carlo Milstein scheme in {time_delta}s : price ${price_milstein}, std {std_milstein}, and Confidence interval [{infinum_milstein},{supremum_milstein}]\n")

Variance has been null 0 times over the 100000 iterations (0.0%)
Monte Carlo Milstein scheme in 0.0577s : price $12.03, std 0.49, and Confidence interval [11.99,12.08]

Price via Fourier Transform

start_time = time.time()
price_FT, error_FT = heston.fourier_transform_price(strike=100, time_to_maturity=1, error_boolean=True)
time_delta = round(time.time() - start_time,4)
infinum = round(price_FT-error_FT, 3)
supremum = round(price_FT+error_FT, 3)
price_FT = round(price_FT, 3)
error_FT = round(error_FT, 8)
print(f"Fourier Transform in {time_delta}s : price ${price_FT}, error ${error_FT} , and Confidence interval [{infinum},{supremum}]\n")

Fourier Transform in 0.1561s : price $11.952, error $0.0 , and Confidence interval [11.952,11.952]

Price via Carr-Madan formula

start_time = time.time()
price_CM, error_CM = heston.carr_madan_price(strike=100, time_to_maturity=1, error_boolean=True)
time_delta = round(time.time() - start_time,4)
infinum = round(price_CM-error_CM, 3)
supremum = round(price_CM+error_CM, 3)
price_CM = round(price_CM, 3)
error_CM = round(error_CM, 14)
print(f"Carr-Madan in {time_delta}s : price ${price_CM}, error ${error_CM} , and Confidence interval [{infinum},{supremum}]\n")

Carr-Madan in 0.0896s : price $11.952, error $1.447e-10 , and Confidence interval [11.952,11.952]

Other

Path simulations

scheme = 'milstein'
S, V = heston.plot_simulation(time_to_maturity=1, scheme="milstein", nbr_points=252)

Characteristic function

One can do a sanity check of the characteristic functions used in the derivation of the Heston original formula

psi1 = heston.characteristic(j=1)
psi2 = heston.characteristic(j=2)

u = np.arange(start=-20, stop=20,step=0.01)

x = np.log(S0)
v = V0
strike = 100
time_to_maturity = 1

Check the real and imaginary parts

fig, axs = plt.subplots(1, 2, figsize=(10,5), sharey=True)

# Plot real part of psi1 and psi2
axs[0].set_title(r'$\mathfrak{Re}(\psi_1)$ and $\mathfrak{Re}(\psi_2)$')
axs[0].plot(u, np.abs(psi1(x, v, time_to_maturity, u)), label=r'$|\psi_1|$', color='orange', linestyle='--')
axs[0].plot(u, psi1(x, v, time_to_maturity, u).real, label=r'$\psi_1$', color='orange')
axs[0].plot(u, psi2(x, v, time_to_maturity, u).real, label=r'$\psi_2$', color='blue')
axs[0].grid(visible=True)
axs[0].set_xlabel(r'$u$')
axs[0].set_ylabel('Real part')
axs[0].legend()

# Plot imaginary part of psi1 and psi2
axs[1].set_title(r'$\mathfrak{Im}(\psi_1)$ and $\mathfrak{Im}(\psi_2)$')
axs[1].plot(u, np.abs(psi1(x, v, time_to_maturity, u)), label=r'$|\psi_1|$', color='orange', linestyle='--')
axs[1].plot(u, psi1(x, v, time_to_maturity, u).imag, label=r'$\psi_1$', color='orange')
axs[1].plot(u, psi2(x, v, time_to_maturity, u).imag, label=r'$\psi_2$', color='blue')
axs[1].grid(visible=True)
axs[1].set_xlabel(r'$u$')
axs[1].set_ylabel('Imaginary part')
axs[1].legend()

plt.tight_layout()

plt.show()

3D plot of the characteristic functions

fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111, projection='3d')
ax.plot(u, psi1(x, v, time_to_maturity, u).real, psi1(x, v, time_to_maturity, u).imag, label=r'$\psi_1$', color='orange')
ax.plot(u, psi2(x, v, time_to_maturity, u).real, psi2(x, v, time_to_maturity, u).imag, label=r'$\psi_2$', color='blue')
ax.set_xticks([-5*i for i in range(6)] + [5*i for i in range(6)])
ax.set_yticks([-1, 0, 1])
ax.set_zticks([-1, 0, 1])
ax.set_xlabel('u')
ax.set_ylabel('Real part')
ax.set_zlabel('Imaginary part')
plt.legend()
plt.show()

Check the integrability over \(R_+\)

psi1 = heston.characteristic(j=1)
integrand1 = lambda u : np.real((np.exp(-u * np.log(strike) * 1j) * psi1(x, v, time_to_maturity, u))/(u*1j))
psi2 = heston.characteristic(j=2)
integrand2 = lambda u : np.real((np.exp(-u * np.log(strike) * 1j) * psi2(x, v, time_to_maturity, u))/(u*1j))

u = np.arange(start=0, stop=40,step=0.01)

plt.figure()
plt.plot(u, integrand1(u) * u**2, label="Integrand 1")
plt.plot(u, integrand2(u) * u**2, label="Integrand 2")
plt.xlabel(r'u')
plt.ylabel(r'Integrand $\times u^2$')
plt.legend()
plt.grid(visible=True)
plt.title(r'Existence of $Q_1$ and $Q_2$')
plt.show()

/home/theo/Documents/packages/hestonpy/src/hestonpy/models/heston.py:258: RuntimeWarning: divide by zero encountered in divide
  gj = lambda u: (rho * sigma * u * 1j - bj - dj(u)) / (
/home/theo/Documents/packages/hestonpy/src/hestonpy/models/heston.py:258: RuntimeWarning: invalid value encountered in divide
  gj = lambda u: (rho * sigma * u * 1j - bj - dj(u)) / (
/home/theo/Documents/packages/hestonpy/src/hestonpy/models/heston.py:264: RuntimeWarning: invalid value encountered in divide
  - 2 * np.log((1 - gj(u) * np.exp(dj(u) * tau)) / (1 - gj(u)))
/home/theo/Documents/packages/hestonpy/src/hestonpy/models/heston.py:267: RuntimeWarning: invalid value encountered in divide
  lambda tau, u: (bj - rho * sigma * u * 1j + dj(u))