spb/wp8_uqo

1"""Volatility modeling for prediction market prices."""23from __future__ import annotations45from dataclasses import dataclass67import numpy as np8import pandas as pd9from arch import arch_model101112@dataclass13class GARCHResult:14    model_type: str15    params: pd.Series16    std_errors: pd.Series17    log_likelihood: float18    aic: float19    bic: float20    conditional_volatility: pd.Series212223@dataclass24class HARResult:25    coefficients: pd.Series26    std_errors: pd.Series27    t_stats: pd.Series28    r_squared: float29    adj_r_squared: float303132def compute_realized_volatility_measures(33    returns: pd.Series,34    windows: list[int] | None = None,35) -> pd.DataFrame:36    """Compute realized volatility, bipower variation, and jump component."""37    if windows is None:38        windows = [1, 5, 22]39    result = pd.DataFrame(index=returns.index)4041    for w in windows:42        result[f"rv_{w}"] = (returns ** 2).rolling(w).sum().apply(np.sqrt)4344    abs_ret = returns.abs()45    result["bv_1"] = (abs_ret * abs_ret.shift(1)).rolling(5).mean() * (np.pi / 2)46    result["jump_1"] = np.maximum(result["rv_5"] ** 2 - result["bv_1"], 0)4748    return result495051def estimate_har_rv(52    rv_series: pd.DataFrame,53    exog: pd.DataFrame | None = None,54) -> HARResult:55    """HAR-RV model (Corsi 2009).5657    RV_{t+1} = c + b_d * RV_t + b_w * RV_t^(w) + b_m * RV_t^(m) + e58    """59    import statsmodels.api as sm6061    if "rv_1" not in rv_series.columns:62        raise ValueError("rv_series must contain 'rv_1' column")6364    rv_d = rv_series["rv_1"]65    rv_w = rv_series.get("rv_5", rv_d.rolling(5).mean())66    rv_m = rv_series.get("rv_22", rv_d.rolling(22).mean())6768    y = rv_d.shift(-1)69    X = pd.DataFrame({"rv_d": rv_d, "rv_w": rv_w, "rv_m": rv_m})70    if exog is not None:71        for col in exog.columns:72            X[col] = exog[col]7374    X = sm.add_constant(X)75    valid = pd.concat([y, X], axis=1).dropna()76    if len(valid) < 10:77        return HARResult(78            coefficients=pd.Series(), std_errors=pd.Series(),79            t_stats=pd.Series(), r_squared=np.nan, adj_r_squared=np.nan,80        )8182    y_clean = valid.iloc[:, 0]83    X_clean = valid.iloc[:, 1:]8485    model = sm.OLS(y_clean, X_clean).fit(cov_type="HAC", cov_kwds={"maxlags": 5})8687    return HARResult(88        coefficients=model.params,89        std_errors=model.bse,90        t_stats=model.tvalues,91        r_squared=model.rsquared,92        adj_r_squared=model.rsquared_adj,93    )949596def estimate_garch(97    returns: pd.Series,98    model_type: str = "GARCH",99    p: int = 1,100    q: int = 1,101    dist: str = "normal",102) -> GARCHResult:103    """Estimate GARCH-family model using the arch package."""104    clean = returns.dropna()105    scaled = clean * 100106107    vol_spec = model_type108    if model_type == "GJR-GARCH":109        vol_spec = "GARCH"110        o = 1111    else:112        o = 0113114    am = arch_model(scaled, vol=vol_spec, p=p, o=o, q=q, dist=dist, mean="Constant")115    result = am.fit(disp="off")116117    cond_vol = result.conditional_volatility / 100118    cond_vol.index = clean.index119120    return GARCHResult(121        model_type=model_type,122        params=result.params,123        std_errors=result.std_err,124        log_likelihood=result.loglikelihood,125        aic=result.aic,126        bic=result.bic,127        conditional_volatility=cond_vol,128    )129