plot_gbm_distribution.py

"""
Figure 4 from the thesis:
2000 geometric Brownian motion paths with log-normal distribution histogram at T=1.

Figure 16: Empirical distribution function of h(X) for different sample sizes.
"""

import numpy as np
import matplotlib.pyplot as plt
from core import sim_brownian_motion, exact_gbm, black_scholes_call

np.random.seed(42)

# --- Figure 4: GBM paths with distribution ---
s0 = 100.0
mu = 0.05
sigma = 0.2
T = 1.0
l = 10
M = 2000

bm = sim_brownian_motion(T, l, M)
t = np.linspace(0, T, 2**l + 1)

gbm = exact_gbm(s0, mu, sigma, t, bm)

fig, axes = plt.subplots(1, 2, figsize=(12, 6), gridspec_kw={'width_ratios': [4, 1]},
                          sharey=True)

# Plot paths
for m in range(M):
    color = 'cyan' if m < M // 2 else 'red'
    axes[0].plot(t, gbm[:, m], color=color, alpha=0.03, linewidth=0.3)

axes[0].set_xlabel('t')
axes[0].set_xlim(0, T)
axes[0].set_title('Geometric Brownian motions')

# Histogram of terminal values
terminal_vals = gbm[-1, :]
axes[1].hist(terminal_vals, bins=50, orientation='horizontal', density=True,
             color='steelblue', edgecolor='none')
axes[1].set_xlabel('Log-normal distribution')
axes[1].set_title('')

fig.tight_layout()
fig.savefig('plot_gbm_distribution.png', dpi=150)
print("Saved: plot_gbm_distribution.png")

# --- Figure 16: Empirical distribution function ---
T0 = 1.0
T1 = 2.0
K = 100.0

fig2, ax2 = plt.subplots(figsize=(8, 5))

for n_samples, color in [(100, 'green'), (1000, 'red'), (10000, 'blue')]:
    # Sample X ~ LogNormal
    X = s0 * np.random.lognormal(mean=(mu - 0.5 * sigma**2) * T0,
                                  sigma=sigma * np.sqrt(T0),
                                  size=n_samples)

    # h(X) = Black-Scholes price conditional on X
    h_X = np.array([black_scholes_call(x, K, T1 - T0, mu, sigma) for x in X])

    # Empirical CDF
    h_sorted = np.sort(h_X)
    ecdf = np.arange(1, n_samples + 1) / n_samples
    ax2.plot(h_sorted, ecdf, color=color, label=f'n = {n_samples}')

ax2.set_xlabel('s')
ax2.set_ylabel('F(s)')
ax2.legend()
ax2.set_title(r'Empirical distribution function of $h(X)$')
fig2.tight_layout()
fig2.savefig('plot_empirical_cdf.png', dpi=150)
print("Saved: plot_empirical_cdf.png")

plt.show()