from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
from cycler import cycler
from quantfinlab.dataio import load_ohlcv, load_option_chain, load_par_yield_curve
from quantfinlab.volatility.realized import log_returns
from quantfinlab.volatility.forecasting import future_realized_variance, rolling_arch_forecasts_weekly, score_forecasts_by_model
from quantfinlab.volatility.har import rolling_har_forecasts
from quantfinlab.volatility.rough import daily_variance, log_variance, fbm_cholesky_paths, moment_scaling, hurst_from_moments, hurst_from_moments_pooled, rough_kernel_forecasts, rough_forecast_frame
from quantfinlab.options.quote_cleaning import attach_spot_from_series, clean_option_quotes, convert_quotes_to_usd_equivalent, surface_ready_quotes
from quantfinlab.options.rates_dividends import attach_rates, add_discount_factors, infer_carry_from_forward, infer_dividend_yield_from_forward
from quantfinlab.options.parity import infer_forwards_from_parity
from quantfinlab.options.iv import implied_vol_table
from quantfinlab.options.surface import surface_grid, fit_surface_panel, surface_iv_grid, surface_residuals
from quantfinlab.calibration.fft_cos import surface_target_grid_quotes, fit_daily_models, compare_fourier_models, residual_by_bucket
from quantfinlab.options.rough import atm_skew_term_structure, skew_power_law, forward_variance_curve, simulate_rbergomi, rbergomi_smile, rbergomi_calibration, fit_rough_heston_dates, compare_heston_rough_heston, rough_heston_residuals, riccati_convergence, rough_delta_grid
from quantfinlab.plotting.volatility import price_and_variance, fbm_roughness, moment_ladder, hurst_linearity, forecast_race, forecast_path, smile_slices, smooth_surface, skew_power, fit_comparison, maturity_error, rbergomi_paths, rbergomi_fit, riccati_stability, parameter_paths, delta_slices, delta_gap
palette = ["#069AF3", "#FE420F", "#00008B", "#008080", "#CC79A7", "#9614fa", "#DC143C", "#7BC8F6", "#0072B2", "#04D8B2", "#800080", "#FF8072"]
plt.rcParams["axes.prop_cycle"] = cycler(color=palette)
data_dir = Path("../data")
if not data_dir.exists():
data_dir = Path("data")
cache_dir = data_dir / "cache" / "project18" / "btc"
cache_dir.mkdir(parents=True, exist_ok=True)
ann_days = 365.0
ann_trading = 365
seed = 7
q_values = (0.5, 1.0, 1.5, 2.0, 3.0)
lags = (1, 2, 4, 8, 16, 32)
horizons = (1, 5, 10, 21, 42, 63)
btc_history_path = data_dir / "btc_usd_ohlcv.csv"
if not btc_history_path.exists():
btc_history_path = data_dir / "btc_usd_ohlcv.csv"
btc_spot = load_ohlcv(btc_history_path, source="yfinance_csv", fields=("close",))
btc_raw = load_option_chain(data_dir / "btc_options_chain.parquet", source="btc_deribit", annualization_days=ann_days)
us_par_btc = load_par_yield_curve(data_dir / "us_treasury_yields.csv", source="us_treasury")
btc_ret = log_returns(btc_spot["close"]).dropna()
btc_var_1d = daily_variance(btc_ret, annualization=ann_trading)
btc_var = (btc_var_1d.rolling(10).sum() / 10.0).dropna().rename("btc_project5_10d_daily_variance")
btc_log_var = log_variance(btc_var)
btc_targets = future_realized_variance(btc_ret, horizons=horizons, annualization=ann_trading)
fbm = fbm_cholesky_paths(h_values=(0.05, 0.10, 0.20, 0.50, 0.80), n_steps=384, n_paths=4, seed=seed)
scaling = moment_scaling(btc_log_var, q_values=q_values, lags=lags)
hurst = hurst_from_moments(scaling, q_col="q", lag_col="lag", moment_col="moment")
hurst_pooled = hurst_from_moments_pooled(scaling)
h = float(np.clip(float(hurst_pooled["H"].iloc[0]) if not hurst_pooled.empty and np.isfinite(float(hurst_pooled["H"].iloc[0])) else hurst["h"].median(), 0.04, 0.35))
train_window = min(180, max(90, len(btc_var) // 2))
rough_fc = rough_kernel_forecasts(btc_var, h=h, horizons=horizons, train_window=train_window, signal_step=7, annualization=ann_trading)
har_fc = rolling_har_forecasts(btc_ret, horizons=horizons, train_window=train_window, refit_every=7, annualization=ann_trading, use_log=True)
btc_forecast_data = rough_forecast_frame(rough_fc=rough_fc, har_fc=har_fc, arch_fc=pd.DataFrame(), rv_targets=btc_targets, horizons=horizons)
forecast_scores = score_forecasts_by_model(btc_forecast_data) if not btc_forecast_data.empty else pd.DataFrame()
btc_quotes = attach_spot_from_series(btc_raw, btc_spot["close"], date_col="date", spot_col="spot", method="previous", overwrite=False)
btc_quotes = convert_quotes_to_usd_equivalent(btc_quotes, unit="base")
btc_unit_audit = pd.DataFrame([{
"input_file": str(btc_history_path.name),
"price_unit": str(btc_quotes["price_unit_detected"].iloc[0]),
"spot_median": float(pd.to_numeric(btc_quotes["spot"], errors="coerce").median()),
"raw_mid_median_btc": float(pd.to_numeric(btc_quotes["mid_raw"], errors="coerce").replace(0.0, np.nan).median()),
"usd_mid_median": float(pd.to_numeric(btc_quotes["mid"], errors="coerce").replace(0.0, np.nan).median()),
}])
clean_quotes, clean_steps = clean_option_quotes(btc_quotes, min_dte=3, max_dte=180, moneyness_range=(0.45, 1.75), max_relative_spread=0.60, closest_atm_pairs=90, min_pairs_per_expiry=0, annualization_days=ann_days)
rate_start = clean_quotes["date"].min() - pd.Timedelta(days=30)
rate_end = clean_quotes["date"].max()
us_par_btc = us_par_btc.loc[(us_par_btc.index >= rate_start) & (us_par_btc.index <= rate_end)].copy()
clean_quotes = attach_rates(clean_quotes, us_par_btc, date_col="date", tau_col="tau", out_col="rate")
clean_quotes = add_discount_factors(clean_quotes, rate_col="rate", tau_col="tau", out_col="discount_factor")
clean_quotes = infer_forwards_from_parity(clean_quotes, date_col="date", expiry_col="expiry", strike_col="strike", option_type_col="option_type", mid_col="mid", discount_col="discount_factor", spot_col="spot", out_col="forward")
clean_quotes = infer_carry_from_forward(clean_quotes, spot_col="spot", forward_col="forward", tau_col="tau", out_col="implied_carry")
clean_quotes = infer_dividend_yield_from_forward(clean_quotes, rate_col="rate", carry_col="implied_carry", out_col="implied_dividend_yield")
iv_quotes = implied_vol_table(clean_quotes, price_cols=("bid", "mid", "ask"), option_type_col="option_type", spot_col="spot", forward_col="forward", strike_col="strike", tau_col="tau", rate_col="rate", discount_col="discount_factor", model="black76", engine="auto")
surface_quotes = surface_ready_quotes(iv_quotes, min_iv=0.05, max_iv=3.50, min_tau=3 / ann_days, max_tau=180 / ann_days, min_weight=0.05, max_weight=20.0)
surface_fit = fit_surface_panel(surface_quotes, date_col="date", k_col="k", tau_col="tau", iv_col="iv_mid", weight_col="surface_weight", visual_params=dict(n_k_basis=10, n_tau_basis=7, degree=3, lambda_k=3.0, lambda_tau=2.0), dupire_params=dict(n_k_basis=9, n_tau_basis=6, degree=3, lambda_k=10.0, lambda_tau=16.0), min_quotes=80, min_expiries=4)
visual_fits_by_date = {pd.Timestamp(k).normalize(): v for k, v in surface_fit["visual_fits"].items()}
btc_tau_values = np.array([7, 10, 14, 21, 30, 45, 60, 90, 120], dtype=float) / ann_days
btc_skew_rows = []
for d, fit in visual_fits_by_date.items():
qd = surface_quotes[surface_quotes["date"].eq(d)].copy()
psi_d = atm_skew_term_structure(qd, fit=fit, k_col="k", tau_col="tau", iv_col="iv_mid", tau_values=btc_tau_values, dk=0.01)
for window, lo, hi in [("short_7_60", 7, 60), ("medium_14_90", 14, 90), ("old_14_150", 14, 150)]:
sub = psi_d[psi_d["tau_days"].between(lo, hi)].copy()
row = skew_power_law(sub, tau_col="tau", skew_col="atm_skew").iloc[0].to_dict()
row["date"] = d
row["window"] = window
row["h_from_skew"] = row.get("h", np.nan)
row["quote_count"] = int(len(qd))
btc_skew_rows.append(row)
btc_skew_history = pd.DataFrame(btc_skew_rows)
surface_summary = surface_fit["fit_summary"]
surface_visual_summary = surface_summary[surface_summary["fit"].eq("visual")].copy() if not surface_summary.empty else pd.DataFrame()
btc_main_candidates = btc_skew_history[btc_skew_history["window"].eq("short_7_60")].copy()
if not surface_visual_summary.empty:
btc_main_candidates = btc_main_candidates.merge(surface_visual_summary[["date", "weighted_rmse", "quote_count", "number_of_maturities"]], on="date", how="left", suffixes=("", "_surface"))
else:
btc_main_candidates["weighted_rmse"] = np.nan
btc_main_candidates["number_of_maturities"] = np.nan
btc_high_quality_skew = btc_main_candidates[btc_main_candidates["r2"].ge(0.80) & btc_main_candidates["n"].ge(6)].copy()
btc_selection_pool = btc_high_quality_skew.copy()
if btc_selection_pool.empty:
btc_selection_pool = btc_main_candidates[btc_main_candidates["n"].ge(4)].copy()
if btc_selection_pool.empty:
btc_selection_pool = btc_main_candidates.copy()
btc_selection_pool["surface_quality_rank"] = btc_selection_pool["weighted_rmse"].rank(method="first", ascending=True, na_option="bottom")
main_date = btc_selection_pool.sort_values(["r2", "surface_quality_rank", "quote_count"], ascending=[False, True, False])["date"].iloc[0]
q_main = surface_quotes[surface_quotes["date"].eq(main_date)].copy()
fit_visual = visual_fits_by_date[pd.Timestamp(main_date).normalize()]
k_grid = surface_grid(q_main, k_col="k", tau_col="tau", k_quantiles=(0.02, 0.98), tau_quantiles=(0.02, 0.98), n_k=65, n_tau=35, annualization_days=ann_days)
sigma_grid = surface_iv_grid(fit_visual, k_grid)
surface_err = surface_residuals(q_main, fits={"smooth": fit_visual}, k_col="k", tau_col="tau", iv_col="iv_mid")
psi = atm_skew_term_structure(q_main, fit=fit_visual, k_col="k", tau_col="tau", iv_col="iv_mid", tau_values=btc_tau_values, dk=0.01)
skew_fit = skew_power_law(psi[psi["tau_days"].between(7, 60)], tau_col="tau", skew_col="atm_skew")
xi = forward_variance_curve(q_main, fit=fit_visual, tau_values=np.array([7, 14, 21, 30, 45, 60, 90, 120], dtype=float) / ann_days)
calibration_quotes = surface_target_grid_quotes(
surface_quotes,
surface_fit["visual_fits"],
min_dte=7,
max_dte=150,
max_abs_log_moneyness=0.45,
min_quotes_per_date=45,
min_source_quotes_per_date=70,
annualization_days=ann_days,
iv_floor=0.05,
iv_cap=3.50,
iv_error_floor=0.020,
iv_error_cap=0.140,
adaptive_k_grid=True,
adaptive_k_min_nodes=7,
adaptive_k_max_nodes=11,
)
date_counts = calibration_quotes.groupby("date").size()
eligible_btc_fit_dates = date_counts[date_counts >= 45].index.sort_values()
if len(eligible_btc_fit_dates) > 12:
fit_dates = eligible_btc_fit_dates[np.linspace(0, len(eligible_btc_fit_dates) - 1, 12).round().astype(int)]
else:
fit_dates = eligible_btc_fit_dates
h_skew_main = float(skew_fit["h"].iloc[0]) if not skew_fit.empty and np.isfinite(skew_fit["h"].iloc[0]) else h
btc_h_source = btc_high_quality_skew.copy()
if btc_h_source.empty:
btc_h_source = btc_main_candidates[btc_main_candidates["r2"].ge(0.70)].copy()
btc_h_source = btc_h_source[btc_h_source["h_from_skew"].between(0.04, 0.45)].copy()
btc_h_by_date = btc_h_source.set_index("date")["h_from_skew"].dropna().sort_index()
btc_heston_cache = cache_dir / "btc_heston_surface_grid_daily.parquet"
btc_heston_fit_cache = cache_dir / "btc_heston_surface_grid_fit.parquet"
btc_rough_cache = cache_dir / "btc_rough_heston_surface_grid_skew_anchor_v2_daily.parquet"
btc_rough_fit_cache = cache_dir / "btc_rough_heston_surface_grid_skew_anchor_v2_fit.parquet"
if btc_heston_cache.exists() and btc_heston_fit_cache.exists():
heston_daily = pd.read_parquet(btc_heston_cache)
heston_prices = pd.read_parquet(btc_heston_fit_cache)
else:
heston_out = fit_daily_models(calibration_quotes, ["heston"], calibration_dates=fit_dates, min_quotes=45, max_nfev=18, engine="cpp", n_terms=112, truncation_width=16.0)
heston_daily = heston_out["params"]
heston_prices = heston_out["fit"]
heston_daily.to_parquet(btc_heston_cache, index=False)
heston_prices.to_parquet(btc_heston_fit_cache, index=False)
if btc_rough_cache.exists() and btc_rough_fit_cache.exists():
rough_daily = pd.read_parquet(btc_rough_cache)
rough_prices = pd.read_parquet(btc_rough_fit_cache)
else:
rough_out = fit_rough_heston_dates(calibration_quotes, calibration_dates=fit_dates, h_start=h_skew_main, h_by_date=btc_h_by_date, h_mode="fixed", min_quotes=45, max_nfev=14, engine="cpp", n_terms=88, truncation_width=16.0, riccati_steps=192, lambda_skew=3.0, surface_fit_map=surface_fit["visual_fits"])
rough_daily = rough_out["params"]
rough_prices = rough_out["fit"]
rough_daily.to_parquet(btc_rough_cache, index=False)
rough_prices.to_parquet(btc_rough_fit_cache, index=False)
model_table = compare_heston_rough_heston(heston_daily=heston_daily, rough_daily=rough_daily)
heston_table = compare_fourier_models(daily=heston_daily)
heston_buckets = residual_by_bucket(heston_prices)
rough_buckets = rough_heston_residuals(rough_prices, scale_col="calib_scale_px")
rough_valid = rough_daily[
rough_daily["success"].astype(bool)
& pd.to_numeric(rough_daily.get("weighted_price_rmse"), errors="coerce").replace([np.inf, -np.inf], np.nan).notna()
& pd.to_numeric(rough_daily.get("finite_price_share", pd.Series(1.0, index=rough_daily.index)), errors="coerce").ge(0.999)
].copy()
if rough_valid.empty:
rough_valid = rough_daily[pd.to_numeric(rough_daily.get("weighted_price_rmse"), errors="coerce").replace([np.inf, -np.inf], np.nan).notna()].copy()
if rough_valid.empty:
raise RuntimeError("No finite BTC rough-Heston parameter row is available for diagnostics.")
display_date = pd.Timestamp(main_date).normalize()
if display_date not in set(pd.to_datetime(rough_valid["date"]).dt.normalize()):
display_date = pd.Timestamp(rough_valid.sort_values("weighted_price_rmse")["date"].iloc[0]).normalize()
rough_fit = rough_prices[rough_prices["date"].eq(display_date)].copy()
heston_fit = heston_prices[heston_prices["date"].eq(display_date)].copy()
q_main_target = calibration_quotes[calibration_quotes["date"].eq(pd.Timestamp(main_date).normalize())].copy()
if q_main_target.empty:
q_main_target = q_main.copy()
rough_param_row = rough_valid[rough_valid["date"].eq(display_date)].sort_values("weighted_price_rmse").iloc[0]
riccati_test = riccati_convergence(q_main_target, params=rough_param_row, n_grid_values=(128, 256, 512), n_terms=96, engine="cpp")
rbergomi = rbergomi_calibration(q_main_target, xi=xi, h_start=h_skew_main, nu_start=2.0, rho_start=-0.70, paths=5200, steps=88, restarts=6, seed=seed, engine="numba", use_sobol=True, lambda_skew=1.0)
rbergomi_sim = simulate_rbergomi(spot=float(q_main["spot"].median()), xi=xi, h=float(rbergomi["params"]["h"]), nu=float(rbergomi["params"]["nu"]), rho=float(rbergomi["params"]["rho"]), tau=float(np.nanmax(q_main["tau"])), paths=1200, steps=80, seed=seed, engine="numba", rate=float(q_main["rate"].median()), dividend_yield=float(q_main["implied_dividend_yield"].median()), antithetic=True)
rbergomi_iv = rbergomi_smile(q_main_target, xi=xi, params=rbergomi["params"], maturity_days=(7, 14, 30, 60, 90), paths=3000, steps=88, seed=seed, engine="numba")
delta_map = rough_delta_grid(q_main, heston_params=heston_daily.dropna().iloc[-1], rough_params=rough_param_row, rbergomi_params=rbergomi["params"], xi=xi, surface_fit=fit_visual, k_values=np.linspace(-0.35, 0.25, 49), tau_days=np.array([7, 14, 30, 60, 90], dtype=float), n_terms=96, riccati_steps=192, bump=0.005, paths=2500, steps=80, seed=seed, engine="cpp", mc_engine="numba")
summary = pd.DataFrame([{"asset": "BTC", "spot_start": float(btc_spot["close"].dropna().iloc[0]), "spot_end": float(btc_spot["close"].dropna().iloc[-1]), "return_days": int(btc_ret.shape[0]), "option_rows_clean": int(clean_quotes.shape[0]), "surface_rows": int(surface_quotes.shape[0]), "surface_target_rows": int(calibration_quotes.shape[0]), "main_date": pd.Timestamp(main_date).date(), "main_date_rule": "best short-skew R2 with surface RMSE tiebreak", "h_from_rv": float(h), "h_from_skew": float(skew_fit["h"].iloc[0]), "h_used_for_rough_heston": float(h_skew_main), "skew_r2": float(skew_fit["r2"].iloc[0]), "skew_high_quality_dates": int(len(btc_high_quality_skew)), "rough_dates": int(rough_daily["date"].nunique()), "heston_dates": int(heston_daily["date"].nunique())}])
for name, table in {"summary": summary, "btc_unit_audit": btc_unit_audit, "hurst": hurst, "forecast_scores": forecast_scores.sort_values(["horizon", "qlike_var"]).head(20) if not forecast_scores.empty else forecast_scores, "btc_skew_history": btc_skew_history.sort_values(["window", "r2"], ascending=[True, False]).head(12), "skew_fit": skew_fit, "model_table": model_table, "heston_table": heston_table, "rbergomi_params": pd.DataFrame([rbergomi["params"]]), "rbergomi_restart_loss": rbergomi["fit"], "riccati_test": riccati_test, "delta_sample": delta_map[np.isclose(delta_map["k"], 0.0, atol=0.015)].copy()}.items():
print(name)
display(table)
fig, axes = plt.subplots(4, 4, figsize=(30, 22), constrained_layout=True)
axes = axes.ravel()
price_and_variance(axes[0], btc_spot, btc_var, title="BTC Price and Project-5 Realized Variance")
fbm_roughness(axes[1], fbm, title="Fractional Brownian Roughness")
moment_ladder(axes[2], scaling, title="BTC Moment Scaling")
hurst_linearity(axes[3], hurst, title="BTC qH Linearity")
forecast_race(axes[4], forecast_scores, metric="qlike_var", title="BTC Forecast Race by QLIKE")
forecast_path(axes[5], btc_forecast_data, model="rough_kernel", horizon=21, title="BTC 21-Day Rough Forecast")
smile_slices(axes[6], q_main, fit=fit_visual, maturities_days=(7, 14, 30, 60, 90), title="BTC IV Smiles")
smooth_surface(axes[7], k_grid, sigma_grid, q_main, title="BTC Smooth IV Surface")
skew_power(axes[8], psi, skew_fit, title="BTC ATM Skew Power Law")
fit_comparison(axes[9], heston_fit, rough_fit, title="BTC Heston vs Rough Heston")
maturity_error(axes[10], heston_buckets, rough_buckets, title="BTC Error by Maturity and Moneyness")
rbergomi_paths(axes[11], rbergomi_sim, title="BTC rBergomi Paths")
rbergomi_fit(axes[12], rbergomi_iv, q_main, title="BTC rBergomi Smile Fit")
riccati_stability(axes[13], riccati_test, title="BTC Fractional Riccati Stability")
delta_slices(axes[14], delta_map, title="BTC Delta Slices")
delta_gap(axes[15], delta_map, base_col="surface_bsm_delta", title="BTC Rough-Model Delta Gap")
fig.suptitle(f"Project 18 BTC Library Implementation | {pd.Timestamp(main_date).date()}", fontsize=22, y=1.01)
plt.show()