平成26年度 ミクロ計量経済学
講義ノート9 マルコフ連鎖モンテカルロ法(MCMC)
このノートでは、マルコフ連鎖モンテカルロ法(Markov chain Monte Carlo, MCMC)の 解説を行う。MCMCは、乱数を発生させる技術の一つであり、直接的に乱数を発生させるこ とが難しい場合に使用される。主な使用法は、ベイズ統計学において、事後分布からMCMC によって乱数を発生させて、その乱数の統計量を調べることで事後分布の性質をみることで ある。さらに、近年では、頻度論に依拠する統計分析においても、MCMCは有用であるこ とが示されている。
まずマルコフ連鎖に関する基本的な知識を概観し、MCMCの代表的な手法であるGibbs samplerとMetropolis and Hastings (MH)法を概観する。そして、頻度論でのMCMCの 有用性を示したChernozhukov and Hong (2003)の手法を紹介する。なお、このノートの最 初の3節はCameron and Trivedi (2005)の13章を参考に書かれている。また、MCMC法 のより包括的な議論と計量経済学におけるベイズ統計への応用については、Chib (2001)を 参照のこと。
9.1 マルコフ連鎖
この節ではマルコフ連鎖に関する基本的な内容を概観する。簡単化のために、離散確率変 数の場合のみを考える。連続確率変数でも同様な議論が成り立つ。重要な結果はマルコフ連 鎖の収束に関する部分であり、マルコフ連鎖が定常分布に収束するという結果がMCMCの 理論的根拠となっている。
マルコフ連鎖とは、次の性質を満たす確率変数列xn, n = 0, 1, . . . ,のことである。xnは 離散変数であり、その取り得る値は1, . . . , mであるとする。xnがマルコフ連鎖であるとは、
Pr(xn+1 = x|xn, xn−1, . . . , x0) = Pr(xn+1= x|xn), ∀x ∈ {1, . . . , m} (1) のように、過去の値に条件づけた確率が、一期前の値にのみ依存する性質を持つということ である。なお、ここでは、Pr(xn+1= y|xn= x)の値はnによらないと考える。Pr(xn+1 = y|xn= x)の値を定める関数K(x, y)をtransition kernelという。
マルコフ連鎖のの議論のためには、行列表記を導入すると便利である。確率txyを
txy = Pr(xn+1 = y|xn= x) = K(x, y) (2) と定義する。行列T を
T =
t11 . . . t1m
... . .. ... tm1 . . . tmm
(3)
と定義する。このように定義されたT を推移行列(transition matrix)と呼ぶ。
推移行列を用いることで、マルコフ連鎖の各時点での分布の移り変わりを、簡便に表現す ることができる。xnのn時点での(周辺)分布をq(n)とする。つまり、q(n)はm × 1のベク トルであり、そのi番目の要素はPr(xn= i)である。すると、(n + 1)時点での分布q(n + 1) は
q(n + 1) = q(n)T (4)
と書ける。このことから、初期分布q(0)がわかると、
q(n) = q(0)Tn (5)
と、その後の分布がすべて、T を用いて表現することができる。
マルコフ連鎖はn → ∞のときに、ある条件のもとで、常に定常分布に収束することが証 明できる。マルコフ連鎖の定常分布とは、
q∗= q∗T (6)
となるような分布q∗のことである。つまり、ある条件の下で、
n→∞lim qT
n= q∗ (7)
が、任意のqについて成り立つ。この収束のための条件は、n時点でどの値をとっていても
(n + 1)時点で任意の値を取り得るということである。
MCMCはマルコフ連鎖の収束に関する性質を利用して、乱数を発生させる方法である。 つまり、分布q∗から乱数を発生させたい場合には、q∗を定常分布として持つマルコフ連鎖 を構築する。そのマルコフ連鎖の条件付き分布が乱数を容易に発生させうる者であれば、適 当な初期分布からマルコフ連鎖を長期間発生させることで、定常分布q∗からの乱数を得る ことができる。
9.2 Gibbs sampler
この節ではGibbs samplerを紹介する。この方法は、同時分布から乱数を発生させること は難しいが、条件付き分布からの乱数発生なら簡単である場合に有用な方法である。
まず、発生させたい乱数が二つの組に分けられる場合のGibbs samplerから紹介する。そ の分布がP (x, y)であるとする。条件付き分布のP (x|y)とP (y|x)からの乱数発生は容易で あるとする。Gibbs samplerのアルゴリズムは次の通りである。
1. n = 1とし、x(0)を適当に定める。
2. y(n)をP (y|x = x(n − 1))という条件付き分布から発生させた乱数の実現値とする。 3. x(n)をP (x|y = y(n))という条件付き分布から発生させた乱数の実現値とする。 4. nの値をひとつ増やす
5. 2-4をN 回繰り返す。
こうして得られる(x(n), y(n))の組はnが十分に大きい場合には、P (x, y)から発生させた 乱数の実現値とみなしてもよいことが、マルコフ連鎖の理論からわかる。
Gibbs samplerによって{(x(n), y(n))}Nn=1というN組の乱数の実現値が得られるわけであ
るが、P (x, y)からの実現値と見なせるのは、nが十分に大きい場合だけであるので、最初の
方の発生させた値は、使用することができない。たとえば、N = 10000とすると、n = 5000 までの数値は捨てて、{(x(n), y(n))}10000n=5001をP (x, y)からの抽出として使用する。このよう に、最初の捨てる部分をburn-in phaseと呼ぶ。どの程度の最初の発生値を捨てるべきかは、 マルコフ連鎖がどれほど早く収束するかによる。収束判定の方法はいくつか提唱されている ので、関連文献を調べるとよい。
なお、Gibbs samplerに限らずMCMCによって乱数を発生させる場合には、通常のi.i.d. のモンテカルロに比べて、多くの回数を必要とする。MCMCによって発生させる乱数は系 列相関を持っている。つまり、(x(n), y(n))と(x(n + 1), y(n + 1))は相関している。この系 列相関のために、たとえば、xの平均を求めようとして{x(n)}Nn=1の平均をとったとしても、 その精度は低くなるおそれがある。そのため、通常のi.i.d.のモンテカルロよりも多くの回 数が必要になるのである。
複数の組の条件付き分布を使用するGibbs samplerもある。たとえば、(x1, x2, . . . , xK)と いうK組の乱数をP (x1, x2, . . . , xK)という分布から発生させたいとする。その際、条件付 きP (xk|x−k)からの発生は可能であるとする。ここでx−kはxkを除いた乱数の組である。 この場合のGibbs samplerのアルゴリズムは、
1. n = 1とし、x2(0), . . . , xK(0)を適当に定める。
2. x1(n)をP (x1|x2(n − 1), . . . , xK(n − 1))から発生させる。 3. x2(n)をP (x2|x1(n), x3(n − 1), . . . xK(n − 1))から発生させる。
4. 以下同様に、xK(n)をP (xK|x1(n), . . . , xn−1(n))から発生させるところまで続ける。 5. nの値をひとつ増やす。
6. 2-5をN 回繰り返す。 となる。
9.3 MH法
つぎにMetropolic-Hastings法(MH法)を紹介する。Gibbs samplerを使うためには、条 件付き分布からの乱数発生が可能である必要がある。MH法は、同時分布や条件付き分布か らの発生がなくとも、乱数生成を可能にする方法である。ただし、MHを使用するためには、 確率あるいは密度の値が計算できる必要がある。
まず表記を定義する。密度p(x)をもつ分布から乱数を発生させたいとする。この密度の 値は計算可能であるとする。MHを使用する際には、提案密度(proposal density, jumping distribution)と呼ばれる密度関数を決める必要がある。これは、Q(x′|x)と表記する。これ は、xが与えられたときの、条件付き密度であり、Q(x′|x)からの乱数は発生可能であると する。Q(x′|x)は勝手に決めることができるので、Q(x′|x)から発生が可能であるという条 件は、MH法の使用可能な状況に制限を加えない。ただし、Q(x′|x)とp(x′)の値が近い方 がMH法はうまく機能することがわかっている。
MH法のアルゴリズムは以下の通りである。 1. n = 1と置く。x(0)を適当に定める。 2. x∗をQ(x|x(n − 1))から発生させる。 3. 次の比
r = p(x
∗)Q(x(n − 1)|x∗)
p(x(n − 1))Q(x∗|x(n − 1)) (8)
を計算する。
4. x(n)を
x(n) =
x∗ 確率min(r, 1)
x(n − 1) 確率(1 − min(r, 1)) (9)
として、決める。 5. 2-4をN 回繰り返す。
この様にして得られる{x(n)}Nn=1はnが十分に大きい場合にはp(x)からの抽出と見なすこ とができる。
MH法の理論的根拠を簡単に述べる。xaとxbの二つの点をとる。今、p(xa)Q(xb|xa) > p(xb)Q(xa|xb)とする。もしx(n − 1) = xaかつx∗= xbであれば、r = 1となり、x(n) = xb となる。したがって、x(n − 1)がp(x)から抽出されている場合の、(x(n), x(n − 1))の同時 密度の(xb, xa)での値は、f (xb, xa) = Q(xb|xa)p(xa)となる。次に、f (xa, xb)の値を計算 する。もしx(n − 1) = xbかつx∗ = xaであれば、r = p(xa)Q(xb|xa)/(p(xb)Q(xa|xb))と なる。したがって、x(n − 1)が、p(x)から抽出されている場合の(x(n), x(n − 1))の同時 密度の(xa, xb)での値は、f (xa, xb) = rQ(xa|xb)p(xb) = Q(xb|xa)p(xa)となる。これから、 f (xa, xb) = f (xb, xa)となる。そのため、x(n)の周辺密度もp(x)となる。よって、x(n − 1) の周辺密度がp(x)であれば、x(n)の周辺密度もp(x)となることがわかる。つまり、p(x) は、MH法で構築するマルコフ連鎖の定常分布である。
Gibbs samplerの時と同様に、nが大きくないとp(x)から生成した乱数とはいえないの で、nの小さい部分を捨ててしまう「burn-in」を行う必要がある。また、アルゴリズムをみ ればわかるように、MH法によって生成される乱数列は系列相関が強いために、通常のモン テカルロよりも多くの回数、乱数発生を行う必要がある。
9.4 Chernozhukov and Hong (2003)の方法
MCMC法は、ベイズ統計で広く使用されているが、頻度論の基づく統計分析でも有用で ある。ここでは、Chernozhukov and Hong (2003)によって提唱された準ベイズ法(ラプラ ス型推定量とも呼ばれる)を紹介する。極値推定量の計算をMCMC法を使って行うもので あり、推定量の計算のための最適化が難しい場合でも、準ベイズ法は実行可能(ただし時間 はかかる)ため有用である。
推定すべき母数θ0は、ある関数M (θ)を最小化する者とする。つまり、 θ0= arg min
θ M (θ) (10)
である。そして、M (θ)に収束するデータから計算できる関数Ln(θ)があるとする(Ln(θ)は Ln(θ)/n → M(θ)を満たすとする)。θ0の推定量として、次の極値推定量
θˆn= arg min
θ Ln(θ) (11)
を考える。しかし、問題によってはこの最小化問題を解くことが非常に難しい。準ベイズ法 は、そのような場合に非常に有用である。
例としては、GMM推定があげられる。真値θ0は次のモーメント条件
E(g(Xi, θ0)) = 0 (12)
を満たすとする。GMM推定量は、
θˆn= arg min
θ n
(1 n
n
∑
i=1
g(Xiθ) )′
W (1
n
n
∑
i=1
g(Xiθ) )
(13)
として、計算する。ただし、W は重み付け行列である。したがってGMM推定量は、ここ で考える問題の設定に当てはまる。モーメント条件の形状によっては、GMM推定量の計算 が難しくなる場合がある。例としては、操作変数を使用した分位点回帰などがあげられる。
準ベイズ法による推定量は、Ln(θ)をもとにした母数の疑似事後分布から乱数を抽出する ことで求められる。θの疑似事後分布とは、π(θ)を適当に定めた事前分布として、
pn(θ) = exp(−Ln(θ))π(θ)
∫ exp(−Ln(θ))π(θ)dθ (14)
である。pn(θ)から、MCMCを用いてθを抽出する。(θ1, . . . , θB)をそのようにして抽出さ れたB個のθの値とする。(θ1, . . . , θB)の平均、あるいは中央値をθ0の推定量とする。これ が準ベイズ法である。なおMH法を使用する場合は、分母の値を計算する必要はない。
準ベイズ法による推定の理論的根拠を簡単に紹介する。まず、疑似事後分布を真値からの 乖離幅によって書き直す。
h =√n(θ − θ0) − J(θ0)−1∆(θ0)/√n (15) とし、
p∗n(h) =√npn(h/√n + J(θ0)−1∆(θ0)/√n) (16) とする。ただし、J(θ0)はLn(θ)の2次微分の極限、∆n(θ)はLn(θ)の一次微分である。す ると、ある条件のもとで、p∗n(h)はn → ∞のとき、平均0分散J(θ0)−1の正規分布の密度 に収束する。したがって、pn(θ)はnが大きければ、平均が
θ0+ J(θ0)−1∆(θ0)/n (17)
であり、分散が、J(θ0)−1/nの正規分布で近似できる。そのため、pn(θ)の平均や中央値は、 θ0への一致性をもち、またその漸近分布は、極値推定量のものとは同じになる。
準ベイズ法を用いることで信頼区間を構築することもできる。たとえば、95%信頼区間を 得る場合には、(θ1, . . . , θB)の2.5%分位点と、97.5%分位点からなる区間を計算するとよい。 そのため、準ベイス法は、推定量の計算が難しい場合だけでなく、推定量の計算はできるが、 信頼区間や標準誤差の計算が難しい場合でも有用である。そのような例としては、分位点回 帰があげられる。
準ベイズ法の疑似事後分布は、津城のベイズ統計における事後分布と同じような性質をも ち、たとえば、モデル選択などにも使用することができる。Inoue and Shintani (2014)は準 ベイズ法によるモデル選択基準を考えている。Inoue and Shintani (2014)では、マクロの DSGEモデルへの適用を主眼においているが、もちろん他の場合にも使用可能である。
参考文献
[1] A. C. Cameron and P. K. Trivedi. Microeconometrics, Methods and Applications. Cambridge University Press, 2005.
[2] V. Chernozhukov and H. Hong. An MCMC approach to classical estimation. Journal of Econometrics, 115:293–346, 2003.
[3] S. Chib. Markov chain Monte Carlo method: Computation and inference. In J. J. Heckman and E. Leamer, editors, Handbook of Econometrics, volume 5, chapter 37, pages 3569–3649. Elsevier Science B.V., 2001.
[4] A. Inoue and M. Shintani. Quasi-Bayesian model selection. mimeo, 2014.
