ファイナンスに関する数値計算 資料
一橋大学商学研究科 小林健太

2015/4/13 10:56 更新

0

はじめに

0.1

この資料の概要

ファイナンスに関係する数値計算を一通りできるようになることを目的とする．確率論に関する知識はあまり前提とせず，あまり数学的な厳密性は突き詰めない予定．
この資料はhttp://cm.hit-u.ac.jp/~kobayashi/lecture/からダウンロードできる．随時，更新する予定．

0.2

参考文献

プログラミングには C 言語を用いるので，その勉強は各自でやってもらいたい．C 言語の参考書としては
柴田望洋「新版 明解 C 言語 入門編」（ソフトバンククリエイティブ）
皆本晃弥「やさしく学べる C 言語入門−基礎から数値計算入門まで」（サイエンス社）を勧める．初心者の人は，これらの本に載っているサンプルコードを順に動かしてみるとよい．
ファイナンスや数値解析の理論面の話については，この講義資料に全て書いていく方針なので，特に教科書等は用意しなくてよいが，この資料で取り上げるトピックについてはJohn C. Hull, “ファイナンシャルエンジニアリング”, 金融財政事情研究会
Steven E. Shreve, “ファイナンスのための確率解析 (1)”, シュプリンガーフェアラーク東京
Peter E. Kloeden, Eckhard Platen “Numerical Solution of Stochastic Diﬀerential Equations”, Springer
森正武, “数値解析”, 共立出版
杉原正顕, 室田一雄, “数値計算法の数理”, 岩波書店
齊藤宣一, “数値解析入門”, 東京大学出版会などの書籍や，個別の論文などから持ってくる予定である．

1

0.3

コンパイラ

数値解析においては，理論の理解も大切だが，実際にプログラムを組んで理論を確かめることも重要である．この資料では，プログラムは C 言語（C++）を用いる．パソコン上で C 言語（C++）を使うにはコンパイラをインストールする必要があるが，フリーのものでも各種の選択肢がある．
初心者にとっては Borland C++ Compiler がインストールも簡単で良いだろう．ただ，仕様が少し古いので，最新技術を用いたプログラミングには少々支障がある（とはいっても，初心者の間は全く関係がないが）．

http://www.embarcadero.com/jp/products/cbuilder/free-compiler
こちらからダウンロードできる．ダウンロードには登録が必要だが，登録しても宣伝メールが送られてくるとかの不都合は特に生じない．

MinGW は軽くてお勧めだが，若干，インストールに苦労することがあるかもしれない．
Visual Studio は開発環境などが色々と揃っていて，そこそこ本格的なプログラミングも可能である．フリーだが登録が必要．
Cygwin は Windows 上の仮想 Linux 環境だが，インストールすると C++コンパイラが付いてくる．サーバー管理等，Linux に興味のある人はこれも良いだろう．
他にも，Open Watcom C/C++，Digital Mars C/C++ Compilers，なんてのもある．

Mac については，私は詳しくは知らないが，色々とやり方はあるようだ．Windows でな
く，最初から Linux 機でプログラムするのも手かもしれない．

0.4

グラフ表示

数値計算結果をグラフ化して表示するために，gnuplot をインストールする必要がある．

http://www.gnuplot.info/
こちらから Download from SourceForge のリンクを辿るとダウンロードサイトに行けるが，バナー広告で「Download」とか書かれているのが出たりするので間違ってクリック
しないように．

32bit ウィンドウズ機なら「gp***-win32-setup.exe」を，64 ビットウィンドウズ機なら「gp***-win64-setup.exe」を選ぶと良いだろう（ただし，***はバージョン番号で数字が
入る）．

2

0.5

エディタ

プログラムを作成するのにエディタもあった方がいい．Windows に付属のメモ帳でもできるが，Terapad や秀丸エディタなどを使うとプログラミングが格段にやり易くなる．

3

1

中心極限定理中心極限定理は，確率論や統計学で最も重要と言われる定理であり，ファイナンス理論においても各種の議論の根拠となるので，ここで確認しておく．
中心極限定理によれば，平均 µ, 分散 σ 2 の独立同分布な確率変数の列を X1 , X2 , · · · とすると，Sn = X1 + X2 + · · · + Xn は，平均 nµ，分散 nσ 2 の正規分布に漸近する．ただし，
分散 σ 2 の有限性は仮定する．もっと正確に言うと
)
(
∫ a
1
Sn − nµ2√e−x /2 dx≤a = √
lim P
n→∞
nσ
2π −∞
が成り立つ．
中心極限定理の証明を考えるには，独立な確率変数の和がどのような性質を持つか調べる必要がある．そこで，まずは 2 つの確率変数を考える．
互いに独立な連続確率変数を X, Y とし，それぞれの確率密度関数を fX , fY とする．このとき，Z = X + Y の確率密度関数 fZ はどのような形になるだろうか？ 少し考えると，
fZ は fX と fY の畳み込みで求まることがわかる．すなわち∫ ∞fZ (x) =fX (t)fY (x − t)dt−∞となる．
問題 1.1 X と Y が正の数 h の倍数のみを取り得る離散確率変数の場合でも，Z = X + Yの分布が畳み込みになることを示せ．
畳み込みは扱いが面倒なのだが，Fourier 変換を用いると容易に扱えるようになる．確率密度関数 fX を持つ確率変数の Fourier 変換を∫ ∞φX (z) =eizx fX (x)dx−∞と定義する（一般的によく用いられる Fourier 変換とは，係数や izx の符号がちょっと違うが，本質的には変わらない）．Y および Z の Fourier 変換も同様に定義すると

φZ (z) = φX (z)φY (z)が成り立つことがわかる．つまり，Z の確率密度関数は畳み込みを用いなければ求まらなかったが，Z の Fourier 変換は単に X の Fourier 変換と Y の Fourier 変換を掛け合わせる
だけで求まる．
問題 1.2 それを証明せよ．

4

さて，

∫ ∞
−∞

eizx fX (x)dx

は e を X の確率密度で重み付けして積分したものであるから，E[eizX ] と書くこともできる．この表記だと，X が離散確率変数の場合でも大丈夫である．よって，以降は
izx

φX (z) = E[eizX ]と定義する．これを確率変数 X の特性関数という．特性関数をこのように書けば，
E[eizZ ] = E[eiz(X+Y ) ] = E[eizX eizX ] = E[eizX ]E[eizY ] より φZ (z) = φX (z)φY (z) はすぐにわかる．
問題 1.3 P (X = 0) = p, P (X = 1) = q = 1 − p なる Bernoulli 分布に従う確率変数の特性関数を求めよ．
問題 1.4 [0, 1] 上の一様分布に従う確率変数の特性関数を求めよ．
問題 1.5 標準正規分布に従う確率変数の特性関数を求めよ．
特性関数が等しければ確率分布も等しいし，ある特性関数が別の特性関数に近付いて行くなら，対応する確率分布も近付いていくことが知られている．もちろん，これは無制限に成り立つわけではないので，関数の範囲は一定の範囲に制限する必要があるし，関数が関数に近付くという意味も，適切な関数空間上で考えないと意味をなさない．しかし，そこはあまり深入りしないことにする．
それでは中心極限定理の証明（というか説明）に入ろう．独立同分布な確率変数の列をX1 , X2 , · · · とし，その平均を µ, 分散を σ 2 とする．このとき，

Yn =Sn − nµ√nσの特性関数が標準正規分布の特性関数に近付くことを言えばよい．
√

E[eizYn ] = E[eiz(Sn −nµ)/( nσ) ]n∏√=E[eiz(Xk −µ)/( nσ) ]
k=1
(
)n
√
= E[eiz(X1 −µ)/( nσ) ]( [])niz(X1 − µ) z 2 (X1 − µ)2√= E 1+
−
+ Rn
2nσ 2
nσ
)n
(
z2
2
+ ERn
−→ e−z /2= 1−2nより，結論が得られる．ただし，Rn は n の −3/2 乗以下の項からなる．
問題 1.6 上の式変形はかなり大雑把であるので，数学的に厳密にせよ．

5

2

二項モデル

2.1

導入

本章では，金融派生商品の一番基本的なモデルである二項モデルについて説明する．まずは一期間で考え，一期間における無リスク金利を r とする．また，完全市場を仮定する．
さて，原市場 S を考える．一単位あたりの S の価値は，時刻 0 では S0 であり，時刻 1 では
99 ％の確率で上昇して SH になり，1 ％の確率で下落して SL になるものとする．さらに金融派生商品 V を考える．V は，時刻 1 で S が SH になれば VH 払い戻され，時刻 1 で S
が SL になれば VL 払い戻される．さて，V の時刻 0 での価値 V0 はいくらになるだろうか？
普通に考えれば，V は 99 ％という高い確率で VH の払い戻しを受けられるので，時刻 0 に
VH
おける価値は，VH を無リスク金利で割り引いたに近い値になると思うかもしれな
1+r
い．しかし，ファイナンス理論においては少し違った考え方をする．

SH VL − SL VH
VH − VL
単位，無リスク資産を単位所持すれば，
SH − SL(1 + r)(SH − SL )V を複製することができる．よって，時刻 0 において複製に要する資金は，

時刻 0 において S を

VH − VLSH VL − SL VH· S0 +SH − SL(1 + r)(SH − SL )となり，無裁定の原則から，これが V0 に等しくなる．

V0 の価格が実際の S の変動確率によらないというのは奇妙に思うかもしれないが，複製によって同じものを作るのに要する費用をその価値とする，という決め方だとそうなるのである．
さて，このようにして求まった V0 をもう少しよく見ると)(SH − (1 + r)S0
1
(1 + r)S0 − SL
· VH +
· VL
V0 =
1+r
SH − SLSH − SLと書けることがわかる．つまり，V0 は，上昇確率 p =(1 + r)S0 − SLと下落確率
SH − SLSH − (1 + r)S0
で時刻 1 における期待値を求め，これを金利 r で割り引いたものとしSH − SLて求めることができる．この仮想的な確率をリスク中立確率という．
q=

リスク中立確率においては，あらゆる金融商品の価格の期待値は無リスク資産に等しくな
る．実際，S についても

pSH + qSL =SH − (1 + r)S0
(1 + r)S0 − SL
· SH +· SL = (1 + r)S0
SH − SLSH − SL6

となっている．
リスク中立確率は仮想的な確率なので実際の確率とは異なる．なぜなら，リスク資産の期待値が無リスク資産の期待値と等しいならば，だれもリスク資産には投資しないので，リスク資産の期待値は無リスク資産よりは大きくなっている筈だからである．しかし，独立な動きをするリスク資産の選択肢が増えれば増えるほど，それらを組み合わせることによってリスクを減らしていくことができるので，リスク資産の期待値は無リスク資産の期待値に近付いて行く（詳しくは現代ポートフォリオ理論を勉強すること）．また，リスク資産の期待値が無リスク資産の期待値よりもどの程度大きくなるべきなのかは，投資家のリスク選好度にも依存するので，数学的な定式化が難しい（心理学なんかも関係しそう）．以上のような理由で，リスク中立確率は実際の確率とは異なるが，ファイナンス理論では，理論の組み立てやすさや客観性の点で有利であることを考慮して，リスク中立確率を用いるのである．

2.2

ランダウの記号

今後よく使うので，ランダウの記号について説明しておこう．|x| が十分小さいときに，ある定数 M > 0 が存在して |f (x)| < M |g(x)| が成り立つことをf (x) = O(g(x))と書くことにする．この O(·) をランダウの記号と言う．例えば

cos x = O(1)sin x = O(x)e − 1 − x = O(x2 )xx3 = O(x2 )1= 1 − x + x2 + O(x3 )1+xなどが成り立つ．簡単に言うと，O(xs ) は，x の s 乗以上のオーダーで減少する項を表し
ている．

2.3

二項ツリー

本章では，前章で説明した 1 期間モデルを多期間に拡張し，オプションの数値計算を行う．
原市場のボラティリティーを σ とし，現時点を t = 0，満期を t = T とする．無リスク金利を r とし，原市場の配当は考えない．確率は全てリスク中立確率で考える．

7

T, tm = m∆t とする．また，ある時間ステッM
プにおいて原市場価格が S であるとするとき，次の時間ステップには確率 p で uS ，確率q = 1 − p で dS となるものとする．ここで，u > 1 > d > 0 であるが，計算を簡略化するために d = 1/u を仮定する．このとき，p, u をどのように決めればよいか考えよう．
近似のため，時間方向を M 分割し，∆t =r はいわゆる瞬間利子率と考えるので，∆t 間に無リスク資産は er∆t 倍になる．よって，リスク中立の考え方から
pu + qd = er∆tが成り立つ必要がある．p + q = 1 より，これを解くと

p=er∆t − d,u−dq=u − er∆tu−dとなる．また，log S の ∆t 間における変動の分散は，中心極限定理より σ 2 ∆t にならねばならないのでp(log u − r∆t)2 + q(log d − r∆t)2 = σ 2 ∆tが成り立つ必要がある．上で求めた p, q を用いると(er∆t − d)(log u − r∆t)2 + (u − er∆t )(log d − r∆t)2 = (u − d)σ 2 ∆tとなる．ここでlog d = −εlog u = ε,と置くと

(er∆t − e−ε )(ε − r∆t)2 + (eε − er∆t )(−ε − r∆t)2 = (eε − e−ε )σ 2 ∆tとなり，展開すると2(eε + e−ε − 2er∆t )εr∆t = (eε − e−ε )(σ 2 ∆t − r2 ∆t2 − ε2 )となる．このままだと厳密に解くのは難しいので，ε の近似的な値を求める．
(
)
ε3
2
4
2
5
2(ε − 2r∆t + O(ε ) + O(∆t ))εr∆t = 2ε + + O(ε ) (σ 2 ∆t − ε2 + O(∆t2 ))3より

2ε2 r∆t = 2σ 2 ∆t − (2 + O(∆t))ε2 + O(ε4 ) + O(∆t2 )となる．これよりε2 =σ 2 ∆t + O(ε4 ) + O(∆t2 )= σ 2 ∆t + O(ε4 ) + O(∆t2 )1 + O(∆t)

8

となり

√

√
σ 2 ∆t + O(ε4 ) + O(∆t2 ) = σ ∆t + O(∆t−1/2 ε4 ) + O(∆t3/2 )√が成り立つので，近似的に ε = σ ∆t と取れば，誤差項はε=O(∆t−1/2 ε4 ) + O(∆t3/2 ) = O(∆t3/2 )となり，ε を相対誤差 O(∆t) で近似できることになる．ここで，相対誤差とは，x の近似e−xxe についてxのことを言う．
x
以上より，√u = eσ ∆t ,√d = e−σ ∆t ,p=er∆t − d,u−dq=u − er∆tu−dと取るのが適当であるということがわかった．これは Cox, Ross and Rubinstein(1979) で得られているものと同じである．
ここで行った O(·) を用いる議論は，漸近展開という手法の簡単な例である．Taylor 展開の項をどこまで取るかは，数学的な勘，もしくは試行錯誤が必要となる．
問題 2.1 d = 1/u を仮定せず，p = q =当か？

1を仮定すると，u と d はどのように取るのが適2

それでは，以上の結果を用いてヨーロピアンオプションの価格を計算してみよう．まず，√
原市場の離散的な価格を Sn = S0 enσ ∆t と設定し，時刻 tm = m∆t において原市場価格がSn であるときのオプション価格を Ym,n とする．
リスク中立の考え方から，時刻 t = tm におけるオプション価格は，t = tm+1 における期待値を無リスク金利で割り引いたものになるので，Ym,n = e−r∆t (pYm+1,n+1 + qYm+1,n−1 )が成り立つ．t = tM のときにはオプション価格はペイオフ関数と等しいとし，順に時間を遡りながら計算すると，t = 0 での Ym,n の値が求められる．

t = 0 において S−N から SN の範囲でオプション価格を求めたい場合には，時刻 t = M においては価格刻みを S−(M +N ) から SM +N の範囲に取る必要がある．
問題 2.2 適当に状況を設定してヨーロピアンオプションの価格を計算し，時刻 t = 0 から満期まで価格がどのように変化していくかをグラフに描け．

9

問題 2.3 適当に状況を設定して，ヨーロピアンオプションの時刻 t = 0 における一つの値が，どのように時間ステップ幅に依存して変化するか調べよ．時間ステップ幅を半分にしていきながら計算するとよい．
プログラミングにおける注意：
時刻 t = 0 のときの価格のみ求めたいなら，無リスク金利による割り引きは最後にまとめて行ってもよい．しかし，アメリカンオプションの計算では，毎回割り引く必要がある．
データの保存について，全ての時間における価格を求める必要がなければ，1 ステップ前に価格を求めた配列に上書きすればよいので，2 次元配列を用いる必要はない．時間変化のグラフを表示させる場合でも，1 ステップごとにファイルに出力すればよい．

t = 0 において，一つの原市場価格に対するオプション価格を求めたいだけなら，2 項ツリーの構造から，計算の手間を約半分に減らすことができる．例えば t1 においては Y1,−1と Y1,1 のみ求まればよく，Y1,0 は計算しなくてもよい．
C 言語の配列では負の添え字は取れないが，ポインタを用いてシフトしてやることにより，添え字を負にすることができる．具体的には，例えばdouble yy[1001];
double *y;y = yy + 500;などとすれば，y[-500] から y[500] が使える．
いつでも権利行使できるオプションを，アメリカンオプションという．ヨーロピアンオプ
ションの場合，時刻 t = tm におけるオプション価格は t = tm+1 における期待値を無リス
ク金利で割り引いたものになったが，もしその値がペイオフ関数より小さければ，即座に権利行使した方が得なので，ペイオフ関数を f (S) とすると()
−r∆tYm,n = max e(pYm+1,n+1 + qYm+1,n−1 ), f (Sn )となる．
問題 2.4 適当に状況を設定してアメリカンオプションの価格を計算し，時刻 t = 0 から満期まで価格がどのように変化していくかをグラフに描け．
問題 2.5 適当に状況を設定して，アメリカンオプションの時刻 t = 0 における一つの値が，どのように時間ステップ幅に依存して変化するか調べよ．時間ステップ幅を半分にしていきながら計算するとよい．

10

3

加速

数列 an がan = a + cλ−n + o(λ−n )なる挙動で収束していくとき，a′n =λan+1 − anλ−1なる数列を作ると，

a′n = a + o(λ−n )となり，収束が加速される．このような計算の効率化を Richardson 加速と呼ぶ．Richardson
加速は何度も繰り返し適用することができる．

Richardson 加速は λ がわからないと適用できないが，λ がわからなくても適用できる Aitken
加速という手法もある．
問題 3.1 ライプニッツ級数

π
1 1 1 1= 1 − + − + − ···4
3 5 7 9は極めて収束の遅い級数として知られている．
2
∑
(−1)k
n

an = 4k=0

2k + 1として a0 , · · · a16 を求め，Richardson 加速を適用して円周率 π を計算せよ．
問題 3.2 前章で求めた，ヨーロピアンオプションおよびアメリカンオプションの計算値に Richardson 加速を適用せよ．ただし，時間方向の分割を 2n 分割して計算した結果を anとする．λ は 2 や 4 などの整数値になる（つまり，分割を倍にしたら誤差は 2 分の 1 とか4 分の 1 とかになる）ので，計算結果からその数値を推定して用いよ．Richardson 加速を繰り返し適用する場合，アメリカンオプションでは，途中からあまり機能しなくなる．その理由も考えてみよ．

4Black-Scholes 方程式

前に用いたヨーロピアンオプションの価格を計算する式

Ym,n = e−r∆t (pYm+1,n+1 + qYm+1,n−1 )11について，時間刻みを細かくしていくとどうなるかを考えてみよう．

Ym+1,n − Ym,n = Ym+1,n − e−r∆t (pYm+1,n+1 + qYm+1,n−1 )
(
)
p−q
1
−r∆t
= Ym+1,n − e(Ym+1,n+1 + Ym+1,n−1 ) +(Ym+1,n+1 − Ym+1,n−1 )2
2= (1 − e−r∆t )Ym+1,n
(
)
1
p−q
−r∆t
−e
(Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1 ) +(Ym+1,n+1 − Ym+1,n−1 )22より，両辺を ∆t で割って

1 − e−r∆t
Ym+1,n − Ym,n=
Ym+1,n
∆t
∆t
(
)
1
p−q
−r∆t
−e
(Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1 ) +(Ym+1,n+1 − Ym+1,n−1 )2∆t2∆tとなる．ここで，s = log S と置き，数値計算スキームの収束性と解の適度な滑らかさを仮定すると，∆t → 0 のとき∂Y
Ym+1,n − Ym,n−→∆t∂t
1 − e−r∆t
Ym+1,n −→ rY∆t1
σ 2 Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1√(Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1 ) =·2∆t
2(σ ∆t)2
σ2S 2 ∂ 2Y
σ 2 S ∂Yσ2 ∂ 2Y
·
=
·
+
·
2 ∂s2
2
∂S 2
2
∂S
p−q
(p − q)σ Ym+1,n+1 − Ym+1,n−1√·(Ym+1,n+1 − Ym+1,n−1 ) = √2∆t∆t2(σ ∆t)(2er∆t − u − d)σ Ym+1,n+1 − Ym+1,n−1√√=·(u − d) ∆t2(σ ∆t)()()σ 2 ∂Y
σ2
∂Y
−→
r−
= r−
S·
2
∂s
2
∂S
−→

が成り立つので，2 項ツリーによる数値計算スキームは，偏微分方程式∂Y
σ2S 2 ∂ 2Y
∂Y=−·
− rS ·
+ rY
2
∂t
2
∂S
∂S
を離散化したものであると考えられる．この方程式は，Black-Scholes 方程式として知られている．Black-Scholes 方程式は，伊藤の補題からも導出することができるが，それについては後の章で述べる．

12

5

経路依存型オプション

5.1

最大値オプション

ペイオフが価格の動いてきた経路に依存するようなオプションを考える．まずは簡単のため，ペイオフが原資産の満期での値だけではなく，満期までの最大値にも依存するとしよう．このようなオプション価格を計算するにはどうすればよいだろうか？
√

前と同様に原市場の離散的な価格を Sn = S0 enσ ∆t と設定する．その上で，時刻 tm = m∆tにおいて，原市場価格が Sn1 ，それまでの最大値が Sn2 であるときのオプション価格をYm,n1 ,n2 とし，関係式を構築すれば計算が可能となる．
具体的には，n1 ≤ n2 に対して{
e−r∆t (pYm+1,n1 +1,n1 +1 + qYm+1,n1 −1,n1 )Ym,n1 ,n2 =
e−r∆t (pYm+1,n1 +1,n2 + qYm+1,n1 −1,n2 )(n1 = n2 )(n1 < n2 )を計算すればよい．

t = 0 において S−N から SN の範囲でオプション価格を求めたい場合には，時刻 t = M においては，原市場の価格刻みを S−(M +N ) から SM +N の範囲に，時刻 0 からの最大値の価格刻みを S−N から SM +N の範囲に取る必要がある．
問題 5.1 適当に状況を設定し，ペーオフが満期までの最大値に依存するようなヨーロピ
アンオプションおよびアメリカンオプションの価格を計算し，時刻 t = 0 における価格をグラフに描け．
プログラミングにおける注意：
関数内で普通に配列を宣言すると，関数ごとに用意されたスタック領域という領域にメモリが確保される．今回のプログラムでは 2 次元配列を用いることになるが，スタック領域には大きな容量の配列を確保することができないので，何らかの方法を考える必要があ
る．一つは，グローバル変数として宣言する方法である．つまり，通常ならint main(){double yy[2000][2000];....
}と変数を宣言するところを，13

double yy[2000][2000];
int main(){....
}と，関数の外で宣言すると，ヒープ領域という領域に大きな容量の配列を確保できる．もしくは，関数内でもint main(){static double yy[2000][2000];....
}というように static を付けるとヒープ領域に確保されるので，大きな容量の配列を使うことができる．
メモリの動的確保を用いるという方法もある．C 言語だと，

int m;
double *a;m = 1000;
a = (double*)malloc(m*sizeof(double));// ここで a[0] 〜 a[m-1] が使える
free(a); // 配列を使い終わったらメモリを解放するとして 1 次元配列を，

int i;
int m, n;double **a;m = 1000;
n = 1000;
a = (double**)malloc(m*sizeof(double*));a[0] = (double*)malloc(m*n*sizeof(double));for (i=1; i<m; i++){ a[i] = a[i-1] + n;}// ここで a[0][0] 〜 a[m-1][n-1] が使える
free(a[0]); // 配列を使い終わったらメモリを解放する
free(a);14とすると 2 次元配列を確保することができる．

C++の場合には，それぞれ，1 次元配列を確保するには
int m;
double *a;m = 1000;
a = new double[m];// ここで a[0] 〜 a[m-1] が使える
delete[] a; // 配列を使い終わったらメモリを解放するとし，2 次元配列を確保するには

int i;
int m, n;double **a;m = 1000;
n = 1000;
a = new double*[m];a[0] = new double[m*n];for (i=1; i<m; i++){ a[i] = a[i-1] + n;}// ここで a[0][0] 〜 a[m-1][n-1] が使える
delete[] a[0]; // 配列を使い終わったらメモリを解放する
delete[] a;とすればよい．ただし，C++の場合には，ファイルの名前は 〜.cpp とする必要がある．
メモリを動的確保した場合にはメモリの解放が必要である．プログラムを終了するまでに使ったメモリを解放しておかないと，プログラム終了後にも使用可能メモリが減ったままになり，ウィンドウズが不安定になることがある．一応，プログラムが正常に終了した場
合には自動的にメモリが解放されることになっているが，エラーで終了したような場合には解法されずに残ってしまうことがあるので，ちゃんと自分で解放するくせを付けておいた方がよい．

15

5.2

平均値オプション

経路依存型オプションとして，最大値に依存するオプションを考えたが，次に，ペイオフ
が原資産の満期での値および，満期までの平均値に依存するようなオプションを考えよう．この場合，最大値の場合と違って，平均値は原資産の価格刻みでは表現できないので，少し工夫が必要になってくる．
√

今までと同様に原市場の離散的な価格を Sn = S0 enσ ∆t と設定する．更に，時刻 0 からの平均値を表すための離散的な価格を用意する．こちらは細かい方が精度は良くなるが，ここでは簡単のため，原市場と同じ離散価格を用いることにしよう．その上で，時刻 tm = m∆tにおいて，原市場価格が Sn1 ，それまでの平均値が Sn2 であるときのオプション価格をYm,n1 ,n2 とし，関係式を構築すれば計算が可能となる．
時刻 tm = m∆t において，原市場価格が Sn1 ，それまでの平均値が Sn2 とすると，時刻
tm+1 = (m + 1)∆t には，原市場価格は，確率 p で Sn1 +1 ，確率 q で Sn1 −1 になる．またその(m + 1)Sn2 + Sn1 −1(m + 1)Sn2 + Sn1 +1，確率 q で B =とき，平均値は，確率 p で A =m+2
m+2になる．A や B は，現市場の離散価格では表されないので，そのままでは tm におけるオ
プション価格を tm+1 での価格からは計算することができない．そこで，tm+1 における，オプション価格のわかっている点から線形補間を持ちいて求めることにする．具体的には
以下のようになる．
){ (A − Sna
Sna +1 − A−r∆t
Ym+1,n1 +1,na +Ym+1,n1 +1,na +1
pYm,n1 ,n2 = e
Sna +1 − Sna
Sna +1 − Sna()}
Snb +1 − B
B − Snb+q
Ym+1,n1 −1,nb +Ym+1,n1 −1,nb +1,
Snb +1 − Snb
Snb +1 − Snb(m + 1)Sn2 + Sn1 +1A=,Sna ≤ A < Sna +1 ,m+2
(m + 1)Sn2 + Sn1 −1
B=,
Snb ≤ B < Snb +1 .
m+2

t = 0 において S−N から SN の範囲でオプション価格を求めたい場合には，時刻 t = M においては，原市場の価格刻みを S−(M +N ) から SM +N の範囲に，時刻 0 からの平均値の価
格刻みは S−M +N から SM +N の範囲に取れば十分である（平均値は原市場そのものほどは大きく変動しないので，実際はもっと減らすこともできるが）．
問題 5.2 適当に状況を設定し，ペイオフが満期までの平均値に依存するようなヨーロピ
アンオプションおよびアメリカンオプションの価格を計算し，時刻 t = 0 における価格をグラフに描け．

16

6

原市場のアセットが複数あるオプション，ノックアウトオ
プション，他

今まで扱ったオプションは原市場のアセットが一つだけであったが，原市場のアセットが複数あるようなオプションを考えることもできる．ここでは簡単のため，互いに独立なアセッ
ト U と V の満期での価格に依存するようなオプションを考える．U と√V のボラティリティー
√をそれぞれ σ1 , σ2 とし，原市場 U, V の離散的な価格を Un = U0 enσ1 ∆t , Vn = V0 enσ2 ∆t と設定する．また，時刻 tm = m∆t において原市場価格がそれぞれ Un1 , Vn2 であるときのオ
プション価格を Ym,n1 ,n2 とする．
後は，U の上昇確率 p1 ，下落確率 q1 ，V の上昇確率 p2 ，下落確率 q2 を√

u1 = eσ1 ∆t ,√u2 = eσ2 ∆t ,√d1 = e−σ1 ∆t ,√d2 = e−σ1 ∆t ,er∆t − d1,
u1 − d1er∆t − d2,p2 =u2 − d2p1 =u1 − er∆t
u1 − d1u2 − er∆tq2 =u2 − d2q1 =によって定め，(Ym,n1 ,n2 = e−r∆t p1 p2 Ym+1,n1 +1,n2 +1 + p1 q2 Ym+1,n1 +1,n2 −1)+ q1 p2 Ym+1,n1 −1,n2 +1 + q1 q2 Ym+1,n1 −1,n2 −1
によってツリーを計算すればよい．

t = 0 において U−N1 から UN1 かつ V−N2 から VN2 の範囲でオプション価格を求めたい場合には，時刻 t = M においては，原市場の価格刻みを U−(M +N1 ) から UM +N1 かつ V−(M +N2 )から VM +N2 の範囲に取ればよい．
問題 6.1 適当に状況を設定し，ペイオフが 2 種類のアセットに依存するようなヨーロピ
アンオプションおよびアメリカンオプションの価格を計算し，時刻 t = 0 における価格をグラフに描け．
その他のオプションについても，ツリー法を用いて同様の方法で計算することができる．
例えば，以下のようなオプションの計算なども簡単にできるだろう．
問題 6.2 適当なノックアウトオプションの価格を計算し，時刻 t = 0 における価格をグラフに描け．
問題 6.3 原市場の，満期での価格および時刻 t = T /2 での価格に依存してペイオフが決まるようなオプションを何か考え，時刻 t = 0 における価格をグラフに描け．
しかし，アセットの数が増えたり，経路に依存する量が増えると，計算量が爆発的に増えていくので，現実的な方法としては適当ではなくなる．また，価格方向の刻みをある程度小さくしようと考えると，時間方向の刻みを相当細かくしなければならないのもツリー法の欠点である．

17

7

差分法

ツリー法で価格方向の値を細かく見ようとすると，時間方向には価格刻みの 2 乗の細かさが必要となり，計算時間という点で問題が生じる．そこで，Black-Scholes 方程式

σ2S 2 ∂ 2Y
∂Y
∂Y
=−
·
− rS ·
+ rY
2
∂t
2
∂S
∂S
から数値計算を考えてみよう．価格は比で動くので，価格は指数的に取るのが自然である．そこで，s = log S とおいて，Y を s の関数として求めよう．

∂Y ds
1 ∂Y
∂Y
=
·
= ·
∂S
∂s ( dS
S) ∂s
2
∂ Y
∂
1 ∂Y
1 ∂Y
1 ∂ 2 Y ds
1 ∂Y
1 ∂ 2Y
=
·
=
−
·
+
·
·
=
−
·
+
·
∂S 2
∂S S ∂sS 2 ∂sS ∂s2 dSS 2 ∂s
S 2 ∂s2より，代入すると()∂Y
σ2 ∂ 2Yσ 2 ∂Y
=− ·
− r−
+ rY
∂t
2 ∂s2
2
∂s

が得られる．これを解くために，時間刻み ∆t，価格刻み ∆s を決め，Y (m∆t, n∆s) の近似を Ym,n とする．これから微分を差分で近似するわけだが，まずは陽解法と呼ばれる手法を示そう．差分法による近似として

Y ≈ Ym+1,n ,Ym+1,n − Ym,n
∂Y
≈
,
∂t
∆t
∂Y
Ym+1,n+1 − Ym+1,n−1
≈
,
∂s
2∆s
Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1
∂ 2Y
≈
,
2
∂s
∆s2
を導入し，Black-Scholes 方程式に代入すると

Ym+1,n − Ym,nσ 2 Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1
=− ·
2
∆t
2
(
) ∆s
2
σ
Ym+1,n+1 − Ym+1,n−1− r−+ rYm+1,n22∆sとなる．上の式を見ると，Ym,n が Ym+1,n−1 , Ym+1,n , Ym+1,n+1 で書けるので，満期から時間を逆に計算して行くことで時刻 0 での値を求めることができる．
問題 7.1 数値計算を行う上では，価格刻みの両端での境界条件が問題となるが，ヨーロ
ピアンオプションの場合およびアメリカンオプションの場合，それぞれどのように設定するのが適当か考えよ．

18

問題 7.2 適当に状況を設定し，陽解法を用いてヨーロピアンオプションおよびアメリカ
ンオプションの価格を計算し，時刻 t = 0 から満期まで価格がどのように変化していくかをグラフに描け．陽解法においては，時間刻みと価格刻みの関係に制限がある．どのようなときに計算が上手く行くかを調べよ．
問題 7.3 ∆t および ∆s の大きさがどの程度，誤差に影響するか予想せよ．また，実際に数値計算の結果から確かめよ．
陽解法においては，最後の rYm+1,n は rYm,n としてもよい（陽解法というのは，次の時間ステップの値が前の時間ステップの値で書き下せるという意味なので）．その場合に，√∆s = σ ∆t と取ると
(
(
√ (
√ (
))
))
∆t
σ2
1
∆t
σ2
1
Ym+1,n+1 +Ym+1,n+1(1 + r∆t)Ym,n =
+
r−
−
r−
2
2σ
2
2
2σ
2
となり，近似的にツリー法と一致する．つまり，ツリー法は陽解法の特別な場合であるといえる．
陽解法が安定に計算できるための条件（安定性条件）を考えてみよう．差分法の安定性を調べるには，厳密にはフーリエ級数や行列の固有値を考える必要があるのだが，数学的に少し難しいので，ここではもう少し大雑把な議論を行おう．数値計算において安定して計算を行うためには，誤差が増幅されないことが必要である．今回は，差分化した方程式が線形で重ね合わせの原理が成り立つので，誤差の種類ごとに分けて考えることができる．
実は差分法においては，交互に符号が変わるような誤差が最も速く増大するということが知られている．そこで，Ym+1,n に (−1)n ε の誤差が生じたときに Ym,n にどの程度の誤差が混入するかを調べてみよう．

(Ym+1,n + (−1)n ε) − (Ym,n + η)∆t
σ 2 (Ym+1,n+1 + (−1)n+1 ε) − 2(Ym+1,n + (−1)n ε) + (Ym+1,n−1 + (−1)n−1 ε)=− ·∆s2( 2 2)n+1
σ(Ym+1,n+1 + (−1) ε) − (Ym+1,n+1 + (−1)n−1 ε)− r−+ r(Ym+1,n + (−1)n ε)22∆sを η について解くと(
η=

2σ 2 ∆t
1 − r∆t −∆s2)となるので，誤差が増幅されないためには

2σ 2 ∆t< 2 − r∆t∆s219(−1)n εが成り立つ必要があるということがわかる．ここで，右辺の ∆t を無視すると，安定に計算できるためには，おおよそ
∆t
1
<
∆s2
σ2
が満たされる必要がある．これを陽解法の安定性条件という．最後の rYm+1,n を rYm,n とする陽解法の場合は)(
2σ 2 ∆t1(−1)n ε1−η=2
1 + r∆t∆sとなるので，安定性条件は

2σ 2 ∆t< 2 + r∆t∆s2となり，右辺の ∆t を無視すると先に述べた条件と同じになる．
安定性条件を満たすためには，時間刻みは空間刻みの 2 乗のオーダーにしなければならず，計算量という点では，ツリー法の難点は解消されていない．そこで登場するのが陰解法である。陰解法では，時刻 tm において差分化した関係式が成り立つように Ym,∗ を決める．具体的にはσ 2 Ym,n+1 − 2Ym,n + Ym,n−1
Ym+1,n − Ym,n
=− ·
2
∆t
2
(
)∆s
2
σ
Ym,n+1 − Ym,n−1− r−+ rYm,n22∆sと差分化することになる．証明には行列の固有値解析が必要なので詳しくは述べないが，陰解法は無条件安定，つまりどんな ∆t や ∆s に対しても安定であることが知られている．
問題 7.4 適当に状況を設定し，陰解法を用いてヨーロピアンオプションおよびアメリカ
ンオプションの価格を計算し，陰解法スキームが無条件安定であることを確かめよ．連立一次方程式の解法には，まずは通常のガウスの消去法を使用せよ．
今回の陰解法に現れる行列は 3 重対角行列になっており，ガウスの消去法を行う際にも，バンド幅 3 のバンド領域にしかアクセスしない．そこで，バンド領域のみを配列に持つことで，メモリを節約し，計算速度を上げることができる．
問題 7.5 陰解法スキームのプログラムを，バンド行列を用いたものに改良せよ．
今まで述べた解法やツリー法において，時間方向の誤差は O(∆t) となっている．これは，時間方向に片側差分を用いているからである．今まで用いた時間差分はYm+1,n − Ym,n∆tである．これは，t = tm や t = tm+1 の値を表現するには片側差分だが，t = (tm + tm+1 )/2の時点での値を表すと考えれば中心差分とみなせる．よって，t = (tm + tm+1 )/2 の時点で20

差分化を行ってみる．右辺は，t = tm での値と t = tm+1 での値の平均で t = (tm + tm+1 )/2での値を表現する．具体的に書くと{Ym+1,n − Ym,n1
σ 2 Ym,n+1 − 2Ym,n + Ym,n−1
=
−
·
∆t
2
2
∆s2
(
)
2
σ
Ym,n+1 − Ym,n−1− r−+ rYm,n22∆s
σ 2 Ym+1,n+1 − 2Ym+1,n + Ym+1,n−1
·
−
2
∆s2
}
(
)
σ 2 Ym+1,n+1 − Ym+1,n−1− r−+ rYm+1,n2
2∆s
このスキームは，クランク・ニコルソン法と呼ばれている．
問題 7.6 クランク・ニコルソン法の精度が，時間、空間、ともに 2 次精度であることを数値計算によって確かめよ．
ツリー法の極限を考えることで，経路依存型オプションも偏微分方程式に落とし込むことができる．ここではアベレージオプションについて考えてみよう．前に考えたときは価格刻みを整数に取るために線形補間を用いたが，価格刻みを整数と限らなければ}{
Ym,n1 ,n2 = e−r∆t pYm+1,n1 +1,na + qYm+1,n1 −1,nb ,(())(m + 1)Sn2 + Sn1 +1(m + 1)Sn2 + Sn1 −1
11
na = √ log, nb = √ log(m + 2)S0(m + 2)S0σ ∆tσ ∆tと書くことができる．ここで

1 − e−r∆t
Ym+1,n1 ,n2 − Ym,n1 ,n2=
Ym+1,n1 ,n2
∆t (
∆t
1
− e−r∆t(Ym+1,n1 +1,n2 − 2Ym+1,n1 ,n2 + Ym+1,n1 −1,n2 )
2∆t
)
p−q
(Ym+1,n1 +1,n2 − Ym+1,n1 −1,n2 )
+
2∆t
(
)
Ym+1,n1 +1,na − Ym+1,n1 +1,n2
Ym+1,n1 −1,nb − Ym+1,n1 −1,n2
−r∆t
−e
p·
+q·
,
∆t
∆t

となり，前半 3 行の極限は通常の Black-Scholes 方程式と同じなので，最後の 1 行についてのみ考える．U を時刻 t までの平均値を表す変数とし，u = log U とすると

Ym+1,n1 +1,na − Ym+1,n1 +1,n2 log Sna − log Sn2
Ym+1,n1 +1,na − Ym+1,n1 +1,n2·=∆tlog Sna − log Sn2
∆t

21

∂Y
Ym+1,n1 +1,na − Ym+1,n1 +1,n2−→log Sna − log Sn2
∂u
( (
)
)
log Sna − log Sn2
1(m + 1)Sn2 + Sn1 +1=log− log Sn2
∆t
∆t
m+2
(
)
1
Sn1 +1 − Sn2
S−U=
log 1 +−→∆t(m + 2)Sn2
tU
これより

Ym+1,n1 +1,na − Ym+1,n1 +1,n2 ∂Y S − U
·
∆t
∂u
tU

が成り立ち，同様にして

Ym+1,n1 −1,nb − Ym+1,n1 −1,n2 ∂Y S − U
·
∆t
∂u
tU
も成り立つ．これらを代入して，∂Y
σ2S 2 ∂ 2Y
S − U ∂Y∂Y
− rS ·
=−
·
−
·
+ rY
2
∂t
2
∂S
∂S
t
∂U
∂Y
にか
∂U
かっている分母の t の扱いが難しくなるので，V = tU なる変数 V を導入する．V は時刻 tまでの価格の積分値を表す．ここでオプション価格を，t, S, V を変数とする関数 Z(t, S, V )で表し，
Y (t, S, U ) = Z(t, S, tU )なる偏微分方程式が得られる．しかし，これを差分法で離散化しようとするとなる関係が成り立っているものとして，Z が満たす方程式を求める．

∂Z
∂Z
∂Y
=
+U ·
,
∂t
∂t
∂V

∂Y
∂Z
=t·
∂U
∂V

が成り立つので，Y に関する偏微分方程式に代入することでσ2S 2 ∂ 2Z
∂Z
∂Z
∂Z
=−
·
− rS ·
−S·
+ rZ
2
∂t
2
∂S
∂S
∂V
が得られる．今までと同様，s = log S とおいて価格を対数化すると()
σ2 ∂ 2Zσ 2 ∂Z∂Z∂Z=− · 2 − r−
− es ·
+ rY
∂t
2 ∂s
2 ∂s
∂V
が得られる．
問題 7.7 以前と同様の考え方により，アベレージオプションにおける陽解法の安定性条件を考察せよ．

22

問題 7.8 陽解法を用いてアベレージオプションの解を計算し，ツリー法の結果と比較せよ．
問題 7.9 陰解法およびクランク・ニコルソン法を用いてアベレージオプションの解を計算し，ツリー法の結果と比較せよ．連立一次方程式の解法はガウスの消去法を用いよ（そのため，空間刻みはそれほど大きくは取れない）．
陰解法やクランク・ニコルソン法を用いる場合には連立一次方程式を解く必要があるが，今回のように空間方向の独立変数が二つ以上ある場合には，空間刻みが大きくなると連立一次方程式の次数が急速に大きくなってしまうので，ガウスの消去法を用いるのは適当ではない．このような場合には，行列の疎行列性を生かした反復法が用いられる．これについては次章で詳しく述べる．
今回はアベレージオプションが満たす偏微分方程式をツリー法の極限として求めたが，最終的にはこれを離散化して解くわけで，ツリー法（離散）→偏微分方程式（連続）→差分化（離散）、というように離散と連続を行ったり来たりするのはあまりにも見通しが良くない．数理ファイナンスの理論では，このような経路依存型オプションが満たす偏微分方
程式を求める際には，確率微分方程式から伊藤の補題を用いて導出する．これについては後ほど解説する．

8

乱数

今後，モンテカルロ法などで乱数が必要になるので，まず作成しておく．まずは [0, 1] の一様乱数を作成する．rand() 関数を用いると，0 以上 RAND_MAX 以下の整数が返るので，( rand() + 0.5 ) / ( RAND_MAX + 1 )などとすれば [0, 1] の一様乱数が得られる．しかし，通常 RAND_MAX はそれほど大きな値ではないので，このままでは細かさがやや不足気味である．そこで，( rand() + ( rand() + 0.5 ) / ( RAND_MAX + 1 ) ) / ( RAND_MAX + 1 )とすれば十分に細かくなる．
問題 8.1 一様乱数を発生させ，[0, 1] を 100 分割して頻度分布を計算し，確率密度関数のグラフを表示せよ．
次に正規乱数を作成する．正規乱数の発生によく用いられるのが Box-Muller 法である．
U1 , U2 を，[0, 1] 上の一様分布に従うような互いに独立な確率変数とするとき，√
Z1 = −2 log U1 cos 2πU2√Z2 = −2 log U1 sin 2πU2で定義される Z1 , Z2 は，平均 0，分散 1 の標準正規分布に従う互いに独立な確率変数となる．

23

問題 8.2 Box-Muller 法により正規乱数を発生させ，(−∞, +∞) を適当に分割して頻度分布を計算し，確率密度関数のグラフを表示せよ．標準正規分布の確率密度関数も同時に表示せよ．
問題 8.3 Box-Muller 法が正しいことを証明せよ．

24