アービトラージ入門〜HJB方程式を用いた最適投資比率の算出〜

#本記事の内容
初めまして、bitcoinの自動売買を趣味で行っているwillowと申します。
Qiitaにはいつも色々教えていただいているので、お返しの意味で最近勉強した内容で記事を作りました。初めての記事となるため、無作法などあるかも知れませんが、ご容赦ください。
内容としては、bitcoinの取引所間スプレッドの変動過程をモデル化して、HJB方程式を解くことで最適な投資比率で導きます。HJB方程式は、ヘッジファンドなどではよく使われている（と言われている）式ですが、実用性の割にブラックショールズ方程式などと比べて知名度が低いような気がします。よって今回は「bitcoinの取引所間スプレッド」を題材に、HJB方程式を投資でどう役立てるのかご紹介しながら、一緒に学べたらと思います。なお、実データは扱わず、実際に取引ができるプログラムなどはつけませんので、その点ご承知ください。
私にとって最も身近な例だったのでbitcoinアービトラージを題材に今回の記事を書いてはおりますが、今回の結果は「平均回帰性を持つと期待される金融商品一般（ロング・ショートのペアトレードポートフォリオなど）」に対して利用でき、また自分で変動過程をモデル化することが出来れば、当記事の考え方をなぞる事で金融商品一般への投資行動に応用することが出来るはずです。

#対象読者
・定量的な手法を用いて投資をしてみたい、またはしている方々に向けて書いているつもりです。
・数学（微積分、確率）と金融工学の基礎知識があると理解しやすいと思います。
・数式が多少分からなくても理解していけるように、出来るだけ日本語の丁寧な説明を心がけます。そのため冗長に感じる所があるかもしれません。

#用語説明
・アービトラージ：当記事では、「同一の本源的価値を持つとされる金融商品（今回の例ではbitcoin現物）について異なる価格がついていた場合に、割高な方を売り、割安な方を買うことで、市場リスクを抑えたままリターンを得る投資手法」を指すものとします。
・取引所間スプレッド：二つの取引所で扱われているbitcoin現物の値段の差を表します。二つの取引所では、全く同じ商品（bitcoin現物）が扱われているので、理想的には同じ値段がつくはずではありますが、実際には様々な要因で異なる値段がついている状態で放置されていることもあります。アービトラージャーは、この取引所間スプレッドの絶対値が大きくなった時に、スプレッドが縮小することに期待して割高な方を売り、割安な方を買う投資行動を行います。
・HJB方程式；Hamilton–Jacobi–Bellman equation。最適制御理論の根幹をなす偏微分方程式で、当記事では、最適なポートフォリオを求めるために利用します。ヘッジファンドでよく使われているという話を聞いたり読んだりします。

#取引所間スプレッドのモデル化
これから、以下のような手順で最適な投資比率を導いていきます。
1.スプレッドと効用関数のモデル化
2.HJB方程式の導出
3.実際に解く！！

まずは、第一ステップとして、スプレッドの変動過程をモデル化します。スプレッドの変動過程には、次のような１次元のOU過程（Ornstein–Uhlenbeck process）を用います。平均回帰性を持っていて、かつ非常にシンプルなモデルです。
$$ dS_t = \kappa_s(F - S_t)dt + \sigma_sdB_t $$
$ S_t$ : 時点tにおけるスプレッドの値
$ \kappa_s$ : スプレッドの平均回帰性の強さ
$ F$ : スプレッドの回帰する先
$ \sigma_s$ : スプレッドのボラティリティ
$ dB_t$ : 時点tのブラウン運動
ここでは、取引所間スプレッドが従う確率微分方程式を仮定しています。ポイントとしては「①スプレッドの大きさ(F - S)に比例する（スプレッドが大きくなると強くなる）平均回帰の力」と「②ランダムな動き」の合算でスプレッドの変動が決まることです。

ここで、bitcoinの本源的価値の差は0であるはずなので、$F=0$となるはずです。よって、式は以下のようになります。
$$ dS_t = -\kappa_sS_tdt + \sigma_sdB_t $$
投資家は、時点tで資産を$X_t$持っていて、その内リスク資産（ここでは、スプレッド）：無リスク資産（ここでは現金）=$\pi_t : (1-\pi_t)$で投資するとします。この時、投資家の持つ資産の変動過程は以下のように表すことが出来ます。投資比率$\pi_t$は、スプレッドの符号により正にも負にもなります。
$$ dX_t = \pi_tX_t dS_t + (1-\pi_t)X_trdt $$
$ r$ : リスクフリーレート
今回は、簡単のため$r=0$とします。すると、スプレッドの変動過程を代入して以下のように式を変形できます。
$$ dX_t = -\pi_tX_t\kappa_sS_tdt + \pi_tX_t\sigma_sdB_t $$

この設定のもとで、投資家は投資満期Tにおける資産の効用が最大になるように、以下のような期待効用最大化の問題を解くこととなります。ここで、投資家の初期資産$X_t = x$とします。
$$\sup_{\pi\in \Lambda(x)}E_t\biggl [D_{tT}U(X_T)\biggl]$$
$$X_t = x$$
$ \Lambda(x)$ : 最初の資金xにおいて、許容可能な投資$\pi$の集合
$ D_{ts} = e^{-\delta (s - t)}$ : 時間選好に関する割引率
$U(X_T)$ : 資産$X_T$についての効用関数
効用関数にはいくつか例がありますが、今回はリスク回避度を$\gamma$においた以下のべき型効用関数を使うこととします。
$U(x) = \frac{x^{1-\gamma}}{1-\gamma},　\gamma > 0 \neq 1$
上の式のように投資家は、「「目標時点Tにおいて保有している資産の効用」の現時点における期待値」を最大化するように投資率$\pi$をコントロールしながら、満期Tに近づいていきます。ここで、再度簡単のために、時間割引率を$\delta=0, D = 1$としましょう。あまり長い式にはしたくありませんし、パラメータが少ない方が投資は楽です。Dに代入して、以下の式を得ます。
$$\sup_{\pi\in  \Lambda(x)}E_t\biggl [U(X_T)\biggl]$$
$$X_t = x$$

#HJB方程式の導出
ここでは、先の章で出てきた「期待効用最大化の式」と、定義した「スプレッドの変動モデル」を合体させて、解くべき方程式をハッキリさせます。
期待効用最大化の式について、時間を離散化して、動的計画原理を適用します。t~Tとしていたのを、$[0 = t_0 < t_1 < , , , < t_N = T]$とします。この時、i = 0, , , N-1に対して求める期待効用は以下の式で表せます。
$$V_{t_{i}} = \sup_{\pi \in A(x)}E_{t_i}\biggl [ U(X_T)\biggl]$$
Vという関数が定義されましたが、これは「初期時点tにおいて最適な投資比率をとった場合の、満期T時点での効用の期待値」です。一期先の$V_{t_{i+1}}$は同様に以下のように表せます。
$$V_{t_{i+1}} = \sup_{\pi \in  \Lambda(x)}E_{t_{i+1}}\biggl [U(X_T)\biggl]$$
$\Delta t = t_{i+1} - t_i$とすると、$V_{t_{i}}$を以下のように書き直せます。
$$V_{t_{i}} = \sup_{\pi \in  \Lambda(x)}E_{t_i}\biggl [U(X_T)\biggl]$$
$$= \sup_{\pi \in  \Lambda(x)}E_{t_i}\biggl [E_{t_{i+1}}\biggl[U(X_T)\biggl]\biggl]$$
$$= \sup_{\pi \in  \Lambda(x)}E_{t_i}\biggl [\sup_{\pi \in  \Lambda(x)}E_{t_{i+1}}\biggl[U(X_T)\biggl]\biggl]$$
結果、$V_{t_i}$の中に$V_{t_{i+1}}$を入れる少しトリッキーな変形をして、以下のようになります。
$$V_{t_{i}} = \sup_{\pi \in  \Lambda(x)}E_{t_i}\biggl [V_{t_{i+1}}\biggl]$$
ここで、左辺を右辺に移項して$\Delta t\rightarrow 0$の極限を取ります。
$$ 0 = \sup_{\pi \in  \Lambda(x)}\frac{1}{dt}E_t\biggl [dV(t,x,s)\biggl]$$
Vは時間t, 初期資産x, スプレッド価格$S_t$（表記の簡略化のためsとします）の変数となり、右辺の期待値の中身は伊藤の公式を使うと、以下のように変形されます。
$$dV(t,x,s) = \frac{\partial V(t,x,s)}{\partial t}dt + \frac{\partial V(t,x,s)}{\partial x}dx + \frac{\partial V(t,x,s)}{\partial s}ds + \frac{1}{2}\frac{\partial^2 V(t,x,s)}{\partial x^2}(dx)^2 + \frac{1}{2}\frac{\partial^2 V(t,x,s)}{\partial s^2}(ds)^2 + \frac{\partial^2 V(t,x,s)}{\partial x \partial s}(dsdx)$$
ここで前章の結果より、以下の式が導かれます。
$dx＝-\pi_t x \kappa_s sdt + \pi_t x\sigma_sdB_t$
$ds = -\kappa_s sdt + \sigma_sdB_t$　
$(dx)^2＝(\pi_t x \sigma_s )^2dt$
$(ds)^2＝\sigma_s^2dt$
$dxds = \pi_tx\sigma_s^2dt$
上の式に代入すると、以下のようになります。$dB_t$の項は期待値をとると0になります。
$$\sup_{\pi \in  \Lambda(x)}E_t\biggl [\frac{\partial V(t,x,s)}{\partial t} - 
\pi_tx\kappa_s s\frac{\partial V(t,x,s)}{\partial x} 
- \kappa_s s \frac{\partial V(t,x,s)}{\partial s} 
+ \frac{1}{2}(\pi_tx\sigma_s)^2\frac{\partial^2 V(t,x,s)}{\partial x^2}
+ \frac{1}{2}\sigma_s^2\frac{\partial^2 V(t,x,s)}{\partial s^2}
+ \pi_tx\sigma_s^2\frac{\partial^2 V(t,x,s)}{\partial x \partial s}\biggl ] = 0$$
$$V(T,x,s) = U(x)$$
これにて、解く対象のHJB方程式が導かれた。ふう、一息。

#解いてみる
先ほど求められた式について、解きましょう。
方程式は、各時点のt、各状態のxについて成り立っている必要があります。よって、この方程式を各時点t各状態xについて解くことで、動的な最適投資比率$ \pi_t^* $を求めることが出来ます。
求められたHJB方程式の$\pi$についての一階の最適性条件より、最適投資比率$ \pi^* $では以下の式を満たします。（先の式の期待値の中身を$\pi_t$で微分して、$ \pi_t = \pi_t^* $を代入しています。）
$$-x\kappa_s s \partial_xV + \pi^* x^2 \sigma_s^2  \partial_{xx}V + \sigma_s^2 x \partial_{xs}V =0$$
$$\pi^*  = \frac{-1}{\sigma_s^2 x^2 \partial_{xx}V}(-x\kappa_s s \partial_xV +  \sigma_s^2 x \partial_{xs}V)$$
Vの偏微分については、簡易表記を使いました。

あとは、Vに形を与えればこの計算を進められます。 V というのがそもそも何だったのか振り返っておくと、「最適な投資比率を選んだ場合の、投資満期における効用の時点tにおける期待値」でした。$V(t,x,s)$について、次のような関数を想定してみましょう。
$V(t,x,s) = U(x)\exp(A(t) +\frac{C(t)}{2}s^2),　A(T) = C(T) = 0$
このような関数系であれば解析的に解を出すことが可能となりますが、仮に解析的に解けない関数を使うのであれば、数値計算で解を求めて行くことになります。A,Cは満期Tで0になるので、tの関数というよりは満期までの時間T-tの関数としてこれが0となる時C(t)も0の値を取る、最適効用の時間割引き的なものと解釈してます。

これを代入すると先の式中の偏微分は、それぞれ以下のような形になります。ここで定義より、$U(x) = \frac{x^{1-\gamma}}{1-\gamma},　\partial_x U(x) = x^{- \gamma}=\frac{1-\gamma}{x}U(x),　\partial_{xx} U(x) = -\gamma x^{- (\gamma + 1) }=-\frac{\gamma(1-\gamma)}{x^2}U(x)$
$\partial_tV = (A'(t) + \frac{C'(t)}{2}s^2)V$
$\partial_sV = C(t)sV$
$\partial_xV = \frac{1-\gamma}{x}V$
$\partial_{xx}V = -\frac{\gamma(1-\gamma)}{x^2}V$
$\partial_{xs}V = \frac{1-\gamma}{x}C(t)sV$
$\partial_{ss}V = (C(t)^2 s^2 + C(t))V$

$ \pi^* $の式に代入すると以下のようになります。（やばい疲れて息切れしてきた）
$$\pi^*  = \frac{-s}{\gamma} (\frac{\kappa_s}{\sigma_s^2 } -  C(t)) $$
まだ解き切ってはいないのですが、この時点で大体の結論が出ています。
最適投資比率$ \pi^* $には、二つの項が存在して、「スプレッドの期待リターン（$-\kappa_s s$）に比例してスプレッドの分散（$\sigma_s^2$）に反比例する第一項（リスクに対するリターンの大きさを表す指標的なものでしょう）」と「時間変化しながら第一項を抑制する第二項（$\gamma > 1 $の領域ではこれは投資満期Tが近づくにつれて投資比率を減少させ、リスクを減らそうとする項です、$\gamma > 1 $では異なる挙動を取ります。）」の合計でできていることがわかります。正直、頑張った割には当たり前のことが出てきた気がします笑。あとは、HJB方程式に代入して、Cの関数を特定出来れば終わりです。（Aは分からなくてもいい）
解き方としては、これに先程の$\pi_t^* $をそのままHJB方程式に代入すると長すぎて沼にハマります。HJB方程式が$D\pi_t^2 + E\pi_t + F$となるようにD,E,Fを置いたのちに式を平方完成し、これが最大値を取る時(D<0に注意)について考えれば、式を幾分かシンプルにできます。
ここで、
$D = \frac{1}{2}(x\sigma_s)^2\partial_{xx} V$
$E = - x \kappa_s s \partial_x V + x \sigma_s^2 \partial_{xs} V$
$F = \partial_tV - \kappa_s s \partial_s V + \frac{1}{2}\sigma_s^2 \partial_{ss} V$

$$D \pi_t^{2} + E \pi_t + F = D(\pi_t + \frac{E}{2D})^2 + F - \frac{E^2}{4D}$$
$\pi_t = -\frac{E}{2D}$の時（当然ですがこれが先の$\pi_t^*$と一致します。）、上式は$F - \frac{E^2}{4D}$となります。次からはこの部分が0と等しくなること（$F = \frac{E^2}{4D}$）について考えましょう。
$$\partial_tV - \kappa_s s \partial_s V + \frac{1}{2}\sigma_s^2 \partial_{ss} V = \frac{1}{2(x\sigma_s)^2\partial_{xx} V} (- x \kappa_s s \partial_x V + x \sigma_s^2 \partial_{xs}V)^2$$
Vの偏微分を代入し、両辺をVで割ります。
$$(A'(t) + \frac{C'(t)}{2}s^2)
 - \kappa_s s C(t)s
 + \frac{1}{2}\sigma_s^2 ((C(t)^2 s^2 + C(t)) = 
\frac{-(1-\gamma)}{2 \sigma_s^2 \gamma} (\sigma_s^2 C(t)- \kappa_s)^2 s^2 $$

この式が、時点tに置いて、あらゆるsで成り立っている必要があります。sの二次式とし、係数を見ていくと、全ての係数が両辺で一致している必要があります。二次の項について、両辺の係数を比べ、整理しました。
$s^2$の項 : $C'(t) + \frac{1}{\gamma}\sigma_s^2 C(t)^2 - \frac{2}{\gamma}\kappa_sC(t) + \frac{1-\gamma}{\sigma_s^2 \gamma} \kappa_s^2 = 0 $

あとは、これを解くだけですね！
解くと、C(T) = 0より、以下のようにC(t)が表せます。
$$C(t) = \frac{k}{\sigma^2} \frac{(1-\gamma) (1 - \exp(-\frac{2k}{\sqrt{\gamma}}(T-t)))}{(1+\sqrt{\gamma}) - (1-\sqrt{\gamma})\exp(-\frac{2k}{\sqrt{\gamma}}(T-t)))}$$
これは、$\gamma > 1$の範囲で、tの増加に対して減少し、t = Tで0になる関数となっています。
先ほどの最適投資比率の式を見返してみて、以下のようにまとめましょう。
$$\pi^*  = \frac{-s}{\gamma} (\frac{\kappa_s}{\sigma_s^2 } -  C(t)) $$

参考までに、sを固定し、満期T=1でtを0から1に動かした際にテキトーなパラメータで計算した$C(t), \pi_t^*$のグラフを載せておきます。なんとなく、満期が近づくにつれて投資比率が小さくなっているのがわかると思います。
$\sigma = 1, s = 1, \gamma = 2, k = 1$
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/567595/0d0c636b-37a0-9bf6-4865-1148492b78db.png)
$s=-1$の場合
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/567595/4569260a-ec7e-a122-321d-5f87f8d02e63.png)

#結論
・bitcoinの取引所間スプレッドの変動過程、投資家の効用関数をモデル化することで、最適投資比率を解析的に求めることが出来ました。
・最適投資配分は「スプレッドの期待リターン/スプレッドの分散に比例する項」と「満期が近づくにつれてリスクを減らそうとする項」のバランスで、各時点各スプレッド価格について決まります。
・よって投資家はその時々でスプレッド価格を観測し、求めた式に入れることで最適投資比率を求め、連続的にリバランスしていくことが求められます。

#今後の改良点
・今回は、売買手数料およびスプレッドが存在しない系で計算を行いましたが、実際にはこれらは残念ながら存在し、エントリーとクローズのタイミングでコストとして資産を圧迫します。「固定のスプレッドがかかる系」、「リターンの分散に比例してスプレッドが変動する系」などでの計算ができたらより実用的になりそうです。

読んでくださりありがとうございました！
これからもたまに、勉強した中でためになりそうなものを記事にまとめるつもりです。よろしくお願いいたします。