上級編 ７． Local Volatility Model とStochastic Volatility Model
7.3 Stochastic Volatility Model
7.3.2 Displaced Log-Normal Heston Model ("DLN-H")
7.3.2.4 DLN-H モデルからヨーロピアンオプションの価格式の導出
さて、7.40、7.41式の SDE は DLN モデルと Heston Model の複合型であり、確率分布が複雑な形状になるため、分布関数や密度関数を解析的に求めるのは困難です。そこで Andersen-Piterbarg本では、特性関数をフーリエ変換するテクニックを使って、高速で精確な近似解を求める方法を紹介しています。フーリエ変換の方法は、もともと最初に Heston が、Heston Modelからオプション価格式を導出する方法として提示したものですが、それを若干修正したものです。ただ、相当難解な数学のテクニックを使っており、理解するのは簡単ではありません。以下にそのテクニックを紹介しますが、理解を容易にする為に、まず解析のプロセスの概略を示します。
- まず S(t)を変数変換し、よりシンプルな SDE の形にします。変換後の確率変数をX(t)とし、その対数 ln X(t) の特性関数を求めます。
- その特性関数ですが、直接は求まりません。ただ、解析していくと、Stochastic Varianceである z(t)の特性関数から求まります（厳密には、z(t)をｔ時まで時間軸で積分した \(\bar z(t)\) の、さらにそれを測度変換した \(\bar z^*(t)\) の特性関数）。\(\bar z^*(t)\) の特性関数は、ファインマン・カッツの定理を使って、２階の微分方程式から導出できます。\(\bar z^*(t)\) の特性関数が解析的に求まれば、そこから Ln X(t)の特性関数も導出できます。
- 特性関数は、確率密度関数と、フーリエ変換により双対関係にあります。それをうまく利用すれば、ヨーロピアンオプションの Payoffの期待値が導出できます。その解析プロセスでは、フーリエ変換をうまく行えるようにする為に、変数変換や、減衰ファクターや、畳み込み積分、などの数学テックニックを使います。
- 最終的に、上記のプロセスによる解（期待値）は、フーリエ積分の形で求まります。その積分形は、それ以上解析不可能で、数値積分で解く事になります。それを、高速で精確に行うテクニックがいくつか紹介されていますが、Andersen-Piterbarg本では、Black モデルでの解析解を制御変数として使う事で、かなり高速で精度の高い数値解の導出プロセスが紹介されています。
では、上記のステップに従って、解析プロセスの解説をします。かなり長く、難解なプロセスです。実務家は、結果だけをありがたく使わせて頂けばいいので、結論部分に飛んで頂いて結構です。この手のテクニックを理解されたい方は、読み進んで下さい。
1. 変数変換
7.40式の SDE はやや複雑なので、変数変換により、より解析しやすい形に変換します。まず、変換後の確率変数を X(t)とし、下記式を使って、S(t)を X(t)に変換します。 \[ \begin{align} & S(t)=\frac 1 b \left[(bS(0)+(1-b)L)X(t)-(1-b)L \right] \\ & ~~or~~~~X(t)=\frac {bS(t)+(1-b)L}{bS(0)+(1-b)L} ~~~~ \tag{7.46} \end{align} \]
これに伊藤の公式を使って 7.40式を変換すると、下記式のように、X(t)を確率変数とする、よりシンプルな SDE が導出できます。 \[ \begin{align} & \frac {dX(t)}{X(t) }=λ~b \sqrt {z(t)}~dW(t) \\ & ~~~ or~~~~~ \ln X(t)=λ~b∫_0^t \sqrt {z(s)}~dW(s) - \frac 1 2 λ^2 b^2 ∫_0^t \sqrt{z(s)}~ds, \tag{7.47} \end{align} \] （1行目の式から、X(t)をさらに\(\ln X(t)\)に変数変換し、伊藤の公式を使って導出したSDEを、さらに積分形式で表したのが 2行目の式になります）
ここで、式中にある \(∫_0^t \sqrt {z(s)}ds\) を、\(∫_0^t \sqrt {z(s)}ds :=\bar z(t)\) と置き、上式の右辺を下記式のように書き換えます。 \[ \ln X(t)=λ~b∫_0^t \sqrt {z(s)}dW(s)-\frac 1 2 λ^2 b^2 \bar z(t), \tag{7.48} \]
2. \(\ln X(t)\) の特性関数を求める。
この式を使って、\(\ln X(t) \) の特性関数を求めます。\(\ln X(t) \) の特性関数を \(Ψ_X (u;t),~~u∈\mathcal C \) と表記すると、特性関数は \[ Ψ_X (u;t):= E \left[e^{u \lnX(t)}\right] \] と定義されます。この式の \(\ln X(t) \) に、7.48式の右辺を代入して、下記のように変形していきます。
\[ \begin{align} E \left(e^{u \ln X(t) }\right)& =E\left[e^{u \left(λb∫_0^t \sqrt {z(s)} dW(s)- \frac 1 2 λ^2 b^2 \bar z(s)\right)} \right] \\ & =E\left[e^{u~λ~b∫_0^t \sqrt{z(s)}dW(s)} e^{- \frac 1 2 u~λ^2 b^2 \bar z(t)}\right] \\ & =E\left[e^{u~λ~b∫_0^t \sqrt {z(s)}dW(s)} \left[ e^{\left(-\frac {1}{2} u^2 λ^2 b^2 \bar z(t) + \frac 1 2 u^2 λ^2 b^2 \bar z(t) \right) } \right] e^{\left(-\frac {1}{2} u~λ^2 b^2 \bar z(t)\right)} \right] \\ & =E\left[ e^{\left(u~λ~b∫_0^t \sqrt{z(s)}dW(s)-\frac 1 2 u^2 λ^2 b^2 \bar z(t)\right)} ×e^{\left(\frac 1 2 u^2 λ^2 b^2 \bar z(t)- \frac 1 2 u~λ^2 b^2 \bar z(t)\right)} \right] \\ & =E\left[ς(t)∙e^{\left(\frac 1 2 u(u-1) λ^2 b^2 \bar z(t)\right)}\right] \tag{7.49} \end{align} \] \[ 但し、~~~~~ς(t)=e^{\left(u~λ~b∫_0^t\sqrt{z(s)}dW(s)- \frac 1 2 (uλb)^2 \bar z(t)\right)} \]なぜ、このように変形したか、よく見てみます。7.49式の右辺 3行目で、期待値の演算内において、\(1=\left[ e^{(-\frac {1}{2} u^2 λ^2 b^2 \bar z(t) + \frac 1 2 u^2 λ^2 b^2 \bar z(t)) } \right]\) を掛けた上で、4行目にある通り、2つの指数関数の積に分解します。その片方を ς(t) と置いて、5行目の形に変形しています。最後の ς(t) の式をよく見てください。ς(t)が指数マルチンゲールとなる確率変数の形になっている事が分かります。指数マルチンゲールの確率変数は、実はラドンニコディム微分と看做せます。さらに、7.49式の右辺の最後の行をよく見てください。\(\bar z(t)\) を確率変数と看做し、\(\frac 1 2 u(u-1) λ^2 b^2\) を特性関数のモーメントと看做せば、式全体として、ς(t) で測度変換した、\(\bar z(t)\) の特性関数の形になっています。すなわち、ln X の特性関数 \(Ψ_X (u;t)\) が、Stochastic Varianceを時間軸で積分した値 \(\bar z(t)\) の、(さらにそれをラドンニコディム微分で測度変換した)特性関数から導出できるという事です。\(\bar z(t)\) の特性関数を \(Ψ_{\bar z}(v(u);t):= E^{\bar P}\left[e^{v(u) \bar z(t)} \right],\) と表記すると、\(Ψ_X(u;t)\) は下記のように \(Ψ_{\bar z}(v(u);t)\) から導出できます。 \[ Ψ_X (u;t)=Ψ_{\bar z}\left( \frac 1 2 λ^2 b^2 u(u-1);t\right) \tag{7.50} \]
3. \(\bar z\)の特性関数を求める
では、\(\bar z\) を測度変換した特性関数をどうやって求めるかですが、下記のプロセスで、解析的に導出できます。
3.1
ラドンニコディム微分に相当する ζ(t) を使って、\( \bar z(t)\) から測度変換された確率変数を \(\bar z^*(t)\) とおくと、その SDE は、ギルザノフの定理を使って、次のように導出できます。 \[ d\bar z^*(t)=\left(θ(z_0-\bar z^*(t))+ρ η λ b u \bar z^*(t)\right)dt+η \sqrt {\bar z^*(t)} d\bar W(t) \tag{7.51} \] 但し \(\bar W(t)\) は、測度変換後の新たな測度 \(\bar P\) 下のブラウン運動。
3.2
この 7.51式は、やはり、中心回帰する平方根過程（CIR過程と同じ）になっています。ここから、\(\bar z^*(t)\) の特性関数は、ファインマン・カッツの定理を使って、2階の偏微分方程式を解く形で求めます。ファインマン・カッツの定理では、ガウス過程の確率変数の期待値は、その SDE のドリフト項係数と拡散項係数を使った 2階の偏微分方程式の解として求まるというものでした（“寄り道＜期待値の演算方法＞” 参照）。特性関数も期待値演算なので、SDEが特定できていれば、この定理が使えます。
3.3
解析のプロセスは飛ばしますが、その偏微分方程式から、2つの未知関数 \(A(v,u),~B(v,u)\) を持つ、連立常微分方程式が導出でき、\(Ψ_{\bar z}(v(u);t)\) は、その解を使って下記式のように表現できます。 \[ \begin{align} & Ψ_{\bar z}(v(u);t)=A(v,u)+B(v,u)z_0 \\ & ~~ 但し~~A(v,u)=\frac {θz_0}{η^2} \left[2 \ln {\left(\frac {2γ}{θ'+γ-e^{-γT}(θ'-γ) }\right)}+(θ'-γ)T\right] \\ & ~~~~~~~~B(v,u)=\frac {2v(1-e^{-γT})}{(θ'+γ)(1-e^{-γT})+2γe^{-γT} }, \\ & ~~~~~~~~γ=γ(v,u)=\sqrt{(θ')^2-2η^2v} \\ & ~~~~~~~~θ'=θ'(u)=θ-ρ~ η~ λ~ b~ u \tag{7.52} \end{align} \]
4. フーリエ変換を使ったヨーロピアンオプション価格の導出
さてここから、この特性関数を使って、ヨーロピアンCallオプションの価格を導出します。ここでは、特性関数の逆フーリエ変換が確率密度関数になるという関係を使い、Payoffの期待値を計算します。ところが Payoff関数は、そのままでは∞に発散するので、それをフーリエ変換する為には、有界かつ可積分な関数に変換する必要があります。その為に、やや技巧的なテクニックを使っているので、それを念頭に読み進めてください。
4.1
まず、解析をすすめ易くする為、オプション行使時の対象資産価格 S(T)と、ストライク価格 K の両方の対数を取り、それぞれ ξ, k と置きます。そして、それを使って、Payoff関数を書き換えます。 \[ \begin{align} & \lnS(T):=ξ, ~~~~~~~~~~~~\lnK:=k~~~~~~~~~~~~~~~ \\ & Payoff=max(e^ξ-e^k,0) \end{align} \]
4.2
この Payoff 関数の期待値が Option価格になります。 \[ EuropeanCall=E\left[\max(e^ξ-e^k,0)\right]=E\left[e^ξ-e^k\left( \min(e^{ξ-k},1)\right)\right] \tag{7.53} \] ここで、上式右辺にある通り Payoff関数を少し変形していますが、期待値演算の中身を有界にする為です。
4.3
ここで ξ の特性関数を \(Ψ_ξ (u)\) と置くと、上式は、下記のようになります。特に下記式の右辺第１項は、特性関数で u=1 とすれば、それは確率変数の期待値になる関係を使っています。 \[ E\left[e^ξ-e^k\left(\min(e^{ξ-k},1)\right)\right]=Ψ_ξ(1)-e^k E\left[\min(e^{ξ-k},1)\right] \tag{7.54} \]
4.4
ここで、ξ の確率密度関数を p(ξ) と置き、右辺の期待値演算を、p(ξ)を使った積分形式で表現すると下記のようになります。p(ξ)は、解析的に求めるのが困難ですが、後で、フーリエ逆変換して特性関数に変換します。 \[ E \left[\min(e^{ξ-k},1)\right]=∫_{-∞}^∞ \min(e^{x-k},1) p(x)dx \] \(\min(e^{x-k},1)~は、(0,1]\) の範囲で値を取るので有界ですが、\(∞\gt x \gt k\) の範囲でずっと 1になり、上式の積分範囲では発散し、可積分になりません。そこで、Payoff関数に、下記式のように減衰ファクターをかけて、フーリエ変換可能な可積分な関数形に変換します。すると、うまいことに被積分関数は、点ｋを基準に評価した２つの関数のConvolution（畳み込み積分）の形になります。 \[ \begin{align} ∫_{-∞}^∞ \min (e^{x-k},1)p(x)dx & =∫_{-∞}^∞ \min (e^{x-k},1) \left[e^{α(k-x)}~e^αx~ e^{-αk}\right]p(x)dx \\ & =e^{-αk} ∫_{-∞}^∞ \left[\min (e^{-(k-x)},1) e^{α(k-x)}\right] \left[e^{αx} p(x)\right] dx \tag{7.55} \end{align} \]
（注：式の右辺は、被積分関数に \( 1=\left[e^{α(k-x)}~e^{αx}~e^{-αk}\right]\) を掛けて、それを \(\min(e^{x-k},1)~側と、p(x)\) 側と、積分の外側、に分解して掛けています。すると、\(\left[\min(e^{-(k-x)},1)~e^{α(k-x)}\right]~も、[e^{αx}p(x)]\) も、有界でかつ積分範囲の両端で 0 に収束し、周期関数（すなわちフーリエ変換可能）になり、なおかつこの２つの関数の積分はConvolutionの関係になります。）
α は任意の正の実数で、後程説明する数値積分が最も効率的になるような値に設定（実践では0.5が使われる）。
4.5
7.55式にある積分は、うまくｋ点を基準としたConvolutionになりました。ところでConvolution のフーリエ変換は、元の 2つの関数を先にフーリエ変換して掛けたものと同値です。式で表すと、下記式のようになります。 \[ \begin{align} & ∫_{-∞}^∞ \left[ \min(e^{-(k-x)},1)~e^α{k-x}\right] \left[e^αk p(x)\right]dx =(f_1*f_2 )(k) \\ & =\left(F^{-1}(F(f_1*f_2 )(ω))\right) (k) \\ & = \left(F^{-1} (g_1 (ω)*g_2 (ω))\right)(k) \tag{7.56} \end{align} \]
（注 ：\(f_1*f_2\) は、Convolution(k を基準点とした x での積分)なので、k の関数。それをフーリエ変換すると\(F(f_1*f_2)(ω)\) となり、ωの関数になる。さらにそれをフーリエ逆変換すると、\(\left(F^{-1}(g_1(ω)*g_2 (ω))\right)(k)\) となって、元のConvolutionである ｋ の関数になる。また \(g_1(ω),~g_2 (ω)\) は、下記の通り、\(f_1,~f_2\) をそれぞれフーリエ変換したもの）
4.6
上記のConvolutionの関係にある 2つの関数\(f_1,~ f_2\)を、それぞれフーリエ変換します。フーリエ変換された関数を \(g_1,~g_2\)と置くと、下記のように解析できます。 \[ \begin{align} g_1 (ω) & =∫_{-∞}^∞ e^{iωx} \min(e^{-x},1) e^{αx} dx ~~~~~~~~~~\\ & =∫_{-∞}^0 e^{iωx}e^{αx}dx+∫_0^∞ e^{iωx}e^{-x}e^{αx}dx \\ & =∫_{-∞}^0 e^{x(α+iω)}dx+∫_0^∞ e^{x(-1+α+iω)}dx \\ & =\frac {1}{α+iω} -\frac {1}{α-1+iω} \\ & =\frac {1}{(α+iω)(1-(α+iω))} \tag{7.57} \end{align} \] \[ \begin{align} g_2 (ω)& =∫_{-∞}^∞ e^{iωx} e^{αx}p(x)dx ~~~~~~~~~~~~~~~~~~~~~~~~~~~\\ & =∫_{-∞}^∞ e^{α+iωx} p(x)dx \\ & =E \left(e^{(α+iω)ξ} \right)=Ψ_ξ (α+iω) \tag{7.58} \end{align} \]
（注：\(g_2\) は、確率密度関数のフーリエ変換が特性関数になる事を利用）
ここから、この 2つの関数の積が、下記にように求まります。 \[ g_1(ω) g_2 (ω)=\frac {Ψ_ξ (α+iω)}{(α+iω)(1-(α+iω))} \tag{7.59} \]
4.7
さらに 7.59式のフーリエ逆変換を行えば、元の期待値演算に戻ります。 \[ \begin{align} \left(F^{-1} (g_1(ω)*g_2 (ω))\right)(k) & =\frac {1}{2π} ∫_{-∞}^∞ e^{-k(α+iω)} \frac {Ψ_ξ(α+iω)}{(α+iω)(1-(α+iω))}dω\\ & =E \left[\min(e^{ξ-k} ,1) \right] \tag{7.60} \end{align} \]
4.8
上記の 7.60式を、元のPayoffの期待値の 7.54式に代入します。 \[ EuropeanCall=Ψ_ξ(1)-\frac {e^k}{2π}∫_{-∞}^∞ e^{-k(α+iω)} \frac {Ψ_ξ(α+iω)}{(α+iω)(1-(α+iω))}dω \tag{7.61} \]
4.9
さらに、右辺の第2項にある複素積分は、実数部分の積分の 2倍になる事から、下記式のようにできます。（証明のプロセスは省略します。） \[ EuropeanCall=Ψ_ξ(1)- \frac {e^k}{π}∫_0^∞ Re\left[\frac {e^{-k(α+iω)} Ψ_ξ (α+iω)}{(α+iω)(1-(α+iω))} \right]dω \tag{7.62} \]
4.10
最後に右辺にある \(Ψ_ξ(α+iω)\) に、7.52式で求めた特性関数を代入すれば、完成です。
5. 積分経路に関する注意点
かなり難解な数学のテクニックを使って、何とか Option 価格式をフーリエ積分の形まで持ってきました。実際の計算では、7.62式の第2項にあるフーリエ積分を、数値積分で求める必要があります。それを、できるだけ高速で精確に行うテクニックがいくつか紹介されていますが、次のセクションで、Andersen-Piterbarg本で示されている、かなり精確なアルゴリズムを紹介します。また、7.62式の被積分関数にある特性関数は、7.52式にある通り、複素数の対数関数、すなわち多価関数が含まれており、これを積分する際は、積分経路が分岐切断を通るかどうかに注意が必要です。幸い、Andersen-Piterbargが導出した、上記方法の場合、分岐切断を通らないので大丈夫です。（Heston Modelで、Hestonが当初紹介したフーリエ積分では、分岐切断を通るので、その対応が必要です。）