Optimal Hedging — Simulations

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso
import torch
import torch.nn as nn
import torch.optim as optim
from scipy.stats import norm

FIGDIR = "../markdown/figures"
np.random.seed(42)
plt.rcParams.update({'font.size': 10, 'axes.spines.top': False, 'axes.spines.right': False})

# ── BSM ──────────────────────────────────────────────────────────────
def bsm_price(S, K, T, sigma, r=0.0, cp=1):
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    if cp == 1:
        return S * norm.cdf(d1) - K * np.exp(-r*T) * norm.cdf(d2)
    return K * np.exp(-r*T) * norm.cdf(-d2) - S * norm.cdf(-d1)

def bsm_delta(S, K, T, sigma, r=0.0, cp=1):
    if T <= 0:
        return float(cp == 1 and S > K)
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T))
    return norm.cdf(d1) if cp == 1 else norm.cdf(d1) - 1

def bsm_gamma(S, K, T, sigma, r=0.0):
    if T <= 0:
        return 0.0
    d1 = (np.log(S/K) + (r + 0.5*sigma**2)*T) / (sigma * np.sqrt(T))
    return norm.pdf(d1) / (S * sigma * np.sqrt(T))

# ── Yield curve setup ─────────────────────────────────────────────────
D = 20                              # number of maturity buckets
maturities = np.linspace(0.5, 10, D)

def pc_level(t):  return np.ones_like(t)
def pc_slope(t):  return (t - t.mean()) / (t.max() - t.min())
def pc_curve(t):
    s = pc_slope(t)
    return s**2 - (s**2).mean()

def dv01_bond(maturity, coupon=0.05, notional=1e6):
    """Dollar DV01 vector across maturity grid (dollar per unit yield change)."""
    cash_times = np.arange(0.5, maturity + 0.01, 0.5)
    r = 0.04
    d = np.zeros(D)
    for t_cf in cash_times:
        cf = coupon * 0.5 * notional if t_cf < maturity else (coupon * 0.5 + 1) * notional
        pv = cf * np.exp(-r * t_cf)
        idx_m = np.argmin(np.abs(maturities - t_cf))
        d[idx_m] += pv * t_cf
    return d

# ── CVaR ──────────────────────────────────────────────────────────────
def cvar95(x):
    q = np.percentile(x, 5)
    return -np.mean(x[x <= q])

T_obs = 1000
sigma_level, sigma_slope, sigma_curve, sigma_idio = 0.006, 0.003, 0.001, 0.0005

U = np.column_stack([
    pc_level(maturities)  / np.linalg.norm(pc_level(maturities)),
    pc_slope(maturities)  / np.linalg.norm(pc_slope(maturities)),
    pc_curve(maturities)  / np.linalg.norm(pc_curve(maturities)),
])

np.random.seed(42)
F = np.column_stack([
    np.random.randn(T_obs) * sigma_level,
    np.random.randn(T_obs) * sigma_slope,
    np.random.randn(T_obs) * sigma_curve,
])
E  = np.random.randn(T_obs, D) * sigma_idio
dR = F @ U.T + E

# PCA from realised covariance
Sigma = np.cov(dR.T)
eigvals, eigvecs = np.linalg.eigh(Sigma)
idx = np.argsort(eigvals)[::-1]
eigvals, eigvecs = eigvals[idx], eigvecs[:, idx]

# DV01 vectors
d_7y, d_2y, d_5y, d_10y = dv01_bond(7.), dv01_bond(2.), dv01_bond(5.), dv01_bond(10.)

pnl_7y  = dR @ d_7y
pnl_2y  = dR @ d_2y
pnl_5y  = dR @ d_5y
pnl_10y = dR @ d_10y

# DV01 hedge (5y)
h_dv01 = np.sum(d_7y) / np.sum(d_5y)
pnl_dv01_hedged = pnl_7y - h_dv01 * pnl_5y

# 3-factor PCA hedge
hedge_pnls  = np.column_stack([pnl_2y, pnl_5y, pnl_10y])
hedge_pvecs = np.column_stack([d_2y,  d_5y,   d_10y])
d7_factors     = eigvecs[:, :3].T @ d_7y
hedge_factors  = eigvecs[:, :3].T @ hedge_pvecs
h_pca          = np.linalg.solve(hedge_factors, d7_factors)
pnl_pca_hedged = pnl_7y - hedge_pnls @ h_pca

cumvar = np.cumsum(eigvals) / np.sum(eigvals)

sd_nohedge, sd_dv01, sd_pca = np.std(pnl_7y), np.std(pnl_dv01_hedged), np.std(pnl_pca_hedged)
print(f"σ: no hedge=${sd_nohedge:,.0f}  DV01=${sd_dv01:,.0f}  PCA=${sd_pca:,.0f}")

fig, axes = plt.subplots(2, 2, figsize=(12, 8))

ax = axes[0, 0]
r_base = np.linspace(0.02, 0.05, D)
for i in range(0, T_obs, 100):
    ax.plot(maturities, (r_base + np.cumsum(dR[i]) * 5) * 100,
            color='steelblue', alpha=0.2, linewidth=0.7)
ax.set(xlabel='Maturity (years)', ylabel='Yield (%)', title='Sample Yield Curves')

ax = axes[0, 1]
for k, (lbl, col) in enumerate(zip(
        ['Level (PC1)', 'Slope (PC2)', 'Curvature (PC3)'],
        ['steelblue', 'darkorange', 'seagreen'])):
    v = eigvecs[:, k]
    ax.plot(maturities, v / np.abs(v).max(), label=lbl, color=col, linewidth=2)
ax.axhline(0, color='grey', linewidth=0.5, linestyle='--')
ax.set(xlabel='Maturity (years)', ylabel='Loading (normalised)', title='PCA Eigenvectors')
ax.legend(fontsize=9)

ax = axes[1, 0]
ax.hist(pnl_7y,         bins=60, alpha=0.5, color='grey',       label='No hedge',        density=True)
ax.hist(pnl_dv01_hedged,bins=60, alpha=0.6, color='steelblue',  label='DV01 hedge (5y)', density=True)
ax.hist(pnl_pca_hedged, bins=60, alpha=0.6, color='darkorange', label='PCA hedge (3f)',  density=True)
ax.set(xlabel='Daily P&L (USD)', ylabel='Density', title='Hedge Performance')
ax.legend(fontsize=9)
ax.text(0.98, 0.97,
        f'σ no hedge: ${sd_nohedge:,.0f}\nσ DV01:     ${sd_dv01:,.0f}\nσ PCA:      ${sd_pca:,.0f}',
        transform=ax.transAxes, ha='right', va='top', fontsize=8, family='monospace',
        bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

ax = axes[1, 1]
ax.plot(np.arange(1, D+1), cumvar * 100, color='steelblue', linewidth=2, marker='o', markersize=4)
ax.axhline(95, color='grey',       linestyle='--', linewidth=0.8, label='95% threshold')
ax.axvline(3,  color='darkorange', linestyle='--', linewidth=0.8, label='3 factors')
ax.set(xlabel='Number of PC factors', ylabel='Cumulative explained variance (%)',
       title='PCA Explained Variance')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/hedg_pca_factors.png', dpi=150, bbox_inches='tight')
plt.show()

hedge_mats = list(np.arange(1.0, 6.5, 0.5)) + [8.0, 9.0, 10.0]
K_hedge    = len(hedge_mats)
pvecs_lasso = np.column_stack([dv01_bond(m) for m in hedge_mats])
pnl_hedges  = dR @ pvecs_lasso

Sigma_HH = np.cov(pnl_hedges.T)
sigma_XH = np.cov(pnl_7y, pnl_hedges.T)[0, 1:]

h_mv   = np.linalg.solve(Sigma_HH, sigma_XH)
pnl_mv = pnl_7y - pnl_hedges @ h_mv

lambdas         = np.logspace(-2, 1, 80) * np.std(pnl_7y)
lasso_coefs     = []
residual_vars   = []
gross_notionals = []
for lam in lambdas:
    L = Lasso(alpha=lam, fit_intercept=False, max_iter=10000)
    L.fit(pnl_hedges, pnl_7y)
    lasso_coefs.append(L.coef_.copy())
    residual_vars.append(np.var(pnl_7y - pnl_hedges @ L.coef_))
    gross_notionals.append(np.sum(np.abs(L.coef_)))

lasso_coefs   = np.array(lasso_coefs)
var_nohedge   = np.var(pnl_7y)
var_mv        = np.var(pnl_mv)
nz_counts     = np.sum(np.abs(lasso_coefs) > 1e-6, axis=1)

fig, axes = plt.subplots(1, 3, figsize=(14, 4))

ax = axes[0]
colors_lasso = plt.cm.tab10(np.linspace(0, 1, K_hedge))
for j in range(K_hedge):
    ax.plot(lambdas / np.std(pnl_7y), lasso_coefs[:, j],
            color=colors_lasso[j], linewidth=1.5)
ax.axvline(1.0, color='grey', linestyle='--', linewidth=0.8)
ax.set_xscale('log')
ax.set(xlabel='λ / σ_X', ylabel='Hedge coefficient', title='Lasso Hedge Path')

ax = axes[1]
ax.plot(lambdas / np.std(pnl_7y),
        np.array(residual_vars) / var_nohedge * 100, color='steelblue', linewidth=2)
ax.axhline(var_mv / var_nohedge * 100, color='darkorange', linestyle='--', linewidth=1.5,
           label=f'Min-variance ({var_mv/var_nohedge*100:.0f}%)')
ax.axhline(100, color='grey', linestyle=':', linewidth=1, label='No hedge (100%)')
ax.set_xscale('log')
ax.set(xlabel='λ / σ_X', ylabel='Residual variance (% of unhedged)',
       title='Residual Variance vs Sparsity')
ax.legend(fontsize=9)

ax = axes[2]
sc = ax.scatter(gross_notionals, np.array(residual_vars) / var_nohedge * 100,
                c=nz_counts, cmap='plasma', s=25, zorder=3)
ax.axhline(var_mv / var_nohedge * 100, color='darkorange', linestyle='--',
           linewidth=1.2, label='Min-variance')
plt.colorbar(sc, ax=ax, label='Active instruments')
ax.set(xlabel='Gross notional of hedge', ylabel='Residual variance (%)',
       title='Efficiency Frontier')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/hedg_lasso.png', dpi=150, bbox_inches='tight')
plt.show()

def phi_star(p, cost_to_risk):
    return p / (1 + cost_to_risk)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

probs = np.linspace(0, 1, 200)
ax = axes[0]
for ratio, lbl, col in [(0.1, 'c/risk = 0.1', 'steelblue'),
                         (0.5, 'c/risk = 0.5', 'darkorange'),
                         (1.0, 'c/risk = 1.0', 'seagreen'),
                         (2.0, 'c/risk = 2.0', 'firebrick')]:
    ax.plot(probs, phi_star(probs, ratio), label=lbl, color=col, linewidth=2)
ax.plot(probs, probs, color='grey', linestyle='--', linewidth=1, label='φ* = p (no cost)')
ax.set(xlabel='Probability of order materialising (p)',
       ylabel='Optimal pre-hedge fraction φ*',
       title='Pre-hedge Fraction vs Order Probability',
       xlim=(0, 1), ylim=(0, 1))
ax.legend(fontsize=9)

cost_ratios = np.linspace(0, 3, 200)
ax = axes[1]
for prob, lbl, col in [(1.0,  'p = 1.0',  'steelblue'),
                        (0.75, 'p = 0.75', 'darkorange'),
                        (0.5,  'p = 0.50', 'seagreen'),
                        (0.25, 'p = 0.25', 'firebrick')]:
    ax.plot(cost_ratios, phi_star(prob, cost_ratios), label=lbl, color=col, linewidth=2)
ax.set(xlabel='Cost-to-risk ratio $c / (\\gamma \\sigma^2 T)$',
       ylabel='Optimal pre-hedge fraction φ*',
       title='Pre-hedge Fraction vs Transaction Cost',
       xlim=(0, 3), ylim=(0, 1))
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/hedg_prehedge.png', dpi=150, bbox_inches='tight')
plt.show()

S0, K_opt, T_opt, sigma_opt, r_opt = 100., 100., 1.0, 0.2, 0.0
dt = 1 / 252
n_steps = int(T_opt / dt)
n_paths = 5000
np.random.seed(42)

dW = np.random.randn(n_paths, n_steps) * np.sqrt(dt)
S_paths = S0 * np.exp(np.cumsum((r_opt - 0.5*sigma_opt**2)*dt + sigma_opt*dW, axis=1))
S_paths = np.column_stack([np.full(n_paths, S0), S_paths])

option_price0 = bsm_price(S0, K_opt, T_opt, sigma_opt, r_opt)

def sim_dh(S_paths, rebal_freq, tc=0.0):
    n_paths, n_full = S_paths.shape
    n_steps_full = n_full - 1
    pnls = np.zeros(n_paths)
    for i in range(n_paths):
        S      = S_paths[i]
        t_rebal = np.arange(0, n_steps_full + 1, rebal_freq)
        cash   = option_price0
        delta  = 0.0
        for step in t_rebal[:-1]:
            t = step * dt
            T_rem = T_opt - t
            S_t   = S[step]
            dn    = bsm_delta(S_t, K_opt, T_rem, sigma_opt, r_opt)
            cash -= tc * abs(dn - delta) * S_t
            cash -= (dn - delta) * S_t
            cash *= np.exp(r_opt * rebal_freq * dt)
            delta = dn
        cash  += delta * S[-1]
        pnls[i] = cash - max(S[-1] - K_opt, 0)
    return pnls

pnl_daily   = sim_dh(S_paths, 1,  tc=0.0)
pnl_weekly  = sim_dh(S_paths, 5,  tc=0.0)
pnl_monthly = sim_dh(S_paths, 21, tc=0.0)
pnl_ww      = sim_dh(S_paths, 1,  tc=0.001)

# Gamma–theta decomposition along a single path
S_path = S_paths[42]
pnl_gt = []
for step in range(n_steps):
    T_rem = T_opt - step * dt - dt
    if T_rem <= 0:
        break
    S_t   = S_path[step]
    sig_r = abs(np.log(S_path[step+1] / S_t)) / np.sqrt(dt)
    g     = bsm_gamma(S_t, K_opt, T_rem + dt, sigma_opt, r_opt)
    pnl_gt.append(0.5 * g * S_t**2 * (sig_r**2 - sigma_opt**2) * dt)
pnl_gt = np.array(pnl_gt)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
bins = np.linspace(-10, 10, 80)
ax.hist(pnl_daily,   bins=bins, alpha=0.5, density=True, color='steelblue',
        label=f'Daily (σ={np.std(pnl_daily):.2f})')
ax.hist(pnl_weekly,  bins=bins, alpha=0.5, density=True, color='darkorange',
        label=f'Weekly (σ={np.std(pnl_weekly):.2f})')
ax.hist(pnl_monthly, bins=bins, alpha=0.5, density=True, color='seagreen',
        label=f'Monthly (σ={np.std(pnl_monthly):.2f})')
ax.hist(pnl_ww,      bins=bins, alpha=0.5, density=True, color='firebrick',
        label=f'WW band tc=0.1% (σ={np.std(pnl_ww):.2f})')
ax.axvline(0, color='black', linewidth=0.5)
ax.set(xlabel='Hedging P&L ($)', ylabel='Density', title='Delta Hedging P&L Distributions')
ax.legend(fontsize=8)

ax = axes[1]
t_ax = np.arange(len(pnl_gt)) * dt
ax.bar(t_ax, pnl_gt, width=dt,
       color=np.where(pnl_gt > 0, 'steelblue', 'firebrick'), alpha=0.7,
       label='Gamma-theta P&L per step')
ax.plot(t_ax, np.cumsum(pnl_gt), color='black', linewidth=1.5, label='Cumulative P&L')
ax.axhline(0, color='grey', linewidth=0.5)
ax.set(xlabel='Time (years)', ylabel='P&L ($)',
       title='Gamma–Theta P&L Decomposition (Single Path)')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/hedg_delta_hedging.png', dpi=150, bbox_inches='tight')
plt.show()

torch.manual_seed(42)
np.random.seed(42)

n_steps_dh = 30
dt_dh      = T_opt / n_steps_dh
tc_dh      = 0.001
n_train    = 4000
n_test     = 2000
gamma_ra   = 0.1

def sim_gbm_torch(n_paths, n_steps, S0, sigma, r, dt, seed=None):
    if seed is not None:
        torch.manual_seed(seed)
    dW  = torch.randn(n_paths, n_steps) * np.sqrt(dt)
    log_S = (r - 0.5*sigma**2)*dt + sigma*dW
    S   = S0 * torch.exp(torch.cumsum(log_S, dim=1))
    return torch.cat([torch.full((n_paths, 1), S0), S], dim=1)

class HedgingPolicy(nn.Module):
    def __init__(self, hidden=64):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(3, hidden), nn.Tanh(),
            nn.Linear(hidden, hidden), nn.Tanh(),
            nn.Linear(hidden, 1), nn.Sigmoid()
        )
    def forward(self, moneyness, t_rem, prev_delta):
        x = torch.stack([moneyness, t_rem, prev_delta], dim=1)
        return self.net(x).squeeze(-1)

def compute_pnl_dh(S, policy, premium, n_steps, dt, tc):
    n_paths = S.shape[0]
    cash    = torch.full((n_paths,), premium)
    delta   = torch.zeros(n_paths)
    for t in range(n_steps):
        t_rem     = torch.tensor((n_steps - t) / n_steps, dtype=torch.float32).expand(n_paths)
        delta_new = policy(torch.log(S[:, t] / K_opt), t_rem, delta.detach())
        dd        = delta_new - delta
        cash      = cash - tc * torch.abs(dd) * S[:, t] - dd * S[:, t]
        cash      = cash * np.exp(r_opt * dt)
        delta     = delta_new
    cash += delta * S[:, -1]
    return cash - torch.clamp(S[:, -1] - K_opt, min=0)

premium0 = float(bsm_price(S0, K_opt, T_opt, sigma_opt, r_opt))
S_train  = sim_gbm_torch(n_train, n_steps_dh, S0, sigma_opt, r_opt, dt_dh, seed=42)

policy    = HedgingPolicy(hidden=64)
optimizer = optim.Adam(policy.parameters(), lr=1e-3)
scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=100, gamma=0.8)
losses    = []

for epoch in range(300):
    idx   = torch.randperm(n_train)[:512]
    pnl   = compute_pnl_dh(S_train[idx], policy, premium0, n_steps_dh, dt_dh, tc_dh)
    loss  = -pnl.mean() + gamma_ra * pnl.var()
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
    scheduler.step()
    if epoch % 30 == 0:
        losses.append(loss.item())

S_test = sim_gbm_torch(n_test, n_steps_dh, S0, sigma_opt, r_opt, dt_dh, seed=99)
with torch.no_grad():
    pnl_dh = compute_pnl_dh(S_test, policy, premium0, n_steps_dh, dt_dh, tc_dh).numpy()

def bsm_hedge_pnl(S_np, premium, n_steps, dt, tc):
    n_paths = S_np.shape[0]
    pnls    = np.zeros(n_paths)
    for i in range(n_paths):
        cash  = premium
        delta = 0.0
        for t in range(n_steps):
            T_rem    = (n_steps - t) * dt
            dn       = bsm_delta(S_np[i, t], K_opt, T_rem, sigma_opt, r_opt)
            cash    -= tc * abs(dn - delta) * S_np[i, t] + (dn - delta) * S_np[i, t]
            cash    *= np.exp(r_opt * dt)
            delta    = dn
        cash    += delta * S_np[i, -1]
        pnls[i]  = cash - max(S_np[i, -1] - K_opt, 0)
    return pnls

pnl_bsm = bsm_hedge_pnl(S_test.numpy(), premium0, n_steps_dh, dt_dh, tc_dh)
print(f"Deep hedge : mean={pnl_dh.mean():.3f}  CVaR95={cvar95(pnl_dh):.3f}")
print(f"BSM delta  : mean={pnl_bsm.mean():.3f}  CVaR95={cvar95(pnl_bsm):.3f}")

m_range = np.linspace(-0.5, 0.5, 100)
bsm_d   = [bsm_delta(S0 * np.exp(m), K_opt, T_opt/2, sigma_opt, r_opt) for m in m_range]
with torch.no_grad():
    mt = torch.tensor(m_range, dtype=torch.float32)
    dh_d = policy(mt, torch.full_like(mt, 0.5), torch.zeros_like(mt)).numpy()

fig, axes = plt.subplots(1, 3, figsize=(14, 5))

ax = axes[0]
bins = np.linspace(-15, 10, 70)
ax.hist(pnl_bsm, bins=bins, alpha=0.6, density=True, color='steelblue',
        label=f'BSM Δ (CVaR95={cvar95(pnl_bsm):.2f})')
ax.hist(pnl_dh,  bins=bins, alpha=0.6, density=True, color='darkorange',
        label=f'Deep Hedge (CVaR95={cvar95(pnl_dh):.2f})')
ax.axvline(0, color='black', linewidth=0.5)
ax.set(xlabel='Hedging P&L ($)', ylabel='Density', title='P&L Distribution (tc = 0.1%)')
ax.legend(fontsize=9)

ax = axes[1]
ax.plot(np.arange(0, 300, 30), losses, color='steelblue', linewidth=2, marker='o', markersize=4)
ax.set(xlabel='Training epoch', ylabel='Risk-adjusted loss', title='Deep Hedging Training Curve')

ax = axes[2]
ax.plot(m_range, bsm_d, color='steelblue',  linewidth=2, linestyle='--', label='BSM delta')
ax.plot(m_range, dh_d,  color='darkorange', linewidth=2,
        label='Learned policy (t = T/2)')
ax.set(xlabel='Log-moneyness ln(S/K)', ylabel='Delta',
       title='Hedge Ratio vs Moneyness', ylim=(-0.05, 1.05))
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/hedg_deep_hedging.png', dpi=150, bbox_inches='tight')
plt.show()

torch.manual_seed(0)
np.random.seed(0)

T_ae = 2000
sigma_ae = np.array([0.006, 0.003, 0.001])
U_ae = np.column_stack([
    pc_level(maturities)  / np.linalg.norm(pc_level(maturities)),
    pc_slope(maturities)  / np.linalg.norm(pc_slope(maturities)),
    pc_curve(maturities)  / np.linalg.norm(pc_curve(maturities)),
])
F_ae   = np.random.randn(T_ae, 3) * sigma_ae
E_ae   = np.random.randn(T_ae, D) * 0.0004
curves = F_ae @ U_ae.T + E_ae

X_ae   = torch.tensor(curves, dtype=torch.float32)
X_mean = X_ae.mean(0)
X_std  = X_ae.std(0)
X_norm = (X_ae - X_mean) / (X_std + 1e-8)

class Autoencoder(nn.Module):
    def __init__(self, D, M, hidden=32):
        super().__init__()
        self.encoder = nn.Sequential(nn.Linear(D, hidden), nn.Tanh(), nn.Linear(hidden, M))
        self.decoder = nn.Sequential(nn.Linear(M, hidden), nn.Tanh(), nn.Linear(hidden, D))
    def forward(self, x):
        z = self.encoder(x)
        return self.decoder(z), z

ae     = Autoencoder(D, 3, hidden=32)
opt_ae = optim.Adam(ae.parameters(), lr=1e-3)
for epoch in range(600):
    idx  = torch.randperm(T_ae)[:256]
    xhat, _ = ae(X_norm[idx])
    loss_ae = ((X_norm[idx] - xhat)**2).mean()
    opt_ae.zero_grad()
    loss_ae.backward()
    opt_ae.step()

ae.eval()
with torch.no_grad():
    x_sample  = X_norm[1500:1510]
    xhat_s, _ = ae(x_sample)
    _, z_all  = ae(X_norm)

x_np    = (x_sample.numpy() * (X_std.numpy() + 1e-8) + X_mean.numpy())
xhat_np = (xhat_s.numpy()   * (X_std.numpy() + 1e-8) + X_mean.numpy())
z_np    = z_all.numpy()

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
for i in range(3):
    ax.plot(maturities, x_np[i] * 100,    color='steelblue',  linewidth=1.5, alpha=0.8,
            label='True' if i == 0 else None)
    ax.plot(maturities, xhat_np[i] * 100, color='darkorange', linewidth=1.5, alpha=0.8,
            linestyle='--', label='Reconstructed' if i == 0 else None)
ax.set(xlabel='Maturity (years)', ylabel='Yield change (×100)',
       title='Autoencoder Reconstruction (M=3 latent dims)')
ax.legend(fontsize=9)

ax = axes[1]
t_ax = np.arange(T_ae) / 252
for k, (lbl, col) in enumerate(zip(
        ['Latent z₁ (level)', 'Latent z₂ (slope)', 'Latent z₃ (curvature)'],
        ['steelblue', 'darkorange', 'seagreen'])):
    ax.plot(t_ax, z_np[:, k], color=col, linewidth=0.8, alpha=0.7, label=lbl)
ax.set(xlabel='Time (years)', ylabel='Latent code value', title='Latent Factor Time Series')
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/ddm_autoencoder.png', dpi=150, bbox_inches='tight')
plt.show()

torch.manual_seed(3)
np.random.seed(3)

K_mn, sig_mn, r_mn = 100.0, 0.2, 0.0
S_vals = np.linspace(60, 140, 30)
T_vals = np.array([0.1, 0.25, 0.5, 1.0, 1.5, 2.0])

rows = [[s / K_mn, t, bsm_price(s, K_mn, t, sig_mn, r_mn) / K_mn]
        for s in S_vals for t in T_vals]
rows = np.array(rows)
X_mn = torch.tensor(rows[:, :2], dtype=torch.float32)
y_mn = torch.tensor(rows[:,  2], dtype=torch.float32)

class StandardMLP(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(nn.Linear(2, 32), nn.Tanh(),
                                 nn.Linear(32, 32), nn.Tanh(),
                                 nn.Linear(32, 1), nn.ReLU())
    def forward(self, x):
        return self.net(x).squeeze(-1)

class MonotoneMLP(nn.Module):
    def __init__(self, hidden=32):
        super().__init__()
        self.W1c = nn.Parameter(torch.randn(hidden, 1) * 0.3)
        self.W1f = nn.Parameter(torch.randn(hidden, 1) * 0.3)
        self.b1  = nn.Parameter(torch.zeros(hidden))
        self.W2  = nn.Parameter(torch.randn(hidden, hidden) * 0.1)
        self.b2  = nn.Parameter(torch.zeros(hidden))
        self.W3  = nn.Parameter(torch.randn(1, hidden) * 0.1)
        self.b3  = nn.Parameter(torch.zeros(1))
    def forward(self, x):
        S = x[:, :1]
        T = x[:, 1:]
        h = torch.tanh(S @ torch.exp(self.W1c).T + T @ self.W1f.T + self.b1)
        h = torch.tanh(h @ torch.exp(self.W2).T  + self.b2)
        return torch.relu(h @ self.W3.T + self.b3).squeeze(-1)

for model, name in [(StandardMLP(), 'standard'), (MonotoneMLP(), 'monotone')]:
    opt = optim.Adam(model.parameters(), lr=5e-3)
    for ep in range(2000):
        loss = ((model(X_mn) - y_mn)**2).mean()
        opt.zero_grad(); loss.backward(); opt.step()
    if name == 'standard': mlp_std  = model
    else:                   mlp_mono = model
    print(f"{name} MSE = {loss.item():.6f}")

S_dense = np.linspace(50, 160, 200)
T_fixed = 0.5
X_dense = torch.tensor(np.column_stack([S_dense / K_mn, np.full(200, T_fixed)]),
                        dtype=torch.float32)

y_bsm   = np.array([bsm_price(s, K_mn, T_fixed, sig_mn, r_mn) / K_mn for s in S_dense])
d_bsm   = np.array([bsm_delta(s, K_mn, T_fixed, sig_mn, r_mn)       for s in S_dense])

def get_delta(model, X):
    Xg = X.clone().detach().requires_grad_(True)
    model(Xg).sum().backward()
    return Xg.grad[:, 0].detach().numpy() / (1 / K_mn)

with torch.no_grad():
    y_std  = mlp_std(X_dense).numpy()
    y_mono = mlp_mono(X_dense).numpy()
d_std  = get_delta(mlp_std,  X_dense)
d_mono = get_delta(mlp_mono, X_dense)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
ax.plot(S_dense, y_bsm  * K_mn, color='black',      linewidth=2,   label='BSM (ground truth)')
ax.plot(S_dense, y_std  * K_mn, color='steelblue',  linewidth=1.5, linestyle='--', label='Standard MLP')
ax.plot(S_dense, y_mono * K_mn, color='darkorange', linewidth=1.5, linestyle=':',  label='Monotone MLP')
ax.set(xlabel='Spot price $S$', ylabel='Call price $C$ ($)', xlim=(50, 160),
       title='Call Price Approximation ($T=0.5$)')
ax.legend(fontsize=9)

ax = axes[1]
ax.plot(S_dense, d_bsm,  color='black',      linewidth=2,   label='BSM delta')
ax.plot(S_dense, d_std,  color='steelblue',  linewidth=1.5, linestyle='--', label='Standard MLP delta')
ax.plot(S_dense, d_mono, color='darkorange', linewidth=1.5, linestyle=':',  label='Monotone MLP delta')
ax.axhline(0, color='grey', linewidth=0.5)
ax.axhline(1, color='grey', linewidth=0.5)
ax.fill_between(S_dense, 0, 1, alpha=0.05, color='green', label='Valid range [0, 1]')
ax.set(xlabel='Spot price $S$', ylabel='Delta $\\partial C / \\partial S$',
       title='Delta Comparison (Standard vs Monotone MLP)',
       xlim=(50, 160), ylim=(-0.2, 1.3))
ax.legend(fontsize=9)

plt.tight_layout()
plt.savefig(f'{FIGDIR}/ddm_monotone_nn.png', dpi=150, bbox_inches='tight')
plt.show()