def merton_jump_paths(S0, r, sigma, lmbda, mu_J, sigma_J, T, n_sims, n_steps):
    dt = T / n_steps
    kappa = np.exp(mu_J + 0.5 * sigma_J**2) - 1
    mu_adj = r - lmbda * kappa  # risk-neutral drift

    noise = np.random.normal(0, 1, size=(n_sims, n_steps))
    jumps = np.random.poisson(lmbda * dt, size=(n_sims, n_steps))
    jump_sizes = np.random.normal(mu_J, sigma_J, size=(n_sims, n_steps))
    J = np.exp(jump_sizes)

    log_returns = (mu_adj - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * noise
    log_returns += jumps * np.log(J)
    exponent = np.cumsum(log_returns, axis=1)
    paths = S0 * np.exp(exponent)
    return np.column_stack((np.full(n_sims, S0), paths))

gbm_paths    = GBM_paths(S0, sigma, T, r, mu, n_sims, n_steps)
merton_paths = merton_jump_paths(S0, r, sigma, lmbda, mu_J, sigma_J, T, n_sims, n_steps)

t = np.linspace(0, T, n_steps + 1)

plt.figure(figsize=(10, 6))
paths_to_plot = 10

# GBM paths (blue shades)
for i in range(paths_to_plot):
    plt.plot(t, gbm_paths[i], color=plt.cm.Blues(0.3 + 0.6 * i / paths_to_plot), lw=1)

# Merton jump paths (red shades)
for i in range(paths_to_plot):
    plt.plot(t, merton_paths[i], color=plt.cm.Reds(0.3 + 0.6 * i / paths_to_plot), lw=1)

#plt.plot(t, gbm_paths.mean(axis=0), 'b', ls='--', alpha=0.5, label='GBM mean')
#plt.plot(t, merton_paths.mean(axis=0),'r', ls='--', alpha=0.5, label='Merton mean')
#plt.legend()

plt.title("GBM (blue) vs Merton JumpDiffusion (red) Paths")
plt.xlabel("Time (years)")
plt.ylabel("Price")
plt.show()

def mc_option_price_from_paths(paths, K, r, T, option_type='call'):
    S_T = paths[:, -1]
    if option_type == 'call':
        payoffs = np.maximum(S_T - K, 0.0)
    elif option_type == 'put':
        payoffs = np.maximum(K - S_T, 0.0)
    else:
        raise ValueError('option_type must be "call" or "put"')
    discounted = np.exp(-r * T) * payoffs
    price = np.mean(discounted)
    stderr = np.std(discounted, ddof=0) / np.sqrt(len(discounted))
    return price, stderr

mc_call, stderr_call = mc_option_price_from_paths(merton_paths, K, r, T, option_type='call')
mc_put, stderr_put 	 = mc_option_price_from_paths(merton_paths, K, r, T, option_type='put')
bs_call_price 		 = bs_call(S0, K, sigma, T, r)
bs_put_price  		 = bs_put(S0, K, sigma, T, r)

print(f"Black-Scholes call price: \t${bs_call_price:.2f}")
print(f"Merton call price: \t\t${bs_call_price:.2f} ± ${stderr_call:.2f}")

print(f"Black-Scholes put price: \t${bs_put_price:.2f}")
print(f"Merton put price: \t\t${mc_put:.2f} ± ${stderr_put:.2f}")

data_calib = yf.download(ticker, start=calib_start, end=calib_end, progress=False)
prices_calib = data_calib['Close'].dropna()
prices_calib = prices_calib.iloc[:, 0] # This must be changed for multiple tickers...

log_returns_daily = np.log(prices_calib / prices_calib.shift(1)).dropna().values
n = len(log_returns_daily)
dt = 1.0 / TRADING_DAYS

# There is a closed-form MLE for GBM:
# Model: x_i = ln(S_{t+dt}/S_t) ~ N((mu - 0.5*sigma^2)*dt, sigma^2 * dt)
# Let m = mean(x_i), s2 = var(x_i) (MLE uses ddof=0)
m  = np.mean(log_returns_daily)
s2 = np.var(log_returns_daily, ddof=0)

# annualized volatility estimate
sigma_hat = np.sqrt(s2 / dt)

# annualized drift estimate
mu_hat = m / dt + 0.5 * sigma_hat**2

print("Closed-form MLE estimates (annualized):")
print(f"  mu_hat (drift) = {mu_hat:.4f}")
print(f"  sigma_hat      = {sigma_hat:.4f}")
print(f"  sample size n  = {n}")

# Negative log-likelihood (for numeric verification)
# Use parameterization theta = [mu, sigma], sigma > 0
def neg_loglike_theta(theta, returns, dt):
    mu, sigma = theta
    if sigma <= 0:
        return 1e10
    # mean and var of returns under model
    mean_r = (mu - 0.5 * sigma**2) * dt
    var_r = (sigma**2) * dt
    # Gaussian log-likelihood
    ll = -0.5 * np.sum(((returns - mean_r)**2) / var_r + np.log(2 * np.pi * var_r))
    return -ll  # negative log-likelihood

# numeric optimization starting from closed form
init = np.array([mu_hat, sigma_hat])
bounds = [(None, None), (1e-8, None)]  # sigma > 0

# This particular optimization algorithm was suggested in similar problems I found online
res = minimize(neg_loglike_theta, init, args=(log_returns_daily, dt), bounds=bounds, method='L-BFGS-B')
if res.success:
    mu_mle_num, sigma_mle_num = res.x
    print("\nNumeric MLE verification (annualized):")
    print(f"  mu_mle_num    = {mu_mle_num:.4f}")
    print(f"  sigma_mle_num = {sigma_mle_num:.4f}")
else:
    mu_mle_num, sigma_mle_num = mu_hat, sigma_hat
    print("\nNumeric optimization failed; falling back to closed-form estimates.")
    print(res.message)

# prefer closed-form
mu_gbm_mle = mu_hat
sigma_gbm_mle = sigma_hat

# plot histogram of daily log-returns with fitted normal pdf overlay
plt.figure(figsize=(10, 5))
counts, bins, patches = plt.hist(log_returns_daily, bins=60, density=True, alpha=0.6, label='Empirical daily log-returns')

# fitted normal using daily sample mean/var
from scipy.stats import norm
mu_daily_hat = (mu_gbm_mle - 0.5 * sigma_gbm_mle**2) * dt
sigma_daily_hat = sigma_gbm_mle * np.sqrt(dt)

x = np.linspace(bins[0], bins[-1], 200)
pdf = norm.pdf(x, loc=mu_daily_hat, scale=sigma_daily_hat)
plt.plot(x, pdf, lw=2, ls='--', label='Fitted Normal (GBM MLE)')

plt.title(f"{ticker} Daily Log-Returns Histogram ({calib_start} to {calib_end})")
plt.xlabel("Daily log-return")
plt.ylabel("Density")
plt.legend()
plt.show()

# this seems to fix some errors
mu_gbm_mle    = float(mu_gbm_mle)
sigma_gbm_mle = float(sigma_gbm_mle)

def merton_log_likelihood(params, returns, dt, n_max=5):
    """
    Compute the negative log-likelihood for Merton jump-diffusion daily returns.

    params = [mu, sigma, lambda, mu_J, sigma_J]
    """
    mu, sigma, lmbda, mu_J, sigma_J = params
    if sigma <= 0 or lmbda < 0 or sigma_J <= 0:
        return np.inf

    # Precompute Poisson probabilities up to n_max jumps
    pois_probs = np.array([
        np.exp(-lmbda * dt) * (lmbda * dt)**n / factorial(n)
        for n in range(n_max + 1)
    ])

    # For each return evaluate the truncated mixture PDF
    total_pdf = np.zeros_like(returns)
    for n in range(n_max + 1):
        mean_n = (mu - 0.5 * sigma**2) * dt + n * mu_J
        var_n = sigma**2 * dt + n * sigma_J**2
        total_pdf += pois_probs[n] * norm.pdf(returns, loc=mean_n, scale=np.sqrt(var_n))

    # Avoid log(0) issues?
    total_pdf = np.maximum(total_pdf, 1e-300)
    log_likelihood = np.sum(np.log(total_pdf))
    return -log_likelihood  # minimize negative log-likelihood

def calibrate_merton_params(log_returns_daily, dt):
    """
    Calibrate Merton jump-diffusion parameters by MLE.
    Returns a dict with the estimated parameters.
    """
    # Initial guesses based on GBM calibration
    m = np.mean(log_returns_daily)
    s = np.std(log_returns_daily, ddof=0)
    sigma0 = s / np.sqrt(dt)
    mu0 = m / dt + 0.5 * sigma0**2
    lmbda0, mu_J0, sigma_J0 = 0.3, -0.1, 0.2
    x0 = np.array([mu0, sigma0, lmbda0, mu_J0, sigma_J0])

    bounds = [
      (None, None),      # mu
      (1e-6, None),      # sigma > 0
      (0.0, 2.0),        # lambda >= 0, upper bound helps stability
      (-1.0, 1.0),       # mu_J (jump mean in log-space)
      (1e-6, 2.0)        # sigma_J > 0
    ]

    res = minimize(
        merton_log_likelihood, x0,
        args=(log_returns_daily, dt),
        method='L-BFGS-B',
        bounds=bounds,
        options={'maxiter': 300}
    )

    if not res.success:
        print("Merton MLE did not fully converge:", res.message)

    mu_mle, sigma_mle, lmbda_mle, mu_J_mle, sigma_J_mle = res.x
    print("\nMerton Jump-Diffusion MLE estimates (annualized):")
    print(f"  mu_mle      = {mu_mle:.4f}")
    print(f"  sigma_mle   = {sigma_mle:.4f}")
    print(f"  lambda_mle  = {lmbda_mle:.4f}")
    print(f"  mu_J_mle    = {mu_J_mle:.4f}")
    print(f"  sigma_J_mle = {sigma_J_mle:.4f}")

    return {
        'mu': mu_mle,
        'sigma': sigma_mle,
        'lambda': lmbda_mle,
        'mu_J': mu_J_mle,
        'sigma_J': sigma_J_mle
    }

# Run calibration on historical daily log-returns
merton_params = calibrate_merton_params(log_returns_daily, dt)

# Save individual values for later use, not really going to need this since it's all below anyway...
mu_merton_mle     = merton_params['mu']
sigma_merton_mle  = merton_params['sigma']
lambda_merton_mle = merton_params['lambda']
mu_J_mle          = merton_params['mu_J']
sigma_J_mle       = merton_params['sigma_J']

x_grid = np.linspace(np.percentile(log_returns_daily, 1), np.percentile(log_returns_daily, 99), 400)

# GBM fitted PDF
mu_daily_gbm = (mu_gbm_mle - 0.5 * sigma_gbm_mle**2) * dt
sigma_daily_gbm = sigma_gbm_mle * np.sqrt(dt)
pdf_gbm = norm.pdf(x_grid, loc=mu_daily_gbm, scale=sigma_daily_gbm)

# Merton fitted PDF
def merton_pdf(x, mu, sigma, lmbda, mu_J, sigma_J, dt, n_max=5):
	pois_probs = np.array([
		np.exp(-lmbda * dt) * (lmbda * dt)**n / factorial(n)
		for n in range(n_max + 1)
	])
	pdf = np.zeros_like(x)
	for n in range(n_max + 1):
		mean_n = (mu - 0.5 * sigma**2) * dt + n * mu_J
		var_n = sigma**2 * dt + n * sigma_J**2
		pdf += pois_probs[n] * norm.pdf(x, loc=mean_n, scale=np.sqrt(var_n))
	return pdf

pdf_merton = merton_pdf(
	x_grid,
	mu_merton_mle,
	sigma_merton_mle,
	lambda_merton_mle,
	mu_J_mle,
	sigma_J_mle,
	dt
)

# Plotting
plt.figure(figsize=(10, 5))
plt.hist(log_returns_daily, bins=60, density=True, alpha=0.6, label="Empirical daily log-returns")
plt.plot(x_grid, pdf_gbm, 'b--', lw=2, label="GBM (Normal fit)")
plt.plot(x_grid, pdf_merton, 'r-', lw=2, label="Merton (Mixture fit)")
plt.title(f"{ticker}: Empirical vs GBM vs Merton fitted PDFs")
plt.xlabel("Daily log-return")
plt.ylabel("Density")
plt.legend()
plt.show()

# Print fitted parameters again for quick reference
print("GBM MLE:")
print(f"  mu_gbm_mle        = {mu_gbm_mle:.4f}")
print(f"  sigma_gbm_mle     = {sigma_gbm_mle:.4f}")
print("\nMerton MLE:")
print(f"  mu_merton_mle     = {mu_merton_mle:.4f}")
print(f"  sigma_merton_mle  = {sigma_merton_mle:.4f}")
print(f"  lambda_merton_mle = {lambda_merton_mle:.4f}")
print(f"  mu_J_mle          = {mu_J_mle:.4f}")
print(f"  sigma_J_mle       = {sigma_J_mle:.4f}")

def merton_risk_neutral_drift(r, lmbda, mu_J, sigma_J):
    """Continuous diffusion drift under risk-neutral measure for Merton model."""
    EY = np.exp(mu_J + 0.5 * sigma_J**2)  # E[e^{J}]
    return r - lmbda * (EY - 1.0)

# Make a single set of unit-starting paths and re-use
rn_drift = merton_risk_neutral_drift(r, lambda_merton_mle, mu_J_mle, sigma_J_mle)

# simulate once with S0 = 1 so I can scale for different S0 cheaply
_unit_paths_key = ("merton_unit_paths", n_sims, n_steps, T, rn_drift, sigma_merton_mle,
                   lambda_merton_mle, mu_J_mle, sigma_J_mle)

unit_paths = merton_jump_paths(1.0, rn_drift, sigma_merton_mle,
                                   lambda_merton_mle, mu_J_mle, sigma_J_mle,
                                   T, n_sims, n_steps)

# price cache keyed by S0 value, makes things faster
_price_cache = {}

def price_merton_mc(S0_query, option_type='call', force_recompute=False):
    """Return MC price (and stderr) for given S0_query using cached unit paths."""
    key = (float(S0_query), option_type)
    if (not force_recompute) and (key in _price_cache):
        return _price_cache[key]
    # scale unit paths
    paths = unit_paths * float(S0_query)
    price, stderr = mc_option_price_from_paths(paths, K, r, T, option_type=option_type)
    _price_cache[key] = (price, stderr)
    return price, stderr

# finite-difference delta routine
def merton_delta_fd(S0_eval, eps=1e-3, option_type='call'):
    """Central finite-difference delta under Merton using risk-neutral MC pricer.
       eps is relative bump (default 0.1%).
    """
    S_up = S0_eval * (1.0 + eps)
    S_dn = S0_eval * (1.0 - eps)
    # compute prices (cached)
    c_up, se_up = price_merton_mc(S_up, option_type=option_type)
    c_dn, se_dn = price_merton_mc(S_dn, option_type=option_type)
    delta = (c_up - c_dn) / (2.0 * eps * S0_eval)
    # estimate stderr of delta via propagation (approx): sqrt(se_up^2 + se_dn^2)/(2 eps S0)
    delta_stderr = np.sqrt(se_up**2 + se_dn**2) / (2.0 * eps * S0_eval)
    return delta, delta_stderr, (c_up, se_up), (c_dn, se_dn)

# evaluate deltas across a grid and compare to Black-Scholes delta (for intuition)
S_grid = np.linspace(0.8 * S0, 1.2 * S0, 9)  # 80% -> 120% around current S0
eps = 1e-3  # 0.1%

merton_deltas = []
merton_errs   = []
bs_deltas     = []

for S_val in S_grid:
    d_m, d_se, up_info, dn_info = merton_delta_fd(S_val, eps=eps, option_type='call')
    merton_deltas.append(d_m)
    merton_errs.append(d_se)
    # use BS delta with diffusion vol = sigma_merton_mle (for reference only)
    try:
        d_bs = bs_call_delta(S_val, K, sigma_merton_mle, T, r)
    except Exception:
        # fallback: if bs_call_delta signature differs (S,K,sigma,T,r) try alternative order
        d_bs = bs_call_delta(S_val, K, sigma_merton_mle, T, r)
    bs_deltas.append(d_bs)

# print results in a neat table
print("\n     S | Merton Delta (MC ± stderr)      | BS Delta (sigma_diffusion)")
print("-------+---------------------------------+---------------------------")
for S_val, d_m, d_se, d_bs in zip(S_grid, merton_deltas, merton_errs, bs_deltas):
    print(f"{S_val:6.2f} | {d_m:6.4f} ± {d_se:6.4f}                 | {d_bs:6.4f}")

# plot deltas with error bars
plt.figure(figsize=(8,4.5))
plt.errorbar(S_grid, merton_deltas, yerr=merton_errs, fmt='o-', label='Merton MC Δ (central FD)', capsize=3)
plt.plot(S_grid, bs_deltas, 's--', label=f'BS Δ (σ = {sigma_merton_mle:.3f})')
plt.title("Delta: Merton (MC FD) vs Black-Scholes (diffusion vol)")
plt.xlabel("Spot price S")
plt.ylabel("Delta (∂Price/∂S)")
plt.grid(alpha=0.3)
plt.legend()
plt.show()

# single-value print at notebook S0
d0, d0_se, up_info0, dn_info0 = merton_delta_fd(S0, eps=eps, option_type='call')
print(f"\nAt S0 = {S0:.2f} (eps={eps*100:.2f}%):")
print(f"  Merton MC delta        = {d0:.4f} ± {d0_se:.4f}")
print(f"  Up price (S0*(1+eps))  = {up_info0[0]:.4f} (stderr {up_info0[1]:.4f})")
print(f"  Dn price (S0*(1-eps))  = {dn_info0[0]:.4f} (stderr {dn_info0[1]:.4f})")
print(f"  BS delta (diffusion σ) = {bs_call_delta(S0, K, sigma_merton_mle, T, r):.4f}")

# backtest parameters
n_mc_price     = 3000   # MC sims for pricing during backtest (reduce if slow)
n_steps_price  = 252    # used only for MC pricer path granularity (daily hedging => 252 steps/year)
tc_pct         = 0.0005 # transaction cost proportion (0.05% per trade) -- change as needed
loss_threshold = -0.02  # fraction / absolute? I'll treat as absolute P&L < loss_threshold (in currency)

_merton_tau_cache = {} # I think I could have used functools cache decorator...

def _rn_drift(r, lmbda, mu_J, sigma_J):
    return r - lmbda * (np.exp(mu_J + 0.5 * sigma_J**2) - 1.0)

def price_and_delta_merton_at_tau(S0_val, tau, eps=1e-3, option_type='call'):
    """
    Return (price, price_stderr, delta, delta_stderr) for Merton model at spot S0_val and time-to-maturity tau.
    Caches unit-paths keyed by tau for speed.
    """
    # cache key based on tau and model params + MC settings
    key = ("merton_tau", float(tau), n_mc_price, n_steps_price,
           float(sigma_merton_mle), float(lambda_merton_mle),
           float(mu_J_mle), float(sigma_J_mle))
    if key not in _merton_tau_cache:
        # create unit paths for this tau under risk-neutral drift
        rn = _rn_drift(r, lambda_merton_mle, mu_J_mle, sigma_J_mle)
        unit_paths = merton_jump_paths(1.0, rn, sigma_merton_mle,
                                       lambda_merton_mle, mu_J_mle, sigma_J_mle,
                                       tau, n_mc_price, n_steps_price)
        _merton_tau_cache[key] = unit_paths
    else:
        unit_paths = _merton_tau_cache[key]

    # price
    paths = unit_paths * float(S0_val)
    price, price_stderr = mc_option_price_from_paths(paths, K, r, tau, option_type=option_type)

    # delta via central finite difference using same cached unit_paths (scale up/down)
    S_up = S0_val * (1.0 + eps)
    S_dn = S0_val * (1.0 - eps)
    # reuse unit_paths to compute up/down price
    price_up, se_up = mc_option_price_from_paths(unit_paths * S_up, K, r, tau, option_type=option_type)
    price_dn, se_dn = mc_option_price_from_paths(unit_paths * S_dn, K, r, tau, option_type=option_type)
    delta = (price_up - price_dn) / (2.0 * eps * S0_val)
    delta_stderr = np.sqrt(se_up**2 + se_dn**2) / (2.0 * eps * S0_val)

    return price, price_stderr, delta, delta_stderr

# GBM (BS) price & delta
def price_and_delta_gbm(S0_val, tau, sigma_use, option_type='call'):
    # Black-Scholes closed-form price and delta
    price = bs_call(S0_val, K, sigma_use, tau, r)
    delta = bs_call_delta(S0_val, K, sigma_use, tau, r)
    # No MC stderr for BS (set to 0)
    return price, 0.0, delta, 0.0

# Load backtest price series
prices_bt = yf.download(ticker, start=backtest_start, end=backtest_end, progress=False)
prices_bt = prices_bt['Close'].dropna()
prices_bt = prices_bt.iloc[:, 0]
prices_bt = prices_bt.copy()  # not sure if necessary but it seemed to fix a bug
dates     = prices_bt.index
n_days    = len(dates)

# time grid: for daily rebalancing I use the calendar of observed prices
# Map each date to time-to-maturity (in years)
date_to_idx = {d: i for i, d in enumerate(dates)}

# maturity time in calendar: assume option sold at first date with maturity T (in years)
start_date = dates[0]
# simply use T as the remaining time from start; for each day i, tau = T * (1 - i/(n_days-1)) approximating trading calendar to maturity
# but more robust: tau = max(T - i/TRADING_DAYS, 0)
# I'll compute tau daily as T - day_index / TRADING_DAYS
taus = np.maximum(T - np.arange(n_days) / float(TRADING_DAYS), 0.0)

# Backtester (single model)
def run_delta_hedge(prices_series, model='gbm', sigma_for_gbm=None, tc=tc_pct, n_mc_price_local=n_mc_price):
    """
    Run discrete daily delta hedging for selling 1 European call at t0 under chosen model.
    model: 'gbm' or 'merton'
    Returns: dict with P&L, intermediate histories, metrics
    """
    S_arr = np.asarray(prices_series).astype(float).ravel()
    dates_local = prices_series.index
    n = len(S_arr)
    assert S_arr.ndim == 1, "S_arr must be 1-D"
    assert np.issubdtype(S_arr.dtype, np.floating)

    # state
    cash              = 0.0
    stock_holding 	  = 0.0
    shares_history 	  = []
    cash_history 	  = []
    portfolio_history = []  # cash + stock value - option liability
    
    # Sell one call at t0: receive model price
    tau0 = taus[0]
    if model == 'gbm':
        price0, _, delta0, _ = price_and_delta_gbm(S_arr[0], tau0, sigma_for_gbm)
    else:
        # for merton compute MC price/delta at tau0
        price0, _, delta0, _ = price_and_delta_merton_at_tau(S_arr[0], tau0)

    # seller is short 1 call; to delta-hedge short call I HOLD +delta shares (long delta).
    # receive option premium into cash
    cash += price0
    # initial buy of delta0 shares
    trade_shares = float(delta0) - stock_holding
    trade_cost = np.abs(trade_shares) * S_arr[0] * tc
    cash -= trade_shares * S_arr[0] + trade_cost
    stock_holding += trade_shares

    # record initial
    shares_history.append(stock_holding)
    cash_history.append(cash)
    portfolio_history.append(cash + stock_holding * S_arr[0])  # does not subtract option yet

    # iterate daily rebalancing from day 1 to n-1 (rebalancing ON each observed day)
    for i in range(1, n):
        # accrue interest on cash from previous day (simple daily)
        cash *= (1.0 + r / TRADING_DAYS)

        tau_i = taus[i]
        S_i = S_arr[i]

        # compute model delta at this day (for remaining time tau_i)
        if tau_i <= 0:
            delta_i = 0.0  # at expiry delta should be 0/1 but rebalancing at expiry handled after loop
            delta_se = 0.0
        else:
            if model == 'gbm':
                _, _, delta_i, delta_se = price_and_delta_gbm(S_i, tau_i, sigma_for_gbm)
            else:
                # use cached MC pricing for this tau
                _, _, delta_i, delta_se = price_and_delta_merton_at_tau(S_i, tau_i)

        # trade to adjust holdings to delta_i (buy/sell shares)
        trade_shares = float(delta_i) - stock_holding
        trade_cost = np.abs(trade_shares) * S_i * tc
        cash -= trade_shares * S_i + trade_cost
        stock_holding += trade_shares

        # record
        shares_history.append(stock_holding)
        cash_history.append(cash)
        portfolio_history.append(cash + stock_holding * S_i)

    # at expiry (last observed price) compute final P&L:
    S_T = S_arr[-1]
    
    # option payoff
    payoff 	   		  = max(S_T - K, 0.0)
    final_cash 		  = cash * (1.0 + r / TRADING_DAYS)  # accrue on last day to settle
    final_stock_value = stock_holding * S_T
    final_pnl 		  = final_cash + final_stock_value - payoff
    assert isinstance(final_pnl, float)

    # produce time series of mark-to-market P&L excluding option (for trajectories): cash + stock (liability not subtracted)
    mm_series = np.array(portfolio_history) - 0.0
    # for clarity compute realized P&L if I liquidated at each day (not subtracting option liability):
    # but I'll focus on final_pnl across strategies.

    out = {
        'model'				: model,
        'final_pnl'			: final_pnl,
        'final_cash'		: final_cash,
        'final_stock_value ': final_stock_value,
        'payoff'		    : payoff,
        'shares_hist'	    : np.array(shares_history),
        'cash_hist'			: np.array(cash_history),
        'mtm_hist'			: np.array(portfolio_history),
        'dates'				: dates_local,
        'S'					: S_arr
    }
    return out

# run both strategies
print("Running hedging backtests. This may take a moment (MC pricing cached per tau).")
res_gbm = run_delta_hedge(prices_bt, model='gbm', sigma_for_gbm=sigma_gbm_mle, tc=tc_pct)
res_merton = run_delta_hedge(prices_bt, model='merton', sigma_for_gbm=None, tc=tc_pct)

# collect all P&L results
pnls = {
    'GBM': res_gbm['final_pnl'],
    'Merton': res_merton['final_pnl']
}

# metrics function
def pnl_metrics(pnl_vals, loss_thresh=loss_threshold):
    # accept scalar or array; convert to array
    arr = np.atleast_1d(pnl_vals).astype(float)
    metrics = {
        'mean'	  : np.mean(arr),
        'std' 	  : np.std(arr, ddof=1) if arr.size > 1 else 0.0,
        'median'  : np.median(arr),
        '5%'	  : np.quantile(arr, 0.05),
        '95%'	  : np.quantile(arr, 0.95),
        'skewness': float(skew(arr)) if arr.size > 1 else 0.0,
        'frac_below_threshold': float((arr < loss_thresh).sum() / arr.size)
    }
    return metrics

# I only have single-run P&Ls for each model here...

print("\nFinal P&L results (selling 1 call and hedging daily):")
for k, v in pnls.items():
    print(f"  {k:6s} final P&L = {v:10.6f}")

# print richer metrics by doing a small rep-sample to get distribution due to MC noise
# I'll do small n_rep to estimate the Monte-Carlo variability of the final P&Lfor distribution
n_rep = 8
print(f"\nEstimating P&L distribution by re-running Merton MC pricing {n_rep} times (this will take extra time)...")
merton_pnls_rep = []
for rep in range(n_rep):
    # clear the tau cache to force fresh unit path draws (introduces MC randomness)
    _merton_tau_cache.clear()
    rres = run_delta_hedge(prices_bt, model='merton', tc=tc_pct)
    merton_pnls_rep.append(rres['final_pnl'])
merton_pnls_rep = np.array(merton_pnls_rep)

gbm_pnls_rep = []
# GBM price/delta are deterministic so no need to resample; but for parity I can simulate small noise via not doing so.
gbm_pnls_rep = np.array([res_gbm['final_pnl']])

print("\nGBM strategy (single run) metrics:")
print(pnl_metrics(gbm_pnls_rep))
print("\nMerton strategy (MC resampling) metrics:")
print(pnl_metrics(merton_pnls_rep))

# histogram and cumulative P&L trajectories
plt.figure(figsize=(12,5))

plt.subplot(1,2,1)
plt.hist(merton_pnls_rep, bins=15, alpha=0.8, label='Merton (resampled)')
plt.axvline(res_gbm['final_pnl'], color='C1', linestyle='--', label='GBM final P&L (single run)')
plt.title("P&L distribution (Merton resampling)")
plt.xlabel("Final P&L")
plt.legend()

plt.subplot(1,2,2)

# cumulative P&L trajectories: compute mark-to-market (cash + stock - liability at each date) for both models
# I derive liability = model price at each date; for Merton we use cached pricing (potentially expensive) but reuse cached tau keys
def mtm_series_for_model(res, model):
    # compute liability (model price) series across dates to subtract
    S_arr = res['S']
    n = len(S_arr)
    liabilities = np.zeros(n)
    for i in range(n):
        tau_i = taus[i]
        if tau_i <= 0:
            liabilities[i] = max(S_arr[i] - K, 0.0)
        else:
            if model == 'gbm':
                liabilities[i] = bs_call(S_arr[i], K, sigma_gbm_mle, tau_i, r)
            else:
                p, _, _, _ = price_and_delta_merton_at_tau(S_arr[i], tau_i)
                liabilities[i] = p
    mtm = res['cash_hist'] + res['shares_hist'] * S_arr - liabilities
    # cumulative realized P&L path
    return mtm

mtm_gbm = mtm_series_for_model(res_gbm, 'gbm')
mtm_merton = mtm_series_for_model(res_merton, 'merton')

plt.plot(dates, np.cumsum(np.diff(np.concatenate([[0.0], mtm_gbm]))), label='GBM cum P&L (mtm)', alpha=0.9)
plt.plot(dates, np.cumsum(np.diff(np.concatenate([[0.0], mtm_merton]))), label='Merton cum P&L (mtm)', alpha=0.9)
plt.title("Cumulative P&L trajectories (mtm)")
plt.xlabel("Date")
plt.ylabel("Cumulative P&L (mtm)")
plt.legend()
plt.tight_layout()
plt.show()

# I won't need this but just in case
_backtest_results = {
    'res_gbm'        : res_gbm,
    'res_merton'	 : res_merton,
    'merton_pnls_rep': merton_pnls_rep,
    'gbm_pnls_rep'	 : gbm_pnls_rep
}

# candidate tickers suggested by chatgpt
candidates = ["TSLA", "NVDA", "AMD", "AAPL", "AMZN", "INTC", "IBM", "GME", "PLTR", "NFLX", "BA", "UAL", "F"]

# settings for quick scan
scan_mc_price = 1000 # fewer MC sims for speed during scanning
scan_n_rep    = 5    # small rep for a tiny distribution

results_scan = []

print("Scanning tickers for return kurtosis (calibration period) and running fast backtests...")
t0 = time.time()
for tick in candidates:
    try:
        df = yf.download(tick, start=calib_start, end=calib_end, progress=False)
        if 'Close' not in df.columns or df['Close'].dropna().shape[0] < 30:
            # skip thin series, although I think this doesn't happen with these samples
            continue
        prices = df['Close'].dropna()
        # compute daily log returns and kurtosis
        lr = np.log(prices / prices.shift(1)).dropna().values
        if lr.size < 20:
            continue
        kurt_excess = kurtosis(lr, fisher=True, bias=False)  # excess kurtosis
        # quick reduce: run a single GBM and Merton backtest with smaller MC for speed
        # I need to calibrate quickly for this ticker: estimate GBM MLE and use default Merton guesses (lightweight)
        # Use the same calibration function if available; otherwise do quick estimates
        # For speed I just re-use the global sigma_gbm_mle if calibration functions not present; but try to calibrate basic sigma:

        mu_sample = np.mean(lr)
        s2_sample = np.var(lr, ddof=0)
        dt = 1.0 / TRADING_DAYS
        sigma_hint = np.sqrt(s2_sample / dt)

        # For quick Merton params use previously calibrated global if available; otherwise tiny default
        sigma_merton_quick = sigma_merton_mle if 'sigma_merton_mle' in globals() else sigma_hint
        lambda_quick = lambda_merton_mle if 'lambda_merton_mle' in globals() else 0.2
        muJ_quick = mu_J_mle if 'mu_J_mle' in globals() else -0.1
        sigmaJ_quick = sigma_J_mle if 'sigma_J_mle' in globals() else 0.15

        # Quick backtest wrapper: I will call run_delta_hedge but with faster MC; to do this I temporarily adjust the global MC cache behavior
        # I can monkey-patch _merton_tau_cache and run_delta_hedge uses cached unit paths per tau; to force smaller MC I set n_mc_price globally if it exists
        # Simpler: call run_delta_hedge but with existing deterministic GBM and Merton settings — if your run_delta_hedge doesn't accept MC size param,
        # then I rely on the previously set global n_mc_price. I'll temporarily lower it if present.
        old_n_mc = globals().get('n_mc_price', None)
        globals()['n_mc_price'] = scan_mc_price

        # download backtest range prices for this ticker
        df_bt = yf.download(tick, start=backtest_start, end=backtest_end, progress=False)
        if 'Close' not in df_bt.columns or df_bt['Close'].dropna().shape[0] < 10:
            # not enough backtest data
            if old_n_mc is not None:
                globals()['n_mc_price'] = old_n_mc
            continue
        prices_bt_local = df_bt['Close'].dropna()
        # Ensure series 1-D
        if isinstance(prices_bt_local, pd.DataFrame) and prices_bt_local.shape[1] == 1:
            prices_bt_local = prices_bt_local.iloc[:,0]

        # Run backtests using your run_delta_hedge function (it expects the series)
        # Here I assume run_delta_hedge uses sigma_gbm_mle and merton params from globals.
        gbm_res = run_delta_hedge(prices_bt_local, model='gbm', sigma_for_gbm=sigma_hint, tc=tc_pct)
        if '_merton_tau_cache' in globals():
            _merton_tau_cache.clear()
        merton_res = run_delta_hedge(prices_bt_local, model='merton', tc=tc_pct)

        # Restore MC setting
        if old_n_mc is not None:
            globals()['n_mc_price'] = old_n_mc

        # ensure all numeric values are scalars, not 0-D arrays otherwise this will BREAK the code
        def to_scalar(x):
            """Convert anything (scalar, 0D/1D array, pd.Series) to float safely."""
            if isinstance(x, (np.ndarray, pd.Series, list)):
                return float(np.ravel(x)[0])
            return float(x)

        gbm_pnl = to_scalar(gbm_res['final_pnl'])
        merton_pnl = to_scalar(merton_res['final_pnl'])
        diff = merton_pnl - gbm_pnl
        kurt_excess_val = to_scalar(kurt_excess)
        results_scan.append({
			'ticker': tick,
			'kurt_excess': kurt_excess_val,
			'gbm_pnl': gbm_pnl,
			'merton_pnl': merton_pnl,
			'diff_merton_minus_gbm': diff
		})

        print(f"{tick:6s}  kurt_excess={kurt_excess_val:6.2f}  GBM_pnl={gbm_pnl:8.3f}  "f"Merton_pnl={merton_pnl:8.3f}  Δ={diff:8.3f}")

    except Exception as e:
        print(f"Skipped {tick}: {e}")
        # restore n_mc_price if changed
        if old_n_mc is not None:
            globals()['n_mc_price'] = old_n_mc
        continue

t1 = time.time()
print(f"\nScan completed in {(t1-t0):.1f}s for {len(results_scan)} tickers.")

# summarize and pick top candidates (largest positive improvement: Merton better than GBM, i.e., diff < 0 means Merton loss smaller)
df_scan = pd.DataFrame(results_scan).sort_values(by='diff_merton_minus_gbm')
display(df_scan)

better = df_scan[df_scan['diff_merton_minus_gbm'] < 0].copy()
if better.empty:
    print("No candidate showed Merton beating GBM in this quick scan.")
else:
    print("\nTickers where Merton produced better final P&L (quick scan):")
    display(better)