上級編 ６． Libor Market Model
6.6 モンテカルロシミュレーション
6.6.3 乱数の生成
6.6.3.2 一様乱数を標準正規乱数に変換
次に、生成された一様乱数を、標準正規乱数に変換するアルゴリズムにかけます。そのアルゴリズムは、いくつか提示されており、どれを使うかはケースバイケースですが、金融実務でよく使われる方法は、Box-Muller 法と、累積分布関数の逆関数を使う方法です。
< Box-Muller 法 >
Box-Muller 法は、２つの独立した一様乱数から、簡単な数式で、２つの独立した標準正規乱数を生成する方法です。具体的には、２つの一様乱数 \(u_1,u_2\) が与えられたとすると、下記式で、それらを２つの標準正規乱数 \(x_1,~ x_2\) に変換できます。
\[ \begin{align} & x_1=\sqrt {-2 \ln {u_1}} \cos {2πu_2} \\ & x_2=\sqrt {-2 \ln {u_1}} \sin{2πu_2} \tag{6.81} \end{align} \]この式の導出プロセスは、Wikipedia などネットで多く出回っているので、それを参照して下さい。簡単に、その考え方だけを述べておきます。
- 指数分布する乱数は、一様乱数から解析的に変換できる。
- 指数分布する乱数は、標準正規乱数からも解析的に変換できる。
- この２つの変換の解析式を使って、標準正規乱数を未知数とした連立方程式を組み立てる。
- それを解けば上式が導出できます。
6.81 式は、一見すると簡単な式に見えますが、式に対数関数と三角関数が含まれています。これらの計算は、単純な四則演算の組み合わせに比べると、CPU内部での計算に比較的時間がかかる事が知られています。ナノ秒の単位の話ですが、塵も積もれば山となるで、実際に計算してみて、他の方法と計算時間を比べながら使う事になるでしょう。
< 標準正規分布関数の逆関数を使う方法 >
標準正規分布関数や、その逆関数が解析的に求まるのであれば、上記の Box-Muller 法のように遠回りする必要はなく、一様分布する乱数を、直接標準正規分布する乱数に変換できます。しかし、標準正規分布については、確率密度関数は解析的に求まっているものの、累積分布関数や、その逆関数は解析的に求まりません。では、どうやってこの方法を使うのかという事ですが、標準正規分布の分布関数やその逆関数については、極めて誤差の小さい近似式があり、それを使います。
標準正規分布する乱数 X が、ある点 x より小さくなる確率、すなわち累積分布関数を下記のように表記します。
\[ P(X \leq x)=Φ(x) ~~~~-∞ \lt x \lt ∞ \]この逆関数を次のように表記します。
\[ Φ^{-1}(u)=x ~~~~ 0 \lt u \lt 1 \]この式に、累積分布確率 u（例えば 0.5）を与えてやれば、それに対応する確率変数 ｘ（例えば 0）を返します。これを使い、一様分布する乱数 u をこの式に代入すると、それに対応する標準正規乱数 ｘ を返します。
ではその近似式ですが、下記の式が良く使われています。ちなみに、\(\Phi^{-1}(u)\) は定義域が \( 0\leq u \lt 0.5~と~ 0.5\leq u \lt 1 \) で対称なので、どちらかの範囲で定義できれば、反対側の領域は、入力値を １－u と変換する事で、同じ関数形を使えます。ちなみに、下記式は、定義域が \(0.5\leq u \lt 1\) の場合を、２つの領域に分けて定義されています。
< \( 0.5 \leq u \leq 0.92 \) の領域 >
\[ x=Φ^{-1}(u) ≈ \frac {\sum_{i=0}^3 a_i \left(u-1/2\right)^{2i+1}}{1+\sum_{i=0}^3 b_i \left(u-1/2\right)^{2i}},~~~~~0.5\leq u \leq 0.92 \tag{6.82} \] \[ \begin{align} where~~~~~& a_0=2.50662823884, \\ & a_1=-18.61500062529, \\ & a_2=41.39119773534,\\ & a_3=-25.44106049637,\\ & b_0=-8.47351093090,\\ & b_1=23.08336743743,\\ & b_2= -21.06224101826,\\ & b_3=3.13082909833 \\ \end{align} \]< \( 0.92 \leq u \lt 1 \) の領域 >
\[ x=Φ^{-1}(u) ≈ \sum_{i=0}^8 c_i \left[ \ln (-\ln (1-u) ) \right]^i,~~~~~ 0.92 \leq u\lt 1 \tag{6.83} \] \[ \begin{align} where~~~~~& c_0=0.3374754822726147,\\ & c_1=0.97616901917186,\\ & c_2=0.1607979714918209,\\ & c_3=0.0276438810333863,\\ & c_4=0.0038405719373609,\\ & c_5=0.0003951896511919,\\ & c_6=0.0000321767881768,\\ & c_7=0.0000002888167364,\\ & c_8=0.0000003960315187 \end{align} \]< \( 0 \leq u \lt 0.5 \) の領域 >
\(\Phi^{-1}(u)\) は 0.5 で対象なので、u を 1-u と置いて、\(-\Phi^{-1}(1-u)\) を計算すればよい。
この近似式は、極めて正確で、生成される x が ±7 標準偏差の範囲内で、最大誤差が \(3×10^{-9}\) のオーダーになる事が知られています。金融工学の分野で使うには、問題ないでしょう。この方法が、Box-Muller 法より劣っている訳ではありません。Box-Muller 法は、演算の中に指数関数と三角関数が含まれており、計算時間がややかかります。また Box-Muller 法は、解析的な式ではあるものの、コンピューター上で扱える数字は、桁数の制約があり、所詮どこかで切り捨てられるので、その差は近似のレベルの違いだけかと思います。
この他にも、一様乱数から標準正規乱数に変換するアルゴリズムはいくつか紹介されています。以下に、それらを簡単にリストアップしますが、詳細はWiki等で確認して下さい。
- Marsagria Polar 法 Zuggratアルゴリズム （いずれも Box-Muller 法のレベルアップ版）
- 受理棄却法（Acceptance-Rejection Method）
- 中心極限定理を使った方法：一様乱数を N 個発生させ、その合計から N/2 を引いた数は、標準正規分布に分布収束する。
この内、受理棄却法と、中心極限定理を使った方法は、ひとつの正規乱数を生成するのに、多数の一様乱数を必要とし、計算効率はいまひとつです。
< 互いに相関のある複数の標準正規乱数の生成 >
LMM で相関行列を分解した行列 \(\bf C（ρ=CC^{\top}）\) の導出方法を説明しました。この C を互いに独立した標準正規乱数のベクトルに作用させれば、生成される乱数のベクトルは、互いに ρ の相関を持ちます。C の導出方法については、解説済なので、そちらをご覧ください。