"""Volatility risk premium study: how large are measured forward spreads, how often are they positive,
and does today's trailing spread predict them? No implemented investment returns are measured.

Inputs (shared with the Yang-Zhang study, ../yang-zhang-article/data):
  benchmark/{SYM}_inputs.csv  FlashAlpha archive: iv30 (constant-maturity 30-day ATM implied vol from daily SVI fits),
                              total_rv (squared 5-minute midpoint returns 09:30-16:00 ET plus the squared overnight return)
  yahoo/{SYM}.csv             daily regular-session OHLC, split-adjusted, not dividend-adjusted
Definitions (annualized, percent; variance in vol points squared):
  trailing observable spread vrp_ante = iv30_t - RV20_t         RV20: close-to-close over the 20 sessions ending t
  forward realized spread    vrp_post = iv30_t - RVfwd_t        RVfwd: benchmark over sessions t+1..t+21
  variance spread             var_post = iv30_t^2 - RVfwd_t^2    an ATM-implied variance measurement proxy, not a traded swap
                                                                  or an option-strategy return (vol points squared, before costs)
Run: python study.py -> data/vrp_results.json
"""
import hashlib, json, math, os, sys
import numpy as np
import pandas as pd

HERE = os.path.dirname(os.path.abspath(__file__))
YZ = os.path.normpath(os.path.join(HERE, "..", "yang-zhang-article"))
sys.path.insert(0, YZ)
import estimators as E  # noqa: E402

D = os.path.join(YZ, "data")
OUT = os.path.join(HERE, "data")
os.makedirs(OUT, exist_ok=True)
SYMS = ["SPY", "QQQ", "IWM", "GLD", "TLT", "XLF", "XLK", "EEM", "AAPL", "NVDA", "TSLA"]
ETFS = SYMS[:8]
H = 21
START = "2017-03-01"


def hac_mean(values, lag=25):
    """Pointwise normal-approximation interval for a mean, on the ORIGINAL session grid.

    Missing observations contribute no centered score, but are not treated as observed zeroes.
    Keeping their positions prevents a gap in clean observations from becoming a one-day lag.
    Bartlett weights give a Newey-West long-run variance; these are descriptive intervals,
    not formal tests of nested forecast-model superiority or simultaneous intervals.
    """
    x = np.asarray(values, float)
    keep = np.isfinite(x)
    n = int(keep.sum())
    assert n > lag + 1
    mean = float(x[keep].mean())
    z = np.where(keep, x - mean, 0.0)
    score_var = z @ z
    for k in range(1, lag + 1):
        score_var += 2 * (1 - k / (lag + 1)) * (z[k:] @ z[:-k])
    se = math.sqrt(max(float(score_var), 0.0)) / n
    return {"mean": mean, "se": se, "ci95": [mean - 1.95996398454 * se, mean + 1.95996398454 * se],
            "observations": n, "lag_sessions": lag}


def forecast_comparison(g, trailing="cc20", forward="rvf", end_col="end21"):
    """Same-date expanding OLS comparisons; forecast outcomes must end strictly before origin.

    A single clean training/evaluation sample is shared by all four models and both targets.
    The spread model restricts the coefficients on IV and RV to be equal and opposite.
    All models include an intercept; forecasts are neither clipped nor tuned.
    """
    iv, rv, out = (g[col].to_numpy(float) for col in ("iv", trailing, forward))
    ends = g[end_col].to_numpy(int)
    assert np.all(np.diff(ends) >= 0)
    valid = np.isfinite(iv) & np.isfinite(rv) & np.isfinite(out)
    targets = np.column_stack([out, iv - out])
    inputs = {"iv_only": iv[:, None], "rv_only": rv[:, None],
              "iv_plus_rv": np.column_stack([iv, rv]), "iv_minus_rv": (iv - rv)[:, None]}
    count = np.r_[0, np.cumsum(valid)]
    y0 = np.where(valid[:, None], targets, 0.0)
    sy = np.vstack([np.zeros(2), np.cumsum(y0, axis=0)])
    means = np.full_like(targets, np.nan)
    predictions = {name: np.full_like(targets, np.nan) for name in inputs}
    origins, train_counts, train_last, coefficients = [], [], [], {}
    # Sufficient statistics preserve expanding estimation without a fit on future observations.
    moments = {}
    for name, x in inputs.items():
        z = np.where(valid[:, None], x, 0.0)
        moments[name] = (np.vstack([np.zeros(x.shape[1]), np.cumsum(z, axis=0)]),
                         np.concatenate([np.zeros((1, x.shape[1], x.shape[1])),
                                         np.cumsum(z[:, :, None] * z[:, None, :], axis=0)]),
                         np.concatenate([np.zeros((1, x.shape[1], 2)),
                                         np.cumsum(z[:, :, None] * y0[:, None, :], axis=0)]))
    for i in np.flatnonzero(valid):
        stop = min(int(np.searchsorted(ends, i, side="left")), max(i - H, 0))
        n = int(count[stop])
        if n < 252:
            continue
        eligible = np.flatnonzero(valid[:stop])
        assert len(eligible) == n and np.all(ends[eligible] < i)
        assert i - eligible[-1] >= H  # at least a 21-session origin embargo
        origins.append(int(i))
        train_counts.append(n)
        train_last.append(int(eligible[-1]))
        means[i] = sy[stop] / n
        for name, x in inputs.items():
            sx, sxx, sxy = (m[stop] for m in moments[name])
            beta = np.linalg.lstsq(sxx - np.outer(sx, sx) / n,
                                   sxy - np.outer(sx, sy[stop]) / n, rcond=1e-12)[0]
            intercept = (sy[stop] - sx @ beta) / n
            predictions[name][i] = intercept + x[i] @ beta
            coefficients[name] = np.vstack([intercept, beta]).tolist()
            if n == 252 or i == np.flatnonzero(valid)[-1]:
                # Independent direct least-squares check of the cumulative implementation.
                design = np.column_stack([np.ones(n), x[eligible]])
                direct = np.r_[1.0, x[i]] @ np.linalg.lstsq(design, targets[eligible], rcond=None)[0]
                assert np.allclose(predictions[name][i], direct, rtol=1e-8, atol=1e-7)
    origins = np.asarray(origins, int)
    assert len(origins) > 100
    chosen, previous_end = [], -1
    for i in origins:
        if i >= previous_end:
            chosen.append(i)
            previous_end = ends[i]
    chosen = np.asarray(chosen, int)
    assert np.all(chosen[1:] >= ends[chosen[:-1]])
    report = {"sample": {"eligible_origins": int(valid.sum()), "forecast_days": int(len(origins)),
                          "first_forecast": str(g.index[origins[0]]), "last_forecast": str(g.index[origins[-1]]),
                          "first_training_days": train_counts[0], "last_training_days": train_counts[-1],
                          "minimum_origin_gap_sessions": int(min(origins - np.asarray(train_last))),
                          "nonoverlap_days": int(len(chosen))}, "last_fit_coefficients": coefficients}
    for j, name in enumerate(("forward_realized_volatility", "forward_volatility_premium")):
        losses = {model: (targets[:, j] - values[:, j]) ** 2
                  for model, values in {"historical_mean": means, **predictions}.items()}
        mse = {model: float(np.nanmean(v)) for model, v in losses.items()}
        models = {model: {"mse": value, "rmse": math.sqrt(value),
                           "oos_r2_vs_mean": 1 - value / mse["historical_mean"],
                           "mse_gain_vs_iv_pct": 100 * (1 - value / mse["iv_only"]),
                           "nonoverlap_mse": float(np.mean(losses[model][chosen]))}
                  for model, value in mse.items()}
        comparisons = {}
        for model in ("iv_plus_rv", "iv_minus_rv", "rv_only"):
            delta = losses["iv_only"] - losses[model]
            comparisons[model + "_vs_iv"] = {"hac25": hac_mean(delta, 25), "hac63": hac_mean(delta, 63)}
        report[name] = {"models": models, "paired_loss_improvement": comparisons}
    return report


def nw_t(x, lag=H + 4):
    x = np.asarray(x, float)
    x = x[np.isfinite(x)]
    n = len(x)
    m = x.mean()
    z = x - m
    v = z @ z / n
    for L in range(1, lag + 1):
        v += 2 * (1 - L / (lag + 1)) * (z[L:] @ z[:-L]) / n
    return m / math.sqrt(v / n), n


def p2(t):
    return math.erfc(abs(t) / math.sqrt(2))


def holm(ps, alpha=0.05):
    order = sorted(range(len(ps)), key=lambda i: ps[i])
    rej = [False] * len(ps)
    for rank, i in enumerate(order):
        if ps[i] <= alpha / (len(ps) - rank):
            rej[i] = True
        else:
            break
    return rej


frames, input_audit = {}, {}
for s in SYMS:
    yp, bp = os.path.join(D, "yahoo", f"{s}.csv"), os.path.join(D, "benchmark", f"{s}_inputs.csv")
    y = pd.read_csv(yp, index_col="date")
    raw = pd.read_csv(bp, index_col="date")
    assert y.index.is_unique and raw.index.is_unique and y.index.is_monotonic_increasing and raw.index.is_monotonic_increasing
    assert len(raw.index.difference(y.index)) == 0
    assert (y[["open", "high", "low", "close"]] > 0).all().all()
    assert (y.high >= y[["open", "close", "low"]].max(axis=1)).all()
    assert (y.low <= y[["open", "close", "high"]].min(axis=1)).all()
    assert (raw.total_rv.dropna() >= 0).all() and (raw.iv30.dropna() > 0).all()
    b = raw
    b = b[b.quality.isna()]
    comp = E.components(y.open.values, y.high.values, y.low.values, y.close.values, y.close.shift(1).values)
    c = pd.DataFrame(comp, index=y.index).iloc[1:]
    df = pd.DataFrame(index=c.index)
    for n in (5, 10, 20, 30):
        df[f"cc{n}"] = E.annualize(c.r.rolling(n).var(ddof=1))
    df["yz20"] = E.annualize(E.rolling_variances(c, 20)["yang_zhang"])
    bt = b.total_rv.reindex(c.index)
    fwd = bt[::-1].rolling(H, min_periods=H).mean()[::-1].shift(-1)          # mean over t+1..t+H, NaN if any day missing
    df["rvf"] = E.annualize(fwd)
    df["rvf_r2"] = E.annualize((c.r ** 2)[::-1].rolling(H, min_periods=H).mean()[::-1].shift(-1))
    df["rvf_cc"] = E.annualize(c.r[::-1].rolling(H, min_periods=H).var(ddof=1)[::-1].shift(-1))
    df["bench20"] = E.annualize(bt.rolling(20, min_periods=20).mean())
    df["iv"] = 100 * b.iv30.reindex(c.index)
    df = df.loc[START:]
    dates = pd.to_datetime(df.index)
    positions = np.arange(len(df))
    end30 = np.searchsorted(dates, dates + pd.Timedelta(days=30), side="right") - 1
    sessions30 = end30 - positions
    daily = bt.reindex(df.index)
    sums = np.r_[0.0, np.cumsum(daily.fillna(0.0))]
    counts = np.r_[0, np.cumsum(daily.notna())]
    complete30 = ((dates + pd.Timedelta(days=30) <= dates[-1]) & (sessions30 > 0)
                  & (counts[end30 + 1] - counts[positions + 1] == sessions30))
    mean30 = np.divide(sums[end30 + 1] - sums[positions + 1], sessions30,
                       out=np.full(len(df), np.nan), where=sessions30 > 0)
    df["rvf_cal30"] = E.annualize(np.where(complete30, mean30, np.nan))
    df["end21"], df["end_cal30"], df["cal30_sessions"] = positions + H, end30, sessions30
    df["rvf_cal"] = 100 * np.sqrt(np.where(complete30, np.maximum(sums[end30 + 1] - sums[positions + 1], 0) * 365 / 30, np.nan))
    df["n_cal"], df["has_bench"] = sessions30, daily.notna()
    df["fwd_end"] = [str(df.index[i + H]) if i + H < len(df) else None for i in positions]
    df["vrp_ante"] = df.iv - df.cc20
    df["vrp_post"] = df.iv - df.rvf
    df["var_post"] = df.iv ** 2 - df.rvf ** 2
    df["sym"] = s
    frames[s] = df
    clean_iv = df.iv.notna()
    clean_trail = clean_iv & df.cc20.notna()
    eligible = clean_trail & df.rvf.notna()
    yy = y.loc[raw.index.min():raw.index.max()].index
    gaps = yy.difference(raw.index)
    close_returns = np.log(y.close / y.close.shift(1))
    close_diff = (raw.ret - close_returns.reindex(raw.index)).dropna()
    input_audit[s] = {"files_sha256": {"benchmark": hashlib.sha256(open(bp, "rb").read()).hexdigest(),
                                       "yahoo": hashlib.sha256(open(yp, "rb").read()).hexdigest()},
                      "benchmark_rows": int(len(raw)), "benchmark_start": str(raw.index.min()),
                      "benchmark_end": str(raw.index.max()), "benchmark_missing_calendar_dates": list(gaps),
                      "quality_flags": {str(k): int(v) for k, v in raw.quality.dropna().value_counts().items()},
                      "missing_iv_rows": int(raw.iv30.isna().sum()), "clean_missing_iv_rows": int(b.iv30.isna().sum()),
                      "source_return_median_abs_difference": float(close_diff.abs().median()),
                      "source_return_max_abs_difference": float(close_diff.abs().max()),
                      "source_return_max_abs_difference_date": str(close_diff.abs().idxmax()),
                      "source_return_abs_difference_above_half_pct_point_rows": int((close_diff.abs() > 0.005).sum()),
                      "stages": {"candidate_session_origins": int(len(df)), "excluded_no_clean_iv": int((~clean_iv).sum()),
                                 "excluded_no_trailing20_after_iv": int((clean_iv & ~clean_trail).sum()),
                                 "excluded_no_complete_forward21_after_inputs": int((clean_trail & ~eligible).sum()),
                                 "eligible_main": int(eligible.sum())},
                      "eligible_start": str(df.index[eligible][0]), "eligible_end": str(df.index[eligible][-1])}
    assert sum(input_audit[s]["stages"][k] for k in ("excluded_no_clean_iv", "excluded_no_trailing20_after_iv",
                                                   "excluded_no_complete_forward21_after_inputs", "eligible_main")) == len(df)
    # Direct horizon/window checks at both boundaries and the middle of the retained sample.
    for i in np.flatnonzero(eligible)[[0, int(eligible.sum()) // 2, -1]]:
        assert i + H < len(df)
        assert np.isclose(df.rvf.iloc[i], E.annualize(daily.iloc[i + 1:i + H + 1].mean()))
        assert np.isclose(df.rvf_cc.iloc[i], E.annualize(c.r.reindex(df.index).iloc[i + 1:i + H + 1].var(ddof=1)))
    for i in np.flatnonzero(complete30)[[0, int(complete30.sum()) // 2, -1]]:
        window = (dates > dates[i]) & (dates <= dates[i] + pd.Timedelta(days=30))
        assert int(window.sum()) == sessions30[i] and int(np.flatnonzero(window)[-1]) == end30[i]
        assert np.isclose(df.rvf_cal30.iloc[i], E.annualize(daily[window].mean()))

all_ = pd.concat(frames.values())
ok = all_.dropna(subset=["iv", "cc20", "rvf"])
res = {"symbols": SYMS, "horizon": H, "start": str(ok.index.min()), "end": str(ok.index.max()), "days": int(len(ok))}

# ---- 1. size and frequency -----------------------------------------------------------------------------------
per = {}
for s, g in ok.groupby("sym"):
    vp = g.var_post
    pos, neg = vp[vp > 0].sum(), -vp[vp < 0].sum()
    t_m, _ = nw_t(g.vrp_post)
    per[s] = {"days": int(len(g)), "iv_mean": float(g.iv.mean()), "rvf_mean": float(g.rvf.mean()),
              "vrp_post_mean": float(g.vrp_post.mean()), "vrp_post_median": float(g.vrp_post.median()),
              "vrp_post_pos": float((g.vrp_post > 0).mean()), "vrp_post_t": float(t_m),
              "var_post_mean": float(vp.mean()), "var_post_p5": float(vp.quantile(0.05)), "var_post_min": float(vp.min()),
              "var_post_min_date": str(vp.idxmin()), "loss_share": float(neg / pos) if pos > 0 else None,
              "vrp_ante_mean": float(g.vrp_ante.mean()), "vrp_ante_pos": float((g.vrp_ante > 0).mean())}
res["per_symbol"] = per
for grp, syms in (("etf", ETFS), ("stock", SYMS[8:]), ("all", SYMS)):
    g = ok[ok.sym.isin(syms)]
    res[f"pooled_{grp}"] = {"vrp_post_mean": float(g.vrp_post.mean()), "vrp_post_median": float(g.vrp_post.median()),
                            "vrp_post_pos": float((g.vrp_post > 0).mean()), "var_post_mean": float(g.var_post.mean()),
                            "vrp_ante_mean": float(g.vrp_ante.mean()), "vrp_ante_pos": float((g.vrp_ante > 0).mean()),
                            "iv_over_rvf": float((g.iv > g.rvf).mean())}

# ---- 2. by year --------------------------------------------------------------------------------------------------
ok_y = ok.copy()
ok_y["year"] = ok_y.index.str[:4]
res["by_year"] = {yr: {"vrp_post_mean": float(g.vrp_post.mean()), "vrp_post_pos": float((g.vrp_post > 0).mean()),
                       "var_post_mean": float(g.var_post.mean()), "days": int(len(g))} for yr, g in ok_y.groupby("year")}

# ---- 3. does the trailing premium predict the forward one? ------------------------------------------------------
pred = {}
for s, g in ok.groupby("sym"):
    x, yv = g.vrp_ante.values, g.vrp_post.values
    beta = np.polyfit(x, yv, 1)
    resid = yv - (beta[1] + beta[0] * x)
    # Newey-West t for the slope via the moment (x - xbar) * resid
    xc = x - x.mean()
    tb, _ = nw_t(xc * resid / np.mean(xc * xc) + beta[0])  # t of beta computed on the influence function
    # Mincer-Zarnowitz in variance: rvf^2 = a + b iv^2
    mz = np.polyfit(g.iv.values ** 2, g.rvf.values ** 2, 1)
    # out-of-sample R^2: forecast vrp_post with an expanding regression on vrp_ante, using only targets known by t
    oos_num = oos_den = 0.0
    xs, ys = g.vrp_ante.values, g.vrp_post.values
    ends = g.fwd_end.values
    for i in range(252 + H, len(g)):
        xi, yi = xs[: i - H], ys[: i - H]
        assert ends[i - H - 1] < g.index[i], (s, g.index[i])
        bb = np.polyfit(xi, yi, 1)
        f_reg = bb[1] + bb[0] * xs[i]
        f_mean = yi.mean()
        oos_num += (ys[i] - f_reg) ** 2
        oos_den += (ys[i] - f_mean) ** 2
    pred[s] = {"corr": float(np.corrcoef(x, yv)[0, 1]), "slope": float(beta[0]), "slope_t": float(tb),
               "r2": float(np.corrcoef(x, yv)[0, 1] ** 2), "oos_r2": float(1 - oos_num / oos_den),
               "mz_intercept": float(mz[1]), "mz_slope": float(mz[0]),
               "iv_rv_corr": float(np.corrcoef(g.iv, g.rvf)[0, 1]), "trail_rv_corr": float(np.corrcoef(g.cc20, g.rvf)[0, 1])}
res["predict"] = pred

# binned calibration: deciles of the trailing premium (pooled) -> mean forward premium
q = pd.qcut(ok.vrp_ante, 10, labels=False)
res["deciles"] = [{"decile": int(k) + 1, "ante_mean": float(g.vrp_ante.mean()), "post_mean": float(g.vrp_post.mean()),
                   "post_pos": float((g.vrp_post > 0).mean()), "var_post_mean": float(g.var_post.mean()), "days": int(len(g))}
                  for k, g in ok.groupby(q)]

# ---- 4. signal: today's trailing premium vs its own past year, terciles -------------------------------------------
MINOBS = 189   # published v1.1 policy: at least 189 of the previous 252 sessions (75%); v1.0 required a complete window


def past_rank(w):
    """Share of the previous sessions' values below today's. Today (last element) is excluded from the reference set."""
    if not np.isfinite(w[-1]):
        return np.nan
    p = w[:-1][np.isfinite(w[:-1])]
    return float((p < w[-1]).mean()) if len(p) >= MINOBS else np.nan


def hilo(pr, v):
    keep = np.isfinite(pr) & np.isfinite(v)
    return float(v[keep & (pr >= 2 / 3)].mean() - v[keep & (pr < 1 / 3)].mean())


def signal_uncertainty(g, calendar, lag=25):
    """Centered conditional-mean influence score, retaining excluded calendar positions."""
    hi, lo = g.pr >= 2 / 3, g.pr < 1 / 3
    diff = float(g.var_post[hi].mean() - g.var_post[lo].mean())
    d = (np.where(hi, (g.var_post - g.var_post[hi].mean()) / hi.mean(), 0)
         - np.where(lo, (g.var_post - g.var_post[lo].mean()) / lo.mean(), 0) + diff)
    uncertainty = hac_mean(pd.Series(d, index=g.index).reindex(calendar), lag)
    assert np.isclose(uncertainty["mean"], diff)
    return uncertainty


rng = np.random.default_rng(20261110)
BLOCK, REPS = 2 * H, 2000
sig = {}
ps = []
recon = {}
for s in SYMS:
    g0 = frames[s].dropna(subset=["iv", "cc20", "rvf"]).copy()
    hist = frames[s].vrp_ante
    pr = hist.rolling(253, min_periods=MINOBS + 1).apply(past_rank, raw=True)   # rank vs the previous 252 sessions only
    g0["pr"] = pr.reindex(g0.index)
    g = g0.dropna(subset=["pr"])
    recon[s] = {"t4_days": int(len(g0)), "t4_mean": float(g0.var_post.mean()), "warmup_days": int(len(g0) - len(g)),
                "warmup_mean": float(g0[g0.pr.isna()].var_post.mean()) if len(g0) > len(g) else None,
                "t6_days": int(len(g)), "t6_mean": float(g.var_post.mean()), "t6_first": str(g.index.min())}
    hi, lo = g[g.pr >= 2 / 3], g[g.pr < 1 / 3]
    uncertainty = signal_uncertainty(g, frames[s].index)
    diff, se = uncertainty["mean"], uncertainty["se"]
    t = diff / se
    lag_t = {str(L): diff / signal_uncertainty(g, frames[s].index, L)["se"] for L in (10, 25, 50)}
    # Resample full calendar blocks, with missing outcomes still missing; never compress gaps.
    calendar = frames[s].loc[g.index.min():g.index.max()].index
    prv, vv = g.pr.reindex(calendar).values, g.var_post.reindex(calendar).values
    nb = int(math.ceil(len(calendar) / BLOCK))
    boots = []
    for _ in range(REPS):
        st = rng.integers(0, len(calendar) - BLOCK + 1, nb)
        ii = (st[:, None] + np.arange(BLOCK)[None, :]).ravel()[: len(calendar)]
        a, b_ = prv[ii], vv[ii]
        if (np.isfinite(b_) & (a >= 2 / 3)).any() and (np.isfinite(b_) & (a < 1 / 3)).any():
            boots.append(hilo(a, b_))
    # non-overlapping windows: every 21st session, all 21 starting offsets
    offs = [hilo(prv[o::H], vv[o::H]) for o in range(H) if (prv[o::H] >= 2 / 3).any() and (prv[o::H] < 1 / 3).any()]
    sig[s] = {"days": int(len(g)), "hi_mean": float(hi.var_post.mean()), "mid_mean": float(g[(g.pr >= 1 / 3) & (g.pr < 2 / 3)].var_post.mean()),
              "lo_mean": float(lo.var_post.mean()), "all_mean": float(g.var_post.mean()), "diff": diff,
              "hi_days": int(len(hi)), "lo_days": int(len(lo)),
              "ci_nw": uncertainty["ci95"], "ci95": uncertainty["ci95"], "hac_se": se,
              "ci_boot": [float(np.percentile(boots, 2.5)), float(np.percentile(boots, 97.5))],
              "t_by_lag": lag_t, "nonoverlap_pos": int(sum(x > 0 for x in offs)), "nonoverlap_n": len(offs),
              "nonoverlap_min": float(min(offs)), "nonoverlap_max": float(max(offs)), "nonoverlap_median": float(np.median(offs)),
              "vol_diff": float(hi.vrp_post.mean() - lo.vrp_post.mean()),
              "t": float(t), "p": float(p2(t)),
              "pos_ante_mean": float(g[g.vrp_ante > 0].var_post.mean()), "neg_ante_mean": float(g[g.vrp_ante <= 0].var_post.mean()),
              "pos_ante_share": float((g.vrp_ante > 0).mean())}
    ps.append(sig[s]["p"])
rej = holm(ps)
for s, rj in zip(SYMS, rej):
    sig[s]["holm"] = bool(rj)
res["signal"] = sig
res["signal_summary"] = {"hi_gt_lo": int(sum(sig[s]["diff"] > 0 for s in SYMS)), "sig_raw": int(sum(sig[s]["p"] < 0.05 for s in SYMS)),
                         "sig_holm": int(sum(sig[s]["holm"] for s in SYMS)),
                         "sig_holm_pos": int(sum(sig[s]["holm"] and sig[s]["diff"] > 0 for s in SYMS))}

res["recon"] = recon
res["signal_summary"].update({"boot_excl_zero": int(sum(sig[s]["ci_boot"][0] > 0 or sig[s]["ci_boot"][1] < 0 for s in SYMS)),
                              "days": int(sum(sig[s]["days"] for s in SYMS)), "minobs": MINOBS, "block": BLOCK, "reps": REPS})

# Preserve the earlier complete-window specification as a declared coverage sensitivity.
strict_sig = {}
for s in SYMS:
    hist = frames[s].vrp_ante
    strict_pr = hist.rolling(253, min_periods=253).apply(lambda w: (w[:-1] < w[-1]).mean(), raw=True)
    g = frames[s].dropna(subset=["iv", "cc20", "rvf"]).copy()
    g["pr"] = strict_pr.reindex(g.index)
    g = g.dropna(subset=["pr"])
    u = signal_uncertainty(g, frames[s].index)
    t = u["mean"] / u["se"]
    strict_sig[s] = {"days": len(g), "diff": u["mean"], "ci95": u["ci95"], "t": t, "p": p2(t)}
strict_rej = holm([strict_sig[s]["p"] for s in SYMS])
for s, reject in zip(SYMS, strict_rej):
    strict_sig[s]["holm"] = bool(reject)
res["signal_complete_window"] = {"per_symbol": strict_sig, "days": sum(g["days"] for g in strict_sig.values()),
    "sig_raw": sum(g["p"] < .05 for g in strict_sig.values()), "sig_holm": sum(strict_rej),
    "note": "Coverage sensitivity requiring all 252 prior session readings; canonical signal permits at least 189 of those 252, following published version 1.1."}

# ---- 4b. does trailing realized variance add information beyond implied variance? ----------------------------
# Expanding out-of-sample forecasts of the next 21 sessions' realized variance (vol points squared). Training labels
# are matured (asserted). Models: IV^2 only, trailing variance only, both. Losses: MSE and QLIKE. Clark-West (2007)
# tests the nested pair IV^2 vs IV^2 + trailing variance.
def qlike(y, f):
    r = y / f
    return float(np.mean(r - np.log(r) - 1))


inc = {}
for s in SYMS:
    g = frames[s].dropna(subset=["iv", "cc20", "rvf", "bench20"])
    y = g.rvf.values ** 2
    X = {"iv": np.c_[np.ones(len(g)), g.iv.values ** 2], "rv": np.c_[np.ones(len(g)), g.cc20.values ** 2],
         "both": np.c_[np.ones(len(g)), g.iv.values ** 2, g.cc20.values ** 2],
         "both_intraday": np.c_[np.ones(len(g)), g.iv.values ** 2, g.bench20.values ** 2]}
    ends = g.fwd_end.values
    F = {k: [] for k in X}
    F["mean"], F["iv_raw"], yy, floors, forecast_dates = [], [], [], [], []
    for i in range(252 + H, len(g)):
        assert ends[i - H - 1] < g.index[i]
        tr = slice(0, i - H)
        for k, x in X.items():
            beta = np.linalg.lstsq(x[tr], y[tr], rcond=None)[0]
            F[k].append(float(x[i] @ beta))
        F["mean"].append(float(y[tr].mean()))
        F["iv_raw"].append(float(g.iv.values[i] ** 2))
        yy.append(y[i])
        floors.append(0.05 * float(np.median(y[tr])))
        forecast_dates.append(g.index[i])
    yy = np.array(yy)
    floors = np.asarray(floors)
    assert np.all(floors > 0)
    out = {"days": int(len(yy)), "first_forecast": forecast_dates[0], "last_forecast": forecast_dates[-1],
           "qlike_floor_method": "At each forecast origin, 5% of the median matured training variance; evaluation targets are never used",
           "qlike_floor_min": float(floors.min()), "qlike_floor_max": float(floors.max())}
    for k, f in F.items():
        f = np.array(f)
        out[f"mse_{k}"] = float(np.mean((yy - f) ** 2))
        out[f"qlike_{k}"] = qlike(yy, np.maximum(f, floors))
        out[f"qlike_clipped_{k}"] = int((f < floors).sum())
    f1, f2 = np.array(F["iv"]), np.array(F["both"])
    cw = (yy - f1) ** 2 - ((yy - f2) ** 2 - (f1 - f2) ** 2)
    cw_u = hac_mean(pd.Series(cw, index=forecast_dates).reindex(frames[s].index))
    tcw = cw_u["mean"] / cw_u["se"]
    f2b = np.array(F["both_intraday"])
    cwb = (yy - f1) ** 2 - ((yy - f2b) ** 2 - (f1 - f2b) ** 2)
    cwb_u = hac_mean(pd.Series(cwb, index=forecast_dates).reindex(frames[s].index))
    tcwb = cwb_u["mean"] / cwb_u["se"]
    out.update({"cw_t": float(tcw), "cw_p": float(0.5 * math.erfc(tcw / math.sqrt(2))),
                "cw_t_intraday": float(tcwb), "cw_p_intraday": float(0.5 * math.erfc(tcwb / math.sqrt(2))),
                "mse_ratio_both_iv": out["mse_both"] / out["mse_iv"], "mse_ratio_rv_iv": out["mse_rv"] / out["mse_iv"],
                "mse_ratio_bothintra_iv": out["mse_both_intraday"] / out["mse_iv"]})
    inc[s] = out
cw_rej = holm([inc[s]["cw_p"] for s in SYMS])
cwb_rej = holm([inc[s]["cw_p_intraday"] for s in SYMS])
for s, a, b_ in zip(SYMS, cw_rej, cwb_rej):
    inc[s]["cw_holm"], inc[s]["cw_holm_intraday"] = bool(a), bool(b_)
res["incremental"] = inc
res["incremental_summary"] = {
    "iv_beats_rv_mse": int(sum(inc[s]["mse_iv"] < inc[s]["mse_rv"] for s in SYMS)),
    "iv_beats_rv_qlike": int(sum(inc[s]["qlike_iv"] < inc[s]["qlike_rv"] for s in SYMS)),
    "both_beats_iv_mse": int(sum(inc[s]["mse_both"] < inc[s]["mse_iv"] for s in SYMS)),
    "both_beats_iv_qlike": int(sum(inc[s]["qlike_both"] < inc[s]["qlike_iv"] for s in SYMS)),
    "cw_sig_raw": int(sum(inc[s]["cw_p"] < 0.05 for s in SYMS)), "cw_sig_holm": int(sum(cw_rej)),
    "cwb_sig_raw": int(sum(inc[s]["cw_p_intraday"] < 0.05 for s in SYMS)), "cwb_sig_holm": int(sum(cwb_rej)),
    "bothintra_beats_iv_mse": int(sum(inc[s]["mse_both_intraday"] < inc[s]["mse_iv"] for s in SYMS)),
    "iv_beats_mean_mse": int(sum(inc[s]["mse_iv"] < inc[s]["mse_mean"] for s in SYMS))}

# ---- 4c. matched estimators on common dates --------------------------------------------------------------------
# The headline pairs a close-to-close trailing window with a 5-minute forward window. Repeat with the same estimator
# on both sides, on the dates where all four series exist.
def oos_r2(x, yv, ends, index):
    num = den = 0.0
    for i in range(252 + H, len(x)):
        assert ends[i - H - 1] < index[i]
        bb = np.polyfit(x[: i - H], yv[: i - H], 1)
        num += (yv[i] - (bb[1] + bb[0] * x[i])) ** 2
        den += (yv[i] - yv[: i - H].mean()) ** 2
    return 1 - num / den


match = {}
SETUPS = {"mixed": ("cc20", "rvf"), "close_both": ("cc20", "rvf_cc"), "intraday_both": ("bench20", "rvf")}
common = {s: frames[s].dropna(subset=["iv", "cc20", "rvf", "rvf_cc", "bench20"]) for s in SYMS}
for name, (tr_, fw_) in SETUPS.items():
    rows = {}
    for s in SYMS:
        g = common[s]
        a, p_ = (g.iv - g[tr_]).values, (g.iv - g[fw_]).values
        rows[s] = {"corr": float(np.corrcoef(a, p_)[0, 1]), "oos_r2": float(oos_r2(a, p_, g.fwd_end.values, g.index)),
                   "post_mean": float(p_.mean()), "post_pos": float((p_ > 0).mean())}
    pooled = pd.concat([pd.DataFrame({"p": (common[s].iv - common[s][fw_])}) for s in SYMS])
    match[name] = {"per_symbol": rows, "oos_neg": int(sum(rows[s]["oos_r2"] < 0 for s in SYMS)),
                   "corr_min": float(min(rows[s]["corr"] for s in SYMS)), "corr_max": float(max(rows[s]["corr"] for s in SYMS)),
                   "oos_min": float(min(rows[s]["oos_r2"] for s in SYMS)), "oos_max": float(max(rows[s]["oos_r2"] for s in SYMS)),
                   "post_mean": float(pooled.p.mean()), "post_median": float(pooled.p.median()), "post_pos": float((pooled.p > 0).mean()),
                   "days": int(len(pooled))}
res["matched"] = match

# ---- 4d. horizon: 21 sessions versus the exact 30 calendar days -------------------------------------------------
hz = ok.dropna(subset=["rvf_cal"])
res["horizon"] = {"days": int(len(hz)), "sessions_min": int(hz.n_cal.min()), "sessions_max": int(hz.n_cal.max()),
                  "sessions_mean": float(hz.n_cal.mean()),
                  "s21": {"mean": float(hz.vrp_post.mean()), "median": float(hz.vrp_post.median()), "pos": float((hz.vrp_post > 0).mean())},
                  "cal30": {"mean": float((hz.iv - hz.rvf_cal).mean()), "median": float((hz.iv - hz.rvf_cal).median()),
                            "pos": float((hz.iv > hz.rvf_cal).mean())},
                  "median_abs_diff": float((hz.rvf - hz.rvf_cal).abs().median()), "corr": float(np.corrcoef(hz.rvf, hz.rvf_cal)[0, 1]),
                  "var_mean_s21": float(hz.var_post.mean()), "var_mean_cal30": float((hz.iv ** 2 - hz.rvf_cal ** 2).mean()),
                  "var_neg_syms_cal30": int(sum(((g.iv ** 2 - g.rvf_cal ** 2).mean() < 0) for _, g in hz.groupby("sym")))}

# ---- 4e. inference for the pooled mean: Newey-West by lag, calendar-block bootstrap, equal weights ----------------
pooled_calendar = frames[SYMS[0]].loc[ok.index.min():ok.index.max()].index
wide = ok.pivot_table(index=ok.index, columns="sym", values="vrp_post").reindex(pooled_calendar)
widev = ok.pivot_table(index=ok.index, columns="sym", values="var_post").reindex(pooled_calendar)
nd = len(wide)
nbk = int(math.ceil(nd / BLOCK))
bm, bv, bp = [], [], []
for _ in range(REPS):
    st = rng.integers(0, nd - BLOCK + 1, nbk)
    ii = (st[:, None] + np.arange(BLOCK)[None, :]).ravel()[:nd]
    w_, v_ = wide.values[ii], widev.values[ii]
    bm.append(np.nanmean(w_)); bv.append(np.nanmean(v_)); bp.append(np.nanmean(np.where(np.isfinite(w_), w_ > 0, np.nan)))
daily_mean = wide.mean(axis=1)
eq = {s: per[s]["vrp_post_mean"] for s in SYMS}
loo = {s: float(ok[ok.sym != s].vrp_post.mean()) for s in SYMS}
res["inference"] = {"dates": int(nd),
                    "mean_ci_boot": [float(np.percentile(bm, 2.5)), float(np.percentile(bm, 97.5))],
                    "var_ci_boot": [float(np.percentile(bv, 2.5)), float(np.percentile(bv, 97.5))],
                    "pos_ci_boot": [float(np.percentile(bp, 2.5)), float(np.percentile(bp, 97.5))],
                    "daily_mean_t_by_lag": {str(L): hac_mean(daily_mean.values, L)["mean"] / hac_mean(daily_mean.values, L)["se"] for L in (10, 25, 50)},
                    "equal_weight_mean": float(np.mean(list(eq.values()))), "equal_weight_var_mean": float(np.mean([per[s]["var_post_mean"] for s in SYMS])),
                    "loo_min": float(min(loo.values())), "loo_max": float(max(loo.values())),
                    "loo_min_sym": min(loo, key=loo.get), "loo_max_sym": max(loo, key=loo.get)}

# ---- 4f. within-symbol deciles (timing) versus pooled deciles (mixes symbols) -------------------------------------
wq = ok.groupby("sym").vrp_ante.transform(lambda v: pd.qcut(v, 10, labels=False))
res["deciles_within"] = [{"decile": int(k) + 1, "ante_mean": float(g.vrp_ante.mean()), "post_mean": float(g.vrp_post.mean()),
                          "post_pos": float((g.vrp_post > 0).mean()), "days": int(len(g))} for k, g in ok.groupby(wq)]

# ---- 4g. coverage and episodes ------------------------------------------------------------------------------------
cov = {}
for s in SYMS:
    f = frames[s]
    cov[s] = {"sessions": int(len(f)), "with_iv": int(f.iv.notna().sum()), "with_trailing": int(f.dropna(subset=["iv", "cc20"]).shape[0]),
              "used": int(f.dropna(subset=["iv", "cc20", "rvf"]).shape[0]), "flagged_sessions": int((~f.has_bench).sum())}
res["coverage"] = cov
spy_v = ok[ok.sym == "SPY"].var_post
yr = spy_v.groupby(spy_v.index.str[:4]).sum()
res["episodes"] = {"spy_sum_by_year": {k: float(v) for k, v in yr.items()}, "spy_total": float(spy_v.sum()),
                   "spy_worst_year": str(yr.idxmin()), "spy_worst_year_sum": float(yr.min()),
                   "spy_mean_ex_worst_year": float(spy_v[spy_v.index.str[:4] != yr.idxmin()].mean()),
                   "pooled_mean_ex_2020": float(ok[ok.index.str[:4] != "2020"].vrp_post.mean()),
                   "pooled_var_mean_ex_2020": float(ok[ok.index.str[:4] != "2020"].var_post.mean()),
                   "var_neg_syms_ex_2020": int(sum(g[g.index.str[:4] != "2020"].var_post.mean() < 0 for _, g in ok.groupby("sym")))}

# ---- 5. choices: RV window and estimator for the trailing premium ------------------------------------------------
ch = {}
for col in ("cc5", "cc10", "cc20", "cc30", "yz20"):
    ante = ok.iv - ok[col]
    ch[col] = {"mean": float(ante.mean()), "pos": float((ante > 0).mean()), "corr_post": float(np.corrcoef(ante, ok.vrp_post)[0, 1]),
               "mad_vs_cc20": float((ante - ok.vrp_ante).abs().median())}
ch["fwd_r2_vs_bench"] = {"median_abs_diff": float((ok.rvf_r2 - ok.rvf).abs().median()), "corr": float(np.corrcoef(ok.rvf_r2.dropna(), ok.rvf[ok.rvf_r2.notna()])[0, 1])}
res["choices"] = ch

# ---- SPY distribution of the variance premium (histogram for the chart) ------------------------------------------
spy = ok[ok.sym == "SPY"].var_post
edges = list(range(-1500, 401, 50))
hist, _ = np.histogram(spy.clip(-1499, 399), bins=edges)
res["spy_hist"] = {"edges": edges, "counts": [int(x) for x in hist], "median": float(spy.median()), "mean": float(spy.mean())}

# ---- 6. bounded validation: common forecasting task, input timing and measurement sensitivity -------------------
validation = {"method": {
    "targets": {"forward_realized_volatility": "Annualized realized volatility of the next 21 sessions, in percent volatility",
                "forward_volatility_premium": "Current 30-day ATM IV minus that forward realized volatility, in vol points"},
    "models": {"iv_only": "intercept + beta_IV * IV", "rv_only": "intercept + beta_RV * trailing RV20",
               "iv_plus_rv": "intercept + beta_IV * IV + beta_RV * trailing RV20",
               "iv_minus_rv": "intercept + beta * (IV - trailing RV20): equal and opposite input coefficients",
               "historical_mean": "Expanding average of the same fully observed training outcomes"},
    "estimation": "Separate per-symbol expanding OLS with an intercept; at least 252 eligible training outcomes; no clipping, hyperparameter tuning or extra signals",
    "timing": "After-close inputs only. Train only on origins at least 21 sessions earlier and whose complete forward windows ended strictly before the forecast date. Sessions retain their original calendar positions even across excluded observations.",
    "uncertainty": "Pointwise 95% normal-approximation Newey-West/Bartlett intervals with 25 and 63 session lags on paired squared-error differences. Positive means lower observed MSE than IV-only. Conditional on the fitted forecast path; not formal nested-model superiority tests, simultaneous intervals, or independent confirmation.",
    "nonoverlap": "Greedy first eligible forecast origin, then the first origin at or after its previous outcome window ends; one predetermined phase, descriptive MSE only",
    "signal_ranking": "Canonical version-1.1 policy: compare the current finite trailing spread with available finite readings in the preceding 252-session window, requiring at least 189 prior readings (75%). Missing prior readings are omitted from the rank, but calendar positions are retained. Scoring additionally requires a complete clean forward outcome. The earlier all-252-readings policy is retained separately in signal_complete_window as a declared coverage sensitivity.",
    "signal_uncertainty": "Centered difference-of-conditional-means influence scores with session-grid HAC25 (10/50-lag checks), 2000 moving-block bootstrap replications of 42 original sessions, and all 21 fixed phases of disjoint forward windows. No skipped dates are compressed into adjacent sessions.",
    "variance_forecasts": "The published v1.1 incremental variance comparison is retained: expanding OLS in squared-volatility units on common eligible dates, IV-only, historical daily variance, both, and both with historical intraday variance; every training label must mature before its forecast date. MSE uses unmodified forecasts; QLIKE clips each forecast at 5% of the median matured TRAINING variance at that origin. Clark-West one-sided adjusted-loss tests use session-grid HAC25; Holm is applied separately to the 11 daily-input and 11 intraday-input tests. These are distinct from the new volatility-unit models and their descriptive intervals.",
    "calendar_annualization": "horizon.cal30 retains the published v1.1 calendar-time measure: sum of future daily variances over 30 calendar days, multiplied by 365/30 before taking the square root. validation.sensitivity.calendar30 instead holds the 252-day mean-daily-variance annualization fixed to isolate window membership. Their different numbers are different annualization conventions, not duplicate estimates of an identical measure.",
    "pooled_inference": "Pooled mean/frequency bootstrap resamples 42-session blocks of the entire symbol matrix together, preserving cross-symbol shocks and the observed missingness pattern. Equal-symbol means, leave-one-symbol-out values, within-symbol deciles and episode exclusions are descriptive checks. Daily-mean HAC t statistics refer to equal-date means, not the observation-weighted pooled mean.",
    "sensitivity": "All four measurement specifications use the same complete origins and conservative training endpoints: the later of the 21-session and 30-calendar-day endpoints. Only the stated realized-vol measure/window changes.",
    "calendar30": "Observed sessions strictly after the origin and on or before origin + 30 calendar days; 252-day annualization retained. Exclude horizons ending beyond the supplied data. This checks the 21-session approximation; it is not a traded variance-swap replication.",
    "input_limits": "Checks validate supplied file structure, calculations, quality exclusions and session alignment. They do not independently verify vendor prices, market-close conventions, historical data vintages or redistribution rights. Yahoo split adjustment is retrospective; return-based measures are scale invariant, but vendor revision history is not reconstructed.",
    "source_comparison": "The archive's 16:00 midpoint returns and Yahoo daily closing-price returns are compared in the input audit, in decimal log-return units, not assumed identical. Observed differences are retained and are not automatically labeled errors; matched-estimator/source specifications address this measurement choice.",
    "interpretation": "Exploratory measurement/forecast validation on the existing sample, in volatility units. No option trades, transaction costs, strategy returns, or model-free implied variance are tested.",
    "sources": ["https://www.federalreserve.gov/pubs/feds/2011/201145/index.html",
                "https://doi.org/10.2307/1913610",
                "https://www.kansascityfed.org/documents/5368/pdf-RWP05-05.pdf"]},
    "input_audit": input_audit, "descriptive": {}, "forecast": {}, "sensitivity": {}}

specs = {"original": ("cc20", "rvf"), "matched_daily": ("cc20", "rvf_cc"),
         "matched_intraday": ("bench20", "rvf"), "calendar30": ("cc20", "rvf_cal30")}
for s in SYMS:
    g = frames[s]
    mask = g[["iv", "cc20", "rvf"]].notna().all(axis=1)
    validation["descriptive"][s] = {"days": int(mask.sum()), "start": str(g.index[mask][0]), "end": str(g.index[mask][-1]),
        "forward_premium_hac25": hac_mean(g.vrp_post.where(mask), 25),
        "forward_premium_hac63": hac_mean(g.vrp_post.where(mask), 63),
        "variance_spread_hac25": hac_mean(g.var_post.where(mask), 25),
        "variance_spread_hac63": hac_mean(g.var_post.where(mask), 63)}
    validation["forecast"][s] = forecast_comparison(g)
    input_audit[s]["sample_reconciliation"] = {"main_descriptive": int(mask.sum()),
        "signal": res["signal"][s]["days"],
        "new_forecast": validation["forecast"][s]["sample"]["forecast_days"]}
    common = g[["iv", "cc20", "rvf", "rvf_cc", "bench20", "rvf_cal30"]].notna().all(axis=1)
    common_g = g.copy()
    common_g.loc[~common, "iv"] = np.nan
    common_g["end_common"] = np.maximum(g.end21, g.end_cal30)
    counts30 = g.loc[common, "cal30_sessions"].value_counts().sort_index()
    sensitivity = {"common_days": int(common.sum()), "start": str(g.index[common][0]), "end": str(g.index[common][-1]),
                   "calendar30_session_counts": {str(int(k)): int(v) for k, v in counts30.items()},
                   "calendar30_median_abs_rv_change": float((g.rvf_cal30 - g.rvf)[common].abs().median()),
                   "daily_median_abs_rv_change": float((g.rvf_cc - g.rvf)[common].abs().median()),
                   "specifications": {}}
    forecast_samples = []
    for spec, (trail, fwd_col) in specs.items():
        q = forecast_comparison(common_g, trail, fwd_col, "end_common")
        forecast_samples.append((q["sample"]["forecast_days"], q["sample"]["first_forecast"], q["sample"]["last_forecast"]))
        vp = (g.iv - g[fwd_col])[common]
        sensitivity["specifications"][spec] = {"trailing_estimator": trail, "forward_estimator": fwd_col,
            "forward_premium_mean": float(vp.mean()), "forward_premium_positive_share": float((vp > 0).mean()),
            "variance_spread_mean": float((g.iv ** 2 - g[fwd_col] ** 2)[common].mean()),
            "trailing_forward_premium_corr": float((g.iv - g[trail])[common].corr(vp)),
            "forecast_sample": q["sample"],
            "premium_spread_oos_r2": q["forward_volatility_premium"]["models"]["iv_minus_rv"]["oos_r2_vs_mean"],
            "realized_vol_models": q["forward_realized_volatility"]["models"],
            "joint_vs_iv_hac25": q["forward_realized_volatility"]["paired_loss_improvement"]["iv_plus_rv_vs_iv"]["hac25"]}
    assert len(set(forecast_samples)) == 1
    validation["sensitivity"][s] = sensitivity

validation["sample_reconciliation"] = {"main_descriptive_symbol_days": res["days"],
    "signal_symbol_days": sum(res["signal"][s]["days"] for s in SYMS),
    "forecast_symbol_days": sum(validation["forecast"][s]["sample"]["forecast_days"] for s in SYMS),
    "sensitivity_common_symbol_days": sum(validation["sensitivity"][s]["common_days"] for s in SYMS),
    "signal_prior_sessions": 252, "signal_minimum_prior_readings": MINOBS,
    "signal_rank_window": "Current reading plus the preceding 252-session window; current reading and at least 189 prior readings required",
    "corrections": ["Version 1.0 compared the current reading with 251 prior sessions and required a complete window. Published version 1.1 adopted a 252-prior-session window with at least 189 available prior readings. Version 1.2 preserves that canonical policy and records the stricter complete-window variant separately.",
                    "Version 1.2 uses centered conditional-mean signal influence scores, original-session HAC/bootstrap/offset grids and training-only QLIKE floors. The matched daily estimator now uses centered sample variance on both sides; published v1.1 used uncentered forward squared returns."],
    "note": "Symbol-days are overlapping observations, not independent trials. Forecast samples exclude initial training history; signal samples require a current reading and at least 189 available readings in the preceding 252 sessions. Four symbols begin in 2018, not 2017."}
validation["forecast_summary"] = {}
for target in ("forward_realized_volatility", "forward_volatility_premium"):
    rows = [validation["forecast"][s][target] for s in SYMS]
    validation["forecast_summary"][target] = {
        "joint_lower_mse_than_iv_symbols": sum(r["models"]["iv_plus_rv"]["mse"] < r["models"]["iv_only"]["mse"] for r in rows),
        "joint_vs_iv_ci25_above_zero_symbols": sum(r["paired_loss_improvement"]["iv_plus_rv_vs_iv"]["hac25"]["ci95"][0] > 0 for r in rows),
        "joint_vs_iv_ci63_above_zero_symbols": sum(r["paired_loss_improvement"]["iv_plus_rv_vs_iv"]["hac63"]["ci95"][0] > 0 for r in rows),
        "spread_positive_oos_r2_symbols": sum(r["models"]["iv_minus_rv"]["oos_r2_vs_mean"] > 0 for r in rows),
        "joint_mse_gain_vs_iv_pct_range": [min(r["models"]["iv_plus_rv"]["mse_gain_vs_iv_pct"] for r in rows),
                                           max(r["models"]["iv_plus_rv"]["mse_gain_vs_iv_pct"] for r in rows)]}
res["validation"] = validation

json.dump(res, open(os.path.join(OUT, "vrp_results.json"), "w"), indent=1, allow_nan=False)
print("days", res["days"], res["start"], res["end"])
print("pooled", {k: round(v, 3) for k, v in res["pooled_all"].items()})
for s in SYMS:
    p, pr_, sg = per[s], pred[s], sig[s]
    print(f"{s:5s} post {p['vrp_post_mean']:+5.2f} pos {p['vrp_post_pos']:.2f} var {p['var_post_mean']:+7.1f} p5 {p['var_post_p5']:+8.1f} loss/gain {p['loss_share']:.2f} | corr {pr_['corr']:+.2f} slope {pr_['slope']:+.2f} t {pr_['slope_t']:+.1f} oosR2 {pr_['oos_r2']:+.3f} MZ b {pr_['mz_slope']:.2f} | sig hi {sg['hi_mean']:+7.1f} lo {sg['lo_mean']:+7.1f} t {sg['t']:+.1f} holm {sg['holm']}")
print("years", {k: (round(v['vrp_post_mean'], 2), round(v['vrp_post_pos'], 2)) for k, v in res["by_year"].items()})
print("deciles", [(d['decile'], round(d['ante_mean'], 1), round(d['post_mean'], 2)) for d in res["deciles"]])
print("signal", res["signal_summary"])
print("choices", {k: {kk: round(vv, 3) for kk, vv in v.items()} for k, v in ch.items()})