パリジャン・オプションの評価と陽的有限差分法

【論 説】

パリジャン・オプションの評価と 陽的有限差分法

久 保 徳 次 郎

1 は じ め に

パリジャン・オプションは，バリア・オプションの一種で，強い経路依存性 を有している．これは，通常のバリア・オプションと違って，原資産価格が一 定時間バリア外（内）にあればノックアウト（ノックイン）されるという特徴を もっている．したがって，このパリジャン・オプションと区別するために，通 常のバリア・オプションのことをワンタッチ・バリア・オプションと呼ぶ場合 がある．

パリジャン・オプションの評価に関しては，Haber-Scho¨nbucher-Wilmott（以

下ではHSW）（1999）では有限差分法，Vetzal-Forsyth（1999）では有限要素法，

そしてKwok-Lau（2001）ではForward Shooting Grid 法を使ってその数値計算

がそれぞれ行われている．とりわけ，HSW のアルゴリズムは単純で取扱いが容 易であると言える．このHSW の数値計算法に関しては，Wilmott（2000）で陽 的有限差分法（陽解法）によるプログラム例が紹介されている．ただし，これに は次のような改善・追加されるべき余地が残されている．まず第1 に，Black-

Sholes（以下ではBS）方程式を直接的に陽解法で近似しているため，計算効率

はあまり高くない．第2 に，バリア値が正確にグリッド点上に乗っていないた め，それが計算精度低下の一因となっている．そして第3 に，プログラム例が ダウン・アンド・アウト型のヨーロピアン・プット・オプションのケースに限

っている．

本稿の目的は，HSW あるいはWilmott の数値計算法に改良を加え，様々な パリジャン・オプションを比較検討することにある．まず上記の第1 の問題点 に関しては，HSW でも指摘されているように，原資産価格を対数変換すること によって改善することができる．しかし，本稿では，BS 方程式を熱方程式

（heat equation）に変換することによってさらに計算効率を改善できるということ を示す．また，この熱方程式への変換による方法では，実行価格が陽表的に出 てこなくなるため，もっぱらバリア値をグリッド点上に乗せることにプログラ ムを集中させることができる．したがって，上記の第2 の問題点も同時に改善 することが可能となる．これに加えて，線形補間法と外挿法を用いることによ って，さらに計算精度と効率性を向上させうるということも示す．そして第3 の問題点に関しては，ヨーロピアン型，アメリカン型，ダウン・アンド・アウ ト型，アップ・アンド・アウト型などを含めたすべてのケースのプログラムを 新たに作成し，それら数値計算結果を比較検討することにする．

また，3 つの陽解法の計算効率と安定条件についても再検討を加えることに する． ただし， 陽解法の安定条件に関しては， 従来の教科書と違って，

Gerschgorin の定理を使って導出を行うことにする．周知のように，陽解法は取 扱いが容易で，ツリー・モデルとも共通点が多い．この数値解法に対数変換や 熱方程式への変換を施すことによって，陰解法やCrank-Nicholson 法よりも精 度の高い計算を実現することが可能となる．

本稿の構成は以下のとおりである．第2 節では，陽的有限差分法による3 つ の数値計算法について，その計算効率と安定条件の側面から比較検討を行うこ とにする．第3 節では，熱方程式への変換を施した場合の陽解法を使って，パ リジャン・オプションの評価について考察を行い，第4 節では，様々なパリジ ャン・オプションの数値計算結果の比較検討を行うことにする．そして第5 節 では，本稿の要約と今後の研究課題について簡単に触れることにする．

2 陽的有限差分法

2. 1 直接的な差分近似の場合

任意の時間をt，原資産価格をS，ドリフトをη，ボラティリティをσ，ウィ ーナー過程をdw とそれぞれ表し，原資産価格が次のような確率過程に従うも のと仮定する．

dS＝ηSdt＋σSdw

この原資産に対する派生証券の価格がF(S, t) という関数形で表せるならば，次 のような周知のBS 方程式をえることができる．

（1）

ただし，r は安全利子率，dは原資産が株式ならば連続配当率，外国通貨なら ば外国利子率である．

Sの最大値と最小値をSmax，Smin とそれぞれ表し，またS のステップ幅をΔS と表すと，グリッド点上のS は

S(i)＝iΔS＋Smin，（ただし，i＝0, 1, ……, n）

と表すことができる．この式より，グリッド点上のS の最大値と最小値は，

S(n)＝Smax＝nΔS＋Smin

S(0)＝Smin

とそれぞれ表すことができる．また，満期時点までの長さをT，t のステップ幅

をΔt とそれぞれ表すと，グリッド点上のt は，

t(j)＝jΔt，（ただし，j＝0, 1, ……, m）

と表すことができる．ただし，

t(m)＝T＝mΔt

である．以下では，F(S (i), t(j)) を単にF(i, j)と表すことにする．

BS 方程式（1）を直接的に陽解法で差分近似する場合，次のような方法がと られる．まず，S に関する1 次導関数には中心差分を，t に関する1 次導関数

∂F

∂t

∂F

∂S

∂²F

∂²S 1

＋ σ2 ²S²＋(r−d)S   −rF＝0

には後退差分をそれぞれ用いる．すなわち，グリッド点上では，

となる．そして，Sに関する2 次導関数には次のような差分近似を用いる．

かくして，以上の差分近似式を（1）に適用すると，次のような結果をえること ができる．

F(i, j)＝a(i) F(i＋1, j＋1)＋b(i) F(i, j＋1)＋c(i) F(i−1, j＋1) （2）

ただし，

a(i)＝{σ²i²＋(r−d) i}Δt/ 2 b(i)＝1−(σ²i²−r)Δt c(i)＝{σ²i²−(r−d) i}Δt/ 2

である．ここで注意すべきことは，a(i)，b(i)，c(i) がi の大きさ，すなわち原 資産価格の大きさに依存しているという点である．これは，数値計算の効率性 を悪化させる一因となる．そこで以下で，原資産価格の対数変換，BS 方程式 の熱方程式への変換を通じて，この点の改善が図れるということを示すことに しよう．

2. 2 対数変換を施して差分近似した場合

いま，x＝lnSとおき，G(x, t)＝F(S, t) となるような変数変換を考えた場合，

BS 方程式（1）は，次のように書き換えることができる．

（3）

さらに，G(i, j)＝G(x(i), t(j)) とおき，さきの場合と同様に，上式を陽解法で差

∂F

∂S

F(i＋1, j)−F(i−1, j)

〜  2ΔS

〜 

∂F

∂S

F(i, j)−F(i, j−1) Δt

〜 〜 

∂²F

∂S²

F(i＋1, j)−2F(i, j)＋F(i−1, j) ΔS²

〜 〜 

∂G

∂t

∂²G

∂x²

∂G

∂x 1

＋ σ2 ²   ＋(r−d−σ²/ 2)   −rG＝0

分近似すると，次のような結果をえることができる．

G(i, j)＝aG(i＋1, j＋1)＋bG(i, j＋1)＋cG(i−1, j＋1) （4）

ただし，

a＝{σ²/ (Δx)²＋(r−d−σ²/ 2) /Δx}Δt / 2 b＝1−{r−σ²/ (Δx)²}Δt

c＝{σ²/ (Δx)²−(r−d−σ²/ 2) /Δx}Δt / 2

である．上式より，対数変換を施した場合の陽解法では，a，b，c は定数とな り，x（原資産価格の自然対数値）の大きさに依存することなく，数値計算の効率 性を高めうるということが分かる．

2. 3 熱方程式に変換して差分近似した場合

次に，BS 方程式（1）を熱方程式に変換した場合について考えてみることに しよう．そこで，

y＝ln (S / K) z＝(T−t)σ²/ 2

とおき，次のように置き換えることにする．

H(y, z)＝F(S, t) exp [k1＋(k2−1) y / 2＋(k2−1)²z / 4] / K ここで，

k1＝r / (σ²/ 2) k2＝(r−d) / (σ²/ 2)

である．すると，BS 方程式は以下のような熱方程式に変換することができる^1）．

（5）

かくして，H(i, j)＝H(y(i), z(j)) とおき，（5）を陽解法で差分近似すると，次の ような結果をえることができる．

1）詳細に関しては，例えばWilmott et al. (1993) などを参照．

∂H

∂z

∂²H

∂y²

＝ 

H(i, j＋1)＝α{H(i＋1, j)＋H(i−1, j)}＋βH(i, j) （6）

ただし，

α＝Δz/Δy² β＝1−2α である．

y の最大値と最小値をymax，yminとそれぞれ表し，またy のステップ幅をΔyと 表すと，グリッド点上のy は，

y(i)＝iΔy＋ymin，（ただし，i＝0, 1, ……, n）

と表すことができる．この式より，グリッド点上のy の最大値と最小値は，

y(n)＝ymax＝nΔy＋ymin

y(0)＝ymin

とそれぞれ表すことができる．

また，z とtの関係に関しては，t＝0_{（現時点）}はz＝Tσ²/ 2に，t＝T（満期時

点）はz＝0 にそれぞれ対応する．zのステップ幅Δz は，その分割数をm とす

ると，

Δz＝Tσ²/ (2m)

と決定でき，グリッド点上のz は，

z(j)＝jΔz，（ただし，j＝0, 1, ……, m）

と表すことができる．この場合，j＝0 は満期時点に，j＝m は現時点にそれぞれ 対応している．

さらに，任意の座標点(i, j) においてオプションを実行したときの価値をP(i, j) と表すと，それは，

P(i, j)＝k3{exp[(k2＋1)y(i) / 2]−exp[(k2−1)y(i) / 2]}，（コールの場合）

P(i, j)＝k3{exp[(k2−1)y(i) / 2]−exp[(k2＋1)y(i) / 2]}，（プットの場合）

と求めることができる．ただし，

k3＝exp[((k2−1)²/ 4＋k1) jΔz]

である．したがって，アメリカン型の場合，z＝0_{（満期時点）}以降の後方計算に

おいて，各グリッド点上におけるオプション価値は，

max[H(i, j), P(i, j), 0], (i＝0, 1, ……, n，j＝0, 1, ……, m)

と決定されていくことになる．ただし，P(0, j) とP(n, j) に関しては，境界条件 より決定する．

以上より，BS 方程式を熱方程式に変換して数値計算を行うことの利点は，ま ず第1 に，（6）においてαとβがy の大きさに依存しないということである．

第2 に，実行価格が陽表的に出てこないため，実行価格がグリッド点上に正確 に乗らないことによる誤差を排除できる．そして第3 に，この実行価格をグリ ッド点上に乗せる手順を省略できるため，プログラミング上の自由度を追加的 に1 つ確保できるということである．かくして，陽解法でオプション価格を数 値計算する場合，熱方程式に変換する方が対数変換の場合よりも優れていると 言うことができる．

2. 4 安定条件

上記の3 つのケースの安定条件について考えてみることにしよう．まず，対 数変換を施した場合から始めると，上記と同様に，x に関するインデックスi の カウントをi＝0, 1, 2,……, n と，またt に関するインデックスj のカウントを

j＝0, 1, 2,……, mとそれぞれ行うことにする．いま，次のようなオプション価

値のベクトルを定義する．

G(:, j)＝(Gn (n−1, j), G(n−2, j), G(n−3, j),……

……, G(2, j), G0(1, j))^T

G(:, j＋1)＝(G(n−1, j＋1), G(n−2, j＋1), ……, G(1, j＋1))^T

ここで，

Gn (n−1, j)＝G(n−1, j)−aG(n, j＋1) G0(1, j)＝G(1, j)−cG(0, j＋1)

である．ただし，G(n, j＋1)，G(0, j＋1) に関しては境界条件より決定する．す ると，（4）より，次式をえることができる．

G(:, j)＝AG(:, j＋1)，( j＝0, 1, ……, m−1) _（7）

ただし，Aは(n−1)×(n−1)の3 重対角行列で，

である．（7）より，後方計算を満期時点から現時点まで繰り返すことによって，

オプション価格を求めることができる．

j＋1 時点の真のオプション価値のベクトルをG−(:, j＋1)，この値からの誤差ベ クトルをE(:, j＋1) とそれぞれ表すと，

G(:, j＋1)＝G−(:, j＋1)＋E(:, j＋1) となるので，（7）の右辺は，

AG(:, j＋1)＝AG−(:, j＋1)＋AE(:, j＋1)

と表すことができる．上式の右辺第2 項は誤差の大きさを表しているが，この 誤差が発散しないためには‖A‖<

−1 となる必要がある．いま，行列Aの任意の 固有値をλiと表すと，Aのスペクトル半径は，

ρ(A)＝max

1<−i<−n−1｜λⁱ｜

と定義できるので，したがって，安定条件‖A‖<

−1 は次のように書き換えるこ とができる．

ρ(A)<−1 （8）

行列Aに対するすべての固有値の範囲に関しては，以下のGerschgorin の定 理を利用することによって確定することができる^2）．いま，n×n の正方行列M

A＝ 

⁝ 

b, c, 0, 0, …  a, b, c, 0, 0, …  0, a, b, c, 0, 0, …  0, 0, a, b, c, 0, 0, …   

  ……  0, 0, a, b, c  0, 0,  ……  0, 0, a, b

⁝ 

⁝ ⁝ 

⁝ 

⁝ 

2）Gerschgorin の定理に関しては，Burden-Faires (2001) を参照．

の成分を{mij}と表すと，複素平面C におけるGerschgorin の円盤Di は次のよ うに定義される．

このとき，行列Mに対するすべての固有値は，

の範囲に限定されることになる．

上記の Gerschgorin の定理を行列 Aに適用すると， この行列に対する

Gerschgorin の円盤Diは，

Di＝

{

｜

−｜a｜＋｜c｜

}

Di＝

{

｜

−｜c｜

}

Di＝

{

｜

−｜a｜

}

となるので，行列Aに対するすべての固有値は，

b−(｜a｜＋｜c｜)<−λⁱ<

−b＋｜a｜＋｜c｜，（ただし，i＝1, 2 , ……, n−1）

の範囲に限定されることになる．安定条件は，この結果と（8）より，

−1<−b−(｜a｜＋｜c｜) かつ1>−b＋｜a｜＋｜c｜ となる．さらに，これらの不等式は，

rΔt<

−2 かつ2σ²/Δx²<

−r

と書き換えることができるので，かくして，対数変換を施した場合の安定条件 は，

σ²Δt /Δx²<

−1 _（10）

と求めることができる．

次に，BS 方程式を陽解法で直接近似した場合の安定条件について考えてみる ことにしよう．この場合は，（10）を使うことによって容易に求めることができ る．すなわち，x＝lnS なので，

Dⁱ＝

{

｜

^Σ  

}

n j＝1 j≠i

中心：mⁱⁱ，半径：

Σ  

n j＝1 j≠i

∪Di n i＝1

Δx＝ΔS/ S

となり，これを（10）に代入すると，安定条件は，

σ²S²Δt/ΔS²<

−1

と求めることができる．上式の左辺の値が最大となるのは，S が最大値をとる ときである．S の最大値はグリッド点上ではnΔS＋Sminとなるので，かくして，

BS 方程式を直接近似した場合の安定条件は，

(nΔS＋Smin)²σ²Δt/ΔS²<

−1 _（11）

と書き換えることができる．とくにSmin＝0 と設定した場合，（11）の安定条件は，

n²σ²Δt<

−1 となる．

最後に，熱方程式に変換した場合の安定条件について考えてみることにしよ う．y に関するインデックスi のカウントをi＝0, 1, 2, ……, n と，またz に関す るインデックスj のカウントをj＝0, 1, 2, ……, m とそれぞれ行うことにする．

いま，次のようなオプション価値のベクトルを定義する．

H(:, j)＝(H(n−1, j ), H(n−2, j ), ……, H(1, j ))^T

H(:, j＋1)＝(Hn(n−1, j＋1), H(n−2, j＋1), H(n−3, j＋1),……

……, H(2, j＋1), H0(1, j＋1))^T ここで，

Hn (n−1, j＋1)＝H(n−1, j＋1)−αH(n, j) H0(1, j＋1)＝H(1, j＋1)−αH(0, j)

である．ただし，H(n, j＋1)，H(0, j＋1)に関しては境界条件より決定する．す ると，（6）より，次式をえることができる^3）．

H(:, j＋1)＝BH(:, j)，(j＝0, 1, ……, m−1) _（12）

ただし，Bは(n−1)×(n−1) の3 重対角行列で，

3）本文のインデックスjに関しては，対数変換を施した場合，現時点がj＝0，満期時点がj＝mと

それぞれなるのに対して，熱方程式に変換した場合は，現時点がj＝m，満期時点がj＝0 と逆にな ることに注意せよ．

である．（12）より，後方計算をz＝0（満期時点）からz＝Tσ²/ 2_{（現時点）}まで 繰り返すことによって，オプション価格を求めることができる．さきの場合と 同様に考えると，この場合の安定条件は，

ρ(B)<−1 （13）

となる．

行列Bに対するGerschgorin の円盤Diは，

Di＝

{

｜

}

Di＝

{

｜

−α

}

となるので，Gerschgorin の定理より，行列Bに対するすべての固有値は，

1−4α<

−λⁱ<

−1，（ただし，i＝1, 2, ……, n−1）

の範囲に限定されることになる．かくして，この結果と（13）より，BS 方程式 を熱方程式に変換した場合の安定条件は，

α<

−1 / 2 _（15）

と求めることができる．

（10），（11）および（15）の各安定条件より明らかなように，陽解法で計算 精度を上げるためにΔS（Δx, Δy）を小さく（n を大きく）しようとする場合，そ れ以上にΔt（Δz）を小さく（m を大きく）する必要がある．精度を上げようとす る場合，直接近似の場合の計算効率の悪さは明らかなので，以下では，パリジ ャン・オプションの評価において，熱方程式への変換と対数変換の場合につい

B＝ 

β,  α,  0,  0,  …  α,  β,  α,  0,  0,  …   0,  α,  β,   α,  0,  0,  …   0,   0,  α,   β,    α,  0,   0,  …   

   ……  0,  0,    α,  β,  α   0,   0,  ……         0,   0,  α,  β   

⁝ 

⁝ 

⁝ 

⁝ 

⁝  ⁝ 

て比較検討を行うことにする．

3 パリジャン・オプションの評価

3. 1 パリジャン・オプションの特徴

ここでは，ノックアウト型のパリジャン・オプションを例に上げながら議論 を進めることにしよう．原資産価格がバリア外にある時間を「バリア時間」，ノ ックアウトされるバリア時間の長さを「トリガー時間」とそれぞれ呼ぶことに する．パリジャン・オプションは，このバリア時間のカウントの仕方によって2 つに分類される．すなわち，バリア時間を累積的にカウントするという場合と，

原資産価格がトリガー時間前にバリア内に戻ればその都度バリア時間をゼロに リセットするという場合とである．前者は累積的パリジャン・オプション，後 者は継続的パリジャン・オプションとそれぞれ呼ばれる．以下では，原資産価 格は連続的にモニターされるものと仮定し，また説明の便宜上，継続的パリジ ャン・オプションの場合に限定して考察を進めることにする^4）．

t時点における原資産価格をS(t)，原資産価格に対する上限のバリア水準をSu， 下限のバリア水準をSdとそれぞれ表すと，原資産価格がt 時点までにバリア外 に連続的にある時間，すなわちバリア時間τは，

τ＝t−sup{t' <−t｜S(t' )>−Su}，（アップ・アンド・アウト型の場合）

τ＝t−sup{t' <−t｜S(t' )<−Sd}，（ダウン・アンド・アウト型の場合）

と表すことができる．継続的パリジャン・オプションの価値がゼロとなるτの 値，すなわちトリガー時間をω（プラスの定数）と表すと，τ＝ωとなる前に原 資産価格がバリア内に戻った場合，τはその時点でゼロにリセットされる．し たがって，τに関しては，τ＝ωとなるまでは次のようにカウントされること になる．

4）累積的パリジャン・オプションのプログラムは，継続的パリジャン・オプションの場合のマイナ ーな修正だけで済むので，以下の本文では，継続的パリジャン・オプションのケースに限って説明 を行うことにする．

dτ S＞Su

dτ＝ −τ- S＝Su，（アップ・アンド・アウト型の場合）

0 S＜Su

0 S＞Sd

dτ＝ −τ- S＝Sd，（ダウン・アンド・アウト型の場合）

dτ S＜Sd

ここで，τ-はτがゼロにリセットされる直前のτの値である．パリジャン・オ プションの価値をF(S, t,τ) と表せるものとすると，BS 方程式は，原資産価格 がバリアの範囲内（S＜SuまたはS＞Sd）にある場合は（1）のままであるが，バリ アの範囲外（S＞SuまたはS＜Sd）にある場合は，

となる．ただし，S＝Su, S＝Sd のときは，次のようなオプション価値の連続性

（跳躍条件）を課すことにする^5）．

F(Su, t,τ)＝F(Su, t, 0)，（アップ・アンド・アウト型の場合）

F(Sd, t,τ)＝F(Sd, t, 0)，（ダウン・アンド・アウト型の場合）

また，満期時点における終期条件は次のようになる．

F(S, T,ω)＝0

F(S, T,τ)＝max[S−K, 0]，（τ＜ω，コール・オプションの場合）

F(S, T,τ)＝max[K−S, 0]，（τ＜ω，プット・オプションの場合）

ただし，Kは実行価格である．

この継続的パリジャン・オプションは，ω＝dt のときは通常のワンタッチ型の バリア・オプションに，ω＝T のときはプレーンバニラ・オプションにそれぞれ 一致する．この点は，バリア時間が累積的にカウントされる累積的パリジャ ン・オプションについても言える．したがって，これら2 種類のパリジャン・

∂F

∂t

∂F

∂S

∂² F

∂² S

∂F

＋    ＋  σ∂τ  1 ²S²     ＋(r−d) S    −rF＝0 2

5）跳躍条件に関しては，Wilmott et al. (1993) などを参照．

オプションの価格は，ともにワンタッチ・バリア・オプションよりも大きく，

プレーンバニラ・オプションよりも小さくなるということが言える^6）．

3. 2 数値計算のためのアルゴリズム

継続的パリジャン・オプションに対する陽的有限差分法の適用例について，

説明の便宜上，ダウン・アンド・アウト型のケースに限って考えることにしよ う．また，対数変換した場合の数値計算法は容易にパリジャン・オプションに 適用できるので，熱方程式に変換した場合についてのみ説明を行うことにする．

まず，グリッド点上に下限のバリアSdが正確に乗るように，y のステップ幅 を決定することにしよう．このバリア値は，熱方程式に変換した場合では，

yd＝ln (Sd/ K)

となる．yの最小値y^minからy^dまでの分割数n^d（自然数）を与え，yのステップ幅 を，

Δy＝(y^d−y^min) / n^d （16）

と決定すれば，グリッド点上にydを正確に乗せることが可能となる．したがっ て，グリッド点上でy^dは，

yd＝ndΔy＋ymin

と表すことができ，さらにy^minからy^maxまでの分割数は，

n＝Round [(ymax−ymin) /Δy＋0.5] （17）

と決定することができる．ただし，Round [・] は変数を四捨五入して整数化する 関数である．

上記のように n とΔy を決定したもとで，現時点（z＝Tσ²/ 2）から満期時点

（z＝0）までの分割数m と Δz を求める場合は，安定条件（15）を満たすように 求めなければならない．すなわち，zのステップ幅はΔz＝σ² T/ (2m) と決定すれ ばよいので，0＜α<

−1 / 2の範囲から任意のα（＝Δz/Δy²）の値を選択すると，

6）継続的パリジャン・オプションと累積的パリジャン・オプションを比べた場合は，バリア時間の リセットの可能性がある分，前者の価格の方が高くなる．

z に関する分割数は，

m＝Round [σ² T/ (2αΔy²)＋0.5]，（ただし，0＜α<

−1 / 2） （18）

と求めることができる．そして，このmを使ってΔzを決定することによって，

安定条件が満たされるということになる^7）．

以上の手順を要約すると次のようになる．まず，nd，yd，ymin，ymaxの4 つの値 をそれぞれ与えることによって，（16），（17）より，Δy とn を決定する．そし て，0＜α<

−1 / 2 の範囲内でαの値を選択し，（18）より，m を求めてΔz を決

定するということになる．

トリガー時間ωに関しては，Δπ＝Δz であることを考慮すると，グリッド点 上では次のように変換されることになる．

mbΔz ただし，

mb＝Round [_σ²ω/ (2Δz)]

である．したがって，バリア時間は，グリッド点上では，次のようにカウント されることになる．

π(k)＝kΔz，（k＝0, 1, 2, ……, mb）

次に，配列の次元を減らす工夫を考えてみることにしよう．周知のように，

配列の次元は低いほどメモリ節約的で，計算速度を速くする．そこでまず，パ リジャン・オプションの評価関数F(S, t,τ) を熱方程式に変換したものをH(y, z,π) と表すことにする．これは，グリッド点上では，H(i, j, k)＝H(y(i), z(j),π

(k)) のように，3 次元配列で表されることになる．しかし，（6）から分かるよう

に，後方計算においては，z(j＋1) 時点のH(・) はz(j) 時点のH(・) のみで決定で きるため，z を省略した2 次元配列を2 つ使うことによって，3 次元配列をな くすことができる．

いま，任意の2 時点を考え，z(j) 時点の評価関数をH1(i, k)＝H1(y(i),π(k)), z

7）本文の（18）より，αの値が小さければ小さいほどmが大きくなる（Δzが小さくなる）という

ことが分かる．

(j＋1) 時点のそれをH2(i, k)＝H2(y(i),π(k))とそれぞれ表すことにする．これら2

つの2 次元配列を使うと，例えば原資産価格がバリア外にある領域では，

H2(i, k)＝α{H1(i＋1, k＋1)＋H1(i−1, k＋1)}＋βH1(i, k＋1)

と計算することができる．さらに，ここで算出されたz(j＋1) 時点のH2(i, k) は，

次のz(j＋2) 時点ではH1(i, k) として計算に利用できるため，後方計算において

は， H2(i, k) を順次H1(i, k) に書き換えていくことにより，効率的な計算を実現

することが可能となる．

満期時点（z＝0）におけるペイオフ値の計算に関しては，Visual Basic .net

（Ver.2003）によるプログラム例で示すと，次のようにして行うことができる．

For k = 0 To m For i = 0 To n

H1(i, k) = 0 y = i * dy + ymin

P = a1 * (Math.Exp(0.5 * (k2 + 1) * y) _ - Math.Exp(0.5 * (k2 - 1) * y)) If k <= mb And P > 0 Then H1(i, k) = P Next i

Next k

For i = 0 To nd H1(i, mb) = 0 Next i

ここで，H1(i, k)＝H1(i, k)，dy＝Δy，ymin＝ymin，k2＝k2，mb＝mb，nd＝ndで ある．また，a1=1 のときはコール・オプションの場合，a1=−1 のときはプッ ト・オプションの場合である．さらに，満期時点（z＝0）から現時点（z＝Tσ² / 2）

までの後方計算に関しては，例えば補論で示すようなプログラムを作成するこ とによって実行することができる．

なお，現時点の原資産価格S(0) に対するオプション価格に関しては，必ずし

もグリッド点上に正確に乗っているわけではない．この点は，以下のような線 形補間法の適用によって極めて正確に補うことができる．

y1 = Math.Log(S0 / K) c1 = (y1 - ymin) / dy c2 = Int(c1)

H3 = H2(c2, 0) + (H2(c2 + 1, 0) - H2(c2, 0)) * (c1 - c2) ただし，S0＝S(0)である．

3. 3 計算結果

まず，パリジャン・オプション価格の計算値の収束率に関して，BS 方程式を 熱方程式に変換した場合と原資産価格を対数変換した場合との比較を行ってみ ることにしよう．第 1 表はヨーロピアン型で，第 2 表はアメリカン型でそれぞ れ計算値が収束する様子を比較したものである．表中のH は熱方程式に変換し た場合を，L は対数変換した場合をそれぞれ表している．ただし，アメリカン 型の場合はプット・オプションの場合だけで比較を行っている．ここでパラメ ータは，S(0)＝100，σ＝0.2，r＝0.02，d＝0.00，T＝0.5，B^d＝90，B^u＝110，

ω＝0.01 である．また，実行価格に関しては，コールの場合はK＝95，プット

の場合はK＝105 である．さらに，熱方程式に変換した場合と対数変換した場

合の収束率を比較するために，x，y のステップ幅をΔx＝Δy＝0.0045 と同一に したもとで，時間分割数（ただし，熱方程式に変換した場合はΔz の数）を漸次増や していくという方法をとっている．ただし，Δx＝Δy とするために，x のステ ップ数は1177（ダウン・アンド・アウトのとき），1178（アップ・アンド・アウトのと き），y のステップ数は848（アップ・アンド・アウトのとき），889（ダウン・アン ド・アウトのとき）とそれぞれ設定している．なお，表中の収束値に関しては，

熱方程式へ変換した場合に外挿法を適用して算出したものである^8）．

8）外挿法に関しては，久保（2003）を参照．

第1 表と第2 表より分かることは，まず第1 に，熱方程式に変換した場合と 対数変換した場合のどちらにおいても，計算値がほぼ一様に収束しているとい うことである．このことは，これら2 つの場合において外挿法を利用すること ができるということを表している．第2 に，これら2 つの場合の同じ時間分割 数に対しては，ほぼ同じ計算値をえているということである．したがって，y の ステップ数の方がxより少ない分，計算効率上，熱方程式に変換した方が優れ ていると言うことができる．

第 1 表 ヨーロピアン型パリジャン・オプションの評価 

第 3 表と第 4 表は，ヨーロピアン型とアメリカン型の場合で，トリガー時間 ωのパリジャン・オプション価格に与える効果について比較したものである．

ただし，これら表中のオプション価格に関しては，熱方程式へ変換した場合に 第 2 表 アメリカン型パリジャン・オプションの評価 

第 3 表 ヨーロピアン型パリジャン・オプションとトリガー時間 

外挿法を適用することにより算出を行っている．2 つの表より明らかなように，

トリガー時間が長くなるほどパリジャン・オプションの価格が大きくなるとい うことが分かる．また，ヨーロピアン型の方がアメリカン型に比べて，トリガ ー時間の長さにより強く影響されるということが分かる．

4 お わ り に

本稿では，陽的有限差分法（陽解法）およびこれを用いたパリジャン・オプシ ョンの評価について考察を行ってきた．具体的には，まず第1 に，3 つのケー スの陽的有限差分法を使ってオプションを評価する場合の安定条件を，

Gerschgorin の定理を使って導出し比較検討を行った．そして第2 に，HSW あ

るいはWilmott の数値計算アルゴリズムの改良と，それを基にした様々な種類

のパリジャン・オプションの評価を試みた．

明らかにした点は次の3 点である．まず第1 に，陽的有限差分法でオプショ ンを評価する場合，BS 方程式を熱方程式に変換して数値計算した方が，原資産 価格を対数変換した場合よりも計算効率を向上させることができる．第2 に，

熱方程式に変換した場合と対数変換した場合のどちらにおいても，パリジャン・

オプション価格の計算値は比較的滑らかに収束値に向かう．したがって，これ らの変換を施した場合は外挿法を適用することができ，効率的な計算を実現す ることが可能となる．そして第3 に，トリガー時間のパリジャン・オプション 価格に与える影響は，ヨーロピアン型の方がアメリカン型よりも大きくなる．

本稿では，原資産価格が連続的にモニターされるパリジャン・オプションの 第 4 表 アメリカン型パリジャン・オプションとトリガー時間 

評価を行ってきた．しかしながら，現実の原資産価格のモニターあるいはバリ ア時間のカウントは，連続的にではなく離散的に行われる^9）．この点も含め，よ り実際的なパリジャン・オプションの評価に関しては他日を期したい．

補論

yn = n * dy + ymin For j = 1 To m

z = j * dz

k3 = Math.Exp((0.25 * (k2 - 1) * (k2 - 1) + k1) * z) For k = 0 To mb

H2(0, k) = a2 * k3 * (Math.Exp(0.5 * (k2 - 1) * ymin) _ - Math.Exp(0.5 * (k2 + 1) * ymin)) H2(n, k) = a3 * k3 * (Math.Exp(0.5 * (k2 + 1) * yn) _

- Math.Exp(0.5 * (k2 - 1) * yn)) Next k

For i = nd To n - 1

H2(i, 0) = (1 - 2 alpha ) * H1(i, 0) _

+ alpha * (H1(i + 1, 0) + H1(i - 1, 0)) y = i * dy + ymin

P = a1 * (k3 * (Math.Exp(0.5 * (k2 + 1) * y) _ - Math.Exp(0.5 * (k2 - 1) * y))) If H2(i, 0) < P Then H2(i, 0) = P

Next i

For k = 0 To mb - 1 For i = 1 To nd - 1

9）この点に関しては，例えばTavella-Randall (2000) などを参照．

H2(i, k) = (1 - 2 alpha ) * H1(i, k + 1) _ + alpha * (H1(i + 1, k + 1) _ + H1(i - 1, k + 1))

y = i * dy + ymin

P = a1 * (k3 * (Math.Exp(0.5 * (k2 + 1) * y) _ - Math.Exp(0.5 * (k2 - 1) * y))) If H2(i, k) < P Then H2(i, k) = P

Next i Next k

For i = 0 To nd - 1 H2(i, mb) = 0 Next i

For k = 1 To mb - 1 H2(nd, k) = H2(nd, 0) Next k

For k = 0 To mb For i = 0 To n

H1(i, k) = H2(i, k) Next i

Next k Next j

（注）アメリカン型パリジャン・オプション（ダウン・アンド・アウト）の評価に おける後方計算のプログラム例．使用言語はVisual Basic .net (Ver.2003)．

ここで，ymin＝ymin，k1＝k1，k2＝k2，k3＝k3，alpha＝α，mb＝mb，nd＝

nd，H1(i, k)＝H1(i, k)，H2(i, k)＝H2(i, k)である．また，a1＝1，a2＝0，

a3＝1 のときはコール・オプションの場合，a1＝−1，a2＝1，a3＝0 のと

きはプット・オプションの場合である．

【参考文献】

Burden, R. L., and J. D. Faires, (2001) Numerical Analysis, Seventh Edition, Brooks/Cole.

Geske, R., and K. Shastri, (1985) Valuation of Approximation: A Comparison of Alternative Approaches, Journal of Financial and Quantitative Analysis, Vol.20, No.1, pp.45-72.

Haber, R. j., P. J. scho¨nbucher, and P. Wilmott, (1999) Pricing Parisian Options, The jour- nal of Derivatives, Vol.6, No.2, pp.71-79.

Kwok, Y. K., and K. W. Lau, (2001) Pricing Algorithms for Options with Exotic Path- Dependence, The journal of Derivatives, Vol.9, No.1, pp.28-38.

Tavella, D., and C. Randall, (2000) Pricing Financial Instruments: The Finite Difference Method, John Wiley ＆Sons.

Vetzal, K. R., and P. A. Forsyth, (1999) Discrete Parisian and Delayed Barrier Options: A General Numerical Approach, Advances in Futures and Options Research, Vol.10, pp.1-15.

Wilmott, P., (2000) Paul Wilmott on Quantitative Finance, John Wiley ＆Sons.

Wilmott, P., J. Dewynne, and S. Howison, (1993) Option Pricing: Mathematical Methods and Computation, Oxford Financial Press.

久保徳次郎，（2003）「ルックバック・オプションの評価と外挿法」『経済学論叢』（同志 社大学）第55 巻第2 巻，pp.1-19.

The Doshisha University Economic Review Vol.57 No.2 Abstract

Tokujiro KUBO, Pricing Parisian Options Using the Explicit Finite Difference Method

Parisian options are the barrier-style options which are ‘knocked out’

(‘knocked in’) only if the price of the underlying asset has been outside (inside) the barrier for more than a certain prescribed time. This paper revises the pricing algorism for the consecutive Parisian options of Haber-Schonbucher-Wilmott (1999). And it presents that the transformation of Black-Sholes equation into heat equation is most efficient in case of the numerical calculation of the Parisian options with the explicit finite difference method, using the revised algorism.