Bybit全銘柄のボラティリティをPythonで推定し
銘柄をスクリーニングする
Yang-Zhang・HAR-RV

"""vol_estimators.py — OHLCを使う高効率ボラティリティ推定量
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

参考:
  Parkinson (1980), Garman & Klass (1980),
  Rogers & Satchell (1991), Yang & Zhang (2000)

すべて「年率化されたボラティリティ」を返す。
bars_per_year は足の種類で変える（1時間足なら 24*365 = 8760）。
"""
from __future__ import annotations

import numpy as np
import pandas as pd

BARS_PER_YEAR = {
    "1": 525_600, "3": 175_200, "5": 105_120, "15": 35_040,
    "30": 17_520, "60": 8_760, "120": 4_380, "240": 2_190,
    "360": 1_460, "720": 730, "D": 365, "W": 52,
}


def _ann(var: pd.Series, bars_per_year: int) -> pd.Series:
    """分散を年率ボラティリティに変換する。"""
    return np.sqrt(var.clip(lower=0)) * np.sqrt(bars_per_year)


def close_to_close(bars: pd.DataFrame, window: int,
                   bars_per_year: int = 8_760) -> pd.Series:
    """終値の対数リターンの標準偏差。最も単純で最も非効率。"""
    r = np.log(bars["close"] / bars["close"].shift(1))
    return r.rolling(window, min_periods=max(window // 2, 5)).std() \
        * np.sqrt(bars_per_year)


def parkinson(bars: pd.DataFrame, window: int,
              bars_per_year: int = 8_760) -> pd.Series:
    """高値・安値のレンジを使う。ギャップとドリフトは無視される。

    σ² = (1 / (4 ln2)) * mean[ (ln(H/L))² ]
    """
    hl = np.log(bars["high"] / bars["low"]) ** 2
    var = hl.rolling(window, min_periods=max(window // 2, 5)).mean() \
        / (4 * np.log(2))
    return _ann(var, bars_per_year)


def garman_klass(bars: pd.DataFrame, window: int,
                 bars_per_year: int = 8_760) -> pd.Series:
    """始値・高値・安値・終値を使う。

    σ² = mean[ 0.5*(ln(H/L))² − (2ln2 − 1)*(ln(C/O))² ]
    """
    hl = np.log(bars["high"] / bars["low"]) ** 2
    co = np.log(bars["close"] / bars["open"]) ** 2
    term = 0.5 * hl - (2 * np.log(2) - 1) * co
    var = term.rolling(window, min_periods=max(window // 2, 5)).mean()
    return _ann(var, bars_per_year)


def rogers_satchell(bars: pd.DataFrame, window: int,
                    bars_per_year: int = 8_760) -> pd.Series:
    """ドリフト（トレンド）があっても偏らない推定量。

    σ² = mean[ ln(H/C)*ln(H/O) + ln(L/C)*ln(L/O) ]
    """
    h, l, o, c = bars["high"], bars["low"], bars["open"], bars["close"]
    rs = (np.log(h / c) * np.log(h / o)
          + np.log(l / c) * np.log(l / o))
    var = rs.rolling(window, min_periods=max(window // 2, 5)).mean()
    return _ann(var, bars_per_year)


def yang_zhang(bars: pd.DataFrame, window: int,
               bars_per_year: int = 8_760) -> pd.Series:
    """ギャップ（オーバーナイト）とドリフトの両方に対応する推定量。

    σ² = σ²_open + k*σ²_close + (1−k)*σ²_rs
      σ²_open  : 前足終値からの始値ギャップの分散
      σ²_close : 足の中の始値→終値の分散
      σ²_rs    : Rogers-Satchell成分
      k = 0.34 / (1.34 + (n+1)/(n−1))
    """
    o, h, l, c = bars["open"], bars["high"], bars["low"], bars["close"]
    c_prev = c.shift(1)

    log_oc_prev = np.log(o / c_prev)          # ギャップ成分
    log_co = np.log(c / o)                    # 足の中の成分
    rs = (np.log(h / c) * np.log(h / o)
          + np.log(l / c) * np.log(l / o))

    min_p = max(window // 2, 5)
    sigma_o = log_oc_prev.rolling(window, min_periods=min_p).var()
    sigma_c = log_co.rolling(window, min_periods=min_p).var()
    sigma_rs = rs.rolling(window, min_periods=min_p).mean()

    k = 0.34 / (1.34 + (window + 1) / (window - 1))
    var = sigma_o + k * sigma_c + (1 - k) * sigma_rs
    return _ann(var, bars_per_year)


ESTIMATORS = {
    "close_to_close": close_to_close,
    "parkinson": parkinson,
    "garman_klass": garman_klass,
    "rogers_satchell": rogers_satchell,
    "yang_zhang": yang_zhang,
}


def all_estimators(bars: pd.DataFrame, window: int = 168,
                   interval: str = "60") -> pd.DataFrame:
    """全推定量を並べたDataFrameを返す。"""
    bpy = BARS_PER_YEAR.get(interval, 8_760)
    return pd.DataFrame({
        name: fn(bars, window, bpy) for name, fn in ESTIMATORS.items()
    })


if __name__ == "__main__":
    from client import make_client
    from kline_fetch import fetch_klines

    c = make_client()
    bars = fetch_klines(c.http, "BTCUSDT", interval="60",
                        start="2026-01-01", verbose=False)

    vols = all_estimators(bars, window=168).dropna()
    print(f"BTCUSDT 1時間足 {len(bars)}本 / 推定窓 168本（1週間）\n")

    print(f"{'推定量':<18}{'直近':>10}{'平均':>10}{'最小':>10}"
          f"{'最大':>10}{'標準偏差':>10}")
    print("-" * 68)
    for name in ESTIMATORS:
        s = vols[name]
        print(f"{name:<18}{s.iloc[-1]:>9.1%}{s.mean():>10.1%}"
              f"{s.min():>10.1%}{s.max():>10.1%}{s.std():>10.1%}")

    print("\n※ 推定量ごとに水準が違う。")
    print("  重要なのは絶対値の一致ではなく、")
    print("  『同じ推定量を一貫して使う』ことと『変化に速く追随する』こと。")

"""vol_efficiency.py — 推定量の効率（分散の小ささ）を実測する
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

方法1: 合成データ（真のボラが既知）で誤差を測る
方法2: 実データで、短窓の推定が長窓の推定をどれだけ再現するか測る
"""
from __future__ import annotations

from dataclasses import dataclass

import numpy as np
import pandas as pd

from vol_estimators import ESTIMATORS


@dataclass
class EfficiencyResult:
    name: str
    rmse: float
    bias: float
    relative_efficiency: float    # close_to_close を1とした比

    def line(self) -> str:
        return (f"{self.name:<18}{self.rmse:>10.5f}{self.bias:>+11.5f}"
                f"{self.relative_efficiency:>12.2f}x")


def simulate_ohlc(n_bars: int = 2000, true_vol: float = 0.60,
                  bars_per_year: int = 8_760, ticks_per_bar: int = 200,
                  drift: float = 0.0, seed: int = 7) -> pd.DataFrame:
    """既知のボラティリティを持つOHLCを合成する。

    足の中を細かいステップで生成し、そこからOHLCを作る。
    真値が分かるので、推定量の誤差を直接測れる。
    """
    rng = np.random.default_rng(seed)
    sigma_per_bar = true_vol / np.sqrt(bars_per_year)
    sigma_per_tick = sigma_per_bar / np.sqrt(ticks_per_bar)
    mu_per_tick = drift / bars_per_year / ticks_per_bar

    price = 100.0
    rows = []
    for _ in range(n_bars):
        steps = rng.normal(mu_per_tick, sigma_per_tick, ticks_per_bar)
        path = price * np.exp(np.cumsum(steps))
        rows.append({
            "open": price, "high": path.max(), "low": path.min(),
            "close": path[-1], "volume": 1.0,
        })
        price = path[-1]

    idx = pd.date_range("2025-01-01", periods=n_bars, freq="1h", tz="UTC")
    return pd.DataFrame(rows, index=idx)


def measure_efficiency(bars: pd.DataFrame, true_vol: float,
                       window: int = 168,
                       bars_per_year: int = 8_760) -> list[EfficiencyResult]:
    """各推定量のRMSEとバイアスを測る。"""
    results = []
    baseline_mse = None

    for name, fn in ESTIMATORS.items():
        est = fn(bars, window, bars_per_year).dropna()
        if est.empty:
            continue
        err = est - true_vol
        mse = float((err ** 2).mean())
        if name == "close_to_close":
            baseline_mse = mse

        results.append(EfficiencyResult(
            name=name, rmse=float(np.sqrt(mse)), bias=float(err.mean()),
            relative_efficiency=(baseline_mse / mse) if (baseline_mse and mse) else 1.0,
        ))

    # baselineが後から確定する場合に備えて再計算
    if baseline_mse:
        for r in results:
            r.relative_efficiency = baseline_mse / (r.rmse ** 2)
    return results


def stability_on_real_data(bars: pd.DataFrame, short_window: int = 48,
                           long_window: int = 336,
                           bars_per_year: int = 8_760) -> pd.DataFrame:
    """実データで、短窓の推定が長窓の推定をどれだけ再現するか。

    真値が分からない実データでは、
    「長窓の推定を近似的な真値とみなす」ことで比較できる。
    """
    rows = []
    for name, fn in ESTIMATORS.items():
        short = fn(bars, short_window, bars_per_year)
        long_ = fn(bars, long_window, bars_per_year)
        df = pd.DataFrame({"short": short, "long": long_}).dropna()
        if len(df) < 50:
            continue
        err = df["short"] - df["long"]
        rows.append({
            "推定量": name,
            "RMSE": float(np.sqrt((err ** 2).mean())),
            "相関": float(df["short"].corr(df["long"])),
            "短窓の標準偏差": float(df["short"].std()),
        })
    out = pd.DataFrame(rows).set_index("推定量")
    if not out.empty:
        base = out.loc["close_to_close", "RMSE"]
        out["相対効率"] = (base / out["RMSE"]) ** 2
    return out


if __name__ == "__main__":
    TRUE_VOL = 0.60

    print(f"=== 合成データ（真のボラ {TRUE_VOL:.0%}）での比較 ===")
    print("  ドリフトなし")
    bars = simulate_ohlc(n_bars=3000, true_vol=TRUE_VOL, drift=0.0)
    print(f"\n{'推定量':<18}{'RMSE':>10}{'バイアス':>11}{'相対効率':>12}")
    print("-" * 52)
    for r in measure_efficiency(bars, TRUE_VOL):
        print(r.line())

    print("\n  強いドリフトあり（年率+200%のトレンド）")
    bars_drift = simulate_ohlc(n_bars=3000, true_vol=TRUE_VOL, drift=2.0)
    print(f"\n{'推定量':<18}{'RMSE':>10}{'バイアス':>11}{'相対効率':>12}")
    print("-" * 52)
    for r in measure_efficiency(bars_drift, TRUE_VOL):
        print(r.line())

    print("\n  -> ドリフトがあると Garman-Klass と Parkinson は過大評価する。")
    print("     Rogers-Satchell と Yang-Zhang はバイアスが小さい。")

    print("\n\n=== 実データでの安定性 ===")
    try:
        from client import make_client
        from kline_fetch import fetch_klines

        c = make_client()
        real = fetch_klines(c.http, "BTCUSDT", interval="60",
                            start="2025-01-01", verbose=False)
        print(stability_on_real_data(real).to_string(
            float_format=lambda v: f"{v:.4f}"))
    except Exception as e:
        print(f"  実データの取得に失敗: {e}")

"""realized_vol.py — 実現ボラティリティとシグネチャプロット
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

参考: Andersen & Bollerslev (1998), Zhang, Mykland & Ait-Sahalia (2005)

RV = Σ r_i²  （r_i は高頻度の対数リターン）

サンプリング間隔を短くするほど理論上は精度が上がるが、
実際にはビッド・アスク・バウンスなどのノイズでRVが膨らむ。
どの間隔が適切かは「シグネチャプロット」で決める。
"""
from __future__ import annotations

import numpy as np
import pandas as pd


def realized_variance(prices: pd.Series, freq: str) -> float:
    """指定間隔にリサンプリングしてRVを計算する。"""
    p = prices.resample(freq).last().dropna()
    if len(p) < 3:
        return float("nan")
    r = np.log(p / p.shift(1)).dropna()
    return float((r ** 2).sum())


def realized_vol_daily(ticks: pd.DataFrame, freq: str = "5min",
                       price_col: str = "price") -> pd.Series:
    """日次の実現ボラティリティ（年率）を計算する。"""
    prices = ticks[price_col]
    out = {}
    for day, group in prices.groupby(prices.index.date):
        rv = realized_variance(group, freq)
        if np.isfinite(rv):
            out[pd.Timestamp(day, tz="UTC")] = np.sqrt(rv * 365)
    return pd.Series(out, name="rv").sort_index()


def signature_plot(ticks: pd.DataFrame,
                   freqs: tuple[str, ...] = ("1s", "5s", "15s", "30s",
                                             "1min", "2min", "5min",
                                             "10min", "15min", "30min"),
                   price_col: str = "price") -> pd.DataFrame:
    """サンプリング間隔ごとのRVを並べる（シグネチャプロット）。

    短い間隔でRVが跳ね上がる＝マイクロストラクチャ・ノイズ。
    プラトー（平坦になる領域）の開始点が適切な間隔。
    """
    prices = ticks[price_col]
    days = prices.groupby(prices.index.date)

    rows = []
    for freq in freqs:
        rvs = []
        for _, g in days:
            rv = realized_variance(g, freq)
            if np.isfinite(rv) and rv > 0:
                rvs.append(np.sqrt(rv * 365))
        if rvs:
            rows.append({"freq": freq, "mean_rv": float(np.mean(rvs)),
                         "std_rv": float(np.std(rvs)), "days": len(rvs)})
    return pd.DataFrame(rows).set_index("freq")


def suggest_sampling(sig: pd.DataFrame, tolerance: float = 0.05) -> str:
    """シグネチャプロットからプラトーの開始点を推定する。

    最も長い間隔のRVを基準に、そこから tolerance 以内に
    収まる最短の間隔を「適切な間隔」とする。
    """
    if sig.empty:
        return "判定できません"
    baseline = sig["mean_rv"].iloc[-1]
    for freq, row in sig.iterrows():
        if abs(row["mean_rv"] - baseline) / baseline <= tolerance:
            return freq
    return str(sig.index[-1])


def two_scale_rv(prices: pd.Series, fast: str = "10s",
                 slow: str = "5min") -> float:
    """Two-Scale RV（ノイズ補正の簡易版）。

    参考: Zhang, Mykland & Ait-Sahalia (2005)
    高頻度のRVからノイズ成分を差し引く考え方。
    """
    rv_fast = realized_variance(prices, fast)
    rv_slow = realized_variance(prices, slow)
    if not (np.isfinite(rv_fast) and np.isfinite(rv_slow)):
        return float("nan")

    n_fast = len(prices.resample(fast).last().dropna())
    n_slow = len(prices.resample(slow).last().dropna())
    if n_fast <= 0 or n_slow <= 0:
        return float("nan")
    # ノイズ分散の推定値を差し引く
    return float(max(rv_slow - (n_slow / n_fast) * rv_fast, 0.0))


if __name__ == "__main__":
    from datetime import date, timedelta

    try:
        from tick_parser import load_range
    except ImportError:
        raise SystemExit("記事08の tick_parser.py が必要です")

    SYMBOL = "BTCUSDT"
    end = date.today() - timedelta(days=1)
    ticks = load_range(SYMBOL, end - timedelta(days=9), end)

    if ticks.empty:
        raise SystemExit("先に記事08でティックデータを取得してください")

    print(f"{len(ticks):,}ティック / "
          f"{(ticks.index[-1] - ticks.index[0]).days}日分\n")

    print("=== シグネチャプロット ===")
    sig = signature_plot(ticks)
    print(f"{'間隔':<8}{'平均RV(年率)':>14}{'標準偏差':>12}{'日数':>7}")
    print("-" * 42)
    for freq, r in sig.iterrows():
        bar = "#" * int(r["mean_rv"] * 30)
        print(f"{freq:<8}{r['mean_rv']:>13.1%}{r['std_rv']:>12.1%}"
              f"{int(r['days']):>7}  {bar}")

    print(f"\n推奨サンプリング間隔: {suggest_sampling(sig)}")
    print("  （これより短いとノイズでRVが膨らむ）")

    print("\n=== 日次の実現ボラティリティ ===")
    rv = realized_vol_daily(ticks, freq=suggest_sampling(sig))
    print(rv.to_string(float_format=lambda v: f"{v:.1%}"))

"""har_rv.py — HAR-RV（Heterogeneous AutoRegressive）モデル
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

参考: Corsi, F. (2009)
      "A Simple Approximate Long-Memory Model of Realized Volatility"

  RV_{t+1} = c + β_d*RV_t + β_w*RV_t^(w) + β_m*RV_t^(m) + ε

  RV^(w) = 直近5期間の平均、RV^(m) = 直近22期間の平均
  （日足なら1週間・1か月に相当）

線形回帰だけで組めるのに、GARCHに匹敵する予測精度が出ることが多い。
"""
from __future__ import annotations

from dataclasses import dataclass, field

import numpy as np
import pandas as pd


@dataclass
class HarModel:
    coefficients: np.ndarray          # [c, β_d, β_w, β_m]
    lags: tuple[int, int, int]
    in_sample_r2: float
    n_train: int
    log_space: bool = True
    residual_std: float = 0.0

    @property
    def names(self) -> list[str]:
        return ["定数", "日次", "週次", "月次"]

    def report(self) -> str:
        lines = [f"HAR-RV（lags={self.lags} / "
                 f"{'対数' if self.log_space else '実数'}空間）"]
        for name, coef in zip(self.names, self.coefficients):
            lines.append(f"  {name:<8}: {coef:>+9.4f}")
        lines.append(f"  In-sample R²: {self.in_sample_r2:.4f} "
                     f"（学習 {self.n_train}件）")
        return "\n".join(lines)


def build_features(rv: pd.Series,
                   lags: tuple[int, int, int] = (1, 5, 22)) -> pd.DataFrame:
    """日次・週次・月次の成分を作る。

    重要: すべて shift(1) してあり、未来を参照しない（記事31）。
    """
    d, w, m = lags
    out = pd.DataFrame(index=rv.index)
    out["rv_d"] = rv.rolling(d).mean().shift(1)
    out["rv_w"] = rv.rolling(w).mean().shift(1)
    out["rv_m"] = rv.rolling(m).mean().shift(1)
    out["target"] = rv
    return out.dropna()


def fit_har(rv: pd.Series, lags: tuple[int, int, int] = (1, 5, 22),
            log_space: bool = True) -> HarModel:
    """最小二乗法でHAR-RVを推定する。

    log_space=True にすると、ボラの非負制約と歪みに対応しやすい。
    """
    series = np.log(rv.clip(lower=1e-8)) if log_space else rv
    feat = build_features(series, lags)
    if len(feat) < 40:
        raise ValueError(f"学習データ不足 ({len(feat)}件)")

    X = np.column_stack([
        np.ones(len(feat)),
        feat["rv_d"].to_numpy(),
        feat["rv_w"].to_numpy(),
        feat["rv_m"].to_numpy(),
    ])
    y = feat["target"].to_numpy()

    coef, *_ = np.linalg.lstsq(X, y, rcond=None)
    pred = X @ coef
    resid = y - pred
    ss_res = float((resid ** 2).sum())
    ss_tot = float(((y - y.mean()) ** 2).sum())

    return HarModel(
        coefficients=coef, lags=lags,
        in_sample_r2=(1 - ss_res / ss_tot) if ss_tot else 0.0,
        n_train=len(feat), log_space=log_space,
        residual_std=float(resid.std(ddof=1)),
    )


def predict_har(model: HarModel, rv: pd.Series) -> float:
    """直近のデータから次期のRVを予測する。"""
    series = np.log(rv.clip(lower=1e-8)) if model.log_space else rv
    d, w, m = model.lags
    x = np.array([
        1.0,
        float(series.tail(d).mean()),
        float(series.tail(w).mean()),
        float(series.tail(m).mean()),
    ])
    pred = float(x @ model.coefficients)
    if model.log_space:
        # 対数正規の期待値補正
        pred = float(np.exp(pred + model.residual_std ** 2 / 2))
    return pred


def walk_forward_har(rv: pd.Series, train_size: int = 250,
                     lags: tuple[int, int, int] = (1, 5, 22),
                     log_space: bool = True) -> pd.DataFrame:
    """Walk-Forwardで予測精度を評価する（記事32の考え方）。"""
    preds, actuals, idx = [], [], []

    for i in range(train_size, len(rv)):
        train = rv.iloc[:i]
        try:
            model = fit_har(train, lags, log_space)
            preds.append(predict_har(model, train))
        except Exception:
            continue
        actuals.append(float(rv.iloc[i]))
        idx.append(rv.index[i])

    if not preds:
        return pd.DataFrame()

    df = pd.DataFrame({"pred": preds, "actual": actuals}, index=idx)
    df["error"] = df["pred"] - df["actual"]
    return df


def evaluate(df: pd.DataFrame, baseline: pd.Series | None = None) -> dict:
    """予測精度を、単純なベースラインと比較する。"""
    if df.empty:
        return {"error": "予測がありません"}

    rmse = float(np.sqrt((df["error"] ** 2).mean()))
    mae = float(df["error"].abs().mean())
    corr = float(df["pred"].corr(df["actual"]))

    # ベースライン: 「前日のRVがそのまま続く」というランダムウォーク予測
    naive = df["actual"].shift(1).dropna()
    aligned = df["actual"].loc[naive.index]
    naive_rmse = float(np.sqrt(((naive - aligned) ** 2).mean()))

    return {
        "件数": len(df),
        "RMSE": round(rmse, 5),
        "MAE": round(mae, 5),
        "相関": round(corr, 4),
        "ベースラインRMSE": round(naive_rmse, 5),
        "改善率": (f"{(1 - rmse / naive_rmse):.1%}"
                   if naive_rmse else "-"),
    }


if __name__ == "__main__":
    from client import make_client
    from kline_fetch import fetch_klines
    from vol_estimators import yang_zhang

    c = make_client()
    bars = fetch_klines(c.http, "BTCUSDT", interval="60",
                        start="2024-01-01", verbose=False)

    # 日次のボラティリティ系列を作る（Yang-Zhangを日次にリサンプル）
    daily_vol = yang_zhang(bars, window=24, bars_per_year=8_760) \
        .resample("D").last().dropna()

    print(f"日次ボラ系列: {len(daily_vol)}件 "
          f"({daily_vol.index[0].date()} 〜 {daily_vol.index[-1].date()})\n")

    model = fit_har(daily_vol)
    print(model.report())

    print(f"\n次期の予測ボラ: {predict_har(model, daily_vol):.1%}")
    print(f"直近の実績    : {daily_vol.iloc[-1]:.1%}")

    print("\n=== Walk-Forward評価 ===")
    wf = walk_forward_har(daily_vol, train_size=200)
    for k, v in evaluate(wf).items():
        print(f"  {k:<18}: {v}")

    print("\n※ 改善率がプラスなら、HAR-RVはランダムウォークより優れている。")
    print("  マイナスなら『前日の値をそのまま使う』ほうが良い（よくある）。")
    print("  予測できないなら、予測せず直近の実績を使えばよい。")

"""garch_compare.py — GARCH(1,1)とHAR-RVを比較する
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

  pip install arch   （無い場合は簡易EWMAで代替する）

実務的な結論: 複雑なモデルが常に優れるわけではない。
              EWMAやHAR-RVで十分なことが多い。
"""
from __future__ import annotations

import numpy as np
import pandas as pd


def ewma_vol(returns: pd.Series, lam: float = 0.94,
             bars_per_year: int = 365) -> pd.Series:
    """指数加重移動平均によるボラ推定（RiskMetrics方式）。

    σ²_t = λ*σ²_{t-1} + (1−λ)*r²_{t-1}
    実装が最も単純で、実務では今でも広く使われる。
    """
    r2 = (returns ** 2).shift(1)
    var = r2.ewm(alpha=1 - lam, adjust=False).mean()
    return np.sqrt(var * bars_per_year)


def fit_garch(returns: pd.Series, bars_per_year: int = 365):
    """GARCH(1,1)を推定する（archパッケージ）。"""
    try:
        from arch import arch_model
    except ImportError:
        return None

    # archはパーセント単位のリターンを想定
    r = returns.dropna() * 100
    model = arch_model(r, vol="Garch", p=1, q=1, mean="Zero", dist="t")
    res = model.fit(disp="off")

    cond_vol = res.conditional_volatility / 100 * np.sqrt(bars_per_year)
    return {
        "result": res,
        "conditional_vol": pd.Series(cond_vol, index=r.index),
        "params": dict(res.params),
        "forecast": float(
            np.sqrt(res.forecast(horizon=1, reindex=False)
                    .variance.iloc[-1, 0]) / 100 * np.sqrt(bars_per_year)),
    }


def compare_models(bars: pd.DataFrame, bars_per_year: int = 365) -> pd.DataFrame:
    """EWMA・GARCH・HAR-RV・Yang-Zhangを並べて比較する。"""
    from har_rv import fit_har, predict_har, walk_forward_har
    from vol_estimators import yang_zhang

    daily = bars["close"].resample("D").last().dropna()
    returns = np.log(daily / daily.shift(1)).dropna()

    yz = yang_zhang(bars, window=24, bars_per_year=8_760) \
        .resample("D").last().dropna()

    rows = []

    # 1) EWMA
    ew = ewma_vol(returns, 0.94, bars_per_year).dropna()
    common = ew.index.intersection(yz.index)
    if len(common) > 30:
        err = ew.loc[common] - yz.loc[common]
        rows.append({"モデル": "EWMA(0.94)",
                     "RMSE": float(np.sqrt((err ** 2).mean())),
                     "相関": float(ew.loc[common].corr(yz.loc[common]))})

    # 2) GARCH
    g = fit_garch(returns, bars_per_year)
    if g is not None:
        cv = g["conditional_vol"]
        common = cv.index.intersection(yz.index)
        if len(common) > 30:
            err = cv.loc[common] - yz.loc[common]
            rows.append({"モデル": "GARCH(1,1)",
                         "RMSE": float(np.sqrt((err ** 2).mean())),
                         "相関": float(cv.loc[common].corr(yz.loc[common]))})

    # 3) HAR-RV（Walk-Forward）
    wf = walk_forward_har(yz, train_size=200)
    if not wf.empty:
        rows.append({"モデル": "HAR-RV (WF)",
                     "RMSE": float(np.sqrt((wf["error"] ** 2).mean())),
                     "相関": float(wf["pred"].corr(wf["actual"]))})

    # 4) ベースライン（前日の値）
    naive_err = (yz.shift(1) - yz).dropna()
    rows.append({"モデル": "前日値（基準）",
                 "RMSE": float(np.sqrt((naive_err ** 2).mean())),
                 "相関": float(yz.shift(1).corr(yz))})

    df = pd.DataFrame(rows).set_index("モデル")
    base = df.loc["前日値（基準）", "RMSE"]
    df["改善率"] = 1 - df["RMSE"] / base
    return df.sort_values("RMSE")


if __name__ == "__main__":
    from client import make_client
    from kline_fetch import fetch_klines

    c = make_client()
    bars = fetch_klines(c.http, "BTCUSDT", interval="60",
                        start="2024-01-01", verbose=False)

    print("=== ボラティリティモデルの比較 ===")
    print("  （Yang-Zhang推定量を近似的な真値として評価）\n")
    df = compare_models(bars)
    print(df.to_string(float_format=lambda v: f"{v:.4f}"))

    g = fit_garch(np.log(bars["close"].resample("D").last().dropna()
                         .pct_change() + 1).dropna())
    if g:
        print(f"\nGARCHパラメータ: "
              f"omega={g['params'].get('omega', 0):.5f} "
              f"alpha={g['params'].get('alpha[1]', 0):.4f} "
              f"beta={g['params'].get('beta[1]', 0):.4f}")
        print(f"  α+β = {g['params'].get('alpha[1]', 0) + g['params'].get('beta[1]', 0):.4f}"
              f"（1に近いほど persistence が高い）")
        print(f"  次期予測: {g['forecast']:.1%}")
    else:
        print("\n（archパッケージが無いためGARCHはスキップ）")

    print("\n※ 『前日値』に勝てないモデルは使う意味がない。")
    print("  複雑なモデルほど良いとは限らない。")

"""atr_tools.py — ATR系指標（外部依存なし）
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

ATRは「価格の単位」で表されるため銘柄間の比較ができない。
NATR（正規化ATR）にすれば比較可能になる。
"""
from __future__ import annotations

import numpy as np
import pandas as pd


def true_range(bars: pd.DataFrame) -> pd.Series:
    """True Range = max(H−L, |H−C_prev|, |L−C_prev|)"""
    h, l, c = bars["high"], bars["low"], bars["close"]
    prev_c = c.shift(1)
    return pd.concat([h - l, (h - prev_c).abs(), (l - prev_c).abs()],
                     axis=1).max(axis=1)


def atr(bars: pd.DataFrame, period: int = 14,
        method: str = "wilder") -> pd.Series:
    """ATR。Wilderの平滑化（既定）と単純移動平均を選べる。

    Wilder方式は TA-Lib の ATR と一致する。
    """
    tr = true_range(bars)
    if method == "sma":
        return tr.rolling(period, min_periods=period).mean()
    # Wilder: alpha = 1/period のEMA
    return tr.ewm(alpha=1 / period, adjust=False,
                  min_periods=period).mean()


def natr(bars: pd.DataFrame, period: int = 14) -> pd.Series:
    """正規化ATR（%）。銘柄間で比較できる。"""
    return atr(bars, period) / bars["close"] * 100


def atr_to_annual_vol(bars: pd.DataFrame, period: int = 14,
                      bars_per_year: int = 8_760) -> pd.Series:
    """ATRを年率ボラティリティ相当に換算する（近似）。

    True Rangeは平均絶対偏差に近い性質を持つため、
    正規分布を仮定すると σ ≈ TR / sqrt(2/π) の関係がある。
    """
    natr_ratio = atr(bars, period) / bars["close"]
    return natr_ratio * np.sqrt(np.pi / 2) * np.sqrt(bars_per_year)


def compare_atr_and_vol(bars: pd.DataFrame, period: int = 14,
                        window: int = 168,
                        bars_per_year: int = 8_760) -> pd.DataFrame:
    """ATR由来のボラと、Yang-Zhangの推定を比較する。"""
    from vol_estimators import yang_zhang

    df = pd.DataFrame({
        "atr": atr(bars, period),
        "natr_pct": natr(bars, period),
        "atr_annual_vol": atr_to_annual_vol(bars, period, bars_per_year),
        "yang_zhang": yang_zhang(bars, window, bars_per_year),
    }).dropna()

    df["ratio"] = df["atr_annual_vol"] / df["yang_zhang"]
    return df


if __name__ == "__main__":
    from client import make_client
    from kline_fetch import fetch_klines

    c = make_client()

    print(f"{'symbol':<12}{'ATR':>14}{'NATR':>9}"
          f"{'ATR年率換算':>13}{'Yang-Zhang':>13}{'比率':>8}")
    print("-" * 70)

    for symbol in ("BTCUSDT", "ETHUSDT", "SOLUSDT", "DOGEUSDT"):
        try:
            bars = fetch_klines(c.http, symbol, interval="60",
                                start="2026-05-01", verbose=False)
            df = compare_atr_and_vol(bars)
            if df.empty:
                continue
            last = df.iloc[-1]
            print(f"{symbol:<12}{last['atr']:>14.4f}"
                  f"{last['natr_pct']:>8.2f}%"
                  f"{last['atr_annual_vol']:>12.1%}"
                  f"{last['yang_zhang']:>12.1%}"
                  f"{last['ratio']:>8.2f}")
        except Exception as e:
            print(f"{symbol:<12} 取得失敗: {e}")

    print("\n※ ATRは価格の単位なので銘柄間で比較できない。")
    print("  NATR（%）なら比較可能。")
    print("  ATRの年率換算はYang-Zhangとおおむね同水準になるが、")
    print("  完全一致はしない（測っている性質が微妙に違う）。")

"""vol_screener.py — 全USDT無期限銘柄のボラティリティ・スクリーナー
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

  ボラが高いだけの銘柄は罠。以下を同時に評価する:
    - 年率ボラティリティ（Yang-Zhang）
    - 24h出来高（流動性）
    - スプレッド（往復コスト）
    - ボラ / スプレッド比（取れる幅がコストに対して十分か）

  記事05の流量ゲートを通し、記事03の取引ルールで最小ロットも確認する。
"""
from __future__ import annotations

import time
from dataclasses import dataclass

import pandas as pd

from client import make_client
from kline_fetch import fetch_klines
from vol_estimators import yang_zhang


@dataclass(frozen=True)
class ScreenConfig:
    category: str = "linear"
    quote: str = "USDT"
    interval: str = "60"          # 1時間足
    lookback_bars: int = 24 * 30  # 30日
    vol_window: int = 24          # ボラ推定の窓（1日）
    min_turnover_24h: float = 5_000_000.0   # 24h売買代金の下限（USD）
    min_vol: float = 0.30         # 年率30%未満は除外
    max_vol: float = 3.00         # 年率300%超は除外（動きすぎ）
    max_spread_bps: float = 5.0   # 実効スプレッドの上限
    top_n: int = 30
    sleep: float = 0.06           # レート制限への配慮（記事05）


def list_candidates(http, cfg: ScreenConfig) -> pd.DataFrame:
    """tickersから候補銘柄を1リクエストで取得し、粗く絞る。"""
    r = http.get_tickers(category=cfg.category)
    if r["retCode"] != 0:
        raise RuntimeError(f"get_tickers失敗: {r['retCode']} {r['retMsg']}")

    rows = []
    for t in r["result"]["list"]:
        sym = t["symbol"]
        if not sym.endswith(cfg.quote):
            continue
        bid, ask = float(t.get("bid1Price") or 0), float(t.get("ask1Price") or 0)
        if bid <= 0 or ask <= 0:
            continue
        mid = (bid + ask) / 2
        rows.append({
            "symbol": sym,
            "last": float(t["lastPrice"]),
            "turnover_24h": float(t.get("turnover24h") or 0),
            "spread_bps": (ask - bid) / mid * 10_000,
            "chg_24h_pct": float(t.get("price24hPcnt") or 0) * 100,
            "funding": float(t.get("fundingRate") or 0),
            "oi_value": float(t.get("openInterestValue") or 0),
        })

    df = pd.DataFrame(rows)
    if df.empty:
        return df

    before = len(df)
    df = df[df["turnover_24h"] >= cfg.min_turnover_24h]
    df = df[df["spread_bps"] <= cfg.max_spread_bps]
    print(f"  一次スクリーニング: {before} → {len(df)}銘柄")
    return df.sort_values("turnover_24h", ascending=False).reset_index(drop=True)


def measure_volatility(http, symbols: list[str],
                       cfg: ScreenConfig) -> pd.DataFrame:
    """候補銘柄それぞれのボラティリティを実測する。"""
    bars_per_year = 8_760 if cfg.interval == "60" else 365
    rows = []

    for i, sym in enumerate(symbols, 1):
        try:
            bars = fetch_klines(http, sym, interval=cfg.interval,
                                limit_bars=cfg.lookback_bars, verbose=False)
            if len(bars) < cfg.vol_window * 3:
                continue

            yz = yang_zhang(bars, cfg.vol_window, bars_per_year).dropna()
            if yz.empty:
                continue

            current = float(yz.iloc[-1])
            median = float(yz.median())
            rows.append({
                "symbol": sym,
                "vol_now": current,
                "vol_median": median,
                "vol_ratio": current / median if median else float("nan"),
                "vol_p10": float(yz.quantile(0.10)),
                "vol_p90": float(yz.quantile(0.90)),
                "bars": len(bars),
            })
        except Exception as e:
            print(f"  [skip] {sym}: {e}")
        finally:
            time.sleep(cfg.sleep)

        if i % 20 == 0:
            print(f"  ... {i}/{len(symbols)}銘柄")

    return pd.DataFrame(rows)


def screen(cfg: ScreenConfig = ScreenConfig()) -> pd.DataFrame:
    """スクリーニング本体。"""
    c = make_client()

    print("=== ボラティリティ・スクリーナー ===")
    cand = list_candidates(c.http, cfg)
    if cand.empty:
        return cand

    targets = cand["symbol"].head(cfg.top_n * 3).tolist()
    print(f"  ボラ実測: {len(targets)}銘柄")
    vol = measure_volatility(c.http, targets, cfg)
    if vol.empty:
        return vol

    df = cand.merge(vol, on="symbol")
    df = df[(df["vol_now"] >= cfg.min_vol) & (df["vol_now"] <= cfg.max_vol)]

    # 1時間あたりの期待変動幅（bps）とスプレッドの比
    # 年率ボラ → 1時間ボラ に換算して bps 表示
    hourly_bps = df["vol_now"] / (8_760 ** 0.5) * 10_000
    df["move_per_hour_bps"] = hourly_bps
    df["edge_ratio"] = hourly_bps / df["spread_bps"].clip(lower=0.1)

    return df.sort_values("edge_ratio", ascending=False) \
             .head(cfg.top_n).reset_index(drop=True)


if __name__ == "__main__":
    result = screen()

    if result.empty:
        print("\n条件に合う銘柄がありません。")
    else:
        print(f"\n{'symbol':<12}{'年率ボラ':>9}{'中央値比':>9}"
              f"{'ｽﾌﾟﾚｯﾄﾞ':>9}{'1h変動':>9}{'比率':>7}{'24h代金':>12}")
        print("-" * 70)
        for _, r in result.iterrows():
            print(f"{r['symbol']:<12}{r['vol_now']:>8.1%}"
                  f"{r['vol_ratio']:>9.2f}"
                  f"{r['spread_bps']:>8.2f}b"
                  f"{r['move_per_hour_bps']:>8.0f}b"
                  f"{r['edge_ratio']:>7.1f}"
                  f"{r['turnover_24h'] / 1e6:>11.1f}M")

        print("\n【読み方】")
        print("  中央値比 > 1.5 : 平常時より明らかにボラが高い（イベント中の可能性）")
        print("  中央値比 < 0.7 : 静かな状態（レンジ戦略向き）")
        print("  比率           : 1時間の期待変動 ÷ スプレッド。")
        print("                   これが小さい銘柄は、動く前にコストで負ける。")
        print("\n※ 銘柄選定は入口にすぎない。実際に発注できるかは")
        print("  記事03の最小ロット・ティックサイズも確認すること。")

"""vol_regime.py — ボラティリティレジームの判定とサイズ調整
本コードは本人利用に限り自由に使用可。再配布・商用利用・販売は禁止。

参考: Moskowitz, Ooi & Pedersen (2012) "Time Series Momentum"
        — ボラの逆数でポジションを調整すると、リスク調整後リターンが改善する
      Harvey et al. (2018) "The Impact of Volatility Targeting"

  ボラターゲティング: position ∝ target_vol / realized_vol
  「ボラが高いときは小さく、低いときは大きく」持つだけの単純な規則が、
  多くの戦略でシャープレシオを改善する。
"""
from __future__ import annotations

from dataclasses import dataclass
from enum import Enum

import numpy as np
import pandas as pd


class Regime(str, Enum):
    CALM = "静穏"
    NORMAL = "通常"
    ELEVATED = "上昇"
    STRESS = "ストレス"


@dataclass(frozen=True)
class RegimeConfig:
    lookback: int = 24 * 60       # レジーム判定に使う履歴（1時間足で60日）
    calm_pct: float = 0.25        # 下位25%以下 → 静穏
    elevated_pct: float = 0.75    # 上位25%以上 → 上昇
    stress_pct: float = 0.95      # 上位5%以上   → ストレス
    target_vol: float = 0.50      # 年率50%を目標
    max_leverage: float = 3.0     # サイズ倍率の上限
    min_leverage: float = 0.0


def classify_regime(vol_series: pd.Series,
                    cfg: RegimeConfig = RegimeConfig()) -> Regime:
    """直近のボラが、過去分布のどこに位置するかで判定する。

    絶対値ではなく分位点で判定するのが要点。
    「年率80%は高い」は銘柄によって成立しない。
    """
    hist = vol_series.tail(cfg.lookback).dropna()
    if len(hist) < 30:
        return Regime.NORMAL

    current = float(hist.iloc[-1])
    rank = float((hist < current).mean())

    if rank >= cfg.stress_pct:
        return Regime.STRESS
    if rank >= cfg.elevated_pct:
        return Regime.ELEVATED
    if rank <= cfg.calm_pct:
        return Regime.CALM
    return Regime.NORMAL


def regime_series(vol_series: pd.Series,
                  cfg: RegimeConfig = RegimeConfig()) -> pd.Series:
    """全期間のレジームを計算する（過去のみ参照、未来を見ない）。"""
    ranks = vol_series.rolling(cfg.lookback, min_periods=30) \
                      .rank(pct=True)

    out = pd.Series(Regime.NORMAL.value, index=vol_series.index, dtype=object)
    out[ranks <= cfg.calm_pct] = Regime.CALM.value
    out[ranks >= cfg.elevated_pct] = Regime.ELEVATED.value
    out[ranks >= cfg.stress_pct] = Regime.STRESS.value
    out[ranks.isna()] = None
    return out


def vol_target_size(current_vol: float,
                    cfg: RegimeConfig = RegimeConfig()) -> float:
    """ボラターゲティングによるサイズ倍率。

    target_vol / current_vol をクリップして返す。
    記事23のロットサイジングに、この倍率を掛ける。
    """
    if not np.isfinite(current_vol) or current_vol <= 0:
        return 0.0
    return float(np.clip(cfg.target_vol / current_vol,
                         cfg.min_leverage, cfg.max_leverage))


def regime_statistics(bars: pd.DataFrame, vol_series: pd.Series,
                      cfg: RegimeConfig = RegimeConfig()) -> pd.DataFrame:
    """レジーム別に、次足リターンの性質を集計する。

    「荒れた相場では順張りが効く／効かない」といった仮説を
    実データで確認するための土台。
    """
    reg = regime_series(vol_series, cfg)
    fwd = bars["close"].pct_change().shift(-1)   # 次足リターン

    df = pd.DataFrame({"regime": reg, "fwd": fwd}).dropna()
    if df.empty:
        return df

    g = df.groupby("regime")["fwd"]
    out = pd.DataFrame({
        "件数": g.size(),
        "平均リターン_bps": g.mean() * 10_000,
        "標準偏差_bps": g.std() * 10_000,
        "勝率": g.apply(lambda s: (s > 0).mean()),
        "最大下落_bps": g.min() * 10_000,
        "最大上昇_bps": g.max() * 10_000,
    })
    # 1トレードあたりのシャープ相当（年率化はしない）
    out["平均/標準偏差"] = out["平均リターン_bps"] / out["標準偏差_bps"]
    order = [r.value for r in
             (Regime.CALM, Regime.NORMAL, Regime.ELEVATED, Regime.STRESS)]
    return out.reindex([o for o in order if o in out.index])


def backtest_vol_targeting(bars: pd.DataFrame, vol_series: pd.Series,
                           signal: pd.Series,
                           cfg: RegimeConfig = RegimeConfig()) -> dict:
    """ボラターゲティングの有無で成績を比較する。

    signal: +1 / 0 / -1 のポジション方向（前バーまでの情報で決定済み）
    """
    ret = bars["close"].pct_change()
    sizes = vol_series.apply(lambda v: vol_target_size(v, cfg)).shift(1)

    df = pd.DataFrame({
        "ret": ret, "signal": signal.shift(1), "size": sizes,
    }).dropna()
    if df.empty:
        return {}

    plain = df["signal"] * df["ret"]
    scaled = df["signal"] * df["size"] * df["ret"]

    def summarize(s: pd.Series) -> dict:
        ann = np.sqrt(8_760)
        std = float(s.std())
        equity = (1 + s).cumprod()
        dd = float((equity / equity.cummax() - 1).min())
        return {
            "累積リターン": float(equity.iloc[-1] - 1),
            "年率ボラ": std * ann,
            "シャープ相当": (float(s.mean()) / std * ann) if std else 0.0,
            "最大DD": dd,
        }

    return {"固定サイズ": summarize(plain),
            "ボラターゲット": summarize(scaled),
            "平均倍率": float(df["size"].mean())}


if __name__ == "__main__":
    from client import make_client
    from kline_fetch import fetch_klines
    from vol_estimators import yang_zhang

    c = make_client()
    bars = fetch_klines(c.http, "BTCUSDT", interval="60",
                        start="2024-01-01", verbose=False)
    vol = yang_zhang(bars, window=24, bars_per_year=8_760)

    cfg = RegimeConfig()
    print(f"現在のレジーム: {classify_regime(vol, cfg).value}")
    print(f"現在の年率ボラ: {vol.iloc[-1]:.1%}")
    print(f"サイズ倍率    : {vol_target_size(float(vol.iloc[-1]), cfg):.2f}倍"
          f"（目標 {cfg.target_vol:.0%}）")

    print("\n=== レジーム別の次足リターン統計 ===")
    stats = regime_statistics(bars, vol, cfg)
    print(stats.to_string(float_format=lambda v: f"{v:.3f}"))

    # 単純な順張りシグナルで、ボラターゲティングの効果を見る
    ma_fast = bars["close"].rolling(24).mean()
    ma_slow = bars["close"].rolling(24 * 7).mean()
    sig = (ma_fast > ma_slow).astype(int) * 2 - 1

    print("\n=== ボラターゲティングの効果（単純な移動平均シグナル） ===")
    res = backtest_vol_targeting(bars, vol, sig, cfg)
    for name in ("固定サイズ", "ボラターゲット"):
        m = res.get(name, {})
        if not m:
            continue
        print(f"\n[{name}]")
        for k, v in m.items():
            print(f"  {k:<12}: {v:>8.2%}" if k != "シャープ相当"
                  else f"  {k:<12}: {v:>8.2f}")
    print(f"\n平均サイズ倍率: {res.get('平均倍率', 0):.2f}倍")

    print("\n※ この比較は手数料・スリッページを含まない。")
    print("  ボラターゲティングは売買回数を増やすため、")
    print("  記事31の枠組みでコストを入れて再評価すること。")