はじめに

この記事はマケデコ Advent Calendar 2024の10日目の記事です。

Realized Volatility (RV) という、金融商品の価格のボラティリティの推定量があります。

RVはティックレベル～分足レベルの高頻度の価格データから計算できます。モデルフリーで精度のよい日次ボラティリティの推定量であることが知られているとともに、連続な価格変動によるボラティリティとジャンプによるボラティリティを分離するように拡張できるのがメリットです。

この記事ではまずRVの理論をざっくり解説したあと、GMOコインのAPIからドル円の分足データを取得してドル円のRVの計算を実装します。

為替介入があった日にジャンプによるボラティリティが大きかったことが確認できました。これは、RVによるボラティリティの推定と、ジャンプ拡散過程を考慮することによるボラティリティの連続成分とジャンプ成分の分離が有用なことを示唆します。

理論

ボラティリティはリスク管理上などで重要な変数ですが真の値は直接観測できない潜在変数であるため、その推定方法は数理ファイナンスの主要な分野となっています。

一定期間の日次収益率のローリング標準偏差とする方法はシンプルながらよく用いられる方法ですが、ローリングウィンドウの間はボラティリティが一定であることを仮定するため、ボラティリティの日々のシャープな変動をとらえることができません。

日々変動するボラティリティを推定するため、日次のボラティリティを数理モデルで表し、日次の収益率からボラティリティを推定する方法があります。有名なGARCHモデルや、状態空間モデルの一種であるStochastic Volatility (SV) モデルはこれに該当します。なお、昨年のアドベントカレンダーの記事ではSVモデルをStanで実装してみました。興味があれば読んでみてください。

一方、高頻度の価格データを用いることでボラティリティの推定量を直接求めるアプローチがあります。特に2000年代～2010年代に注目された分野です。最も代表的な推定量であるAndersen and Bolleslev (1998) が示したRVについて、簡単に解説します。

連続過程

ある時点$s$における資産の対数価格を$p(s)$とします。

いま、$t$日に等間隔に$n$回価格が観測されるとして¹、$t$日の観測時点$i$における日中収益率を

$$ r_{t, i} = p(t-1 + i/n) - p(t-1 + (i-1)/n), \quad i=1, 2, \dots, n $$

と定義します。要するに1時点前との対数収益率です。

このとき、$t$日のRVは以下で求められます。

$$ RV_{t} = \sum_{i=1}^n r_{t, i}^2 $$

日中収益率をそれぞれ2乗して足し合わせただけですね。

RVは真のボラティリティの一致推定量であることが知られています²。簡単な計算方法に見えますが、対数価格がセミマルチンゲールな確率過程に従うとき、RVは真のボラティリティの一致推定量であるという理論的な裏付けがあります。

いま具体的に、$p(s)$が以下の確率過程に従うとします³。$W(s)$は標準ブラウン運動です。

$$ dp(s) = \mu(s) ds + \sigma(s) dW(s) $$

$t$日における真のボラティリティ$IV_{t}$ (Integrated Volatility; IV) は、瞬間的なボラティリティ$\sigma(s)$を一日分積分したものです。

$$ IV_{t} = \int_{t-1}^{t} \sigma(s)^2 ds $$

ただし、積分範囲の$t-1, t$は、それぞれ$t-1$日, $t$日の最後の観測時点であることを示します。

ここで問題になるのは、$\sigma(s)$は直接観測できないパラメータなので$IV_{t}$も直接観測できないということです。そのため、何かしらの方法で推定値を求めることになります。

$n \rightarrow \infty$のとき$RV_{t} \rightarrow IV_{t}$となるため⁴、$n$を大きく取るとき、RVはIVの精度のよい一致推定量となります。

なお、$t$日における対前日の対数収益率を$r_{t}$とすると、以下の$\sigma_{t}$が$t$日における真のボラティリティです。

$$ \begin{aligned} r_{t} &= E_{t-1} [r_{t}] + \epsilon_{t} \\ \epsilon_{t} &= \sigma_{t} z_{t}, \quad \sigma_{t} > 0, \quad z_{t} \sim i.i.d., \quad E[z_{t}] = 0, \quad Var[z_{t}] = 1 \end{aligned} $$

$E_{t-1} [r_{t}]= 0, z_{t} \sim N(0, 1)$とすれば$r_{t} \sim N(0, \sigma_{t}^2)$です。$RV_{t}$は$\sigma_{t}^2$であるというとRVとは何かイメージしやすいのではないでしょうか。

ジャンプ過程

実際の対数価格は連続な確率過程ではなく、しばしば不連続なジャンプを含みます。以上の議論をジャンプ過程に拡張してみましょう。

$p(s)$は以下のジャンプ拡散過程に従うとします。

$$ dp(s) = \mu(s) ds + \sigma(s) dW(s) + \kappa(s) dN(s) $$

$N(s)$は時点$s$でジャンプがある場合$dN(s) = 1$、ない場合$dN(s) = 0$となるポアソン過程です。$\kappa(s)$は時点$s$におけるジャンプの大きさを表します。

連続の場合と同様に、$RV_{t} = \sum_{i=1}^n r_{t, i}^2$とすると、$n \rightarrow \infty$のとき、

$$ RV_{t} \rightarrow \int_{t - 1}^{t} \sigma(s)^2 ds + \sum_{t-1 < s \leq t} \kappa(s)^2 $$

となります。

右辺第1項と右辺第2項はそれぞれ連続 (Continuous) な価格変動によるボラティリティとジャンプ (Jump) の価格変動によるボラティリティです。それぞれの推定量を$C_{t}, J_{t}$と表します。

ジャンプ過程: ジャンプ部分に由来するボラティリティの分離

$C_{t}, J_{t}$は異なる性質を持つリスクなので、$RV_{t}$を$C_{t}, J_{t}$に分解できると嬉しいことがありそうです。分解の方法については提案されている手法がいくつかありますが、Barndorff-Nielsen and Shephard (2004, 2006) の方法を説明します。

次の式の$BV_{t}$ (Bipower Variation; BV) は$n \rightarrow \infty$で連続部分のボラティリティに確率収束します。

$$ BV_{t} = \mu_{1}^{-2} \sum_{i=2}^{n} |r_{t, i}| |r_{t, i-1}| $$

ただし、$\mu_{1} = 2^{1/2} \Gamma(1) \Gamma(1/2)^{-1}$です。

以上より、$n \rightarrow \infty$のとき、$C_{t} = BV_{t}, J_{t} = \max(RV_{t} - BV_{t}, 0)$となります。

ただし、ジャンプの存在が有意かどうかを検定することが実証上は多いです。ジャンプが存在しないという帰無仮説のもとでは、以下の統計量$Z_{t}$は漸近的に標準正規分布に従います。

$$ Z_{t} = \frac{\log RV_{t} - \log BV_{t}}{((\mu_{1}^{-4} + 2\mu_{1}^{-2} - 5) TQ_{t} BV_{t}^{-2} / n)^{1/2}} $$

ただし、$TQ_{t}$ (Tri-power Quarticity; TQ) は

$$ TQ_{t} = n \mu_{4/3}^{-3} \sum_{i=3}^{n} |r_{t, i}|^{4/3} |r_{t, i-1}|^{4/3} |r_{t, i-2}|^{4/3} $$

であり、$\mu_{4/3} = 2^{2/3} \Gamma(7/6) \Gamma(1/2)^{-1}$です。

有意水準を5%のように設定して⁵、帰無仮説が棄却されれば$J_{t} = RV_{t} - BV_{t}$、棄却されなければ$J_{t} = 0$とします。

nの設定

以上の議論は$n \rightarrow \infty$で行ってきたので、実装する際にはティックデータを用いるのが一見最もよさそうに思われます。

しかし、実際のマーケットでは、価格には様々なノイズ（マイクロマーケットストラクチャーノイズといいます）が含まれています。一例を挙げると、価格には呼値が設定されていますが、最も高い買値で約定するか最も安い売値で約定するかによって、真の価格は変わらなくても実際の価格はBidとAskを行ったり来たりします。なのでnを細かくするほどRVにはノイズが多く含まれ、IVの一致推定量としての精度が落ちます。

ではノイズを考慮してnをどの程度に取るのが望ましいかですが、アセットクラスや流動性によって多少異なるものの、RVでは5分間隔の価格データを用いるのがおおむねよいとされます (Liu, Patton and Sheppard (2015))。以前より実証では5分がよく用いられていたのですが、5分のRVは他の時間間隔や他のボラティリティの推定量と比べてだいぶいい推定精度のパフォーマンスを示すことがこの論文で示されました⁶。

実装

RVを実装してみましょう。上の結果の導出は難しいですが、結果の実装は難しくありません。

データの取得

GMOコインが為替の分足データをAPI経由で無料で提供しています。ここからドル円の5分足をJSTの2023/10/30から2024/12/6まで取得します。分足が簡単に手に入りレスポンスも速い、いいAPIです。

環境はpython=3.11.5, polars=1.16.0, plotnine=0.14.3です。

import datetime
import json
import math
import time

import plotnine as pn
import polars as pl
import requests
import scipy as sp

endpoint = "https://forex-api.coin.z.com/public/v1/klines"
dates = [i.strftime("%Y%m%d") for i in pl.date_range(
    start=datetime.date(2023, 10, 28),
    end=datetime.date(2024, 12, 6),
    interval="1d",
    eager=True
)]

res = []
for date in dates:
    params = {
        "symbol": "USD_JPY",
        "priceType": "ASK",
        "interval": "5min",
        "date": date,
    }
    resp = requests.get(endpoint, params=params)
    # データが存在しない日（市場が開いていない日）は空のリスト("[]")が返る
    res.append(pl.DataFrame(json.loads(resp.text)["data"]))
    time.sleep(1)

df = pl.concat([i for i in res if not i.is_empty()])
df = (
    df
    .with_columns(
        open=pl.col("open").cast(float),
        high=pl.col("high").cast(float),
        low=pl.col("low").cast(float),
        close=pl.col("close").cast(float),
    )
)

取得したデータをpolars.DataFrameで持ちます。

shape: (82_248, 5)
┌───────────────┬─────────┬─────────┬─────────┬─────────┐
│ openTime      ┆ open    ┆ high    ┆ low     ┆ close   │
│ ---           ┆ ---     ┆ ---     ┆ ---     ┆ ---     │
│ str           ┆ f64     ┆ f64     ┆ f64     ┆ f64     │
╞═══════════════╪═════════╪═════════╪═════════╪═════════╡
│ 1698616800000 ┆ 149.641 ┆ 149.685 ┆ 149.641 ┆ 149.669 │
│ 1698617100000 ┆ 149.669 ┆ 149.669 ┆ 149.636 ┆ 149.644 │
│ 1698617400000 ┆ 149.644 ┆ 149.677 ┆ 149.638 ┆ 149.67  │
│ …             ┆ …       ┆ …       ┆ …       ┆ …       │
│ 1733517900000 ┆ 150.103 ┆ 150.123 ┆ 150.096 ┆ 150.105 │
│ 1733518200000 ┆ 150.105 ┆ 150.109 ┆ 150.075 ┆ 150.105 │
│ 1733518500000 ┆ 150.106 ┆ 150.111 ┆ 150.091 ┆ 150.096 │
└───────────────┴─────────┴─────────┴─────────┴─────────┘

openTime列はUTCのエポックミリ秒なので、JSTに直します。

次に、RVなどの計算において何時から何時までを1日として扱うかですが、APIのパスパラメータでdate=yyyymmddを指定すると、JSTのyyyymmddの6:00-翌日5:59（月曜日はJST7:00-5:59）のデータが得られるので、JSTの6:00-5:59を1日とすることにします。

df = (
    df
    .rename({"openTime": "openTimeUtc"})
    .with_columns(
        timestamp_utc=pl.from_epoch(pl.col("openTimeUtc").cast(int), time_unit="ms")
    )
    .with_columns(
        timestamp_jst=pl.col("timestamp_utc")+datetime.timedelta(hours=9),
        timestamp=pl.col("timestamp_utc")+datetime.timedelta(hours=9)-datetime.timedelta(hours=6)
    )
    .with_columns(date=pl.col("timestamp").dt.date())
)

こんなデータができました。表示上の都合で一部の列に絞っています。

JSTの2023/10/30の7:00-7:05は始値が149.641, 終値が149.669ということを示します。

shape: (82_248, 6)
┌─────────┬─────────┬─────────┬─────────┬─────────────────────┬────────────┐
│ open    ┆ high    ┆ low     ┆ close   ┆ timestamp_jst       ┆ date       │
│ ---     ┆ ---     ┆ ---     ┆ ---     ┆ ---                 ┆ ---        │
│ f64     ┆ f64     ┆ f64     ┆ f64     ┆ datetime[ms]        ┆ date       │
╞═════════╪═════════╪═════════╪═════════╪═════════════════════╪════════════╡
│ 149.641 ┆ 149.685 ┆ 149.641 ┆ 149.669 ┆ 2023-10-30 07:00:00 ┆ 2023-10-30 │
│ 149.669 ┆ 149.669 ┆ 149.636 ┆ 149.644 ┆ 2023-10-30 07:05:00 ┆ 2023-10-30 │
│ 149.644 ┆ 149.677 ┆ 149.638 ┆ 149.67  ┆ 2023-10-30 07:10:00 ┆ 2023-10-30 │
│ …       ┆ …       ┆ …       ┆ …       ┆ …                   ┆ …          │
│ 150.103 ┆ 150.123 ┆ 150.096 ┆ 150.105 ┆ 2024-12-07 05:45:00 ┆ 2024-12-06 │
│ 150.105 ┆ 150.109 ┆ 150.075 ┆ 150.105 ┆ 2024-12-07 05:50:00 ┆ 2024-12-06 │
│ 150.106 ┆ 150.111 ┆ 150.091 ┆ 150.096 ┆ 2024-12-07 05:55:00 ┆ 2024-12-06 │
└─────────┴─────────┴─────────┴─────────┴─────────────────────┴────────────┘

RVの実装

前のセクションで述べた結果を実装します。ジャンプの検定に使う有意水準は5%としています。

mu_1 = 2**(1/2) * math.gamma(1) * math.gamma(1/2)**(-1)
mu_4over3 = 2**(2/3) * math.gamma(7/6) * math.gamma(1/2)**(-1)
alpha = 0.95

df_volatility = (
    df
    # 収益率は100倍して%表記にする
    .with_columns(ret=(pl.col("close").log() - pl.col("close").shift(1).log()) * 100)
    .group_by("date")
    .agg(
        n=pl.len(),
        rv=(pl.col("ret")**2).sum(),
        bv=mu_1**(-2) * (pl.col("ret").abs() * pl.col("ret").shift(1).abs()).sum(),
        tq=pl.len() * mu_4over3**(-3) * (pl.col("ret").abs()**(4/3) * pl.col("ret").shift(1).abs()**(4/3) * pl.col("ret").shift(2).abs()**(4/3)).sum(),
    )
    .sort("date")
    .with_columns(
        z=(pl.col("rv").log() - pl.col("bv").log()) / ((mu_1**(-4) + 2 * mu_1**(-2) - 5) * pl.col("tq") * pl.col("bv")**(-2) / pl.col("n"))**(1/2)
    )
    .with_columns(
        j=pl.when(pl.col("z") > sp.stats.norm.ppf(alpha)).then(pl.col("rv") - pl.col("bv")).otherwise(pl.lit(0))
    )
    .with_columns(
        c=pl.col("rv") - pl.col("j")
    )
)

以下のとおり求められました。

shape: (288, 8)
┌────────────┬─────┬──────────┬──────────┬──────────┬───────────┬──────────┬──────────┐
│ date       ┆ n   ┆ rv       ┆ bv       ┆ tq       ┆ z         ┆ j        ┆ c        │
│ ---        ┆ --- ┆ ---      ┆ ---      ┆ ---      ┆ ---       ┆ ---      ┆ ---      │
│ date       ┆ u32 ┆ f64      ┆ f64      ┆ f64      ┆ f64       ┆ f64      ┆ f64      │
╞════════════╪═════╪══════════╪══════════╪══════════╪═══════════╪══════════╪══════════╡
│ 2023-10-30 ┆ 276 ┆ 0.199465 ┆ 0.094773 ┆ 0.012525 ┆ 13.415503 ┆ 0.104692 ┆ 0.094773 │
│ 2023-10-31 ┆ 288 ┆ 0.374825 ┆ 0.30842  ┆ 0.333002 ┆ 2.266396  ┆ 0.066405 ┆ 0.30842  │
│ 2023-11-01 ┆ 288 ┆ 0.194893 ┆ 0.19418  ┆ 0.086115 ┆ 0.052693  ┆ 0.0      ┆ 0.194893 │
│ …          ┆ …   ┆ …        ┆ …        ┆ …        ┆ …         ┆ …        ┆ …        │
│ 2024-12-04 ┆ 288 ┆ 0.487489 ┆ 0.446522 ┆ 0.22751  ┆ 1.786992  ┆ 0.040967 ┆ 0.446522 │
│ 2024-12-05 ┆ 288 ┆ 0.471903 ┆ 0.458031 ┆ 0.306132 ┆ 0.537139  ┆ 0.0      ┆ 0.471903 │
│ 2024-12-06 ┆ 288 ┆ 0.486899 ┆ 0.455146 ┆ 0.425198 ┆ 1.023669  ┆ 0.0      ┆ 0.486899 │
└────────────┴─────┴──────────┴──────────┴──────────┴───────────┴──────────┴──────────┘

polarsはめちゃくちゃ読みやすいですね。polarsは大きなDataFrameを高速に処理できることがメリットとしてよく言われますが、認知負荷にやさしい構文なのが一番好きなポイントです。Rのtidyverseに近い書き方なのでtidyverseで育ったわたしにとっても使いやすいです。

プロット

RVは%表記です。積み上げ面グラフで連続部分$C_{t}$とジャンプ部分$J_{t}$に分けています。

p1 = (
    pn.ggplot(
        df_volatility
        .select("date", "j", "c")
        .unpivot(index="date", variable_name="variable", value_name="value"),
    )
    +pn.geom_area(pn.aes("date", "value", fill="variable"), color="gray", size=0.2, alpha=0.8)
    +pn.scale_fill_brewer(type="qual", palette="Set2")
    +pn.scale_x_date(date_labels="%Y/%m")
    +pn.scale_y_continuous(breaks=range(0, 100, 1))
    +pn.theme_minimal()
    +pn.labs(x="date", y="volatility", title="realized volatility", subtitle="decomposed into continuous component (c) and jump component (j)")
    +pn.theme(figure_size=(8, 3), dpi=200, legend_position="right")
)

cとjは必ずしも同じような動きをしていないのが面白いです。

さて、ジャンプ部分jの値が一番大きい日は2024/7/11、次に大きい日は2024/4/29でした。

shape: (288, 8)
┌────────────┬─────┬──────────┬──────────┬───────────┬──────────┬──────────┬──────────┐
│ date       ┆ n   ┆ rv       ┆ bv       ┆ tq        ┆ z        ┆ j        ┆ c        │
│ ---        ┆ --- ┆ ---      ┆ ---      ┆ ---       ┆ ---      ┆ ---      ┆ ---      │
│ date       ┆ u32 ┆ f64      ┆ f64      ┆ f64       ┆ f64      ┆ f64      ┆ f64      │
╞════════════╪═════╪══════════╪══════════╪═══════════╪══════════╪══════════╪══════════╡
│ 2024-07-11 ┆ 288 ┆ 2.826537 ┆ 1.161507 ┆ 26.76854  ┆ 4.341744 ┆ 1.66503  ┆ 1.161507 │
│ 2024-04-29 ┆ 276 ┆ 5.346737 ┆ 3.789614 ┆ 85.518574 ┆ 3.002976 ┆ 1.557123 ┆ 3.789614 │
│ 2024-09-27 ┆ 288 ┆ 2.879915 ┆ 1.608273 ┆ 19.819779 ┆ 4.57689  ┆ 1.271642 ┆ 1.608273 │
│ …          ┆ …   ┆ …        ┆ …        ┆ …         ┆ …        ┆ …        ┆ …        │
│ 2024-12-03 ┆ 288 ┆ 0.588244 ┆ 0.561047 ┆ 0.635917  ┆ 0.724247 ┆ 0.0      ┆ 0.588244 │
│ 2024-12-05 ┆ 288 ┆ 0.471903 ┆ 0.458031 ┆ 0.306132  ┆ 0.537139 ┆ 0.0      ┆ 0.471903 │
│ 2024-12-06 ┆ 288 ┆ 0.486899 ┆ 0.455146 ┆ 0.425198  ┆ 1.023669 ┆ 0.0      ┆ 0.486899 │
└────────────┴─────┴──────────┴──────────┴───────────┴──────────┴──────────┴──────────┘

この2日は何があったのでしょうか？

2024/7/11のドル円チャートはこちらです。

2024/4/29はこちらです。

気付きましたか？この2日は為替介入が行われた日なんですね。為替介入であれば当然、ジャンプのような非連続な価格変動が起こります。

連続な過程ではなくジャンプ拡散過程でモデリングするのがより適切であること、そして一般にイレギュラーな要因で生まれるジャンプ部分のボラティリティを分離して把握することが大切なことが分かるよい例でした。

おわりに

RVの理論と実装を簡単に解説してみました。

金融市場の解析でも機械学習や深層学習の話題が多いですが、数理ファイナンスのトピックや論文も知っておくといろいろ使えることが多いと思います。

参考文献

• Andersen, T. and Bollerslev, T. (1998). Answering the Skeptics: Yes, Standard Volatility Models do Provide Accurate Forecasts. International Economic Review, 39(4), 885-905.
• Barndorff-Nielsen, O. E., and Shephard, N. (2004). Power and Bipower Variation with Stochastic Volatility and Jumps, Journal of Financial Econometrics, 2, 1-48.
• Barndorff-Nielsen, O. E., and Shephard, N. (2006). Econometrics of Testing for Jumps in Financial Economics Using Bipower Variation, Journal of Financial Econometrics, 4, 1-30.
• Hansen, P. R., and Lunde, A. (2005). Realized Variance and Market Microstructure Noise. Journal of Business & Economic Statistics, 24, 127-161.
• Liu, L. Y., Patton, A. J., and Sheppard K. (2015). Does anything beat 5-minute RV? A comparison of realized measures across multiple asset classes. Journal of Econometrics, 187(1), 293-311.