金融工学でのモンテカルロ法(05/23)：一般の分布に従う乱数(2) - 「大人の教養・知識・気付き」を伸ばすブログ

　今回から、金融工学におけるシミュレーションについて学んでいく。テキストとして以下を使う。今回はP.52-P.58まで。

モンテカルロ法の金融工学への応用 (シリーズ現代金融工学)

作者:祥二, 湯前,輝好, 鈴木
朝倉書店

Amazon

power-of-awareness.com

4. 一般の分布に従う乱数
- 4.4 正規分布
次回

4. 一般の分布に従う乱数

　モンテカルロ法は

問題に沿った（同時）分布に従う（多変量）乱数列の生成
その乱数列を使った計算

の2つの部分に分けて考えられる。乱数列の生成は更に

一様分布に従う乱数列 $\{u_1,u_2,\cdots,u_n\}$ の生成
それを元にした必要な同時分布に従う乱数列 $\{\boldsymbol{x}_1,\boldsymbol{x}_2,\cdots,\boldsymbol{x}_n\}$ の生成

に分けられる。

4.4 正規分布

　ここでは特に正規分布に従う乱数の生成方法を説明する。

4.4.1　逆関数法

　標準正規分布の分布関数の逆関数 $\Phi^{-1}(x)$ を用いると、一様分布に従う乱数列から、標準正規分布に従う乱数列が得られる。ただし、 $\Phi^{-1}(x)$ は解析的に求めることはできないため、Moroの方法などで近似的に算出する*1。

4.4.2　極座標法

　2変量の確率変数 $(X,Y)$ の同時分布関数が

$\begin{aligned} \Phi(x,y)=\int_{-\infty}^{x}\int_{-\infty}^{y}\displaystyle{\frac{1}{2\pi}}\exp\left(\displaystyle{-\frac{u^2+v^2}{2}}\right)du dv \end{aligned}$

と書けるとき、 $(X,Y)$ は2次元標準正規分布に従うという。2次元標準正規分布では、周辺分布 $F_X, F_Y$ はいずれも標準正規分布になる。また $X,\ Y$ は独立である。極座標、一様分布の逆関数法、指数関数の逆関数法を組み合わせると、正方形 $(0,1)$ に一様に分布する点列から、2次元標準正規分布に従う点列が得られる。
　2変量標準正規分布の密度関数 $f(r,\theta)$ の確率密度関数は

$\begin{aligned} f(r,\theta)=\displaystyle{\frac{1}{2\pi}}\exp\left(\displaystyle{-\frac{r^2}{2}}\right)=\displaystyle{\frac{1}{2\pi}}\exp\left(\displaystyle{-\frac{s}{2}}\right),\ s=r^2 \end{aligned}$

で書ける。確率変数 $X\sim U(0,2\pi)$ に従うとき、 $g(X)$ の期待値は

$\begin{aligned} E[g(X)]=\int_{0}^{2\pi}g(\theta)\displaystyle{\frac{1}{2\pi}}d\theta=\int_{0}^{1}g(2\pi\theta)d\theta \end{aligned}$

と書ける。また $X\sim Exp(1/2)$ であるとき、 $g(\sqrt{X})$ の期待値は、

$\begin{aligned} E[g(X)]=\int_{0}^{\infty}g(\sqrt{s})\displaystyle{\frac{1}{2}}\exp\left(-\displaystyle{\frac{s}{2}}\right)ds=\int_{0}^{1}g(\sqrt{-2\log{(1-s)}})ds \end{aligned}$

である。
　以上から、確率変数 $X(r,\theta)=X(\sqrt{s},\theta)$ が2次元標準正規分布に従うとき、 $g(X)$ の期待値は

$\begin{aligned} E[g(X)]&=\int_{0}^{2\pi}\int_{0}^{\infty}g(r,\theta)\displaystyle{\frac{1}{2\pi}}\exp\left(-\displaystyle{\frac{s}{2}}\right)rdrd\theta\\ &=\int_{0}^{\infty}\left[\int_{0}^{1}g(\sqrt{s},2\pi\theta)d\theta \right]\displaystyle{\frac{1}{2}}\exp\left(-\displaystyle{\frac{s}{2}}\right)ds\\ &=\int_{0}^{1}\int_{0}^{1}g(\sqrt{-2\log{(1-s)}},2\pi\theta)d\theta ds \end{aligned}$

で計算できる。
　これを活用して2次元標準正規分布に従う乱数を得るにはたとえばBox-Muller法がある*2。
　またBox-Muller-Marsaglia法もある。
　まず単位円の中に一様に分布する点列をつくり、指数分布の逆関数法を用いて2次元標準正規分布に変換する方法である。

(1)	正方形 $(0,1)^2$ に一様に分布する点列 $(u_n^{(1)},u_n^{(2)})$ を生成する。
(2)	点列 $(u_n^{(1)},u_n^{(2)})$ から逆関数法によって正方形 $(-1,1)^2$ の中に一様に分布する点列 $(x_n^{(1)},x_n^{(2)})$ を生成する。
(3)	点列 $(x_n^{(1)},x_n^{(2)})$ から不要な点を取り除いて単位円の内側に一様に分布する点列 $(y_n^{(1)},y_n^{(2)})$ を生成する。点 $(y_n^{(1)},y_n^{(2)})$ を極座標 $(r_n,\theta_n)$ $\begin{aligned}r_n&=\sqrt{s_n},\ s_n=(y_n^{(1)})^2+(y_n^{(2)})^2,\\\theta_n&=\tan^{-1}\left(\displaystyle{\frac{y_n^{(1)}}{y_n^{(2)}}}\right)\end{aligned}$ で表示すると、 $\{r_n\}$ は動径方向 $(0,1)$ で一様に分布し $\{\theta_n\}$ は角方向で $(0,2\pi)$ で一様に分布する。
(4)	動径方向 $r_n$ に指数分布の逆関数法を用いると、2次元標準正規分布に従う点列 $\{t_n,\theta_n\}$ $\begin{aligned}t_n=\sqrt{-2\log{(1-s_n)}}\end{aligned}$ が得られる。　 $(t_n,\theta_n)$ を再びデカルト座標系 $(z^{(1)},z^{(2)})$ で表示すると $\begin{aligned}z_n^{(1)}&=y_n^{(1)}\sqrt{\displaystyle{\frac{-2\log{(1-s_n)}}{s_n}}},\\z_n^{(2)}&=y_n^{(2)}\sqrt{\displaystyle{\frac{-2\log{(1-s_n)}}{s_n}}}\end{aligned}$ 点列 $\{(z_n^{(1)},z_n^{(2)})\}$ は2次元正規分布に従う。

4.4.3　一般の多次元正規分布

　 $k$ 次元確率ベクトル $\boldsymbol{X}=(X^{(1)},\cdots,X^{(k)})$ が

$\begin{aligned} f(\boldsymbol{x})=\displaystyle{\frac{1}{(2\pi)^{n/2}\sqrt{|\Sigma|}}}\exp{\left(-\displaystyle{\frac{^{t}(\boldsymbol{x}-\boldsymbol{\mu})\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu})}{2}}\right)} \end{aligned}$

にしたがうとき、 $\boldsymbol{X}\sim N_k(\boldsymbol{\mu},\Sigma)$ である。 $\boldsymbol{\mu}=\boldsymbol{0},\ \Sigma=I$ であるとき、 $N_k(\boldsymbol{0},\Sigma)$ を $k$ 次元標準正規分布という。
　 $X\sim N_k(\boldsymbol{0},I)$ で、 $Y\sim N_k(\boldsymbol{0},\Sigma)$ とする。さらに正則行列 $C$ を用いて $Y=CX$ と表されたとする。このとき、 $\Sigma=C ^{t}C$ である。この $C$ は $Y$ をCholesky分解することで得られる。

　 $k$ 次元分散共分散行列を $\Sigma=(\sigma_{ij})$ 、下三角行列 $C=(c_{ij})$ の各要素を
$\begin{aligned} c_{11}&=\sqrt{\sigma_{11}},\\ c_{i1}&=\displaystyle{\frac{\sigma_{i1}}{\sqrt{\sigma_{11}}}},\ j=2,3,\cdots,k,\\ c_{jj}&=\displaystyle{\sqrt{\sigma_{jj}-\sum_{l=1}^{j-1}c_{jl}^2}},\ j=2,3,\cdots,k,\\ c_{ij}&=\displaystyle{\frac{1}{c_{jj}}\left(\sigma_{ij}-\sum_{l=1}^{j-1}c_{il}c_{jl}\right)},\ j\lt i,\ i=2,3,\cdots,k-1 \end{aligned}$
とすれば、
$\begin{aligned} \Sigma=C ^{t}C \end{aligned}$
となる。

　逆関数法およびCholesky分解を順に用いると $k$ 次元超立方体 $(0,1)^k$ に一様に分布する点列 $U$ から、多次元標準正規分布に従う確率変数 $\boldsymbol{X}$ を介して、一般の多次元正規分布に従う確率変数 $\boldsymbol{Y}$ が得られる：

(1)	$k$ 次元超立方体 $(0,1)^k$ に一様に分布する点列 $\{\boldsymbol{u}_n\}=\{(u_n^{(1)},\cdots,u_n^{(k)})\}$ を生成する。
(2)	$\{\boldsymbol{u}_n\}$ から逆関数法又は極座標法により $k$ 次元標準正規分布に従う点列 $\{\boldsymbol{x}_n\}=\{(x_n^{(1)},\cdots,x_n^{(k)})\}$ を生成する。
(3)	欲しい $k$ 次元正規分布の分散共分散行列 $\Sigma$ をCholesky分解して $C$ を求める。
(4)	点列 $\{\boldsymbol{y}_n\}=\{(y_n^{(1)},\cdots,y_n^{(k)})\}$ を $\boldsymbol{y}_n=C\boldsymbol{x}_n$ で得る。