From d484f9c2364a358a695f334c8f01f681cf562ad9 Mon Sep 17 00:00:00 2001 From: local Date: Sun, 8 Feb 2026 13:40:00 +0000 Subject: [PATCH] working AMC algorithm tested against quantlib --- python/american-mc/.vscode/settings.json | 5 + python/american-mc/CONTEXTRESUME.md | 19 +++ python/american-mc/main.py | 145 +++++++++++++++++++++++ python/american-mc/main2.py | 133 +++++++++++++++++++++ 4 files changed, 302 insertions(+) create mode 100644 python/american-mc/.vscode/settings.json create mode 100644 python/american-mc/CONTEXTRESUME.md create mode 100644 python/american-mc/main.py create mode 100644 python/american-mc/main2.py diff --git a/python/american-mc/.vscode/settings.json b/python/american-mc/.vscode/settings.json new file mode 100644 index 0000000..a8c2003 --- /dev/null +++ b/python/american-mc/.vscode/settings.json @@ -0,0 +1,5 @@ +{ + "python-envs.defaultEnvManager": "ms-python.python:conda", + "python-envs.defaultPackageManager": "ms-python.python:conda", + "python-envs.pythonProjects": [] +} \ No newline at end of file diff --git a/python/american-mc/CONTEXTRESUME.md b/python/american-mc/CONTEXTRESUME.md new file mode 100644 index 0000000..681aae6 --- /dev/null +++ b/python/american-mc/CONTEXTRESUME.md @@ -0,0 +1,19 @@ +Current Progress Summary: + + The Baseline: We fixed a standard Longstaff-Schwartz (LSM) American Put pricer, correcting the cash-flow propagation logic and regression targets. + + The Evolution: We moved to Bermudan Swaptions using the Hull-White One-Factor Model. + + The "Gold Standard": We implemented a 100% exact simulation from scratch. Instead of Euler discretization, we used Bivariate Normal sampling to jointly simulate the short rate rt​ and the stochastic integral ∫rs​ds. This accounts for the stochastic discount factor (the convexity adjustment) without approximation error. + + The Current Frontier: We were debating the Risk-Neutral Measure (Q) vs. the Terminal Forward Measure (QT). + + We concluded that while QT simplifies European options, it makes Bermudan LSM "messy" because it introduces a time-dependent drift shift: DriftQT​=DriftQ​−aσ2​(1−e−a(T−t)). + +Pending Topics: + + Mathematical Proof: The derivation of the "Drift Shift" via Girsanov’s Theorem. + + Exercise Boundary Impact: How the choice of measure (and the resulting drift) visually shifts the optimal exercise boundary in simulation. + + Beyond One-Factor: Potential move toward Two-Factor models or non-flat initial term structures. diff --git a/python/american-mc/main.py b/python/american-mc/main.py new file mode 100644 index 0000000..033a863 --- /dev/null +++ b/python/american-mc/main.py @@ -0,0 +1,145 @@ + +import numpy as np +import matplotlib.pyplot as plt +from scipy.stats import norm + +# ----------------------------- +# Black-Scholes European Put +# ----------------------------- +def european_put_black_scholes(S0, K, r, sigma, T): + d1 = (np.log(S0/K) + (r + 0.5*sigma**2)*T) / (sigma*np.sqrt(T)) + d2 = d1 - sigma*np.sqrt(T) + return K*np.exp(-r*T)*norm.cdf(-d2) - S0*norm.cdf(-d1) + +# ----------------------------- +# GBM Simulation +# ----------------------------- +def simulate_gbm(S0, r, sigma, T, M, N, seed=None): + if seed is not None: + np.random.seed(seed) + dt = T / N + # Vectorized simulation for speed + S = np.zeros((M, N + 1)) + S[:, 0] = S0 + for t in range(1, N + 1): + z = np.random.standard_normal(M) + S[:, t] = S[:, t-1] * np.exp((r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * z) + return S + +# ----------------------------- +# Longstaff-Schwartz Monte Carlo (Fixed) +# ----------------------------- +def lsm_american_put(S0, K, r, sigma, T, M=100_000, N=50, basis_deg=2, seed=42, return_boundary=False): + S = simulate_gbm(S0, r, sigma, T, M, N, seed) + dt = T / N + df = np.exp(-r * dt) + + # Immediate exercise value (payoff) at each step + payoff = np.maximum(K - S, 0) + + # V stores the value of the option at each path + # Initialize with the payoff at maturity + V = payoff[:, -1] + + boundary = [] + + # Backward induction + for t in range(N - 1, 0, -1): + # Identify In-The-Money paths + itm = payoff[:, t] > 0 + + # Default: option value at t is just the discounted value from t+1 + V = V * df + + if np.any(itm): + X = S[itm, t] + Y = V[itm] # These are already discounted from t+1 + + # Regression: Estimate Continuation Value + coeffs = np.polyfit(X, Y, basis_deg) + continuation_value = np.polyval(coeffs, X) + + # Exercise if payoff > estimated continuation value + exercise = payoff[itm, t] > continuation_value + + # Update V for paths where we exercise + # Get the indices of the ITM paths that should exercise + itm_indices = np.where(itm)[0] + exercise_indices = itm_indices[exercise] + V[exercise_indices] = payoff[exercise_indices, t] + + # Boundary: The highest stock price at which we still exercise + if len(exercise_indices) > 0: + boundary.append((t * dt, np.max(S[exercise_indices, t]))) + else: + boundary.append((t * dt, np.nan)) + else: + boundary.append((t * dt, np.nan)) + + # Final price is the average of discounted values at t=1 + price = np.mean(V * df) + + if return_boundary: + return price, boundary[::-1] + return price + +# ----------------------------- +# QuantLib American Put (v1.18) +# ----------------------------- + +def quantlib_american_put(S0, K, r, sigma, T): + import QuantLib as ql + + today = ql.Date().todaysDate() + ql.Settings.instance().evaluationDate = today + maturity = today + int(T*365) + payoff = ql.PlainVanillaPayoff(ql.Option.Put, K) + exercise = ql.AmericanExercise(today, maturity) + + spot = ql.QuoteHandle(ql.SimpleQuote(S0)) + dividend_curve = ql.YieldTermStructureHandle(ql.FlatForward(today, 0.0, ql.Actual365Fixed())) + riskfree_curve = ql.YieldTermStructureHandle(ql.FlatForward(today, r, ql.Actual365Fixed())) + vol_curve = ql.BlackVolTermStructureHandle(ql.BlackConstantVol(today, ql.NullCalendar(), sigma, ql.Actual365Fixed())) + + process = ql.BlackScholesMertonProcess(spot, dividend_curve, riskfree_curve, vol_curve) + engine = ql.FdBlackScholesVanillaEngine(process, 200, 400) + option = ql.VanillaOption(payoff, exercise) + option.setPricingEngine(engine) + return option.NPV() + +# ----------------------------- +# Main script +# ----------------------------- +if __name__ == "__main__": + # Parameters + S0, K, r, sigma, T = 100, 100, 0.05, 0.2, 1.0 + M, N = 200_000, 50 + + # European Put + eur_bs = european_put_black_scholes(S0, K, r, sigma, T) + print(f"European Put (BS): {eur_bs:.4f}") + + # LSM American Put + lsm_price, boundary = lsm_american_put(S0, K, r, sigma, T, M, N, return_boundary=True) + print(f"LSM American Put (M={M}): {lsm_price:.4f}") + + # QuantLib American Put + ql_price = quantlib_american_put(S0, K, r, sigma, T) + print(f"QuantLib American Put: {ql_price:.4f}") + + print(f"Lower bound (immediate exercise): {max(K-S0,0):.4f}") + + # Plot Exercise Boundary + times = [b[0] for b in boundary] + boundaries = [b[1] for b in boundary] + + plt.figure(figsize=(8,5)) + plt.plot(times, boundaries, color='orange', label="LSM Exercise Boundary") + plt.axhline(K, color='red', linestyle='--', label="Strike Price") + plt.xlabel("Time to Maturity") + plt.ylabel("Stock Price for Exercise") + plt.title("American Put LSM Exercise Boundary") + plt.gca().invert_xaxis() + plt.legend() + plt.grid(True) + plt.show() diff --git a/python/american-mc/main2.py b/python/american-mc/main2.py new file mode 100644 index 0000000..9e49845 --- /dev/null +++ b/python/american-mc/main2.py @@ -0,0 +1,133 @@ +import numpy as np +import matplotlib.pyplot as plt + +# --------------------------------------------------------- +# 1. Analytical Hull-White Bond Pricing +# --------------------------------------------------------- +def bond_price(r_t, t, T, a, sigma, r0): + """Zero-coupon bond price P(t, T) in Hull-White model for flat curve r0.""" + B = (1 - np.exp(-a * (T - t))) / a + A = np.exp((B - (T - t)) * (a**2 * r0 - sigma**2/2) / a**2 - (sigma**2 * B**2 / (4 * a))) + return A * np.exp(-B * r_t) + +def get_swap_npv(r_t, t, T_end, strike, a, sigma, r0): + """NPV of a Payer Swap: Receives Floating, Pays Fixed strike.""" + payment_dates = np.arange(t + 1, T_end + 1, 1.0) + if len(payment_dates) == 0: + return np.zeros_like(r_t) + + # NPV = 1 - P(t, T_end) - strike * Sum[P(t, T_i)] + p_end = bond_price(r_t, t, T_end, a, sigma, r0) + fixed_leg = sum(bond_price(r_t, t, pd, a, sigma, r0) for pd in payment_dates) + + return np.maximum(1 - p_end - strike * fixed_leg, 0) + +# --------------------------------------------------------- +# 2. Joint Exact Simulator (Short Rate & Stochastic Integral) +# --------------------------------------------------------- +def simulate_hw_exact_joint(r0, a, sigma, exercise_dates, M): + """Samples (r_t, Integral[r]) jointly to get exact discount factors.""" + M = int(M) + num_dates = len(exercise_dates) + r_matrix = np.zeros((M, num_dates + 1)) + d_matrix = np.zeros((M, num_dates)) # d[i] = exp(-integral from t_i to t_{i+1}) + + r_matrix[:, 0] = r0 + t_steps = np.insert(exercise_dates, 0, 0.0) + + for i in range(len(t_steps) - 1): + t, T = t_steps[i], t_steps[i+1] + dt = T - t + + # Drift adjustments for flat initial curve r0 + alpha_t = r0 + (sigma**2 / (2 * a**2)) * (1 - np.exp(-a * t))**2 + alpha_T = r0 + (sigma**2 / (2 * a**2)) * (1 - np.exp(-a * T))**2 + + # Means + mean_r = r_matrix[:, i] * np.exp(-a * dt) + alpha_T - alpha_t * np.exp(-a * dt) + # The expected value of the integral is derived from the bond price: E[exp(-I)] = P(t,T) + mean_I = -np.log(bond_price(r_matrix[:, i], t, T, a, sigma, r0)) + + # Covariance Matrix Components + var_r = (sigma**2 / (2 * a)) * (1 - np.exp(-2 * a * dt)) + B = (1 - np.exp(-a * dt)) / a + var_I = (sigma**2 / a**2) * (dt - 2*B + (1 - np.exp(-2*a*dt))/(2*a)) + cov_rI = (sigma**2 / (2 * a**2)) * (1 - np.exp(-a * dt))**2 + + cov_matrix = [[var_r, cov_rI], [cov_rI, var_I]] + + # Sample joint normal innovations + Z = np.random.multivariate_normal([0, 0], cov_matrix, M) + + r_matrix[:, i+1] = mean_r + Z[:, 0] + # Important: The variance of the integral affects the mean of the exponent (convexity) + # mean_I here is already the log of the bond price (risk-neutral expectation) + d_matrix[:, i] = np.exp(-(mean_I + Z[:, 1] - 0.5 * var_I)) + + return r_matrix[:, 1:], d_matrix + +# --------------------------------------------------------- +# 3. LSM Pricing Logic +# --------------------------------------------------------- +def price_bermudan_swaption(r0, a, sigma, strike, exercise_dates, T_end, M): + # 1. Generate exact paths + r_at_ex, discounts = simulate_hw_exact_joint(r0, a, sigma, exercise_dates, M) + + # 2. Final exercise date payoff + T_last = exercise_dates[-1] + cash_flows = get_swap_npv(r_at_ex[:, -1], T_last, T_end, strike, a, sigma, r0) + + # 3. Backward Induction + # exercise_dates[:-1] because we already handled the last date + for i in reversed(range(len(exercise_dates) - 1)): + t_current = exercise_dates[i] + + # Pull cash flows back to current time using stochastic discount + # Note: If path was exercised later, this is the discounted value of that exercise. + cash_flows = cash_flows * discounts[:, i] + + # Current intrinsic value + X = r_at_ex[:, i] + immediate_payoff = get_swap_npv(X, t_current, T_end, strike, a, sigma, r0) + + # Only regress In-The-Money paths + itm = immediate_payoff > 0 + if np.any(itm): + # Regression: Basis functions [1, r, r^2] + A = np.column_stack([np.ones_like(X[itm]), X[itm], X[itm]**2]) + coeffs = np.linalg.lstsq(A, cash_flows[itm], rcond=None)[0] + continuation_value = A @ coeffs + + # Exercise decision + exercise = immediate_payoff[itm] > continuation_value + + itm_indices = np.where(itm)[0] + exercise_indices = itm_indices[exercise] + + # Update cash flows for paths where we exercise early + cash_flows[exercise_indices] = immediate_payoff[exercise_indices] + + # Final discount to t=0 + # The first discount factor in 'discounts' is from t1 to t0 + # But wait, our 'discounts' matrix is (M, num_dates). + # Let's just use the analytical P(0, t1) for the very last step to t=0. + final_price = np.mean(cash_flows * bond_price(r0, 0, exercise_dates[0], a, sigma, r0)) + return final_price + +# --------------------------------------------------------- +# Execution +# --------------------------------------------------------- +if __name__ == "__main__": + # Params + r0_val, a_val, sigma_val = 0.05, 0.1, 0.01 + strike_val = 0.05 + ex_dates = np.array([1.0, 2.0, 3.0, 4.0]) + maturity = 5.0 + num_paths = 100_000 + + price = price_bermudan_swaption(r0_val, a_val, sigma_val, strike_val, ex_dates, maturity, num_paths) + + print(f"--- 100% Exact LSM Bermudan Swaption ---") + print(f"Parameters: a={a_val}, sigma={sigma_val}, strike={strike_val}") + print(f"Exercise Dates: {ex_dates}") + print(f"Calculated Price: {price:.6f}") \ No newline at end of file