定数ボラティリティのHJMモデルを最尤法で解く

（2024/4/24 節分けと文言を微修正）

# 0. はじめに

## 0-1. HJMモデルって何？🤔

金融工学におけるHeath-Jarrow-Mortonモデル（以下HJMモデル）とは金利の期間構造を確率微分方程式で記述したモデルである．ここでの金利は観測時点$t$における満期$T$のフォワードレート$ f(t,T) $を指す．イメージとして**金利の時間発展はトレンドとホワイトノイズで表現される**というものであり，数式では
```math
\begin{align}
     df(t,T) = 
    \begin{pmatrix}
        \alpha(t,T) \quad  \sigma (t,T)
    \end{pmatrix}
    \begin{pmatrix}
        dt \\ dW^{P}(t)
    \end{pmatrix}
\end{align}
```

と表現される．ここで$\alpha$はドリフト項であり，$\sigma$はボラティリティ項である．上式は観測時点$t$におけるさまざま満期$T$の金利を包括的に表現しており，したがってHJMモデルは$t$時点の**イールドカーブ全体をモデリングしている**と解釈できる．

## 0-2. 前提知識🏫
リスク管理クオンツで1-3年目の方を読者として想定している．具体的には以下の知識を前提とする．参考書籍についても記す．

- 金融工学 [[村上] マルチンゲールアプローチ入門](https://www.kindaikagaku.co.jp/book_list/detail/9784764904804/)
    - 資産価格の基本定理
    - 伊藤の公式
    - 測度変換（Girsanovの定理）
- 統計学 [[ビショップ] パターン認識と機械学習 上](https://www.maruzen-publishing.co.jp/item/b294524.html)
    - 正規分布のパラメータに対する最尤法

## 0-3. この記事を読めば📖

- わかる😁
    - HJMモデルにおけるドリフトのボラティリティによる表現
    - 定数ボラティリティを仮定したときの最尤推定量
- わからない😢
    - さまざまな関数形のボラティリティ
    - マルチファクターHJMモデル

# 1. ドリフトをボラティリティで表す：無裁定

無裁定条件を仮定すればHJMモデルのドリフト$\alpha(t,T)$はボラティリティ$\sigma(t,T)$で表される：
$$ \alpha(t,T) \approx \sigma(t,T) \int^T_t \sigma (t,u) \ du. $$
よってボラティリティを学習させればドリフトは自然に定まり，これがパラメータ学習の実装をシンプルにしている．$\approx$となっている理由は確率微分方程式を駆動する**ブラウン運動が異なる世界**の間で同一視しているからである．しかしリスク管理の上では**この違いを意識する必要はない．** なぜなら無裁定条件から一物一価の法則が得られ，確率測度$P$と確率測度$Q$で同じ時価が得られるからである．

以下，証明を載せる．

## 1-1. 証明

主張は以下となる：
> 瞬間フォワードレート$f(t,T)$が確率微分方程式
```math
\begin{align}
     df(t,T) = 
    \begin{pmatrix}
        \alpha(t,T) \quad  \sigma (t,T)
    \end{pmatrix}
    \begin{pmatrix}
        dt \\ dW^{P}(t)
    \end{pmatrix}
\end{align}
```
に従っているとする．市場が無裁定であるならば，実確率$P$と同値な測度$Q$が存在し，測度$Q$のもとで$f(t,T)$は
```math
\begin{align}
     df(t,T) = 
    \begin{pmatrix}
        \sigma(t,T) \int^T_t \sigma (t,u) \ du \quad \sigma (t,T)
    \end{pmatrix}
    \begin{pmatrix}
        dt \\ dW^{Q}(t)
    \end{pmatrix}
\end{align}
```
という確率微分方程式に従う．

<details>
<summary>（クリックで証明を展開）</summary>

資産として満期$T$の債券$P(t,T)$全体と銀行預金$B(t)$を考える．無裁定という仮定から，ある同値マルチンゲール測度$Q$のもとで**すべての$T$に対して**相対価格$P(t,T)/B(t)$がマルチンゲールとなる．$P(t,T)/B(t)$がマルチンゲールである事実を利用するために，まずこれを$\alpha$と$\sigma$で表現する．今確率微分法手式の形で瞬間フォワードレート（金利の）時間発展に関する情報が手元にある．ここで金利の議論に持ち込むため，$P(t,T)/B(t)$の対数を取って伊藤の公式で計算していく：
```math
\begin{align}
    \frac{d\frac{P(t,T)}{B(t)}}{\frac{P(t,T)}{B(t)}}
    & = d \log \frac{P(t,T)}{B(t)} + \frac{1}{2}\ d\left\langle\log \frac{P(t,T)}{B(t)}\right\rangle.
\end{align}
```
確率微分の線形性から$d \log \left(P(t,T)/B(t)\right) =d \log P(t,T) - d \log B(t) $となるため
それぞれを計算していく．第一項については
```math
\begin{align}
    d \log P(t,T)
    & = - d \left(\int^T_t f(t,u)\ du\right) \\
    & = f(t,t)\ dt - \int^T_t df(t,u)\ du \\
    & = r(t) \ dt - \int^T_t \left(\alpha(t,T) \ dt +\sigma (t,T) \ dW^{P}_t\right) \ du \\
    & = r(t) \ dt - \left(\int^T_t\alpha(t,u)\ du\right)\ dt - \left(\int^T_t\sigma(t,u)\ du\right)\ dW^P_t \\
\end{align}
```
となり，第二項について$d \log B(t)=r(t)dt$となるので
```math
\begin{align}
    d \log \frac{P(t,T)}{B(t)}
    &= - \left(\int^T_t\alpha(t,u)\ du\right)\ dt - \left(\int^T_t\sigma(t,u)\ du\right)\ dW^P_t \\
    & = - \left(\int^T_t\alpha(t,u)\ du\right)\ dt - \Sigma(t,T)\ dW^P_t
\end{align}
```
となる．ただしボラティリティを$\Sigma$とおいた：
$$\Sigma(t,T) := \int^T_t \sigma(t,u) \ du. $$
さらに二次変分については
```math
d \left \langle \log \frac{P(t,T)}{B(t)} \right\rangle = \left(\Sigma(t,T)\right) ^2 \ dt
```
と計算される．最終的に相対価格$P(t,T)/B(t)$が$\alpha$と$\sigma$で表現される：
```math
\begin{align}
    \frac{d\frac{P(t,T)}{B(t)}}{\frac{P(t,T)}{B(t)}}
    &= - \left(\int^T_t\alpha(t,u)\ du + \frac{1}{2}\left(\Sigma(t,T)\right) ^2\right)\ dt - \Sigma(t,T)\ dW^P_t \\
    & = -\Sigma(t,T)\left\{\left(\frac{-\int^T_t\alpha(t,u)\ du+\frac{1}{2}\left(\Sigma(t,T)\right)^2}{-\Sigma(t,T)}\right)\ dt + dW_t^P\right\} \\
    & = -\Sigma(t,T)\left(\lambda(t,T)\ dt + dW_t^P\right).
\end{align}
```
ただしリスクの市場価格を$\lambda(t,T)$とおいた：
```math
\lambda(t,T) = \frac{-\int^T_t\alpha(t,u)\ du+\frac{1}{2}\left(\Sigma(t,T)\right)^2}{-\Sigma(t,T)}.
```
満期$T$の債券$P(t,T)$全体と銀行預金$B(t)$が無裁定という仮定から**リスクの市場価格$ \lambda(t,T) $は満期$T$によらず定まっている**はずである．よって単に$\lambda(t)$と表すことができ，これによって微分が計算される：
```math
\begin{align}
    \int^T_t\alpha(t,u)\ du
    &= \lambda(t) \Sigma(t,T) + \frac{1}{2}\left(\Sigma(t,T)\right)^2, \\
    \alpha(t,T)
    & = \lambda(t)\frac{\partial}{\partial T}\int^T_t \sigma(t,u) \ du+ \frac{1}{2} \frac{\partial}{\partial T} \left(\int^T_t \sigma(t,u) \ du\right)^2 \\
    & = \lambda(t)\sigma(t,T)+\sigma(t,T)\Sigma(t,T) \\
    & = \sigma(t,T)\left(\lambda(t)+\Sigma(t,T)\right).
\end{align}
```
$\Sigma$の定義に戻れば無裁定条件による表現
```math
\begin{align}
     df(t,T) = 
    \begin{pmatrix}
        \sigma(t,T) \int^T_t \sigma (t,u) \ du \quad \sigma (t,T)
    \end{pmatrix}
    \begin{pmatrix}
        dt \\ dW^{Q}(t)
    \end{pmatrix}
\end{align}
```
が得られる．ドリフトがボラティリティで表された．

</details>

# 2. ボラティリティをキャリブレーションする

## 2-1. 関数形を決める：Absolute Volatility Model

次は金利のヒストリカルデータからボラティリティを学習させる．まずボラティリティの最も簡単な形として定数$\sigma$が考えられる．他には年限$T-t$を反映するため$ \sigma e^{-a(T-t)} $という形なども考えられる．両方のボラティリティ形を導入するために，駆動するブラウン運動も複数用意したマルチファクターHJMモデルを導入することもできるが今回は前者のみの形，すなわちフォワードレートが
```math
\begin{align}
     df(t,T) = 
    \begin{pmatrix}
        \sigma^2 (T-t) \quad \sigma
    \end{pmatrix}
    \begin{pmatrix}
        dt \\ dW^{Q}(t)
    \end{pmatrix}
\end{align}
```
という形で表されたとする．

以下，$ \sigma $をキャリブレーションしていく．

## 2-2. キャリブレーションする：最尤法

フォワードレートの変動データ
```math
\begin{align}
    \mathcal{D} = \left\{\Delta f(t,T)\right\}_{t,T}
\end{align}
```
が**何かしらの方法で手に入ったとする．**（一般にフォワードレートはマーケットにないため，何かしらの方法で算出する必要があるが今回は割愛する．）満期$T$のデータだけをまとめて
```math
\begin{align}
    \mathcal{D}_T = \left\{\Delta f(t,T)\right\}_{t \leq T}
\end{align}
```
とするとその尤度は
```math
\begin{align}
    L(\mathcal{D}_ T \mid \sigma)
    & = \Pi_{\mathcal{D}_T} \mathcal{N}\left( \Delta f(t,T) \middle | \sigma^2(T-t), \sigma \right) \\
    & = \Pi_{\mathcal{D}_T} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left(-\frac{1}{2} \left(\frac{\Delta f(t,T)-\sigma^2(T-t)}{\sigma}\right)^2\right)
\end{align}
```
となる．ここで尤度$L(\mathcal{D}_T \mid \sigma)$とはパラメータの値が$\sigma$であるときにデータ$\mathcal{D}_T$が観測される確率であり，最尤推定法とは$L(\mathcal{D}_T \mid \sigma)$を最大化する$\sigma$，すなわち観測データ$\mathcal{D}$が最頻値となっているような$\sigma$を推定量とする手法を指す．$L(\mathcal{D}_T \mid \sigma)$を最大化する$ \sigma $は対数尤度$\log L(\mathcal{D}_T \mid \sigma)$をも最大化するため，
```math
\begin{align}
    \log L(\mathcal{D}_T \mid \sigma)
    & = \Sigma_{\mathcal{D}_T} \left\{ -\frac{1}{2} \log ( 2 \pi \sigma^2 ) - \frac{1}{2} \left(\frac{\Delta f(t,T)-\sigma^2(T-t)}{\sigma}\right)^2 \right\} \\
    & = \Sigma_{\mathcal{D}_T} \left\{-\frac{1}{2} \log (2 \pi \sigma^2) - \frac{1}{2} \frac{\Delta f(t,T)^2}{\sigma^2} + \Delta f(t,T)(T-t) - \frac{1}{2}\sigma^2(T-t)^2 \right\}
\end{align}
```
を計算する．最大値を与える$\sigma$を求めるため微分する．ここでテクニカルだが$\sigma^2$で微分して計算を簡単にする：
```math
\begin{align}
    \frac{d\log L(\mathcal{D}_T \mid \sigma)}{d\sigma^2}
    & = \Sigma_{\mathcal{D}_T} \left(-\frac{1}{2\sigma^2} + \frac{1}{2} \frac{\Delta f(t,T)^2}{\sigma^4} - \frac{1}{2}(T-t)^2 \right)
\end{align}
```
ここで求めたのはあくまで$\mathcal{D}_T$の対数尤度であるので，$\mathcal{D}$の尤度を求めるために$T$について対数尤度を足していく：
```math
\begin{align}
    \Sigma_T\frac{d\log L(\mathcal{D}_T \mid \sigma)}{d\sigma^2}
    & = \Sigma_T\left\{\Sigma_{\mathcal{D}_T} \left(-\frac{1}{2\sigma^2} + \frac{1}{2} \frac{\Delta f(t,T)^2}{\sigma^4} - \frac{1}{2}(T-t)^2\right)\right\}.
\end{align}
```
微分が$0$となるような$\sigma$を求めるために$\sigma^2$に対する二次方程式を解いていく．（なおこのような$\sigma$は最大で2つ存在し，その場合はそれぞれについて尤度を判定していくことになる．また$\sigma$が存在しなかった場合はほかの手法が必要となってくる．）

以上で最尤法によりキャリブレーションすることができた．

# 3. ここまででわかったこと😁
- HJMモデルはフォワードレートをドリフトとホワイトノイズで表現するモデルである
- 適当な測度変換によりドリフトをボラティリティで記述することができる
- ボラティリティの関数形として最もシンプルな形を仮定すると最尤法で解くことができる


# -1. おわりに🙇

私は金融機関でリスク管理業務をしており，この度社内の勉強会で運用チームやシステム調達チーム向けに金利モデルについて発表することになった．しかし社内のモデルは実用に耐えうる複雑なモデルとなっており，それをそのまま発表しては参加者の理解が進まないと感じた．そこで定数ボラティリティのみの1ファクターモデルについて簡単な解法を発表することで，まずはHJMモデルに慣れてもらおうと考えている．この記事はそのための資料であり，社内に限るものでもないため一般に公開した．本記事が何かの助けとなれば幸いである．