講義ノート
7
：
DSGE
モデルのパラメータ推定
蓮見 亮
2013
年
3
月
7
日
目次
1 ベイズ統計学入門 118 1.1 ベイズ統計学の枠組み . . . 118 1.2 確率分布 . . . 123 1.3 ベイズ推定の手順 . . . 125 1.4 ベイズ推定の例：ベルヌーイ過程♠ . . . 127 1.5 ベイズ統計学でのモデル比較♠♠ . . . 128 1.6 ベイズ統計学は頻度論的統計学とどう違うか♠♠ . . . 130 2 マルコフ連鎖モンテカルロ（MCMC）法♠ 131 2.1 モンテカルロ法♠ . . . 132 2.2 マルコフ連鎖♠ . . . 132 2.3 メトロポリス・ヘイスティングス（M-H）アルゴリズム♠ . . . 134 2.4 MCMC法の収束判定♠♠ . . . 136 3 DSGEモデルのパラメータのベイズ推定 137 3.1 状態空間モデル♠ . . . 138 3.2 カルマン・フィルタ♠ . . . 138 3.3 DSGEモデルのパラメータのベイズ推定の手順 . . . 140 4 例：ニューケインジアン・モデルのパラメータのベイズ推定 142 4.1 モデル、データ、事前分布 . . . 142 4.2 事後分布、歴史分解 . . . 145最適成長モデルやそれを改良したモデルに誤差項を付加して、実証分析・予測を可能
にしたモデルをDSGEモデル（Dynamic Stochastic General Equilibrium Model：動学
的・確率的一般均衡モデル）という。誤差項には、現実データとモデルの理論値の差を埋 める役割がある。前章まで各モデルのパラメータは“手置き”していたが、誤差項を含む DSGEモデルはデータを用いたパラメータ推定（データに合うようなパラメータ値を推 測すること）が可能であることが知られている。例えば、ノート5で紹介したニューケイ ンジアン・モデルには技術ショックと金融政策ショックという2つの誤差項を含むため、 DSGEモデルとしてほぼそのまま実証分析に用いることができる。 DSGEモデルのパラメータは、多くの場合とりうる値の範囲が決まっている。例えば、 ニューケインジアン・モデルの価格非改定確率ωは、0 < ω < 1のいずれかの値でなけれ ばならない。このような事前の制約がある場合に、ベイズ推定という方法を用いると、そ のような事前の制約を事前情報という形で自然に織り込むことができる。ベイズ推定を含 むベイズ統計学と呼ばれる一連の理論体系について、1節で紹介する。統計学には事前情 報を用いない頻度論的統計学という理論体系もあり、それとベイズ統計学の違いについて も併せて述べる。 知られている方法によりデータを用いてDSGE モデルのパラメータを推定するには、 対数線形近似したモデルを用意する必要がある。対数線形近似したモデルのBlanchard
and Kahnの方法などによる行列を用いた線型VAR表現は、3節で説明するように状態 空間モデルと呼ばれる時系列モデルの遷移方程式とみなすことができる。状態空間モデル のパラメータは、マルコフ連鎖モンテカルロ法（MCMC法、2節）とカルマン・フィル タというアルゴリズムを組み合わせることで推定できる。 ノート5の補論で既に対数線形近似したニューケインジアン・モデルを用意してあるた め、データを用いたニューケインジアン・モデルのパラメータの推定例を4節で紹介する。
1
ベイズ統計学入門
まず、ベイズ統計学に関する基本的な事項について述べる。“入門”と銘打っていると はいえ、ベイズ統計学では理論的一貫性を重視することから、実際のベイズ推定の例は、 全て1節で説明する基本原理の応用に過ぎない。1.1
ベイズ統計学の枠組み
ベイズ統計学とは、以下のような特徴を有する統計学の理論体系である。 1. ベイズ統計学でいう確率は、主観確率である。 2. 統計的推論に必ずベイズの定理を用いる。3. 統計的推論に何らかの形で事前情報を用いる。 4. 未知パラメータ、未観測データを確率変数とみなす。 5. 観測されたデータは固定された値（非確率変数）とみなす。 ここでいう統計的推論とは、観測データから未知パラメータのもっともらしい値を推測す ることをいう。 確率の定義と主観確率 まず、確率の定義について、主観確率よりも狭い概念である頻度論的確率から説明して いく。一般に、確率とは物事の起こりやすさ、起こりにくさを数値化したものであり、事 象A が起こる確率P (A)は[0, 1]のいずれかの値として表される。確率が大きいほど起 こりやすい事象だと我々は考える。頻度論的確率とは、事象Aが起こる確率P (A)を P (A) = 実際にAが起こった回数 Aが起こりうる回数 (1) と定義する。例えば、事象A′ をある地点で雨が観測された日数とする。去年一年間を観 測期間とすると、その期間に事象A′ が起こった確率P (A′)は、 P (A′) = 雨が観測された日数 365 で与えられる。 では、明日雨が観測される確率（事象A′′ とする）についてはどのように考えたらよい だろうか。期間を明日に限定すれば、頻度論的確率の定義からは分母のA′′が起こりうる 回数は1以外にありえない。分子は雨が観測されるか（= 1）雨が観測されないか（= 0） なので、P (A′′)は1か0かの値しかとり得ない。これに対して、ベイズ流の主観確率で は、確率をもっと自由に考える。例えば、明日の天気についていえば、自分がどう思うか を[0, 1]のいずれかの値で表すことが許される。したがって、天気予報から判断してある 雨が降る確率は0.3である、といった表現が許容される。なお、ベイズ流の主観確率は頻 度論的確率も包含する概念であり、(1)式で定義されるようなP (A)も確率として取り扱 う。確率の定義は他にもいくつかあり、以下で説明するように例えば数学では、確率とい う概念をより抽象化して、全事象という集合に属する事象AをP という物差し（関数） で測った大きさとして定義する。ベイズ流の主観確率は確率の数学的定義とも矛盾しない と考えてよい。 条件付確率とベイズの定理 ベイズの定理とは、条件付確率 P (B|A) = P (A∩ B) P (A) (2)
に関連する定理である。(2)式は条件付確率を定義する式であり、右辺分母の P (A∩ B) は事象 A が起こりかつ事象 B が起こる確率を意味する。それを事象A が起こる確率 P (A)で割ったものが、左辺の条件付確率P (B|A)であり、事象Aが起こったという条件 のもとで事象Bが起こる確率をいう。 例えば、トランプ（ジョーカーを含まない）のカードを一枚引くというゲームを考えよ う。事象Aと事象Bは、それぞれ以下のような集合の部分集合とする。 A∈ {♣, ♢, ♡, ♠} B∈ {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13}. この定義によると、例えばP ({♣} ∩ {1})はクラブのエースを引く確率であるから、1/52 である。では、P ({1}|{♣})は、クラブを引いたという条件のもとでそれがエースである 確率であるから1/13である。(2)式を用いると P ({1}|{♣}) = P ({♣} ∩ {1}) P ({♣}) = 1 52 1 4 = 1 13 (3) だから、確かに矛盾なく定義できている。 条件付確率の定義(2)式のAとBを入れ替えると P (A|B) = P (A∩ B) P (B) (4) であるが、P (A∩ B) = P (B|A)P (A)を代入すると（これは(2)式の変形である） P (A|B) = P (B|A)P (A) P (B) (5) と表すことができ、この(5)式をベイズの定理という。これを統計的推論に用いるのがベ イズ統計学である（具体例については1.4節参照）。 事前情報、事前分布 事前情報とは、統計的推論の対象となる未知パラメータについて我々が知っている情報 をいう。事前情報は事前分布という形で確率を用いて表現する。事前分布は、我々が主観 的に決めなければならない*1。例えば、ニューケインジアン・モデルの価格非改定確率ω をベイズ推定しようとする場合、ωは定義上確率を意味することから0 < ω < 1のいず れかの値でなければならないという事前情報を我々は用いる。さらに、事前情報は確率と して表現しなればならないが、その際に必ず任意性が残る（ω = 0, 1はどちらも不都合な *1事前分布を主観的に決めなければならないことは、ベイズ統計学でいう確率は主観確率でなければならな いことの1つの理由である。
ので除外する）。ω は0に近い可能性が高いのか、1に近いのか可能性が高いのか、それ とも0と1で満遍なく可能性があるのか、などである。他の推定例などの情報を参考にす る場合も多々あるが、たとえそのようなものがなくとも、我々は自己の信念に基づき事前 分布を決めなければならない。この例では推定しようとする未知パラメータについて事前 の情報があったが、たとえなんら情報がない場合でも、情報がないという形で事前情報を 用いる。 確率の数学的定義 上記の 4.に確率変数という用語があるため、その説明に必要な確率の数学的定義を説 明しておく。数学では、確率という概念をより抽象化して、全事象という集合に含まれる 集合（つまり全事象の部分集合）である事象E を確率測度P という物差しで測った大き さとして定義する。全ての可能な結果を根元事象といい、根元事象を全て集めた集合が全 事象である（つまり、根元事象は全事象の要素）。 例として、ジョーカーを含まないトランプの山から一枚のカードを引くという試行を考 えよう。この場合、根元事象は、クラブの1に対してω1、クラブの2に対してω2といっ た順に番号を割り当てていくことにする。全事象は、定義によりΩ = {ω1, ω2, . . . , ω52} と表される。クラブの1がでるという事象はE1 ={ω1}と表されるが、これは全事象の 部分集合である。これを図示すると図1のようになる。可能な限り分割した後のタイルの 一枚一枚の上の“点”が根元事象、それらを集めたタイル全体が全事象である。我々は、 ある1回の試行において、根元事象のどれか1つが実現する（あるいは実現している）と 考える。 例えば、1のカードがでるという事象は E2 = {ω1, ω14, ω27, ω40}と表すことができ、 クラブのカードがでるという事象はE3 = {ω1, ω2, . . . , ω13}と表すことができる。これ らはいずれも全事象（タイル全体）の部分集合である。それをタイルで表現したのが図2 である。タイルの組み合わせが事象であることがわかる。 このタイルの“大きさ”を確率測度P で測ったものが確率である。トランプの山から一 枚カードを引くという試行がフェアなゲームならば、タイル一枚一枚の大きさを確率測度 P で測ると全て同じ大きさで、P (Ω) = 1よりタイル一枚一枚の大きさは 1/52である。 したがって、上記の事象の例では、 P (E1) = 1/52 P (E2) = 4/52 = 1/13 P (E3) = 13/52 = 1/4 といった具合に、ごく当たり前の結果が得られる。この事例だとタイルをこれ以上区切っ ても意味がないが、それは根元事象に一般の性質である。例えば、クラブの1のタイルを
! "! #! $! %! & ' ( ) * +, ++ +-! +.! /"! /&#! 0! /#! 図1 全事象Ωと根元事象ωi ! !"#$%&' !"#$%&'( )*+,-./01 23 !"#$%&4(1 図2 事象の例（2つ） 2つに分割しても、どちらの根元事象が実現したのか区別がつかない。 タイル一枚一枚の可能な組み合わせを全て集めたもの、すなわちΩの部分集合の全体 F を事象の族という*2。事象は全事象Ωの部分集合だから、事象の族F の要素である。 (Ω,F, P )を確率空間という。確率測度P は関数で、P :F → [0, 1]と書ける。 確率変数 確率の数学的定義が説明できたところで、確率変数の定義について述べる。ある確率 空間(Ω,F, P )が定義されていたとする。確率変数X とは、ΩからRへの関数をいう。 すなわち、X : Ω → Rである。よく用いられる関数の表記法を使うと、確率測度P は P (E)と書けるのに対し、確率変数X はX(ω)と書ける。以下に確率変数の例を1つ挙 げておく。 【例】トランプを用いた賭けの賞金 一定の参加料を支払うことで、引いたトランプのカードに応じて賞金がもらえるという 賭けを考える。賭けの内容は、ある一定金額yを単位とし、引いたカードの数字がaとす *2部分集合全体はべき集合とも呼ばれる。
ると、そのカードがクラブならばay、ダイアならば2ay、ハートならば3ay、スペードな らば4ayの賞金がもらえるものとする。このとき、賞金は根元事象の関数となるから確率 変数である。 ここでは、確率変数の説明にとどめ、上記の4.と5.の具体的な意味については1.4節 で実例に即して説明する。
1.2
確率分布
連続的な値をとる確率変数の、そのとりうる各々の値に対する起こりやすさは、確率密 度関数により記述できる。p(x)を確率変数X(ω)の確率密度関数とし、X(ω)がある実数 zに対してX(ω)≤ z という値をとるという事象をEX≤z とすると、 P (EX≤z) = ∫ z −∞p(x)dx (6) という関係がある*3。この関係からわかるように、確率変数はより高い確率で確率密度の “厚い”区間の値をとる。 全事象をR、事象の族をボレル集合族B（おおざっぱに、実数線上の全ての開集合、閉 集合およびその和集合を全て集めたもの）としたとき、その上で定義される確率測度P を確率分布という。上記の議論から明らかなように、確率変数の値域と確率密度関数が与 えられれば、確率分布が定義できる。 ベイズ推定の際によく用いられる確率分布はいくつかあるが、確率変数のとりうる値が 連続か離散かによって大まかに分類される。連続型の一次元確率分布の例としては、正規 分布、ガンマ分布、ベータ分布、一様分布などがある。それぞれの確率分布の定義域、お よび確率密度関数（離散型については確率関数）は表1の通りである。このうち、正規分 布、ガンマ分布、ベータ分布、一様分布の確率密度関数をグラフに描くと、図3のように なる。離散型の一次元確率分布の例としては、ベルヌーイ分布、二項分布などがある*4。 *3確率変数が多次元の場合は多重積分になる。 *4 離散型の場合は、確率測度P が積分を用いずに書け、ベルヌーイ分布の場合は、0≤ θ ≤ 1をパラメー タとして P (x = k; θ) = θk(1− θ)1−k, k∈ 0, 1 (7) で、二項分布の場合は、同じく0≤ θ ≤ 1をパラメータとして P (x = k; θ) = n! k!(n− k)!θ k(1− θ)n−k, k = 0, 1, 2, . . . , n (8) である（n! =∏ni=1i）。確率分布 確率変数の値域 確率密度関数p(x) 平均 分散 正規分布 N(µ, σ2) −∞ < X < ∞ √ 1 2πσ2 exp ( − (x−µ)2 σ2 ) µ σ2
ガンマ分布 Ga(a, s) X ≥ 0 saΓ(a)1 xa−1exp(−xs) as as2
ベータ分布 Be(a, b) 0≤ X ≤ 1 Γ(a)Γ(b)Γ(a+b) xa−1(1− x)b−1 a+ba (a+b)2ab(a+b+1)
一様分布 Unif(a, b) a≤ X ≤ b, a, b ∈ R b−a1 b−a2 (b−a)12 2
表1 主な一次元確率分布（Γ(·)はガンマ関数） −4 −2 0 2 4 0.0 0.1 0.2 0.3 0.4 N(µ = 0, σ2 = 1 ) 0 1 2 3 4 0.0 0.2 0.4 0.6 0.8 1.0 Ga(s = 1, a = 1 ) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 1.5 Be(a = 2, b = 2 ) −1.0 −0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Unif(a = −1, b = 1 ) 図3 正規分布、ガンマ分布、ベータ分布、一様分布の確率密度 確率密度関数p(x)をもつ確率変数Xの平均（期待値）は、 E(X) = ∫ ∞ −∞ xp(x)dx (9) で定義され、分散は平均をµとすると Var(X) = ∫ ∞ −∞ (x− µ)2p(x)dx (10) で定義される。分散は、確率分布のばらつきを示す数値であり、Xの分散をσ2 とすると、 例えばチェビシェフの不等式より P ({|X − µ| ≥ kσ}) ≤ 1 k2, k > 1 (11)
がいえる。
1.3
ベイズ推定の手順
ベイズの定理は、確率だけでなく確率密度関数についても成り立つ。確率密度関数につ いて、p(a, b)を確率変数A, Bの同時確率密度関数とすると、Aが与えられたときのBの 条件付確率密度関数p(b|a)は p(a|b) = p(a, b) p(a) (12)で与えられる。ただし、p(a)はAの周辺確率密度関数で、p(a) =∫−∞∞ p(a, b)dbである。
この定義を用いると、確率密度関数についてのベイズの定理 p(a|b) = p(b|a)p(a) p(b) (13) が得られる*5。ベイズ推定は、この(13)式を用いて以下の順序で行う。 *5なお、p(a)とp(b)はそれぞれ関数型が異なるので、正確にはpA(a)とpB(b)と書くべきだが、煩雑に なるので下付き文字は適宜省略する。p(a|b)とp(b|a)についても同様である。その場合には、引数の記 号が違えば違う関数であると解してほしい。
ベイズ推定の手順 1. パラメータが与えられた場合のデータの条件付分布（モデルともいう）を与え る。(13)式のp(a|b)の部分にあたるが、aをデータy、bをパラメータθで置 き換えpy|θ(y|θ)と表記する。 2. 事前分布pθ(θ)を与える。(13)式のp(b)の部分にあたる。 3. 観測データyobs を条件付分布py|θ(y|θ)のyに代入する。データの代入後は、 py|θ(y|θ)の部分は条件付分布ではなく尤度関数l(θ; yobs)となる。 4. ベイズの定理から、事後分布pθ|y(θ|yobs)が pθ|y(θ|yobs) = l(θ; y obs)p θ(θ) py(yobs) (14) と与えられる。分母py(yobs)は分子のθによる積分 py(yobs) = ∫ l(θ; yobs)p(θ)dθであってθ に依存しない定数なので、比例記号 ∝を用いて pθ|y(θ|yobs)∝ l(θ; yobs)p(θ) (15) と書いてもよい。 5. パラメータの事後分布pθ|y(θ|yobs)を何らかの形で評価する。 事後分布には観測データから得られた未知パラメータのもっともらしい値に関する情 報が含まれているが、5.でそれを“取り出す”ということである。具体的な評価の方法は 複数ある。例えば、θ が一次元であれば、事後分布の平均は∫−∞∞ θpθ|y(θ|yobs)dθの積分 計算によって得られる。多次元であれば、MCMC法（2節）を用いることが多い。分母 py(yobs)は周辺尤度といい、データyobsが得られる確率（実際には確率密度関数の厚みを 意味するので、確率に比例する量）を意味することから、モデル比較（1.5節）に用いる。 以上より、ベイズ推定とは、モデルと事前分布とデータから得られた事後分布を評価す ることによって、未知パラメータのもっともらしい値を推測する一連の手続きであること がわかる*6。 事前分布および条件付分布（モデル）には、何らかの確率分布の確率密度関数を指定す る。まず、事前分布にどのような確率分布を使えばよいかについてだが、多くの場合、推 定しようとする未知パラメータのとりうる値を手がかりとする。未知パラメータθのとり うる値が−∞ < θ < ∞ならば正規分布、0 < θ <∞ならばガンマ分布、0≤ θ ≤ 1なら ばベータ分布を選ぶことが多い。a ≤ θ ≤ bで一様に分布しているとするならば、一様分 *6この場合、事前分布は事前確率密度関数、事後分布は事後確率密度関数というほうが正確だが、慣例的に それぞれ事前分布、事後分布と呼ばれることが多い。
布を選ぶ。実用的にも、このような事前分布を用いると事後分布の評価がしやすいことが 多い。未知パラメータが多次元である場合には、それぞれが独立の事前分布に従うとする ことが多いが、もちろん相互間の相関や階層構造を考えてもよい*7。DSGEモデルのパラ メータは多くの場合とりうる値が決まっているため、それに応じて事前分布の確率分布を 選べばよい。 そのように確率分布を選んだ後で、何らかの利用したい情報があればそれを事前分布の パラメータ（ハイパーパラメータと呼ばれることがある）に反映させることができる。な んら情報がないという情報を反映させたいならば、事前分布の分散（ばらつき）を充分に 大きく取ればよい*8。正規分布、ガンマ分布、ベータ分布にはそれぞれパラメータが2つ 存在するため、それらに事前情報を反映させる。 ベイズの定理で条件付分布 py|θ(y|θ)にあたる部分にも、同様に何らかの確率分布の確 率密度関数を仮定することになる。どのような確率分布を想定するかは、どのようなモデ ルのパラメータを推定しようとしているかに応じて自然と決まる。例えば次の1.4節の例 では、ベルヌーイ分布を仮定する。DSGEモデルのパラメータをベイズ推定しようとす る場合、多くはモデルとして多変量正規分布を仮定する（3節参照）。 以下、数学的・技術的な話題が続くため、DSGEモデルのパラメータの推定例のみに関 心のある読者は、3節に進んでいただきたい。
1.4
ベイズ推定の例：ベルヌーイ過程
♠
表になる（つまり裏にならない）確率をパラメータとして任意の一定の値に設定できる 特殊なコインを用いてコイントスを行い、表になった回数と裏になった回数を記録して いくことでパラメータの値をベイズ推定するという実験を考えよう*9。毎回の試行をベル ヌーイ試行といい、その繰返しをベルヌーイ過程という。推論する人は、パラメータの値 を知らないものとする。 コイントスを行った回数をN 回とする。表がでたときには 1、裏がでたときには0と 結果を記録していくことにする。すると、データは0と1の列の集合k ={ki}Ni=1, ki ∈ {0, 1}として得られる。推論する人にとっての未知パラメータである表になる（裏になら *7ここでいう独立とは、大雑把にいうと、互いに無関係であることをいう。 *8情報がない場合の最も客観的な事前分布については、Jeffreysの局所一様事前分布などいくつか提案され ているが、決定的なものはない。Jeffreysの局所一様事前分布は事前分布としてimproperな(つまり積 分すると発散する)ため、周辺尤度の計算ができなくなるという不都合が生じる。 *9コインを折り曲げたりすれば、多少なりとも表になる（裏にならない）確率を変えられるだろう。ない）確率をθとすると、尤度関数は f (θ; k) = N ∏ i=1 θki(1− θ)1−ki (16) で与えられる。θ のとりうる値は[0, 1] なので、事前分布として、ベータ分布Be(a, b)、 すなわち π(θ)∝ θa−1(1− θ)b−1 (17) を仮定する。ベータ分布のパラメータa, bに事前情報を反映させる*10。すぐ下の結果か らわかるが、基準化定数はここでは無視してよいので、比例記号∝を用いている。ここ までは、未知パラメータ、データとも確率変数として取り扱う。 尤度関数と事前分布の掛け算が事後分布なので、事後分布は、 π(θ|k) ∝ θK+a−1(1− θ)N−K+b−1, K =∑ki (18) で与えられる。これは、θがBe(K + a, N − K + b)に従うことを意味するので、いった ん観測データが得られれば（それは固定された値で非確率変数とみなす）、事後分布の平 均や分散は簡単に求まる。 乱数を使ってモンテカルロ法（2.1節参照）によりN を変えてシミュレーションを行 い、事前分布と事後分布の関係を見たのが図4である。ここで、未知パラメータ θ は確 率1でθ = 0.6であるとして固定し、事前分布のパラメータa, bはa = 2, b = 2とおい た。図4の左下にあるように、これは0.5をモード（頂点）とする左右対称の分布である。 N = 5, 20, 100, 500とN を増加させていくと、N が小さいときには事後分布に広がり があるが、N が大きくなるにつれて尤度の裾が狭まり、事後分布が固定したθの値の周辺 の狭い範囲に収束していくことがわかる。つまり、N が大きいほど高い精度でθ の値が 言い当てられるようになる。 この例では、問題の設定自体がとてもシンプルだったため、事後分布の平均や分散が簡 単に求まったが、多くの場合には事後分布を解析的に評価することはできない。その際に は、2節で説明するMCMC法を用いる。
1.5
ベイズ統計学でのモデル比較
♠♠
周辺尤度 py(yobs) は、モデルからデータ yobs が得られる確率に比例する量を意味 するので、これの大小を比較することでモデル比較ができる。例えば、m 個のモデル *10ベータ分布Be(a, b)の平均は a a+b、分散は ab (a+b)2(a+b+1) である。0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 N = 5 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 N = 20 0.0 0.2 0.4 0.6 0.8 1.0 0 2 4 6 8 10 N = 100 0.0 0.2 0.4 0.6 0.8 1.0 0 5 10 15 20 25 30 N = 1000 図4 未知パラメータθの分布（点線が事前分布、実線が事後分布、垂線が真の値） M1, M2, . . . , Mmがあるとし、モデルMi のパラメータをθi と書く。(14)式を Mi を明 示的に使って書き直すと、モデルMi に対するθi の条件付事後分布は、 p(θi|y, Mi) = p(y|θi, Mi)p(θi|Mi) p(y|Mi) (19) と表される*11。 次に、Mi を未知パラメータ、yをデータとしてベイズの定理を適用する。ベイズの公 式(13)式のaにyを、bにMiを代入すると、Miの事後分布 p(Mi|y) = p(y|Mi)p(Mi) p(y) (21) が得られる。ここで、p(y|Mi)は(19)式の分母で、モデルMi の周辺尤度である。p(Mi) *11 条件付確率の定義から、 p(b|a, c) = p(a|b, c)p(b|c) p(a|c) (20) は容易に証明できる。
は、事前分布にあたる*12。なお、分母p(y)は、 p(y) = m ∑ i=1 p(y|Mi)p(Mi) (22) で、定数とみなせる。 モデルiとモデルj のどちらがもっともらしいかは、p(Mi|y)とp(Mj|y)の比を取れば よく、 p(Mi|y) p(Mj|y) = p(y|Mi)p(Mi) p(y|Mj)p(Mj) (23) でみればよい。脚注*12のように事前分布をp(Mi) = p(Mj) = 1/mとすると、 p(Mi|y) p(Mj|y) = p(y|Mi) p(y|Mj) (24) であって、周辺尤度の大小を比較すればよい。なお、周辺尤度の比p(y|Mi)/p(y|Mj)はベ イズファクタと呼ばれる*13。p(y|Mi)の具体的な計算方法については、例えば“Bayesian Econometrics”（Koop[2003]）を参照されたい。
1.6
ベイズ統計学は頻度論的統計学とどう違うか
♠♠
本筋からやや外れるが、まとめとしてベイズ統計学の頻度論的統計学との違いを簡単に 述べておく*14。 理念として • 理論的一貫性を重視する。アドホックな手続きを認めない。例えば、GMMと最尤 法のどちらで推定するか、AICとBICのどちらを用いるかといった問題も生じな いし、推定量の小標本特性を良くするような微調整などもない。 •（ベイズの定理を用いない）検定や有意水準といった概念はない*15。 • 確率は常に主観確率である。 *12 事前分布はp(Mi) = 1/mと置くのが自然である場合が多い。これは、事前には全てのモデルが同じく らいもっともらしいと仮定することを意味する。 *13周辺尤度（正確には、−2 ln (p(y|Mi))）の近似がBICにあたる *14ここでいう頻度論的統計学とは、ネイマン・ピアソン流の仮説検定理論、フィッシャー流の最尤法、およ びそれらを発展させた非ベイズの統計学、計量経済学をいう。 *15 MCMC法の収束判定にはネイマン・ピアソン流の仮説検定理論（例えばGewekeの収束判定統計量、 2.4節）を用いることがあるので、それらの概念を完全に否定するわけではないようだ。もちろん、これ はあくまで補助的な検定のレベルであり、手軽に計算できるという実用面での事情による。理論として • パラメータとデータの区別が曖昧である*16。 • パラメータは確率分布として存在すると仮定する。パラメータに真の値があるとは 考えない。一致性といった概念もない。 • 得られたデータは確率変数とはみなさない。 • 漸近理論（大数の法則、中心極限定理）に依存しない。つまり、無限回の試行（頻 度論的確率）を前提としない。 • 漸近的には最尤法と数値的に同等の結果が得られる。通常、大標本では事前分布が 事後分布に影響を与ない。 • モデルを常にパラメトリックに指定する。ノンパラメトリックな推定はありえ ない。 実用面 • 最尤法に比べて推定できるモデルのクラスが広い。ただし、これはMCMC法の柔 軟性によるものであって、ベイズ統計学の本質に由来するものではない。 • 事前分布が主観的であることをもって、ベイズ統計学を批判する向きもあるが、事 前分布の分散を大きい値にとることで客観性を確保できる。また、ほとんどの場 合、真のモデルを知ることができない以上、モデル選択については客観性を装いた くとも実際上不可能であり、事前分布だけが主観的であるわけではない。 • 事前情報を利用できるという特徴を生かすべきとする考え方もある*17。
2
マルコフ連鎖モンテカルロ（
MCMC
）法
♠
マルコフ連鎖モンテカルロ（Markov Chain Monte Carlo：MCMC）法は、確率密度
関数の平均や分散、中央値・パーセンタイルを求めるための1つの方法である。ベイズ推 定する際の事後分布の評価に用いることができる。事後分布が正規分布やベータ分布と いった知られた分布に従うならば、平均や分散は簡単に計算できるが、そうでない場合は MCMC法を用いて、事後分布に従う擬似乱数の列を生成して平均や分散を計算する。 MCMC法は、モンテカルロ法とマルコフ連鎖という2つのアルゴリズムを組み合わせ た方法である。DSGEモデルのパラメータ推定の際にも、MCMC法のうちメトロポリ *16したがって、例えば欠損値に対応しやすい。 *17例えば、DSGE-VARモデルなど。
ス・ヘイスティングス（M-H）アルゴリズムと呼ばれる方法を用いる。
2.1
モンテカルロ法
♠
モンテカルロ法とは、コンピュータにより生成される擬似乱数を用いる数値計算法の総 称である*18。まず一様乱数（[0, 1]で一様に分布する乱数）を生成し、それを加工するこ とで、主要な確率分布に従う乱数を生成できる*19。 数値計算ソフトウェアには、多くの場合、一様乱数を含む擬似乱数生成のアルゴリズム が組み込まれており、例えば統計・数値計算用のフリーソフトウェアであるRはメルセン ヌ・ツイスタと呼ばれるアルゴリズムで一様乱数を生成する。その一様乱数を用いて、例 えばボックス・ミュラー法と呼ばれる方法で正規分布に従う乱数を生成できる。ガンマ分 布、ベータ分布に従う乱数の生成法も知られているが、数値計算ソフトウェア自体にそれ ら主要な確率分布に従う乱数を返す関数が用意されていることも多い。2.2
マルコフ連鎖
♠
確率変数を時系列順に並べたものを確率過程という。例えば、1.4節の例のようにコイ ントスの結果を表が1、裏が0とし、i番目の結果をkiと書くことにすると(ki ∈ {0, 1})、 ki を数列として並べた{k1, k2, . . .}は一つの確率過程である。マルコフ連鎖とは、未来 の状態が現在の状態だけで決定され、過去の経路と無関係である確率過程をいう。 【例】天気の過程 天気 Wt を{晴,雨}のいずれかのみで記録していくことにしよう。明日の天気は今日 の天気だけで決まり、今日が晴れならば明日晴れの確率は0.9、雨の確率は0.1、今日雨な らば明日晴れの確率は0.5、雨の確率は0.5とする。このとき天気の列{W1, W2, . . .}は マルコフ連鎖である。この過程は Q = [ 0.9 0.5 0.1 0.5 ] (25) という行列で特徴づけられ（推移核という）、π晴+ π雨= 1という条件を用いて [ 0.9 0.5 0.1 0.5 ] [ π晴 π雨 ] = [ π晴 π雨 ] (26) という方程式を解くことで、定常状態での晴れと雨の確率がそれぞれπ晴 = 5/6, π雨 = 1/6と求まる。 *18もちろん、擬似乱数に限らず物理乱数を用いてもよい。 *19結論を先にいってしまうと、後述のM-Hアルゴリズムを用いることで、理論上いかなる確率分布に従う 乱数も生成できる。0 50 100 150 200 −4 −2 0 2 4 at 図5 AR(1)過程のシミュレーション例 上記の例では、状態は離散変数だったが、次に状態が連続変数である場合のマルコフ連 鎖の例を示す。 【例】AR(1)過程 εt を平均ゼロ、分散σ2 の正規分布に従うそれぞれ独立な確率変数とし、−1 < ρ < 1 とする。確率変数at がa0 = ¯aとして at+1 = ρat+ εt+1 (27) という過程に従うものとしよう。このとき、どのようなa0 = ¯aを初期値としても、充分 大きいT に対してat′, t′ ≥ T は(27)式で表されるAR(1)過程の定常分布 at′ ∼ N ( 0, σ 2 1− ρ2 ) (28) に従う。この場合の推移核q(at, at+1)は、正規分布N(ρat, σ2)の確率密度関数である。 この例は、すぐ後に説明するようにMCMC法の基本原理と密接なかかわりがある。 図5は、パラメータをσ2 = 1, ρ = 0.5とし、初期値を0, −100, 100の3通りに設定し てat の実現値を乱数を用いてシミュレーションしたものである。初期値を定常分布の外 れ値にしたとしても、この例では比較的早く定常分布N(0, 4/3)に収束することがわかる。
2.3
メトロポリス・ヘイスティングス（
M-H
）アルゴリズム
♠
MCMC法は、ある関心ある分布（ベイズ推定では事後分布）を定常分布とするような マルコフ過程を作り、適当な初期値を置いて順繰りに乱数を用いて計算することで定常 分布に収束させ、その関心ある分布を評価する（例えば平均、分散を計算する）方法であ る。いったん定常分布に収束させることができれば、その後は定常分布である事後分布か ら乱数サンプリングしているものとみなせる。上記のAR(1)過程の例では、MCMC法 を使わなくとも定常分布から直接サンプリング可能だが、ベイズ推定の場合、事後分布か ら直接乱数をサンプリングできることはまれである。その際に、MCMC法を用いて計算 を行う。 MCMC 法は、いくつかのアルゴリズムの総称だが、そのうちの 1 つとして M-H （Metropolis-Hastings：メトロポリス・ヘイスティングス）アルゴリズムが知られてい る*20。これは、理論上どのような分布に従う乱数もサンプリングできるとても一般的な 方法である。 M-Hアルゴリズムでは、ある推移核qを用いて以下のような手順で任意の確率密度関 数f (θ)から乱数を生成する*21。 M-Hアルゴリズム 初期化：θ(0) を設定 繰り返しi（i ≥ 1）： step 1. ˜θ ← q(θ(i−1), θ) step 2. 採択確率 α(θ(i−1), ˜θ) = min { f (˜θ)q(˜θ, θ(i−1)) f (θ(i−1))q(θ(i−1), ˜θ), 1 } (29) を計算 step 3. 確率α(θ(i−1), ˜θ)で、θ˜を採択して、θ(i) = ˜θとおく それ以外はθ˜を棄却して、θ(i) = θ(i−1) とおく f (θ)の基準化定数は未知でよい。 M-Hアルゴリズムにおいて、q(ϕ, θ) = q(θ, ϕ)が成り立つような推移核qを用いると、 *20ギブスサンプラーと呼ばれるアルゴリズムもMCMC法の一種である。ギブスサンプラーは、M-Hアル ゴリズムと組み合わせることもできる。 *21qはあくまでstep 1.での推移核であって、M-Hアルゴリズム全体としての推移核は別に存在する。上記のα(θ(i−1), ˜θ)が、 α(θ(i−1), ˜θ) = min { f (˜θ) f (θ(i−1)), 1 } (30)
と簡単になる。これを酔歩連鎖（random walk chain）という*22。θが一次元のときに限
らず、多次元であってもM-Hアルゴリズムを用いることができる（以下、そのときには 太字でθと書く）。酔歩連鎖の採択確率は、0.2∼ 0.5が望ましいとされる*23。 以下のような手順により、推移核q にrandom walk（誤差項に平均ゼロの正規分布を 仮定する）を用いた酔歩連鎖により、任意の確率密度関数f (θ)から乱数を生成できる。 酔歩連鎖（θが一次元のとき） 初期化：θ(0) を設定 繰り返しi（i ≥ 1）： step 1. ˜θ ← θ(i−1)+ ϵ, ϵ∼ N(0, σ2) step 2. 採択確率 α(θ(i−1), ˜θ) = min { f (˜θ) f (θ(i−1)), 1 } (31) を計算 step 3. 確率α(θ(i−1), ˜θ)で、θ˜を採択して、θ(i) = ˜θとおく それ以外はθ˜を棄却して、θ(i) = θ(i−1) とおく 【例：酔歩連鎖による乱数生成（ガンマ分布）】 M-Hアルゴリズムの適用例として、酔歩連鎖によりガンマ分布Ga(s = 5, a = 1)に従 う乱数を生成してみよう*24。ガンマ分布の確率密度関数の関数型については既知とする。 酔歩連鎖のパラメータσ2は、0.1, 20, 5000の3通りを仮定してサンプリングを行った。 酔歩連鎖の繰り返し回数は各ケース1,200で、うち最初の200は定常状態に至っていない ものとみなして捨てた。各ケースについて残った1,000サンプルを用いてヒストグラムを 描いたのが図6である。一般的な方法（脚注*24参照）でガンマ分布に従う乱数を生成し た場合の結果も併せて示した（図中における“参考”）。曲線がGa(s = 5, r = 1)の確率 密度、棒グラフがヒストグラムである。 *22酔歩連鎖は、Metropolisアルゴリズムとも呼ばれる。その名前からわかるように、こちらが先に発見さ れ、それを一般化したのがM-Hアルゴリズムである。
*23例えば、“Bayesian Econometrics”（Koop[2003]）の98頁を参照。
*24。実用上、酔歩連鎖を用いることなく、任意のパラメータのガンマ分布に従う独立な乱数をより簡単に生
酔歩連鎖のアルゴリズムstep 3.での採択確率は、各ケースそれぞれ0.944, 0.492, 0.045 であった。図6からわかるように、酔歩連鎖の3ケースの中ではσ2 = 20と設定した場 合に最も満遍なくGa(s = 5, r = 1)から乱数が生成できている。σ2 をどのような値にす べきか決める際には、目的とする分布f の分散がひとつの手がかりになる。この例では、 Ga(s, r)の分散はs/r2であることからσ2 = 5 が目安となる。σ2がそれより大幅に小さ いときには採択確率が高くなりすぎ、逆にそれより大幅に大きいときには採択確率が低く なりすぎる。いずれの場合もサンプリングの効率が低下する。すなわち、サンプリング回 数を非常に多くしないと、定常分布からサンプリングしたことにならない。 σ2 = 0.1 0 5 10 15 0.00 0.10 0.20 σ2 = 10 Density 0 5 10 15 0.00 0.05 0.10 0.15 0.20 σ2 = 5000 0 5 10 15 0.00 0.05 0.10 0.15 0.20 (参考) Density 0 5 10 15 0.00 0.05 0.10 0.15 0.20 図6 酔歩連鎖による乱数生成 θ が多次元のときには、step 1.で多変量正規分布MN(0, Σ)から乱数を生成する。3.3 節で説明するように、DSGEのパラメータ推定にはこの酔歩連鎖を用いる。
2.4
MCMC
法の収束判定
♠♠
MCMC法を実行する場合、どの程度の繰返し回数で定常分布へ収束するかは一概には いえない。推移核次第では、定常分布への収束に多くの繰返し回数を必要とする場合や、 定常分布に収束しない場合もある。したがって、何らかの方法により定常分布へ収束しているかどうかの収束判定を行う必要がある。 ひとつのやり方として、グラフに描くという方法がある。図6のように各チェイン（連 鎖）の軌跡をプロットするか、あるいは図5のようにヒストグラムを描いて、サンプリン グに偏りがないかを目で確認する。 もう少し客観的なやり方として、収束判定の統計量を計算するという方法がある。収束 判定の統計量はいくつか提案されているが、そのうち比較的簡単に計算できるものとし て、Gewekeの収束判定統計量がある。これは、頻度論的統計学の平均の差の検定の統計 量であり、その帰無仮説は、チェインの最初の10%と最後の50%の平均θ¯A, ¯θB が等し いことである。ここで、帰無仮説のもとでは、 Z = ¯ θA− ¯θB √ SθA(0)/NA+ SθB(0)/NB (32) が平均ゼロ、分散1の正規分布に従う*25。t検定同様、大雑把にZ 値は絶対値2以下を 目安とする*26。このとき、Z 値の絶対値が2以下であればそのチェインの最初の10%と 最後の50%の平均が等しいとみなし、定常分布へ収束していると判断する。
3
DSGE
モデルのパラメータのベイズ推定
状態空間モデルとは、観測されない潜在変数の過程を観測可能な観測変数から推測する ことを目的とする時系列分析のモデルをいい、観測方程式と遷移方程式という2つの方程 式からなる。この章の冒頭で述べたように、対数線形近似したDSGEモデルの線型VAR 表現に元のモデルにある外生のショック項ηt を加えた [ ˆ xt+1 ˆ st+1 ] = D [ ˆ xt ˆ st ] + Rηt+1 (33) という方程式は、状態空間モデルの遷移方程式とみなすことができる。状態空間モデルの パラメータは、前節のMCMC法と3.2 節で説明するカルマン・フィルタの漸化式と呼ば れる公式を組み合わせることで推定できる。 *25xt, t = 1, 2, . . . , N が独立同分布の場合、xt の平均の分散の推定量は Var(xt)/N で与えられる （Var(xt)はxt の標本分散）。MCMC法によりある分布に従う乱数xtを生成したときには自己相関が あるが（つまりxtそれぞれは独立ではない）、その場合、xtの平均の分散の推定量はSx(0)/N で与え られる。ただし、Sx(0)はxtの自己共分散関数のフーリエ変換（スペクトル密度）に0を代入したもの であり、自己共分散の総和とほぼ等しい。 *26正確には、有意水準次第で閾となる値は変わる。3.1
状態空間モデル
♠
状態空間モデルのうち最も基本となる線型ガウス型状態空間モデルは、t = 1, 2, . . . に 対して yt = Zαt+ ϵt, ϵt ∼ MNm(0, H), (34) αt+1= Dαt+ Rηt+1, ηt ∼ MNr(0, Q), (35) α1 ∼ MNr(a1, P1) (36) と表現することができる*27。 (34)式は観測方程式で、観測変数yt はm次元ベクトル、潜在変数αt はr次元ベクト ル、Z は(m, r)型の行列である。観測変数yt はデータとして得られるが、潜在変数αt はそのような形では得られないものと考える。(35)式は遷移方程式で、D, Rは(r, r)型 の行列である*28。誤差項ϵt, ηt とα1 が多変量正規分布に従うことから、yt, αt の条件付 分布は全て多変量正規分布に従う。3.2
カルマン・フィルタ
♠
パラメータϑ = {Z, H, D, R, Q, a1, P1}が既知のもとで、データYT ={yi}Ti=1が得ら れる確率は p(YT) = p(y1, y2, . . . , yT)= p(yT|YT−1)p(yT−1|YT−2)· · · p(y2|y1)p(y1)
(37) に比例する（全てϑの条件付であるが省略して書く）。このp(YT)はベイズ推定の手順に おける条件付分布py|ϑ(y|ϑ)の部分にあたる。p(yt|Yt−1)が上記の理由より多変量正規分 布に従うことから、条件付期待値E(yt|Yt−1)、条件付分散Var(yt|Yt−1)と実際のデータ の値YT = YTobs が与えられれば尤度l(ϑ; YTobs)が計算できる*29。すなわち、多変量正規 分布のパラメータは平均と分散の2つだけなので、その両方がわかれば確率密度関数を再 現できる（脚注*30参照）。これらは、カルマン・フィルタの漸化式と呼ばれる公式を用い *27ガウス型とは、正規分布の別名がガウス分布であることから、誤差項が全て正規分布に従うという意味で ある。 *28カルマン・フィルタの初期化(36)については、脚注*36も参照。α1はデータを用いて設定できる。ベ イズ統計学の考え方からすれば、α1も事前分布とみなしてそのパラメータa1, P1の値を外生的に決め てしまってもよい。 *29条件付期待値とは、大雑把に、何らかの情報が与えられた場合の確率変数の平均を意味し、条件付分散・ 共分散も同様である。条件付期待値は付け焼刃の知識で議論できるほど易しくはないのだが、本稿はあく まで統計学ではなくマクロ経済学を主題とすることから、細部には触れず重要な部分のみ説明する。
て求められることが知られている。まず、at = E(αt|Yt−1), Pt = Var(αt|Yt−1)と定義 すると、 at+1 = E(αt+1|Yt) = DE(αt|Yt) (38) Pt+1 = Var(αt+1|Yt) = DVar(αt|Yt)D⊤+ RQR⊤ (39) である。次に、 νt = yt− E(yt|Yt−1) = yt− Zat (40) Ft = Var(yt|Yt−1) = Var(νt|Yt−1) (41) Mt = Cov(αt, νt|Yt−1) (42) と定義すると（Mt はαt とνt の条件付共分散である）、 [ αt νt ] Yt−1 ∼ MNr+m ([ at 0 ] , [ Pt Mt Mt⊤ Ft ]) (43) で、{νt, Yt−1}とYt は情報として同等であること、およびこの節の補足の多変量正規分 布の条件付分布に関する定理より E(αt|Yt) = E(αt|νt, Yt−1) = at+ MtFt−1νt (44) Var(αt|Yt) = Var(αt|νt, Yt−1) = Pt − MtFt−1Mt⊤ (45) である。一方で、 Mt = E ( (αt− at)νt⊤|Yt−1 ) = E((αt− at)(yt− Zat)⊤|Yt−1 ) = E((αt− at)(Zαt+ ϵt− Zat)⊤|Yt−1 ) = E((αt− at)(αt− at)⊤Z⊤|Yt−1 ) = PtZ⊤ (46) Ft = Var(Zαt+ ϵt− Zat|Yt−1) = Var(Zαt+ ϵt|Yt−1) = ZPtZ⊤+ H (47) である。途中でαt とϵtが無相関であることを用いている。 結局、at, Pt が、t = 1, 2, . . . T に対してa1, P1から νt = yt− Zat (48) Ft = ZPtZ⊤+ H (49) Kt = DPtZ⊤Ft−1 (50) Lt = D− KtZ (51) at+1 = Dat + Ktνt (52) Pt+1 = DPtL⊤t + RQR⊤ (53)
という漸化式で計算できる。E(yt|Yt−1) = Zat, Var(yt|Yt−1) = Ft = ZPtZ⊤+ H なの で、いったん観測データとしてYT = YTobs が得られれば、この漸化式を使ってパラメー タϑの関数として尤度l(ϑ; YTobs)が計算できる。 補足：多変量正規分布の条件付分布に関する定理 k 次元ベクトルyが多変量正規分布MN(µ, Σ)に従うとする。y, µ, Σを y = [ y(1) y(2) ] (54) µ = [ µ(1) µ(2) ] (55) Σ = [ Σ(11) Σ(12) Σ(21) Σ(22) ] (56) と分割する。y(1) はk1次元、y(2) はk2 次元とする（k = k1 + k2）。µ, Σについても同 様である。このとき、y(2) = y(2)′ が与えられたときのy(1) の条件付分布は、多変量正規 分布MN(µ(1|2), Σ(1|2))に従う。ただし、 µ(1|2) = µ(1)+ Σ(12)Σ−1(22)(y(2)′ − µ(2)) (57) Σ(1|2) = Σ(11)− Σ(12)Σ−1(22)Σ⊤(12) (58) である*30。
3.3
DSGE
モデルのパラメータのベイズ推定の手順
前節までで、DSGEモデルのパラメータをベイズ推定するための方法が揃った。対数 線形近似したDSGEモデルの線型VAR表現を状態空間モデルの遷移方程式とみなすこ とは既に述べた。観測データとしてYT = YTobs が得られているとする。観測方程式につ いて、モデル変数 αt = [ ˆ xt ˆ st ] (59) との対応関係Z は通常外生的に与え、誤差項については一部の変数について恒等的にゼ ロと仮定し、その他については非対角成分をゼロとし対角成分のみを未知パラメータと する。H と(33)式のD、誤差項の分散Qおよび一般にRも推定しようとする未知パラ メータθ の関数であることから、明示的にH(θ), D(θ), Q(θ), R(θ) と書くことにする。2.3節の酔歩連鎖を利用することで、未知パラメータθ の事前分布pθ(θ)と状態空間モ デル yt = Zαt+ ϵt, ϵt ∼ MNm(0, H(θ)), αt+1= D(θ)αt+ R(θ)ηt+1, ηt ∼ MNr(0, Q(θ)), α1 ∼ MNr(a1, P1) (60) の尤度l(θ; YT)（カルマン・フィルタの漸化式を用いて計算する）から、(31)式の採択 確率の計算に含まれるf (θ)に尤度と事前分布の積l(θ; YT)pθ(θ)を用いることで、パラ メータθの事後分布pθ|y(θ|YT)∝ l(θ; YT)pθ(θ)を求めることができる。 結局、以下の手順により未知パラメータθの事後分布がシミュレートできる。 DSGEモデルの未知パラメータのベイズ推定の手順 初期化：θ(0) を設定 繰り返しi（i ≥ 1）： step 1. ˜θ← MN(θ(i−1), Σ)
step 2. Blanchard and Kahnの方法などによりD(θ), R(θ)を求める step 3. カルマン・フィルタの公式により尤度l( ˜θ; YTobs)を計算 step 4. 事前確率pθ( ˜θ)と採択確率p = min [ l( ˜θ;y)pθ( ˜θ) l(θ(i−1);y)pθ(θ(i−1)), 1 ] を計算 step 5. 確率pでθ˜を採択して、θ(i) = ˜θとおく それ以外はθ˜を棄却して、θ(i) = θ(i−1) とおく 上記については、いくつかの注意点がある。 • step. 1のΣは、事後分布のモード周りのヘッセ行列の逆行列を頼りに決めること が多い。 • 尤度は、実際には対数尤度を計算する。 • 繰り返しの回数に特に決まりはない。2.4節で説明したような収束判定方法を手が かりに、サンプリングが適切にできているかを判断する。 • チェインを複数用意し並列計算すると、計算の効率化や収束の判断に有用である。
• Blanchard and Kahnの方法ではR(θ)が求まらないので、必要に応じて他のアル
ゴリズム（例えばSimsの方法）を用いる。ただし、次節の例ではRはθに依存し
4
例：ニューケインジアン・モデルのパラメータのベイズ
推定
それでは、ノート5の対数線形近似したニューケインジアン・モデルと1980年∼1999 年までの日本のデータを用いてモデルのパラメータを上記の手順によりベイズ推定してみ よう。4.1
モデル、データ、事前分布
モデル 対数線型化したニューケインジアン・モデルは、以下のような6本の方程式 πt = βπt+1+ κˆxt (61) ˆ xt = ˆxt+1− (ˆit − πt+1) + νt (62) ˆit = ϕππt+ ϕyxˆt+ vt (63) vt+1= ρvvt+ ut+1 (64) ˆ at+1 = ρAˆat + εt+1 (65) νt = ˆat+1− ˆat (66) からなる*31。(65)式、(66)式よりν t = ˆat+1− ˆat = (ρA− 1)ˆat+ εt+1であるが、ここ ではt期の変数の値が決まる際にt + 1期の外生変数はゼロであると予期されているもの として、(62)式におけるνt は、(ρA− 1)ˆatで置き換える。また、ϕπ > 1, ϕy > 0が合理 的期待均衡解の存在条件となることから、予めパラメータϕπ をϕ˜π = ϕπ− 1と変数変換 しておく。割引因子βは推定すべきパラメータに含めず、β = 0.99とおく。 したがって、以下のような5本の方程式 πt = βπt+1+ κˆxt (67) ˆ xt = ˆxt+1− (ˆit − πt+1) + (ρA− 1)ˆat (68) ˆit = (1 + ˜ϕπ)πt+ ϕyxˆt+ vt (69) vt+1= ρvvt+ ut+1 (70) ˆ at+1 = ρAˆat + εt+1 (71) か ら な る ニ ュ ー ケ イ ン ジ ア ン・モ デ ル を 推 定 す る 。ut, εt は そ れ ぞ れ 独 立 に 平 均 ゼ ロ 、標 準 偏 差 σu, σε の 正 規 分 布 に 従 う も の と し 、推 定 す べ き パ ラ メ ー タ を θ′ = *31κ = (1−ω)(1−ωβ)(γ+1) ω である。{γ, ω, ˜ϕπ, ϕy, ρA, ρv, σu, σε}と書く*32。ここで、 αt = ˆ xt πt ˆit vt ˆ at , ηt = [ ut εt ] , R = 0 0 0 0 0 0 1 0 0 1 (72)
と定義する。このモデルは、あるθ′ が与えられると、例えばBlanchard and Kahnの方
法により、
αt+1= D(θ′)αt+ Rηt+1, ηt ∼ MN(0, Q(θ′)) (73)
と表現できる。
データ
データとして、GDPギャップ、GDPデフレーター（季節調整値、前期比）、短期金利の
3系列を用いる（図7）。それぞれxobst , πobst , iobst の記号を当てる。いずれも四半期データ
で、期間は政策金利がゼロ金利になる前の1980年第2四半期∼1999年第1四半期、観測
数は76である*33。
観測方程式は、
xobst − ¯xobs = ˆxt (74)
πtobs− ¯πobs = πt+ ϵπ,t (75)
(iobst − ¯iobs)/4 = ˆit (76)
とする*34。ただし、x¯obs, ¯πobs,¯iobsはそれぞれの系列の期間平均である。対数線型化した
際、各変数の定義を定常状態からの乖離にしたことに対応している。ϵπ,tは物価上昇率の 観測誤差を意味する。データ3系列に対して誤差項が3個以上ないと尤度が計算できなく なるため、観測誤差は1つ以上必要である。したがって、推定すべきパラメータは、ϵπ,t の標準偏差も含めてθ ={γ, ω, ˜ϕπ, ϕy, ρA, ρv, σu, σε, σϵπ}とする。(76)式の1/4は、年 利表示されている名目金利を四半期表示に変換するという意味である。 *32標準偏差は分散の平方根をとったものである。
*33データソースは、いずれもOECD Economic Outlook No. 90である。
*34このとき、(60)式におけるZは Z = 10 01 00 00 00 0 0 1 0 0 (77) と書ける。
ＧＤＰギャップ ％ 1980 1985 1990 1995 -2 0 2 4 ＧＤＰデフレーター（前期比伸び率） ％ 1980 1985 1990 1995 -0.5 0.5 1.0 1.5 2.0 2.5 短期金利 ％ 1980 1985 1990 1995 2 4 6 8 10 12 図7 データ 事前分布 θ の要素それぞれに対して、表2のように事前分布を設定した*35。労働供給の弾力性 γ、金融政策ルールの係数ϕ˜π, ϕy はとりうる値がゼロより大きい実数のため、ガンマ分布 を確率分布に選んだ。価格改定できない確率ωは0 < ω < 1のいずれかの値をとること *35 なお、対数線形化した段階で、元の非線形のモデルにあった労働の不効用の相対ウエイトµ、需要の価格 弾性値ηは識別不可能になる。
パラメータ 分布の種類 平均 標準偏差 γ ガンマ分布 1.0 0.5 労働供給の弾力性 ω ベータ分布 0.8 0.1 価格改定できない確率 ˜ ϕπ ガンマ分布 0.5 0.25 金融政策ルールの係数（物価） ϕy ガンマ分布 0.5 0.25 金融政策ルールの係数（GDPギャップ） ρA ベータ分布 0.8 0.05 AR(1)項の係数（技術ショック） ρv ベータ分布 0.8 0.1 AR(1)項の係数（金融政策ショック） σε 逆ガンマ分布 0.5 0.5 εtの標準偏差 σu 逆ガンマ分布 0.5 0.5 utの標準偏差 σϵπ 逆ガンマ分布 0.5 0.5 ϵπ,t の標準偏差 表2 パラメータの事前分布 から、ベータ分布を確率分布に選んだ。AR(1)項の係数は絶対値が1未満でなければな らないが、マイナスの値はとりえないとして、同じくベータ分布を確率分布に選んだ。誤 差項の標準偏差σε, σu, σϵπ は逆数がガンマ分布に従うと仮定した（これが逆ガンマ分布 の定義である）*36。
4.2
事後分布、歴史分解
上記のようなセッティングの元で、3.3節の方法によりニューケインジアン・モデルの パラメータの事後分布を求めた。チェインの本数は2、サンプリング回数は各チェイン 125,000で、うち最初の25,000は定常状態に至っていないものとみなして捨てた。した がって、200,000サンプルを用いて、事後分布の平均と標準偏差を計算した。M-Hアルゴ リズムでの採択確率は、それぞれ0.304, 0.306であった。 事後分布の平均、標準偏差とGeweke の収束判定統計量を表3に示した。Z1, Z2 はそ れぞれ1番目、2番目のチェインのGewekeの収束判定統計量である。いくつかの変数に ついて、1番目のチェインのZ の絶対値がやや大きいが、2番目のチェインのZ の絶対 値は小さいこと、およびチェインの軌跡と分布のヒストグラム（掲載せず）から、事後分 布が適切に推定できていると判断した。 事前分布と事後分布の平均を比較すると、労働供給の弾力性のパラメータγは事後平均 のほうが大きい。これは、家計の労働供給が事前に予想したのと比べてより非弾力的であ*36 推定はDynare（バージョン4.3.1）を用いて行った。Dynareとは、数値計算ソフトであるOctaveも
しくはMATLAB上で動く、DSGEモデルのシミュレーション・パラメータ推定のためのツールであ
る。Dynareはカルマン・フィルタの共分散行列の初期化をリアプノフ方程式を解くことによって行う。
なお、上記においてモデルを例えばBlanchard and Kahnの方法で(73)式のように表現できると書い
パラメータ 平均 標準偏差 Z1 Z2 γ 2.271 0.511 −2.509 1.004 ω 0.348 0.091 −0.222 0.369 ˜ ϕπ 0.430 0.037 0.195 0.834 ϕy 0.875 0.455 −3.268 0.043 ρA 0.926 0.013 −3.648 0.344 ρv 0.394 0.191 −0.776 0.391 σu 0.229 0.095 −0.679 0.452 σε 0.958 0.010 0.411 -0.381 σϵπ 0.879 0.016 −0.186 0.415 表3 パラメータの事後分布とGewekeの収束判定統計量 ることを意味する。価格改定できない確率ωは事後平均のほうが小さく、事前に予想し たよりも物価は弾力的である。金融政策ルールについて、物価の係数ϕ˜π は事前と事後で 平均がそれほど変わらないが、GDPギャップの係数ϕy は事後平均のほうが大きく、中 央銀行はGDPギャップの変動をより重要視して政策金利を上下させていることわかる。 AR(1)項の係数の事後平均は、技術ショックに関するρAは大きいが、金融政策ショック に関するρv は比較的小さい。これは、技術ショックの影響がより持続的であることを意 味する。この差は、次に説明する歴史分解の結果に現れている。 (73)式より、 αt = Dt−1α1+ t−1 ∑ i=1 Dt−i−1Rηi+1 (78) と表すことができるので、あるパラメータθ に対してαt はα1 とηt の積み上げにより 分解できる。yt = Zαt + ϵt なので、yt はα1 とηt の積み上げと観測誤差ϵt により分 解できる。そのように観測方程式(74)∼(76)の左辺を分解したのが図8の歴史分解であ る*37。 まず、GDPギャップについてみると90年前後のバブル景気は、その多くの部分が金 融政策ショックにより説明されている。その後については、技術ショックが常に負の影響 を与えているが、97年以降については金融政策ショックも負の影響を与えている。GDP デフレータの伸び率については、観測誤差としてノイズ扱いされた部分を除けば、GDP *37正確には、事後分布からリサンプリングした各θ(i) について分解を計算し、その平均を使って作図して いる。期間平均からの乖離を分解しているので、例えば短期金利の線グラフはマイナスの値もとる。
1.5 1 0.5 0 0.5 1 1.5 2 2.5 3 1980 81 82 83 84 85 86 87 88 89 1990 91 92 93 94 95 96 97 98 99 % GDP!"#!"!#$%&'()* +%, -./0 12345678 !"#$% GDP&'()*) 6 4 2 0 2 4 6 1980 81 82 83 84 85 86 87 88 89 1990 91 92 93 94 95 96 97 98 99+ % GDP !"# $%& '()*+,"- ./+,"-GDP !"# 1.5 1 0.5 0 0.5 1 1.5 2 2.5 1980 81 82 83 84 85 86 87 88 89 1990 91 92 93 94 95 96 97 98 99$ % !"# $!% "&'()*+, -.)*+, !"# 図8 歴史分解 ギャップと比例した動きをしていることがわかる。短期金利のほとんどは、技術ショック と初期値によって説明されている。金利の水準は趨勢的に下がっているが、今回の推定結 果では、80年代初頭の金利の高さは初期値とみなされ、それが徐々に解消されていくとく 結果になった。
