Vasicek モデル

## 概要

金利は期間に対して一定ではなく，**期間構造**（**金利の期間構造**）を持ちます．金利の動きをモデル化することで，利回り曲線（イールドカーブ）として表される期間構造を表現することができ，債券価格の評価やリスク管理，シミュレーション等の実用につなげることができます．

本記事では，証券アナリスト資格試験にも頻出な金利の期間構造モデルのなかで最も基本的な **Vasicek モデル** について整理していこうと思います．目標はイールドカーブ（スポットレート）を導出することです．

## Visicek モデルの定式化

$t$ 時点における期間 $\Delta t$ の短期金利を $r_t$　とします．Vasciek モデルでは将来時点 $t + \Delta t$ の金利 $r_{t + \Delta t}$ が以下に従います：

$$
r_{t + \Delta t} = r_t + \kappa \left( \theta - r_t \right) \Delta t + \sigma \Delta z.　\tag{1}
$$

ここで，$\Delta z \sim N \left( 0, \Delta t \right)$, つまり $\Delta z$ は平均 0，分散 $(\Delta t)^2$ の正規分布に従うとします．また，$\theta$ は長期的な金利の回帰水準，$\kappa$ は平均回帰への調整速度， $\sigma$  は金利ボラティリティを表すパラメタです．

式 (1) を確率微分方程式の形に直すと，

$$
dr_t = \kappa (\theta - r_t) dt + \sigma dz_t \tag{2}
$$

となります．ここで，$dr_t$ は微小時間 $dt$ における瞬間短期金利 $r_t$ の変動幅を表し， $dz_t$ は平均ゼロ，分散 $dt$ の系列相関のない正規分布に従う確率変数，つまりウィナー過程です．右辺の第 1 項はトレンド項，第 2 項が攪乱項を表しています．

式 (2) は [Ornstein–Uhlenbeck 過程](https://ja.wikipedia.org/wiki/%E3%82%AA%E3%83%AB%E3%83%B3%E3%82%B7%E3%83%A5%E3%82%BF%E3%82%A4%E3%83%B3%EF%BC%9D%E3%82%A6%E3%83%BC%E3%83%AC%E3%83%B3%E3%83%99%E3%83%83%E3%82%AF%E9%81%8E%E7%A8%8B)　と呼ばれる確率過程であり，解析解を求められることが知られています．

次節でこの解析解を書き下していきます．

## 解析解の導出

### 瞬間短期金利の導出

まず式 (2) の両辺に $e^{\kappa t}$ をかけ，式を整理します．

$$
\begin{align}
&&e^{\kappa t} dr_t = \kappa \theta e^{\kappa t} dt - \kappa r_t e^{\kappa t} dt + \sigma e^{\kappa t} dz_t \\\\\
&\Rightarrow& e^{\kappa t} dr_t + \kappa r_t e^{\kappa t} dt = \kappa \theta e^{\kappa t} dt + \sigma e^{\kappa t} dz_t \\\\
&\Rightarrow& d \left( e^{\kappa t} dr_t \right) = \kappa \theta e^{\kappa t} dt + \sigma e^{\kappa t} dz_t \tag{3}
\end{align}
$$

時刻 $t$ から $s$ まで積分することで，

$$
\begin{align}
e^{\kappa s} r_s - e^{\kappa t} r_t =& \kappa \theta \int_{t}^{s} e^{\kappa s} ds + \sigma \int_{t}^{s} e^{\kappa u} d z_{u} \\\\
=& \theta \left( e^{\kappa s} - e^{\kappa t} \right) + \sigma \int_{t}^{s} e^{\kappa u} d z_{u}
\end{align}
$$

と計算できます．したがって，

$$
r_s = e^{- \kappa (s - t)} r_t + \theta \left( 1 - e^{- \kappa (s - t)} \right) + \sigma e^{- \kappa s} \int_{s}^{t} e^{\kappa u} d z_{u} \tag{4}
$$

と導くことができ，$r_s$ も正規分布に従うことが分かります．

では実際に，$r_s$ のモーメント（平均と分散）を確認してみましょう．1 次のモーメントは

$$
\begin{align}
\mathbb{E} \left[ r_s \right] =& e^{- \kappa (s - t)} r_t + \theta \left( 1 - e^{- \kappa (s - t)} \right)+ \sigma e^{- \kappa s} \mathbb{E} \left[ \int_{t}^{s} e^{\kappa u} d z_{u} \right] \\\\
=& e^{- \kappa (s - t)} r_t + \theta \left( 1 - e^{- \kappa (s - t)} \right) 
\end{align}
$$

と求まります．また，

$$
\begin{align}
\mathbb{E} \left[ r_s^2 \right] =& \left\\{ e^{- \kappa (s - t)} r_t + \theta \left( 1 - e^{- \kappa (s - t)} \right) \right\\}^2  + \sigma^2 e^{- 2 \kappa s} \mathbb{E} \left[ \left( \int_{t}^{s} e^{\kappa u} d z_{u} \right)^2 \right] \\\\
=&  {\mathbb{E} \left[ r_s \right]}^2  +  \sigma^2 e^{- 2 \kappa s} \mathbb{E} \left[ \int_{t}^{s} e^{2 \kappa u} du \right] \qquad \text{(伊藤の等長定理より)} \\\\
=&  {\mathbb{E} \left[ r_s \right]}^2  + \frac{\sigma^2 }{2 \kappa} \left( 1 - e^{- 2 \kappa (s - t)} \right)
\end{align}
$$

より，2 次のモーメントは

$$
\mathrm{Var} \left[ r_t \right] =  \frac{\sigma^2 }{2 \kappa} \left( 1 - e^{- 2 \kappa (s - t)} \right)
$$

となります．以上より，

$$
r_t \sim N \left( e^{- \kappa (s - t)} r_t + \theta \left( 1 - e^{- \kappa (s - t)} \right) , \frac{\sigma^2 }{2 \kappa} \left( 1 - e^{- 2 \kappa (s - t)} \right) \right) \tag{5}
$$

となります．

### ゼロクーポン債の価格の導出

$r_t$ を時刻 $t$ から満期 $T$ までの瞬間短期金利の積分を $R \left( t, T \right)$ とします．

$$
\begin{align}
R \left( t, T \right) &= \int^{T}\_{t} r_s ds \\\\
&= \int^{T}\_{t} \left[ e^{- \kappa (s - t)} r_t + \theta \left( 1 - e^{- \kappa (s - t)} \right) + \sigma e^{- \kappa (s - t)} \int_{t}^{s} e^{\kappa (u - t)} d z_{u} \right] ds \\\\
&= \theta \left( T - t \right) + \left( r_t - \theta \right) \frac{1 - e^{- \kappa (T - t)}}{\kappa} + \sigma \int^{T}\_{t} \int_{t}^{s} e^{- \kappa (s - t)}  e^{\kappa (u - t)} d z_{u} ds
\end{align}
$$

であり，第 3 項はFubini の定理より積分順序をいれかえて，

$$
\begin{align}
\sigma \int^{T}\_{t} \int_{t}^{s} e^{- \kappa (s - t)}  e^{\kappa (u - t)} d z_{u} ds &= \sigma \int_{t}^{T} \left[ \int_{u}^{T} e^{- \kappa (s - u)} ds \right] dz_u \\\\
&= \sigma \int_{t}^{T} \frac{1 - e^{- \kappa (T - u)}}{\kappa} dz_u
\end{align}
$$

となり平均ゼロの伊藤積分項になります．これまでをまとめると，

$$
\begin{align}
& R(t, T) = M(t, T) + \sigma \int_{t}^{T} B(u, T) dz_u \tag{6} \\\\
& M(t, T) = \theta \left( T - t \right) + \left( r_t - \theta \right) B(t, T) \tag{7} \\\\
& B(t, T) = \frac{1 - e^{- \kappa (T - t)}}{\kappa} \tag{8}
\end{align}
$$

となり，$R (t, T)$ も正規分布にしたがいます．

満期 $T$ のゼロクーポン債の$t$ 時点における価格を $P (t, T)$ とします：

$$
P (t, T) = \mathbb{E} \left[ e^{- R (t, T)} \right].
$$

$R (t, T)$ は正規分布に従い，$e^{- R (t, T)}$ は対数正規分布にしたがうので

$$
P (t, T)　= \exp \left( - \mathbb{E} \left[ R(t, T) \right] + \frac{1}{2} \mathrm{Var} \left[ R (t, T) \right] \right) \tag{9}
$$

が成立します．

また，$R(t, T)$ の平均と分散はそれぞれ

$$
\mathbb{E} \left[ R (t, T) \right] = M(t, T) = - \theta \left( B(t, T) - (T - t) \right) + B(t, T) r_t \tag{10}
$$

$$
\begin{align}
\mathrm{Var} \left[ R(t, T) \right] &= \sigma^2 \mathbb{E} \left[ \left( \int_{t}^{T} B(u, T) dz_u \right)^2 \right] \\\\
&= \sigma^2 \mathbb{E} \left[ \int_{t}^{T} {B(u, T)}^2 du \right] \qquad \text{(伊藤の等長定理より
)}\\\\
&= \frac{\sigma^2}{\kappa^2} \int_{t}^{T} \left[ 1- 2 e^{- \kappa (T - u)} + e^{ - 2 \kappa (T - u)} \right] du \\\\
&= \frac{\sigma^2}{\kappa^2} \left[ (T - t) - \frac{2}{\kappa} \left( 1 - e^{- \kappa (T - t)} \right) + \frac{1}{2 \kappa} \left( 1 - e^{- 2 \kappa (T - t)} \right) \right] \\\\
&= - \frac{\sigma^2}{\kappa^2} \left( B(t, T) - (T - t) \right) - \frac{\sigma^2}{2 \kappa} {B(t, T)}^2 \tag{11}
\end{align}
$$
と導出できます．

式 (10) と (11) を式 (9) に代入して，

$$
\begin{align}
& P(t, T) = A (t, T) e^{- B(t, T ) r_t} \tag{12} \\\\
& A(t, T) = \exp \left[ \left( \theta - \frac{\sigma^2}{2 \kappa^2} \right) \left( B(t, T) - (T - t) \right) - \frac{\sigma^2}{4 \kappa} {B(t, T)}^2 \right] \tag{13} \\\\
& B(t, T) = \frac{1 - e^{- \kappa (T - t)}}{\kappa} \tag{8}
\end{align}
$$

と計算できます．これがよく見る Vasicek モデルによって導出されるゼロクーポン債の価格式です．

### イールドカーブの導出

ゼロクーポン債の価格が導けたので，イールドカープも明示的に表すことができます．

$$
\begin{align}
Y (t, T) &= - \frac{1}{T - t} \ln P(t, T) \\\\
&= - \frac{1}{T - t} \ln A(t, T) + \frac{1}{T - t} B(t, T) r_t
\end{align} \tag{14}
$$

## イールドカーブのパラメタ依存性

パラメタ $\theta, \kappa, \sigma$ を変化させイールドカーブのパラメタ依存性を確認していきます．

### python コード
イールドカーブは下記のコードより描きました．

```python
from typing import TypedDict

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.figure import Figure


class VasicekParameter(TypedDict):
    kappa: float
    theta: float
    sigma: float


class Vasicek:
    @staticmethod
    def B(*, kappa: float, tau: float) -> float:
        """
        Calculate the Vasicek B function.

        Args:
            kappa (float): The mean reversion speed.
            tau (float): The time to maturity.
        """
        return (1 - np.exp(-kappa * tau)) / kappa

    @staticmethod
    def logA(*, kappa: float, theta: float, sigma: float, tau: float) -> float:
        """
        Calculate the log of the Vasicek A function.

        Args:
            kappa (float): The mean reversion speed.
            theta (float): The long-term mean level.
            sigma (float): The volatility of the process.
            tau (float): The time to maturity.
        """
        b_tau = Vasicek.B(kappa=kappa, tau=tau)
        return (
            (theta - sigma**2 / (2 * kappa**2)) * (b_tau - tau)
            - (sigma**2 / (4 * kappa)) * (b_tau**2)
        )

    @staticmethod
    def A(*, kappa: float, theta: float, sigma: float, tau: float) -> float:
        """
        Calculate the Vasicek A function.

        Args:
            kappa (float): The mean reversion speed.
            theta (float): The long-term mean level.
            sigma (float): The volatility of the process.
            tau (float): The time to maturity.
        """
        return np.exp(Vasicek.logA(kappa=kappa, theta=theta, sigma=sigma, tau=tau))

    @classmethod
    def price(
        cls, *, kappa: float, theta: float, sigma: float, tau: float, r0: float
    ) -> float:
        """
        Calculate the price of a zero-coupon bond using the Vasicek model.

        Args:
            kappa (float): The mean reversion speed.
            theta (float): The long-term mean level.
            sigma (float): The volatility of the process.
            tau (float): Time to maturity.
            r0 (float): Current short rate.
        """
        b_tau = cls.B(kappa=kappa, tau=tau)
        a_tau = cls.A(kappa=kappa, theta=theta, sigma=sigma, tau=tau)
        return a_tau * np.exp(-b_tau * r0)

    @classmethod
    def spot_rate(
        cls, *, kappa: float, theta: float, sigma: float, tau: float, r0: float
    ) -> float:
        """
        Calculate the spot rate using the Vasicek model.

        Args:
            kappa (float): The mean reversion speed.
            theta (float): The long-term mean level.
            sigma (float): The volatility of the process.
            tau (float): Time to maturity.
            r0 (float): Current short rate.
        """
        b_tau = cls.B(kappa=kappa, tau=tau)
        loga_tau = cls.logA(kappa=kappa, theta=theta, sigma=sigma, tau=tau)
        return (b_tau * r0 - loga_tau) / tau

    def get_yield_curve(
        cls,
        *,
        r0: float,
        maturities: np.ndarray,
        conditions: list[VasicekParameter],
        label_template: str = "kappa={kappa}, theta={theta}, sigma={sigma}",
        is_percent: bool = True,
    ) -> Figure:
        """
        Draw the yield curve using the Vasicek model.

        Args:
            r0 (float): Current short rate.
            maturities (np.ndarray): Array of maturities.
            conditions (list[VasicekParameter]): List of conditions for the Vasicek model.
            label_template (str): Template for the label of each curve.
            is_percent (bool): Whether to display the rates in percentage.
        """
        fig, ax = plt.subplots()
        for condition in conditions:
            kappa = condition["kappa"]
            theta = condition["theta"]
            sigma = condition["sigma"]
            rates = [
                cls.spot_rate(kappa=kappa, theta=theta, sigma=sigma, tau=tau, r0=r0)
                for tau in maturities
            ]
            if is_percent:
                rates = [rate * 100 for rate in rates]
            ax.plot(maturities, rates, label=label_template.format(**condition))
        ax.set_xlabel("Maturity (years)")
        ylabel = "Spot Rate (%)" if is_percent else "Spot Rate"
        ax.set_ylabel(ylabel)
        ax.legend()
        return fig
```

### θ 依存性

現在の短期金利 $r_t = 2.00 %$ として 長期回帰水準 $\theta$ 依存性を確認してみましょう．
以降，$\theta = 3.0 \\%, \kappa = 0.3, \sigma = 0.5 \\%$ を基準として各パラメタを変化させていきます．

![vasicek_yield_curve_theta.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3560007/50acb4a1-6aa7-4fbe-a1f2-d7e8a3c3f0f1.png)

図から確認できる通り，$\theta$ は長期金利水準に主に影響し， $\theta$ が高い（低い）ほど長期金利水準が高く（低く）なります．

### κ 依存性

$\kappa$ は短期金利が長期金利水準に収束するためのスピードを表しており，$\theta > rt$ の場合は $\kappa$ が大きいほど中期金利の水準を押しあげます．

![vasicek_yield_curve_kappa.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3560007/5786a381-cf1a-4e93-85f5-e413045a9d27.png)

### σ 依存性

金利ボラティリティの水準 $\sigma$ が上昇すると，長期金利を下方に押し下げようとする力が強くなることがわかります．これは，コンベクシティ効果によるものです．

![vasicek_yield_curve_sigma.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3560007/7313e48e-20c9-4cc9-84a8-c58d6ac91b6b.png)

## まとめ

本記事では Vasicek モデルによりスポットレート（イールドカーブ）の解析解を導出し，その挙動をコーディングにより確認しました．