back testing script for bull put spreads

scholes_back_test.py

import math
from datetime import datetime, timedelta
import numpy as np
import pandas as pd
import yfinance as yf


# --- 1. Black-Scholes Pricing Engine ---
def norm_cdf(x: float) -> float:
  """Standard normal cumulative distribution function using math.erf."""
  return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0)))


def bs_put_price(
    S: float, K: float, T: float, r: float = 0.045, sigma: float = 0.18
) -> float:
  """Calculates European Put Option price via Black-Scholes model.

  Parameters:
      S (float): Spot price
      K (float): Strike price
      T (float): Time to expiration in years (DTE / 365)
      r (float): Risk-free interest rate (default 4.5%)
      sigma (float): Volatility (annualized)
  """
  if T <= 0:
    return max(0.0, K - S)
  if S <= 0 or K <= 0 or sigma <= 0:
    return max(0.0, K - S)

  d1 = (math.log(S / K) + (r + 0.5 * sigma**2) * T) / (sigma * math.sqrt(T))
  d2 = d1 - sigma * math.sqrt(T)

  put_price = K * math.exp(-r * T) * norm_cdf(-d2) - S * norm_cdf(-d1)
  return max(0.0, put_price)


# --- 2. Advanced Path-Dependent Put Risk Scanner ---
def scan_put_risk_classified(
    ticker: str,
    lookback_years: float = 2.0,
    span_days: int = 40,
    drop_pct: float = 0.07,
    profit_target_pct: float = 0.50,
    fixed_iv: float = 0.18,
) -> pd.DataFrame:
  """Evaluates short put positions using daily Black-Scholes mark-to-market

  pricing and early profit-taking logic to categorize trades into Paths A-D.
  """
  end_date = datetime.now()
  start_date = end_date - timedelta(days=int(lookback_years * 365 + span_days))

  df = yf.download(ticker, start=start_date, end=end_date, progress=False)
  if df.empty:
    raise ValueError(f"No data found for {ticker}")

  prices = (
      df["Close"][ticker] if isinstance(df.columns, pd.MultiIndex) else df["Close"]
  )

  results = []

  for i in range(len(prices)):
    entry_date = prices.index[i]
    entry_price = float(prices.iloc[i])
    target_exp_date = entry_date + timedelta(days=span_days)

    window = prices.loc[
        (prices.index > entry_date) & (prices.index <= target_exp_date)
    ]
    if window.empty or window.index[-1] < (target_exp_date - timedelta(days=4)):
      continue

    strike_price = entry_price * (1.0 - drop_pct)
    init_T = span_days / 365.0
    entry_option_val = bs_put_price(
        entry_price, strike_price, init_T, sigma=fixed_iv
    )
    target_option_val = entry_option_val * (1.0 - profit_target_pct)

    early_exit = False
    exit_day = None
    exit_reason = None
    breached = False

    # Track path day-by-day
    for current_date, spot in window.items():
      days_remaining = (target_exp_date - current_date).days
      T_remaining = max(0.0, days_remaining / 365.0)
      curr_spot = float(spot)

      if curr_spot <= strike_price:
        breached = True

      curr_option_val = bs_put_price(
          curr_spot, strike_price, T_remaining, sigma=fixed_iv
      )

      # Check early profit exit condition
      if curr_option_val <= target_option_val and not early_exit:
        early_exit = True
        exit_day = current_date
        exit_reason = "Profit Target Hit"
        break  # Exit position on profit target hit

    exp_price = float(window.iloc[-1])
    expired_itm = exp_price < strike_price

    # Path Classification Logic
    if early_exit and not breached:
      path = "Path A (Clean Win)"
    elif early_exit and breached:
      path = "Path B (Saved Win)"
    elif not early_exit and breached and not expired_itm:
      path = "Path C (Whipsaw/Stress)"
    elif not early_exit and breached and expired_itm:
      path = "Path D (Toxic Loss)"
    else:
      path = "Path A (Clean Win)"  # Expired OTM cleanly without early exit trigger

    results.append({
        "entry_date": entry_date.strftime("%Y-%m-%d"),
        "entry_price": round(entry_price, 2),
        "strike": round(strike_price, 2),
        "exp_price": round(exp_price, 2),
        "breached": breached,
        "early_exit": early_exit,
        "path": path,
    })

  res_df = pd.DataFrame(results)

  # Aggregation
  total = len(res_df)
  path_counts = res_df["path"].value_counts()

  print(f"\n================ PATH-DEPENDENT RISK SCAN: {ticker} ================")
  print(
      f"Params: {lookback_years}Y Lookback | {span_days} DTE | -{drop_pct*100:.1f}%"
      f" OTM | {profit_target_pct*100:.0f}% Profit Target"
  )
  print(f"Total Rolling Windows Evaluated: {total}\n")
  print(
      f"  Path A (Clean Win)       : {path_counts.get('Path A (Clean Win)', 0):>3d} "
      f"({path_counts.get('Path A (Clean Win)', 0)/total*100:5.1f}%)"
  )
  print(
      f"  Path B (Saved Win)       : {path_counts.get('Path B (Saved Win)', 0):>3d} "
      f"({path_counts.get('Path B (Saved Win)', 0)/total*100:5.1f}%)"
  )
  print(
      f"  Path C (Whipsaw/Stress)  : {path_counts.get('Path C (Whipsaw/Stress)', 0):>3d} "
      f"({path_counts.get('Path C (Whipsaw/Stress)', 0)/total*100:5.1f}%)"
  )
  print(
      f"  Path D (Toxic Loss)      : {path_counts.get('Path D (Toxic Loss)', 0):>3d} "
      f"({path_counts.get('Path D (Toxic Loss)', 0)/total*100:5.1f}%)"
  )
  print("===================================================================\n")

  return res_df


if __name__ == "__main__":
  df_res = scan_put_risk_classified(
      ticker="SPY",
      lookback_years=4.0,
      span_days=30,
      drop_pct=0.04,
      profit_target_pct=0.45,
      fixed_iv=0.19,
  )