上級編 ６． Libor Market Model
6.6 モンテカルロシミュレーション
6.6.7 リスク感応度（Greeks）の計測
6.6.7.4 Likelihood Ratio Method
6.6.7.4.6 LMM で LRM を使う方法
フォワード Libor の確率過程が LMM で記述されている場合、価格計算は MCS を使わざるを得ませんが、感応度も LRM によるサンプル平均で導出できるでしょうか？
古典的 LMM は、最大で 80～160 のフォワード Libor ベクトルが、3～5 個のブラウン運動に駆動され、それぞれ相関を持ちながら幾何ブラウン運動をすると仮定します。その時、Arbitrage Free の条件を満たす為にドリフト項で Convexity 調整が行われますが、そのドリフト項が、Liborベクトル（確率変数）に依存する事になります。この事から、LMM の場合、LRM が直ちには使えない理由が２つある事が分かります。
- 古典的な LMM から導出される将来のフォワード Libor の確率分布は、ドリフト項に確率変数が含まれる事から、微小時間の間に変動するので、確率密度関数が解析的に求まりません。
- また、LMMでは、最大で 80～160 個ある確率変数を、せいぜい3～5 個のブラウン運動が駆動する形で記述されているので、Libor 間の共分散行列がランクダウンされています（非正則行列）。その為逆行列を求める事ができません。6.199 式に示した、正規分布する確率変数ベクトルの同時確率密度関数は、式中に共分散行列 Σ の逆行列が含まれており、Σ の逆関数が求まらなければ、同時確率密度関数が定義できません。
という事で、一見すると非常に難しそうなのですが、Glasserman本やGlasserman-Zhaoの論文では、若干の工夫で、これらの問題点を解決し、LLM でも LRM が使えるような方法を提案しています。
まず、1.の問題ですが、LMM の式は、Black-Scholes モデルのように解析的に解く事はできませんが、初期値を与えれば、微小時間 ｈ の間における変化を、離散近似する事は可能です（というか、離散近似しないと LMM は実践で使えません。Section 6.6.5 離散化バイアス 参照）。PWD 法の説明の所で、すでに陽的オイラー法を使った離散近似式（6.186式）を示していましたが、それを再記します。
\[ \hat L_i^{(l)}(t+h) ≈ \hat L_i^{(l)}(t) e^{\left(\hat μ_i^{(l)}(t,{\bf L^{(l)}}(t))-\frac 1 2 ‖λ_i‖^2 \right)h+\sqrt{h} {\bf λ_i~z^{(l)}(t)}} \\ {\bf z ~ \sim N(0,1)}~~~ i=1,…,M,~~~l=1,…,n \tag{6.186} \]この式の両辺の対数を取れば、下記式のようになります。
\[ \ln L_i(t+h) ≈ \ln L_i(t)+ \left(μ_i(t,{\bf L(t)})- \frac 1 2 ‖{\bf λ_i(t)}‖^2 \right)h + \sqrt{h}~{\bf λ_i(t)~z, ~~~z ~\sim \mathscr{N}(0,I)} \tag{6.206} \]（6.186式にあった上付き文字のサンプルインデックス表記ははずしています。）
すると、ドリフト項係数も拡散項係数も、期間 h の間では一定とみなすので、\(\ln L_i(t+h)\) は、t→t+h 間で正規分布します。またその正規分布の平均は、\(\ln L_i(t)+(μ_i (t,{\bf L}(t))-\frac 1 2 ‖{\bf λ_i}(t)‖^2)h\) 、分散は \(‖{\bf λ_i}(t)‖^2 h\)、となります。ここから、\(\ln L_i(t+h)\) の確率密度関数が解析的に求まります。ただし、これは \(\ln L_i(t+h)\) 単独の限界確率密度関数であり、すべての Libor の同時確率密度関数の解析解を求めるには、さらに 2.の問題を解決する必要があります。
次に、2.の問題の解決方法ですが、これから説明する事を理解しやすくするため、表記を少し変えます。まず、時間軸を \(ｔ＝0\) から間隔 ｈ で離散化し、その各時点を次のように表記します。
\[ t_0≡0 ;~ t_1=t_0+h≡h,;~ t_2=t_1+h≡2h; …,t_m=t_{m-1}+h≡mh \]すなわち、時間軸を \(t_0~から~t_m\) まで、等間隔で m 個に分割します。ここで、m の決め方ですが、LMM におけるブラウン運動ベクトルの数 ｄ と、確率変数となるフォワードLiborの数 M を使って、ｄ × ｍ ＝Ｍ になるように ｍ を決めます。なぜそうするかは後で分かります。
こうやって分割された離散時間の間における \(\ln L_i(kh)\) の平均と分散は、先ほどの離散化近似式から、
\[ \begin{align} Average & = \ln L_i((k-1)h)+ \left(μ_i\left[(k-1)h,{\bf L}((k-1)h)\right]-\frac 1 2 ‖{\bf λ_i}((k-1)h)‖^2 \right)h \\ Variance & =‖{\bf λ_i}((k-1)h)‖^2 h \tag{6.207} \end{align} \]となります。この平均と分散を、表記を見やすくする為、次のように簡略表記します。
\[ \begin{align} \left(μ_i(t,{\bf L}(k-1))-\frac 1 2 ‖{\bf λ_i}((k-1)h)‖^2 \right)h & ≡a_i^{(kh)}h, \\ ‖{\bf λ_i}((k-1)h)‖^2 h & ≡‖{\bf λ_i}^{(kh)}‖^2 h \\ \end{align} \\ \ln L_i(kh) ~ \sim \mathscr{N} \left(\ln L_i((k-1)h)+a_i^{(kh)}h,~~~~‖λ_i^{(kh)}‖^2 h\right) \] \[\tag{6.208}\]\(a_i^{(kh)}~ も~(‖{\bf λ_i}^{(kh)}‖^2~ も~k-1~ →~k \) 間で一定の定数になります。
以上の表記を使って、\(L_i\) の拡散過程を離散時間で近似し、再帰的に \(ｔ＝0~から~mh\) まで順番に求めていくと、
\[ \begin{align} \ln L_i(h) & =\ln L_i(0)+a_i^{(h)}h+ \sqrt{h} {\bf λ_i}^{(h)}~{\bf ∙z_1,~~~~~z_1 ~\sim \mathscr{N}(0,1)} \\ \ln L_i(2h) & =\ln L_i(h)+a_i^{(2h)}h+ \sqrt{h} {\bf λ_i}^{(2h)}~{\bf ∙z_2,~~~~~z_2 ~\sim \mathscr{N}(0,1)} \\ \ln L_i(3h) & =\ln L_i(2h)+a_i^{(3h)}h+ \sqrt{h} {\bf λ_i}^{(3h)}~{\bf ∙z_3,~~~~~z_3 ~\sim \mathscr{N}(0,1)} \\ & ~~~~~~~~~~~~\vdots \\ \ln L_i(mh) & =\ln L_i(mh-1)+a_i^{(mh)}h+ \sqrt{h} {\bf λ_i}^{(mh)}~{\bf ∙z_m,~~~~~z_m ~\sim \mathscr{N}(0,1)} \tag{6.209} \end{align} \]但し、\(\bf z_1,z_2,…,z_m\) は、それぞれ独立な、d 次元の標準正規乱数ベクトルになります。そして、これらの式を繋げて一本化し、初期値 \(\ln L_i(0)\) が与えられた時の、 \(\ln L_i(mh)\) を表すと、
\[ \ln L_i(mh)= \ln L_i(0) + h\sum_{j=1}^m a_i^{(jh)} + \sqrt{h} \left[ {\bf λ_i^{(h)} |λ_i^{(2h)} |\dots|λ_i^{(mh)}}\right]∙ \left[ \matrix{ {\bf z_1} \\ {\bf z_2} \\ \vdots \\ {\bf z_m} } \right] \tag{6.210} \]但し、\(\left[ {\bf λ_i}^((h) )|{\bf λ_i}^{(2h)}|\dots|{\bf λ_i}^{(mh)}\right]\) は、d 次元のベクトル \({\bf λ_i}^{(jh)}\) を、横につなげた M(=d×m) 次の行ベクトルで、\(\left[\matrix {{\bf z_1}\\ \vdots\\{\bf z_m}} \right]\) は、d 次元のベクトル \(\bf z_j\) を縦につなげたM次の列ベクトルになります。
さらに、すべてのフォワード Libor に対応する行ベクトル \(\sqrt{h}\left[ {\bf λ_i}^{(h)}|{\bf λ_i}^{(2h)} | \dots | {\bf λ_i}^{(mh)}\right],~~i=1,…,M\) を、縦に積み重ね、Ｍ×Ｍ 次元の行列を生成し、それを C と表記します。すると、\(\bf C∙C^{\top}\)は M×M 次の共分散行列になり、それを下記のように Σ と表記します。
\[ \sqrt{h} \left[\matrix {{\bf λ_1}^{(h)}|{\bf λ_1}^{(2h)} | \dots | {\bf λ_1}^{(mh)} \\ {\bf λ_2}^{(h)}|{\bf λ_2}^{(2h)} | \dots | {\bf λ_2}^{(mh)} \\ ~~~~~~~~ \vdots \\ {\bf λ_M}^{(h)}|{\bf λ_M}^{(2h)} | \dots | {\bf λ_M}^{(mh)} } \right]≡{\bf C}, ~~~~~~ \bf C∙C^{\top}≡Σ \]先ほど、m の決め方を説明しましたが、この Ｍ×Ｍ 次の行列を作りたかったからです。そして、やりたい事は、こうやって無理やり生成した正方行列 Σ を、正則行列にしたいという事です。それは、あまり心配する必要はなく、既に説明した LMM における相関行列の生成方法に従えば、Σ は、ほぼ問題なく正則行列になります。
さらに、分布の平均値の式に含まれる \(a_i^{(jh)},~~i=1,…,M~~ j=1,…,m\) の値ですが、LMM のドリフト項係数から導出されるので、そのままでは確率変数となるLiborベクトルに依存し、計算が非常にややこしくなります。そこで、PWD の所で使ったのと同じ簡略化テクニックを使います。すなわち、Liborを ｔ＝0 時の値で固定し定数化、それを \(\bar a_i^{(jh)}\) と表記します。これにより計算負荷を大幅に軽減できます。
以上の表記を使い、すべてのフォワードLibor \(L_1,~L_2,~…,~L_M\) の離散近似式をひとつにまとめてベクトル表記すると
\[ \begin{align} & \ln{\bf L}(mh)=\ln {\bf L(0)}+ h \sum_{j=1}^m {\bf \bar a^{(jh)}} + {\bf C∙z}, ~~~~~\bf C∙C^{\top}=Σ \\ & {\bf z ~\sim \mathscr{N}(0,1)},~~~~~~~~ \ln {\bf L(mh) }~ \sim \mathscr{N}(\ln {\bf L(0)} + h \sum_{j=1}^m {\bf \bar a^{(jh)},~~~Σ}) \tag{6.211} \end{align} \]となり、ベクトル \(\ln{\bf L}(mh)\) は、平均 \( \ln {\bf L}(0)+ h \sum_{j=1}^m {\bf \bar a^{(jh)}} \) 、共分散\(\bf Σ\) で同時正規分布し、なおかつ Σ は正則行列で、その逆行列を求める事ができます。すると、ベクトル \(\ln {\bf L}(mh)\) の同時確率密度関数は、それらを 6.199 式に代入すれば求まります。
\[ \small \begin{align} & g(\ln {\bf L}(mh)) \\ & = \frac {1}{\sqrt{2π |Σ|}} \exp\left[- \frac 1 2 \left(\ln {\bf L}(mh)-\ln {\bf L}(0) -h\sum_{j=1}^m {\bf \bar a^{(jh)}} \right)^{\top} {\bf Σ^{-1}} \left(\ln {\bf L}(mh)-\ln {\bf L}(0) - h\sum_{j=1}^m {\bf \bar a^{(jh)}} \right)\right] \end{align} \] \[\tag{6.212}\]これで、1. 2.の問題はクリアされました。
この同時確率密度関数には、パラメータとして、Libor の初期値 \(L_i(0),~i=1,…,M\) が含まれています（\(\bar a{(jh)}\) のパラメータになっています）。それらに対する感応度を求めたければ、まず 6.201 式を使って \( \frac {∂g\left(\ln L(mh),\bar a^{(jh)}(L_k(0))\right)/∂L_k(0)}{g\left(\ln L(mh),\bar a^{(jh)}(L_k(0))\right)}\) を求め、それを 6.202 式に代入すればいいだけです。しかしその前に、6.201 式には、平均ベクトルのパラメータでの微分 \(\frac {∂a(θ)}{∂θ_k}\) が含まれており、それを求めなければいけません。以下にそれを示します。
まず 6.211 式の \(\bf a(θ)\) における \(\bf \bar a^{(jh)}\) をもとの表記に戻します。
\[ {\bf a(θ)} ≈ {\bf \ln L(0)} + h \sum_{j=1}^m {\bf \bar a^{(jh)} } =\ln {\bf L(0)}+ h\sum_{j=1}^m (\bf \bar μ(jh,L(0))-1/2 ‖λ(jh)‖^2 ) \]この式の両辺を \(L_k(0)\) で偏微分します。分かり易くするために、\(\bf a(θ)\) を、ベクトルの要素に分解し、それを \(L_k(0)\) で偏微分した式を示します。
\[ \begin{align} & \frac {∂a_i(…,L_k (0))}{∂L_k (0)} =\frac {∂}{∂L_k(0)} \left[\ln L_i(0)+h\sum_{j=1}^m \left(\bar μ_i(jh,{\bf L}(0))- \frac 1 2 ‖{\bf λ_i}(jh)‖^2 \right) \right] \\ & =\frac {∂}{∂L_k(0)} \left[\ln L_i(0)+h\sum_{j=1}^m \left(\sum_{q=η(jh)}^i \frac {δL_q(0) {\bf λ_i}(jh)∙ {\bf λ_q}(jh)}{1+δL_q(jh)}- \frac 1 2 ‖{\bf λ_i}(jh)‖^2 \right) \right] \\ & =\frac {1}{L_k(0)}I_{\{i=k\}} +h \sum_{j=1}^m \sum_{q=η(jh)}^i \frac {d}{d L_k(0)} \left(\frac {δL_q (0){\bf λ_i}(jh)∙{\bf λ_q}(jh)}{1+δL_q(jh)}\right) \\ & =\frac {1}{L_k(0)}I_{\{i=k\}} +h \sum_{j=1}^m \frac {δ{\bf λ_i}(jh)\cdot {\bf λ_k}(jh)}{(1+δL_k(0))^2}I_{\{η(jh)≤k≤i\}},~~~~~ i=1,…,M \end{align} \] \[\tag{6.213}\]上式の最後の行の右辺第 2 項は、PWD の所で、すでに導出済み（6.188式）なので、それを使いました。この式は、すべて、Deterministic な関数である \({\bf λ_i}(jh)~ と~ {\bf λ_k}(jh)\) にしか依存しないので、すべてのサンプル経路で共通です。従って、あらかじめ計算しておけば、サンプルごとの計算は省けるので、計算負荷が大幅に軽減できます。 最終的に、これを 6.201 式に代入すれば、6.202 式から \(L_k(0)\) に対する感応度が求まります。
6.6.7.4.7 ガンマの計測
すでに述べた通り、ＭＣＳにおけるガンマの計測は、ＰＷＤ法では対応不可能で、有限差分近似では、推定誤差がかなり大きくなるので、精確性で劣ります。一方、ＬＲＭでは、下記式の右辺の積分にある\(\frac {g''(x,θ)}{g(x,θ)}\) が求まれば、最後の期待値演算を、サンプル平均で近似する事によって、比較的簡単に導出できます。但し、推定誤差は、有限差分近似よりは小さいものの、無い訳ではありません。
\[ \begin{align} \frac{∂^2}{∂~θ^2} E\left[\frac {Payoff(T,x)}{Numeraire(T,x)}\right] & =\int \frac {Payoff(T,x)}{Numeraire(T,x)}~fb \frac {∂~^2g(x,θ)}{∂~θ^2}d(x) \\ & =\int \frac {Payoff(T,x)}{Numeraire(T,x)}~ \frac {g''(x,θ)}{g(x,θ)} g(x,θ)~d(x) \\ & =E \left[\frac {Payoff(T,x)}{Numeraire(T,x)}~ \frac {g''(x,θ)}{g(x,θ)}\right] \end{align} \]上式の、期待値演算の微分と、微分の期待値演算（積分に相当）の交換が成立する為には、確率密度関数のパラメータでの2階微分が可能である必要がありますが、よく使われる正規分布の確率密度関数は、極めて滑らかで、それが問題になる事はありません。
確率密度関数が、正規分布のそれの場合、すでに \(\frac{g'(x,θ)}{g(x,θ)}\) の導出方法は、説明しましたが(6.200、6.201式)、\(\frac{g''(x,θ)}{g(x,θ)}\) の導出も、同様の手順で解析できます。下記にその結果だけを示します。
\[ \begin{align} & \frac {g''({\bf x,a(θ),Σ})}{g({\bf x,a(θ),Σ})} \\ & =\left[({\bf x-a(θ)})^{\top} {\bf Σ^{-1}} \frac {∂{\bf a(θ)}}{∂θ_k}\right]^2 -\left(\frac {{\bf ∂a(θ)}}{∂θ_k}\right)^{\top} {\bf Σ^{-1}} \frac {∂{\bf a(θ)}}{∂θ_k}+ ({\bf x-a(θ)})^{\top}{\bf Σ^{-1}} \frac {∂^2~{\bf a(θ)}}{∂θ_k^2} \end{align} \]サンプルで生成された上式にPayoff 関数を掛けてたものが、ガンマの不偏推定量になり、そのサンプル平均を取れば、ガンマのシミュレーション値が求まりす。
１点だけ注意点。LRM による感応度の導出は、所詮シミュレーション値なので、サンプルのばらつきによる推定誤差が発生します。その推定誤差は、1 階微分であるデルタやベガに比べ、2 階微分となるガンマは、より大きくなります。上の式の第 1 項は、確率変数 ｘ の2乗に比例するので、推定誤差のオーダーも、それに比例します。
他の感応度計算方法と、ＬＲＭ を組み合わせて、ガンマの計算も可能です。例えば、デルタを ＬＲＭ で２回計算した後、有限差分近似で、ガンマを計算する事ができます。あるいは、ＬＲＭ と ＰＷＤを組み合わせ、どちらか一方で、デルタを計算した後、他方でガンマを計算する事ができます。この方が、ＬＲＭだけでガンマを計算するより推定誤差を小さくできる事が経験的に知られています。
ガンマや、ボルガ、バンナといった、2 階微分に相当する感応度は、日々のリスク管理で絶えず監視しヘッジするリスク量ではなく、市場が大きく動いた時に、デルタやベガでは説明しきれない損益の変動を把握するためのリスク量になります。こういったリスク量は、シナリオ分析やストレステストを、時折行って把握するのが、より実践的かと思います。特に、Deep in the money や Deep out of the money のオプションのポジションや、Barrier Optionのように、バリア近辺でガンマが非常に強いポジションが多くある場合は、そういったリスク管理手法が有効です。