working AMC algorithm tested against quantlib

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": []
+}

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.

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()

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}")