Risk Neutrality P vs Q

play.py

# The script simulates a simplified Heston model under both the physical measure P and risk-neutral measure Q, using the same stock drift so that differences come only from the volatility risk premium.

# It plots:
# - The terminal stock-price distributions S_T under P and Q.
# - Their left and right tail probabilities.
# - The latent volatility distributions sqrt(v_t), both pooled across all simulated times and at maturity only.

# The P-to-Q change is implemented through Heston's volatility-risk parameter lambda, which modifies the variance mean-reversion speed kappa and long-run variance theta.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

rng = np.random.default_rng(12345)

# ------------------------------------------------------------
# Parameters
# ------------------------------------------------------------

S0 = 100.0
T = 1.0
steps = 100
dt = T / steps
N = 200_000

# Same stock drift under P and Q so only volatility-risk
# adjustment affects the terminal distribution
mu = 0.03

v0 = 0.04
sigma = 0.45
rho = 0.0

# Physical variance dynamics
kappa_P = 3.0
theta_P = 0.04

# Heston-style volatility-risk-premium adjustment.
# Convention:
#   kappa_Q = kappa_P + lambda
#   theta_Q = kappa_P * theta_P / kappa_Q
lam = -1.0

kappa_Q = kappa_P + lam
theta_Q = kappa_P * theta_P / kappa_Q

print("Parameters:")
print(f"P: kappa={kappa_P:.3f}, theta={theta_P:.4f}, "
      f"long-run vol={np.sqrt(theta_P):.2%}")
print(f"Q: kappa={kappa_Q:.3f}, theta={theta_Q:.4f}, "
      f"long-run vol={np.sqrt(theta_Q):.2%}")
print()


# ------------------------------------------------------------
# Simulation
# ------------------------------------------------------------

def simulate(kappa, theta):
  v = np.full(N, v0)
  logS = np.full(N, np.log(S0))

  # Store a subsample of v at every time point.
  # This approximates the time-pooled distribution
  #   (1/T) integral f(v_t) dt.
  pooled_v = []
  pooled_sample_size = 3000

  for _ in range(steps):
    vp = np.maximum(v, 0.0)

    zv = rng.standard_normal(N)
    zs = rng.standard_normal(N)

    # rho = 0, so stock and variance innovations are independent.
    # Correct log-price Heston dynamics include the Ito correction.
    logS += (
      (mu - 0.5 * vp) * dt
      + np.sqrt(vp * dt) * zs
    )

    # Full-truncation Euler approximation of CIR variance.
    v += (
      kappa * (theta - vp) * dt
      + sigma * np.sqrt(vp * dt) * zv
    )
    v = np.maximum(v, 0.0)

    idx = rng.choice(N, pooled_sample_size, replace=False)
    pooled_v.append(v[idx].copy())

  return (
    np.exp(logS),
    v.copy(),
    np.concatenate(pooled_v)
  )


ST_P, vT_P, pooled_v_P = simulate(kappa_P, theta_P)
ST_Q, vT_Q, pooled_v_Q = simulate(kappa_Q, theta_Q)


# ------------------------------------------------------------
# Utility for KDE
# ------------------------------------------------------------

def subsample(x, n=80_000):
  if len(x) <= n:
    return x
  return x[rng.choice(len(x), n, replace=False)]


# ------------------------------------------------------------
# Plot 1: terminal price distribution
# ------------------------------------------------------------

lo = min(
  np.quantile(ST_P, 0.001),
  np.quantile(ST_Q, 0.001)
)

hi = max(
  np.quantile(ST_P, 0.999),
  np.quantile(ST_Q, 0.999)
)

x = np.linspace(lo, hi, 700)

kde_P = gaussian_kde(subsample(ST_P))
kde_Q = gaussian_kde(subsample(ST_Q))

plt.figure(figsize=(8, 5))

plt.plot(x, kde_P(x), label="Physical P")
plt.plot(x, kde_Q(x), label="Risk-neutral Q")

plt.axvline(
  S0 * np.exp(mu * T),
  linestyle="--",
  label="Common forward / mean"
)

plt.xlabel(r"$S_T$")
plt.ylabel("Density")
plt.title(
  "Terminal Heston distribution\n"
  "same stock drift, rho = 0"
)
plt.legend()
plt.tight_layout()

plt.savefig(
  "heston_terminal_distribution_P_vs_Q.png",
  dpi=180
)

plt.show()


# ------------------------------------------------------------
# Plot 2: terminal tail probabilities
# ------------------------------------------------------------

left = np.linspace(
  np.quantile(ST_P, 0.002),
  np.quantile(ST_P, 0.40),
  150
)

right = np.linspace(
  np.quantile(ST_P, 0.60),
  max(
    np.quantile(ST_P, 0.998),
    np.quantile(ST_Q, 0.998)
  ),
  150
)

cdf_P = np.array([(ST_P <= z).mean() for z in left])
cdf_Q = np.array([(ST_Q <= z).mean() for z in left])

sf_P = np.array([(ST_P >= z).mean() for z in right])
sf_Q = np.array([(ST_Q >= z).mean() for z in right])

plt.figure(figsize=(8, 5))

plt.semilogy(
  left,
  np.maximum(cdf_P, 1 / N),
  label="P: left tail"
)

plt.semilogy(
  left,
  np.maximum(cdf_Q, 1 / N),
  label="Q: left tail"
)

plt.semilogy(
  right,
  np.maximum(sf_P, 1 / N),
  linestyle="--",
  label="P: right tail"
)

plt.semilogy(
  right,
  np.maximum(sf_Q, 1 / N),
  linestyle="--",
  label="Q: right tail"
)

plt.xlabel(r"$S_T$")
plt.ylabel("Tail probability")
plt.title("Terminal tail probabilities")
plt.legend()
plt.tight_layout()

plt.savefig(
  "heston_terminal_tails_P_vs_Q.png",
  dpi=180
)

plt.show()


# ------------------------------------------------------------
# Plot 3: latent volatility distribution
# ------------------------------------------------------------

# Plot volatility sqrt(v), rather than variance v.
vol_pool_P = np.sqrt(pooled_v_P)
vol_pool_Q = np.sqrt(pooled_v_Q)

vol_T_P = np.sqrt(vT_P)
vol_T_Q = np.sqrt(vT_Q)

xmax = max(
  np.quantile(vol_pool_P, 0.995),
  np.quantile(vol_pool_Q, 0.995),
  np.quantile(vol_T_P, 0.995),
  np.quantile(vol_T_Q, 0.995)
)

xv = np.linspace(0, xmax, 700)

k_pool_P = gaussian_kde(subsample(vol_pool_P))
k_pool_Q = gaussian_kde(subsample(vol_pool_Q))

k_T_P = gaussian_kde(subsample(vol_T_P))
k_T_Q = gaussian_kde(subsample(vol_T_Q))

plt.figure(figsize=(8, 5))

plt.plot(
  xv,
  k_pool_P(xv),
  label="P: pooled over all times"
)

plt.plot(
  xv,
  k_pool_Q(xv),
  label="Q: pooled over all times"
)

plt.plot(
  xv,
  k_T_P(xv),
  linestyle="--",
  label="P: terminal only"
)

plt.plot(
  xv,
  k_T_Q(xv),
  linestyle="--",
  label="Q: terminal only"
)

plt.xlabel(r"Latent volatility $\sqrt{v_t}$")
plt.ylabel("Density")
plt.title("Heston latent-volatility distribution: P vs Q")
plt.legend()
plt.tight_layout()

plt.savefig(
  "heston_latent_volatility_P_vs_Q.png",
  dpi=180
)

plt.show()


# ------------------------------------------------------------
# Summary
# ------------------------------------------------------------

print("Terminal stock:")
print(f"E_P[S_T] = {ST_P.mean():.4f}")
print(f"E_Q[S_T] = {ST_Q.mean():.4f}")
print(f"Theoretical common mean = {S0 * np.exp(mu*T):.4f}")
print()

print("Latent volatility:")
print(f"P pooled mean vol   = {vol_pool_P.mean():.2%}")
print(f"Q pooled mean vol   = {vol_pool_Q.mean():.2%}")
print(f"P terminal mean vol = {vol_T_P.mean():.2%}")
print(f"Q terminal mean vol = {vol_T_Q.mean():.2%}")