こんにちは、毛糸です。
ファイナンス（金融工学）において、正規分布は資産の収益率をモデル化するために頻繁に用いられます。
投資対象となる資産は通常1つだけではなく、複数資産を扱いたいことも多いですから、その場合には多次元の正規分布を考えなければなりません。
本記事では統計プログラミング言語Rで多次元正規分布に従う確率変数ベクトルを生成する方法について説明します。
ライブラリ「MASS」による多次元（多変量）正規分布の乱数発生
多次元（多変量）正規分布に従う確率変数ベクトルは、多変量解析用パッケージのMASSのmvrnorm()を使って発生させることが出来ます（公式リファレンスのPDFはこちら）。
mvrnorm()は、発生させる確率変数ベクトルの個数n、期待値ベクトルmu、分散共分散行列Sigmaを与え、n組の確率変数ベクトルを返す関数です。
例えば、( mu=(1,1))、分散共分散行列( Sigma=left( begin{array}{cc}1&0\0&1end{array}right))の2次元正規分布に従う確率変数( (x_1,x_2))を発生させるには、以下のように記述します。
#MASSライブラリを読み込む library(MASS) #期待値ベクトル mu0<-c(1,1) #分散共分散行列 Sigma0<-rbind( c(1,0), c(0,1) ) #多次元正規分布に従う確率変数ベクトルを1組発生 mvrnorm(1,mu0,Sigma0) #[1] 1.8590647 -0.6548381
参考>> 元データ分析の会社で働いていた人の四方山話_多変量正規分布
例題：年金資産の収益率
（出所：https://www.gpif.go.jp/gpif/portfolio.html）
[ begin{split}
Sigma=diag(S) cdot P cdot diag(S)
end{split} ]で表せます。ただし( diag(S))は( S)を対角成分に持つ対角行列で、「( cdot)」は通常の行列積（Rでは%*%）です。
#各資産クラスの期待リターン mu<-c(2.6/100, 6.0/100, 3.7/100, 6.4/100, 1.1/100) #各資産クラスの分散（標準偏差の２乗） sigma<-c(0.047,0.251,0.126,0.273,0.005) #相関行列 Rho<-rbind( c(1,-0.16,0.25,0.09,0.12), c(-0.16,1,0.04,0.64,-0.1), c(0.25,0.04,1,0.57,0.15), c(0.09,0.64,0.57,1,-0.14), c(0.12,-0.1,-0.15,-0.14,1)) #分散対角行列 sigma_diag<-diag(sigma) #分散共分散行列 Sigma<-sigma_diag%*%Rho%*%sigma_diag #多次元正規分布の発生 X <- mvrnorm(10000, mu, Sigma)
結果の確認
こうして得られた6次元確率変数ベクトルの1万個について、標本平均と標本標準偏差を計算すると、大数の法則により、パラメタとして与えたmuとsigmaに近くなるはずです。
Xは、行にサンプル数n、列に確率変数ベクトルの要素が並んでいます。列に対して平均meanと標準偏差sdを適用するには、apply(X,MARGIN=2,mean)という関数を使います。MARGIN=1はXの「列」方向に関数を適用するという意味です。
参考>>24. apply() ファミリー
#各要素の平均を計算（経済中位ケース） apply(X,2,mean)*100 #[1] 2.591798 5.739046 3.773379 6.390514 1.097378 2.793776
#各要素の標準偏差を計算 apply(X,2,sd)*100 #[1] 4.6795340 25.1717768 12.6962985 27.2644081 0.5005285 1.9114364
いずれも理論値に近い値になっています。