上級編 ６． Libor Market Model
6.6 モンテカルロシミュレーション
6.6.3 乱数の生成方法
クオンツモデルで使われるSDEの大半は、拡散項がブラウン運動によって駆動されると仮定しています。ブラウン運動は、初期値０からスタートすると、Ｔ時間後の平均が０で分散がＴとなる正規分布をする事が知られています。分散をＴで割って正規化すれば、標準正規分布になります。モンテカルロシミュレーションでは、このブラウン運動を、コンピューターアルゴリズムで生成した標準正規乱数で表現します。「はじめに」のセクションで述べたように、ここでの課題は
- 質の高い乱数を
- できるだけ高速で生成する
という事になります。
6.6.3.1 一様乱数の生成
ところで、コンピューターを使った乱数生成アルゴリズムでは、標準正規乱数を直接生成する事はできません。まず、0 から 1 の間で一様分布する乱数を生成し、それを別のアルゴリズムを使って標準正規乱数に変換します。変換のアルゴリズムは、いくつか公表されていますが、それについては後ほど紹介します。ここではまず一様乱数の生成方法について簡単に解説します。
一様乱数とは、一定の区間 \([a,b)\) 間で、ある数 \(r,~~(a \leq r \lt b) \) が発生する確率（すなわち確率密度）がすべて同じになるような確率分布をする乱数です。コンピューターで生成される一様乱数は、\([a,b)~を~[0,1)\) としたものです。
簡単な乱数発生アルゴリズムで、よく使われるものに線形合同法(Linear Congruential Generators)と呼ばれているものがあります。簡単にアルゴリズムを説明すると、
- まず何等かの（整数の）シード数 \(S_0\) を与えます
- そのシード数を、まず 1次式で線形変換し、それを何等かの整数 ｍで割った時の余りを求め、それを \(s_1\) とします。すなわち \(S_1=(aS_0+b)~~mod(m)\) (\( ~(x)mod(m)\) は x を m で割った時の余り）。 a,b,m は整数
- 2.で得られた \(S_1\) を ｍで割った値を \(U_1\) とすると、\(0\leq U_1 \lt 1\) となります。
- 2. 3. と同じ操作を、今度は \(S_1\) をシード数として使って行い、生成された数を \(S_2,~U_2\) とします。
以上の操作を繰り返していけば、\(U_1,~U_2,…\) と一様乱数列が生成できます。
この時、定数 a,b,m の選択に注意する必要があります。というのは、上記のアルゴリズムでは、生成された乱数 \(S_i\) が、たまたま当初のシード数 \(S_0\) と同じになった場合、その後に生成される乱数列は、全く同じ数列が繰り返される事になります。生成される乱数の種類は ｍ－１ 個なので、上記のアルゴリズムを ｍ－１ 回繰り返せば、必ずどこかで当初のシード数に戻ります。これを乱数生成アルゴリズムの周期性（Cycle）と言いますが、望むべくは、同じ数に戻るのは、ｍ－１ 個のすべての数が生成された後がいいでしょう。これを Full Period と呼びます。うまく a,b,m を設定すれば、そのような Full Period の乱数列が生成出来る事が証明されています。具体的には、c を 0 とした上で ｍ を素数とし、a をその原始根（Primitive root）とすれば、そうなりますが、この辺りの数学的な説明は、私の専門外であり、飛ばさせていただきます。（原始根の簡単な説明は、原始根の定義と具体例（高校生向け） | 高校数学の美しい物語 (manabitimes.jp)などが、簡明かつ解りやすい解説をしています）
m と a の具体例として、下記のような例が先人達によって示されています。
\[ \begin{align} & m=2^{31}-1,~~~~ a=16807 \\ & m=2^{31}-1,~~~~ a=39373 \\ & m=2^{31}-1,~~~~ a=742938285 \\ & m=2^{31}-1,~~~~ a=950706376 \\ & m=2147483399,~~~~ a=40692 \\ & m=2147483563,~~~~ a=40014 \end{align} \]これらの a,m を選択すれば、周期は \( 2^{31}-1 \) になります。かなりの数の乱数が、周期を繰り返す前に生成できます。ちなみに、\(2^{31}-1 \) はメルセンヌ素数と呼ばれる素数で、コンピューターの 32 ビットで表現できる最大の素数になります。
さて、上記のような線形合同法は、約 21 億個の乱数が、ランダムに登場し、一見よさそうに見えます。アルゴリズムで Full Period 分だけ乱数を生成し、それらを ｍ で割れば、\([0,1)\) 間の一様乱数がうまく生成できます。しかし、よくよく分析してみると、分布に偏りが出来る事が明らかになっています。というのは、この乱数を、生成された順番に 2 個とり、それを 2 次元空間上の点とみなして、平面にプロットすると、点に偏りが見られます。より高次の空間の点でプロットすれば、よりはっきりします。LMM のように高次元のモデルでは、これは大きな問題です。
幸いな事に、先ほど述べた周期の問題と、高次の空間での乱数の偏りの問題を解決したアルゴリズムがいくつか発表されており、その中でも有名なのは、メルセンヌ・ツイスタという乱数生成アルゴリズムがあります。このアルゴリズムは、周期が \( 2^{19937}-1\) (いくつかあるメルセンヌ・ツイスタのアルゴリズムの中で MT19937 と呼ばれているもの)と、非常に長周期で、かつ高次元の空間へプロットしても、偏りがほとんど見られない事が確認されています。物理の世界ではどうか判りませんが、金融工学の分野では、これで十分です。また、C 言語や C++ 言語で書かれたそのアルゴリズムは、そこそこ高速です。
一様乱数の生成アルゴリズムについては、これ以上詳細な説明は、能力外であり、省略します。実務家にとっては、数学の先人達が考え出したこれらのアルゴリズムを、有難く使わせて頂く事でいいでしょう。実践編で解説している、クオンツモデルのオープンソースである QuantLib でも、モンテカルロシミュレーションには主にメルセンヌ・ツイスタを使っており、そのソースコードは、そこで確認できます。