plot_level_variance.py

"""
Figures 10 & 11 from the thesis:
- Level variance Var(R^(l) - R^(l-1)) for Euler and Milstein MLMC (Figure 10)
- Level mean E(|R^(l) - R^(l-1)|) (Figure 11)

European call option, Black-Scholes model.
"""

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

np.random.seed(42)

s0 = 100.0
mu = 0.05
sigma = 0.2
T0 = 1.0
T1 = 2.0
K = 100.0
L = 8
n_reps = 10000

dt = T1 - T0
level_var_euler = np.zeros(L)
level_var_milstein = np.zeros(L)
level_mean_euler = np.zeros(L)
level_mean_milstein = np.zeros(L)

for l in range(1, L + 1):
    n_fine = 2**l
    dt_fine = dt / n_fine
    n_coarse = 2**(l - 1)
    dt_coarse = dt / n_coarse

    bm = sim_brownian_motion(dt, l, n_reps)
    bm_coarse = bm[::2, :]

    # --- Euler ---
    S_fine = np.full(n_reps, s0, dtype=float)
    for k in range(1, n_fine + 1):
        dW = bm[k, :] - bm[k - 1, :]
        S_fine = S_fine + mu * S_fine * dt_fine + sigma * S_fine * dW

    S_coarse = np.full(n_reps, s0, dtype=float)
    for k in range(1, n_coarse + 1):
        dW_c = bm_coarse[k, :] - bm_coarse[k - 1, :]
        S_coarse = S_coarse + mu * S_coarse * dt_coarse + sigma * S_coarse * dW_c

    R_fine = european_call_payoff(S_fine, K)
    R_coarse = european_call_payoff(S_coarse, K)
    diff_euler = R_fine - R_coarse

    level_var_euler[l - 1] = np.var(diff_euler)
    level_mean_euler[l - 1] = np.mean(np.abs(diff_euler))

    # --- Milstein ---
    S_fine_m = np.full(n_reps, s0, dtype=float)
    for k in range(1, n_fine + 1):
        dW = bm[k, :] - bm[k - 1, :]
        S_fine_m = (S_fine_m + mu * S_fine_m * dt_fine
                    + sigma * S_fine_m * dW
                    + 0.5 * sigma**2 * S_fine_m * (dW**2 - dt_fine))

    S_coarse_m = np.full(n_reps, s0, dtype=float)
    for k in range(1, n_coarse + 1):
        dW_c = bm_coarse[k, :] - bm_coarse[k - 1, :]
        S_coarse_m = (S_coarse_m + mu * S_coarse_m * dt_coarse
                      + sigma * S_coarse_m * dW_c
                      + 0.5 * sigma**2 * S_coarse_m * (dW_c**2 - dt_coarse))

    R_fine_m = european_call_payoff(S_fine_m, K)
    R_coarse_m = european_call_payoff(S_coarse_m, K)
    diff_milstein = R_fine_m - R_coarse_m

    level_var_milstein[l - 1] = np.var(diff_milstein)
    level_mean_milstein[l - 1] = np.mean(np.abs(diff_milstein))

    print(f"l={l}: Var_E={level_var_euler[l-1]:.6e}, Var_M={level_var_milstein[l-1]:.6e}")

levels = np.arange(1, L + 1)

# Compute slopes
slope_var_e = np.mean(np.diff(np.log2(level_var_euler)) / np.diff(levels))
slope_var_m = np.mean(np.diff(np.log2(level_var_milstein)) / np.diff(levels))
slope_mean_e = np.mean(np.diff(np.log2(level_mean_euler)) / np.diff(levels))
slope_mean_m = np.mean(np.diff(np.log2(level_mean_milstein)) / np.diff(levels))

# --- Figure 10: Level Variance ---
fig1, ax1 = plt.subplots(figsize=(8, 6))
ax1.semilogy(levels, level_var_euler, 'b-*', base=2, label=r'$\hat{Z}$ (Euler)')
ax1.semilogy(levels, level_var_milstein, 'r-*', base=2, label=r'$\check{Z}$ (Milstein)')
ax1.text(5, level_var_euler[4], f'slope={slope_var_e:.5f}', fontsize=10, color='blue')
ax1.text(5, level_var_milstein[4], f'slope={slope_var_m:.5f}', fontsize=10, color='red')
ax1.set_xlabel('l')
ax1.set_ylabel(r'Var$(R^{(l)} - R^{(l-1)})$')
ax1.legend()
ax1.set_title('Level Variance')
fig1.tight_layout()
fig1.savefig('plot_level_variance.png', dpi=150)
print(f"\nVariance slopes: Euler={slope_var_e:.4f} (exp: -1), Milstein={slope_var_m:.4f} (exp: -2)")
print("Saved: plot_level_variance.png")

# --- Figure 11: Level Mean ---
fig2, ax2 = plt.subplots(figsize=(8, 6))
ax2.semilogy(levels, level_mean_euler, 'b-*', base=2, label=r'$\hat{Z}$ (Euler)')
ax2.semilogy(levels, level_mean_milstein, 'r-*', base=2, label=r'$\check{Z}$ (Milstein)')
ax2.text(5, level_mean_euler[4], f'slope={slope_mean_e:.5f}', fontsize=10, color='blue')
ax2.text(5, level_mean_milstein[4], f'slope={slope_mean_m:.5f}', fontsize=10, color='red')
ax2.set_xlabel('l')
ax2.set_ylabel(r'$E(|R^{(l)} - R^{(l-1)}|)$')
ax2.legend()
ax2.set_title('Level Mean')
fig2.tight_layout()
fig2.savefig('plot_level_mean.png', dpi=150)
print(f"Mean slopes: Euler={slope_mean_e:.4f} (exp: -0.5), Milstein={slope_mean_m:.4f} (exp: -1)")
print("Saved: plot_level_mean.png")

plt.show()