The following script reproduces the analysis and plots:
# Claude launch vs H100 rental price spike analysis
# Inputs:
#   - claude_launch.csv (columns: model_name, version_date)
#   - h100_price.csv    (columns: date, source_type, rental_price)

import os
from datetime import date
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from scipy.stats import hypergeom

# ----------------------------
# Load data
# ----------------------------
claude = pd.read_csv("claude_launch.csv")
h100 = pd.read_csv("h100_price.csv")

claude["version_date"] = pd.to_datetime(claude["version_date"]).dt.normalize()
h100["date"] = pd.to_datetime(h100["date"]).dt.normalize()
h100 = h100.dropna(subset=["source_type"]).copy()

# Aggregate launches by day (some days have multiple model releases)
launches_by_date = (
    claude.groupby("version_date")["model_name"]
    .apply(lambda x: ", ".join(sorted(x)))
    .to_dict()
)
launch_counts = claude["version_date"].value_counts().sort_index()

event_days_all = sorted(launch_counts.index.tolist())

# ----------------------------
# Build daily "as-of" price series per source_type (forward-fill only)
# ----------------------------
date_range = pd.date_range(h100["date"].min(), h100["date"].max(), freq="D")

price_asof = {}
for st, dfst in h100.groupby("source_type"):
    s = dfst.set_index("date").sort_index()["rental_price"]
    # forward-fill; keep NaN before first observation (no backfill)
    price_asof[st] = s.reindex(date_range).ffill()

price_asof_df = pd.DataFrame(price_asof)

log_price = np.log(price_asof_df)
returns = log_price.diff()

# Valid event days must have data at d-1 and d+7 (for 0..7 metrics)
valid_events = {}
for st in price_asof_df.columns:
    lp = log_price[st]
    evs = []
    for d in event_days_all:
        if (d - pd.Timedelta(days=1) in lp.index) and (d + pd.Timedelta(days=7) in lp.index):
            if pd.notna(lp.loc[d - pd.Timedelta(days=1)]) and pd.notna(lp.loc[d + pd.Timedelta(days=7)]):
                evs.append(d)
    valid_events[st] = evs

# ----------------------------
# Spike enrichment test
# ----------------------------
def spike_enrichment(st: str, window: int, q: float = 0.99) -> dict:
    """Top (1-q) positive spikes concentration within +/- window days of launches."""
    rs = returns[st].dropna()
    idx = rs.index
    evs = valid_events[st]

    # days covered by event windows
    near = set()
    for d in evs:
        for k in range(-window, window + 1):
            dt = d + pd.Timedelta(days=k)
            if dt in idx:
                near.add(dt)

    N = len(idx)         # population size (days)
    K = len(near)        # "success" days (near launch)
    thresh = rs.quantile(q)
    spikes = rs[rs >= thresh].index
    n = len(spikes)      # number of spike days
    x = sum([1 for d in spikes if d in near])

    expected = n * K / N
    p_enrich = hypergeom(N, K, n).sf(x - 1)  # P(X >= x)
    return {
        "source_type": st,
        "window": window,
        "threshold": float(thresh),
        "N_days": N,
        "near_days": K,
        "num_spikes": n,
        "spikes_near": x,
        "expected_near": expected,
        "p_enrichment": float(p_enrich),
    }

# Example: compute for multiple windows
spike_rows = []
for st in price_asof_df.columns:
    for w in [3, 7, 14]:
        spike_rows.append(spike_enrichment(st, w, q=0.99))
spike_df = pd.DataFrame(spike_rows)

# ----------------------------
# Event study (average cumulative % change around launches)
# ----------------------------
def bootstrap_ci(mat: np.ndarray, n_boot: int = 3000, ci: float = 0.95, seed: int = 0):
    rng = np.random.default_rng(seed)
    n = mat.shape[0]
    boot_means = []
    for _ in range(n_boot):
        sample = mat[rng.integers(0, n, size=n)]
        boot_means.append(sample.mean(axis=0))
    boot_means = np.vstack(boot_means)
    alpha = (1 - ci) / 2
    lower = np.quantile(boot_means, alpha, axis=0)
    upper = np.quantile(boot_means, 1 - alpha, axis=0)
    return lower, upper

def event_study_paths(st: str, W: int = 30):
    lp = log_price[st]
    evs = valid_events[st]
    rel_days = np.arange(-W, W + 1)

    mat = []
    for d in evs:
        idx = d + pd.to_timedelta(rel_days, unit="D")
        baseline = lp.loc[d - pd.Timedelta(days=1)]
        mat.append(lp.reindex(idx).values - baseline)
    mat = np.vstack(mat)

    mean_log = mat.mean(axis=0)
    lo, hi = bootstrap_ci(mat, n_boot=3000, seed=7)

    mean_pct = (np.exp(mean_log) - 1) * 100
    lo_pct = (np.exp(lo) - 1) * 100
    hi_pct = (np.exp(hi) - 1) * 100
    return rel_days, mean_pct, lo_pct, hi_pct

# ----------------------------
# Plotting helpers
# ----------------------------
OUTDIR = "claude_h100_report_assets"
os.makedirs(OUTDIR, exist_ok=True)

def save_price_plot(st: str):
    s = price_asof_df[st].dropna()
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(s.index, s.values)
    for d in valid_events[st]:
        ax.axvline(d, linestyle="--", linewidth=1, alpha=0.7)
    ax.set_title(f"H100 rental price - {st} (as-of daily)")
    ax.set_xlabel("Date")
    ax.set_ylabel("Rental price")
    ax.xaxis.set_major_locator(mdates.MonthLocator(interval=3))
    ax.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m"))
    fig.autofmt_xdate()
    ax.grid(True, alpha=0.3)
    path = os.path.join(OUTDIR, f"price_{st}.png")
    fig.tight_layout()
    fig.savefig(path, dpi=200)
    plt.close(fig)

def save_event_study_plot(st: str, W: int = 30):
    x, mean_pct, lo_pct, hi_pct = event_study_paths(st, W=W)
    fig, ax = plt.subplots(figsize=(10, 4))
    ax.plot(x, mean_pct, label="Mean")
    ax.fill_between(x, lo_pct, hi_pct, alpha=0.2, label="95% CI")
    ax.axvline(0, linestyle="--", linewidth=1)
    ax.axhline(0, linewidth=1, alpha=0.5)
    ax.set_title(f"Event study: cumulative price change around Claude launches ({st})")
    ax.set_xlabel("Days relative to launch (day 0 = launch date)")
    ax.set_ylabel("Cumulative change vs day -1 (%)")
    ax.grid(True, alpha=0.3)
    path = os.path.join(OUTDIR, f"event_study_{st}.png")
    fig.tight_layout()
    fig.savefig(path, dpi=200)
    plt.close(fig)

for st in price_asof_df.columns:
    save_price_plot(st)
    save_event_study_plot(st)

print("Done. Plots saved to:", OUTDIR)
print("Spike enrichment summary (top 1% positive spikes):")
print(spike_df[["source_type", "window", "num_spikes", "spikes_near", "p_enrichment"]])