403 頁∼ 424 頁
非対称
t
分布を用いた多変量確率的ボラティリティモデルの
推定
中島 上智
∗Estimation of Multivariate Stochastic Volatility Model with Skew-t
Distribution
Jouchi Nakajima∗
複数の株価や為替レート等の資産価格収益率について, 多変量のボラティリティの変動を捉え るモデルとして, 多変量確率的ボラティリティ (Multivariate Stochastic Volatility, MSV) モデル があり, 近年, 様々な拡張が試みられている. 資産価格の収益率分布は, 状況によって左右非対称 となることが知られており, 本稿では, MSV モデルの誤差項に非対称 t 分布を適用したモデルと そのベイズ推定法について解説する. また, 分布の非対称性の程度を表すパラメータに縮小推定 法を適用する方法についても紹介する. 米国株式市場 S&P500 の日次収益率データを用いた実証 分析によると, 非対称 t 分布を用い, かつ縮小推定法を適用すると, MSV モデルの予測精度が向 上することが示される.
Multivariate stochastic volatility (MSV) models and their extensions have been widely de-veloped to capture multivariate volatility of returns of asset prices, such as stock prices and foreign exchange rates. As a return distribution of asset prices tends to be skewed, this pa-per introduces the MSV models with a skew-t distribution. A shrinkage method for skewness parameters is also introduced. An empirical analysis using daily stock returns in the S&P500 is provided to show that the skew-t distribution and the shrinkage method improve predictive ability of the MSV model.
キーワード: 多変量確率的ボラティリティモデル, 非対称 t 分布, マルコフ連鎖モンテカルロ法, 株価収益率
1. はじめに
金融市場で観察される株価や為替レート等の資産価格収益率のボラティリティは, 時々刻々 と変化しているという考え方が一般的であり, その変動を捉えるモデルとして, ARCH モデ ル (Engle (1982), Bollerslev (1986)) 等と並び, 確率的ボラティリティ (Stochastic Volatility,
∗ 日本銀行調査統計局：〒 103-0021 東京都中央区日本橋本石町 2-1-1 (E-mail: [email protected]).
SV) モデルがよく用いられている. SV モデルは, ファイナンス理論との親和性や予測精度 の良さから, モデルの拡張や実証分析が幅広く行われている (例えば, Black (1976), Ghysels et al. (2002), Shephard (2005), 渡部 (2000) を参照). その発展形として, 複数の資産価格 の収益率について, その相関とボラティリティを同時に定式化したものとして, 多変量 SV (Multivariate SV, MSV) モデルがある (例えば, Chib et al. (2009) を参照). この多変量 モデルについても様々な拡張が試みられているが, そのほとんどが, 誤差項に正規分布やス チューデントの t 分布等, 左右対称な分布を仮定している. もっとも, 資産価格の収益率分 布は, 状況によって左右非対称となることが知られており, 特にリスク管理等の観点からは モデル化に際して重要な側面である. Nakajima (2017, 2020) は, MSV モデルに非対称な 誤差分布を導入したモデルとその推定方法を提案しており, 本稿はそれらについて解説を 行う. MSV モデルには, 変数間の相関関係をどのように表現するかによって, いくつかの定式 化の方法がある. その中で, 効率的なベイズ推定方法が可能であるという点から, 2 種類の 定式化が有効であるとされている. 1 つは共分散行列をコレスキー分解の要領で定式化す るコレスキー MSV モデル (Asai et al. (2006), Gouri´eroux et al. (2009) 等) であり, もう 1 つは共通変動部分をファクターとして抽出するファクター SV モデル (Aguilar and West (2000), Chib et al. (2006), Lopes and Carvalho (2007) 等) である. Nakajima (2017) は前 者, Nakajima (2020) は後者の定式化を用いて, それぞれ非対称 t 分布付きの MSV モデル を提案している. 非対称性の程度を表すパラメータについて, k 個の資産価格の収益率分布を表現するため に, k 個の非対称性パラメータを設定することが可能であるが, その全てが必要かどうかは 実証上の問題となる. 例えば, 資産価格間で相関が高く, 分布の非対称性についても同様の 傾向がみられるのであれば, k 個より少ない非対称性パラメータのセットで k 個の資産価格 の収益率分布を十分に表現できる可能性がある. 特に, 分析対象の資産価格の系列数が非常 に多い場合に, 非対称性のパラメータの数を減らすことができれば, 推定の効率性が向上し, モデルの予測精度が高まることが期待される. Nakajima (2017, 2020) は, ベイズ推定法の 枠組みにおいて変数選択のために用いられるスパース事前分布 (sparsity prior; George and McCulloch (1993, 1997), Clyde and George (2004) 等) を, 非対称性パラメータに適用する ことにより, 必要のない非対称性パラメータを探索する方法を提案している. 本稿でも, こ の方法について解説を行い, 米国の業種別の日次収益率データを用いて, このスパース事前 分布が実証分析上, モデルの予測精度の向上につながることを示す. 本稿の以下の構成は次の通りである. まず, 第 2 節で MSV モデルを解説し, 第 3 節でそ の推定方法について紹介する. 第 4 節で S&P500 の業種別収益率データを用いた実証分析 を行い, 第 5 節で今後の課題について述べる. 補論では, 推定方法のアルゴリズムの詳細部
分を解説する.
2. 非対称 t 分布と MSV モデル
2.1 一般化双曲型非対称 t 分布
左右非対称で裾が正規分布より厚い収益率分布を表現できる分布形として, 一般化双曲型 (Generalized Hyperbolic, GH) 分布がしばしば用いられる (Barndorff-Nielsen (1977), 増 田 (2002) を参照). GH 分布は, 正規分布や t 分布, 逆ガウス分布等を包含する分布族であ り, 様々な裾の厚さおよび左右非対称性を表現することができるため, しばしば金融データ の分析に用いられてきた. もっとも, Prause (1999), Aas and Haff (2006), Nakajima and Omori (2012) 等で議論されているように, GH 分布は非常にフレキシブルな分布であり, GH 分布のパラメータのうちいくつかは識別が難しい場合がある. 特に, ボラティリティが 潜在変数となっている SV モデルに, GH 分布をそのまま適応することは難しいと考えられ る. そこで, 本稿では, Nakajima and Omori (2012) に倣って, GH 分布族の特別なケース であり, いくつかのパラメータ空間を限定することによって得られる, GH 非対称 t 分布を 用いる.
GH 分布は一般に, 正規尺度平均混合 (normal variance-mean mixture) で表現される. こ
のことを利用して, GH 非対称 t 分布に従う確率変数 witは次式で定義される. wit = µi+ βizit+ √ zitεit, εit ∼ N(0, 1), zit ∼ IG(νi/2, νi/2). (2.1) ここで, IG は逆ガンマ (Inverse Gamma) 分布である. さらに, SV モデルに組み込む場合 には, E[wit] = 0 と仮定し, µi =−βiciとおく. ただし, ci = E[zit] = νi/(νi− 2) である. また, witの分散が有限であると仮定するため, νi> 4 とする. 中島・大森 (2011) が示すよ うに, βiは分布の左右の非対称性を規定するパラメータであり, νiを一定とした場合, βiの 絶対値が大きくなるほど, 分布の非対称性 (歪み) が大きくなる. 一方, νiはスチューデント の t 分布の自由度に対応するパラメータであり, 分布の裾の厚さを規定する. βiを一定とし た場合, νiが小さくなるほど裾の厚さが大きくなる. 2.2 MSV モデル k× 1 の変数ベクトルを yt = (y1t, . . . , ykt)0とおき, yt ∼ N(0, Σt) となる時系列過程 (t = 1, 2, . . .) を考える. 標準的な MSV モデルでは, ytが多変量正規分布にしたがってお り, 分散共分散行列 Σtが時刻 t によって変動すると仮定する. 1 節で述べた 2 種類の MSV モデルのうち, 1 つ目のコレスキー MSV モデルは, Σtを次
式のようにコレスキー分解する. Σt = A−1Λ2t(A0)−1. ただし, A は次式のような下三角行列, Λtは対角行列である. A = 1 0 · · · 0 −a21 . .. . .. ... .. . . .. . .. 0 −ak1 · · · −ak,k−1 1 , Λt= λ1t 0 · · · 0 0 . .. . .. ... .. . . .. . .. 0 0 · · · 0 λkt . ytについては, 次式のように書ける. yt= A−1Λtet. (2.2) ただし, et∼ N(0, I) である. ここで, A は (y1t, . . . , ykt) の変数間の相関を捉える行列であ り, Λtはボラティリティの変動, すなわち確率的ボラティリティを捉える行列である.
こうしたコレスキー分解を用いた分散共分散行列の分解は, Pinheiro and Bates (1996), Pourahmadi (1999), Smith and Kohn (2002), George et al. (2008) 等で使われているほか, MSV モデルとしては, Cogley and Sargent (2005), Primiceri (2005), Lopes et al. (2012) 等で扱われている. Nakajima (2017) では, 行列 A が時間を通じて変動するように, Atと 定式化しているが, 本稿では簡単化のために時刻 t に依存しないと仮定する. Nakajima (2017) は, 式 (2.2) の誤差項 etを, GH 非対称 t 分布に従う誤差項に置き換え ることにより, 非対称 t 分布付きの MSV モデルを提案している. 具体的には, 式 (2.1) の確 率変数 witを用いて, 誤差項ベクトルを wt= (w1t, . . . , wkt)0と定義し, 式 (2.2) を次式に拡 張する. yt= A−1Λtwt. ここでの留意点として, A−1は下三角行列であるため, 1 番目の変数である y1tは GH 非対 称 t 分布にしたがっているが, 2 番目以降の (y2t, . . . , ykt) は, それ自体は必ずしも GH 非対 称 t 分布にしたがうわけではない. これは, GH 非対称 t 分布にしたがう確率変数の和は必 ずしも GH 非対称 t 分布とならないためである. Λtの対角成分について SV モデルの定式化を行うために, hit= log(λ2it) と変換し, k× 1 のボラティリティのベクトルを ht= (h1t, . . . , hkt)0と定義する. このボラティリティの変
動については, 次式のように定式化する. ht+1 = µ + Φ(ht− µ) + ηt. (2.3) ただし, Φ は対角行列であり, i 番目の対角成分を φiと書く. 式 (2.3) が定常過程となるよ うに, |φi| < 1 と仮定する. また, µ = (µ1, . . . , µk)0, ηt= (η1t, . . . , ηkt)0と定義する. さら に, 式 (2.1) の誤差項 εitを並べたベクトルを εt= (ε1t, . . . , εkt)0と定義し, εtと ηtの分布 を次式のように仮定する. εt ηt ∼ N 0, I Q Q S . ただし, S と Q は対角行列であり, 対角成分をそれぞれ σ2 i, ρiσiと書く. ρiは観測方程式 の誤差項 εitとボラティリティの AR(1) 過程の誤差項 ηitの相関係数である. 株価の収益率 を分析する際には, この相関係数がレバレッジ効果 (leverage effect) を捉える. これは, 収 益率が低下した翌日には, 上昇した翌日よりも, ボラティリティの上昇を引き起こすことが 多いという非対称性を捉えている. このレバレッジ効果が生じる場合, 相関係数 ρiは負の
値となる (Yu (2005), Omori et al. (2007)).
以上の定式化から, ˜yt= Aytと定義したベクトル ˜yt= (˜y1t, . . . , ˜ykt)0の要素 ˜yitについ て, 次式が得られる. ˜ yit = {βi(zit− ci) +√zitεit} exp(hit/2), hi,t+1 = µi+ φi(hit− µi) + ηit, εit ηit ∼ N(0, Ωi), Ωi= 1 ρiσi ρiσi σ2i . すなわち, Φ, S, Q を対角行列と仮定したために, MSV モデルは, k 個の単変量の SV モデ ルに分解できる (Chib et al. (2006), Asai et al. (2006), Gouri´eroux et al. (2009) 等).
コレスキー MSV モデルを使う場合に留意しなければならないのは, 推定結果が ytの中
の変数の順序に依存することである. このため, 実証分析では ytの中の変数をどのような
順序に並べるかに気を付けなければならない. Nakajima and Watanabe (2011) は変数の 順序に不確実性があると仮定し, リバーシブルジャンプ MCMC 法を用いて, モデルのパラ メータの推定と同時に変数の順序も推定する方法を提案している. 簡易的な方法について は, 後述のファクター MSV モデルの実証分析および Nakajima (2017) を参照されたい.
次に, MSV モデルの 2 つ目の定式化の方法である, ファクター MSV モデルを紹介する. k 個の変量よりも少ない数 (p 個) の共通因子 (ファクター) を ft= (f1t, . . . , fpt)0と定義し,
ytが次式のように ftの線形結合と誤差項で表されると仮定する. yt = Bft+ et, et∼ N(0, Λt), (2.4) ft ∼ N(0, Ψt). (2.5) ただし, B は k×p のファクター係数 (factor loading) 行列, Λtは上述と同じ対角成分 λ2itをも つ対角行列と仮定する. Ψtも対角行列と仮定し, その対角成分を ψ2jtと書く. yt∼ N(0, Σt) に関して, k× k の分散共分散行列 Σtは次式のように分解される. Σt= BΨtB0+ Λt. (2.6) ここで Λtと Ψtの対角成分に, SV モデルを導入する. hit= log λ2it, (i = 1, . . . , k), およ び, hk+j,t= log ψ2jt, (j = 1, . . . , p) と定義し, 全ての h`t (` = 1, . . . , q; q = k + p) が単変量 の SV モデルに従うと仮定する. Nakajima (2020) は, etと ftが GH 非対称 t 分布にした がうと考え, 以下のような定式化を提案した. yt = Bft+ Λ1/2t wt, wt= (w1t, . . . wkt)0, ft = Ψ1/2t vt, vt= (wk+1,t, . . . wk+p,t)0, wit = βi(zit− ci) +√zitit, i = 1, . . . , q. ここで, ht= (h1t, . . . , hqt)0, t= (1t, . . . , qt)0と定義して, SV モデルを次式で定義する. ht+1 = µ + Φ(ht− µ) + ηt. ただし, µ = (µ1, . . . , µq)0, Φ = diag(φ1, . . . , φq), ηt= (η1t, . . . , ηqt)0であり, t ηt ∼ N 0, I R R W , と仮定する. W と R は, それぞれ σ2 i と ρiσiを対角成分とする対角行列である. コレス キー MSV モデルと同様に, ファクター MSV モデルも q 個の単変量 SV モデルに分解でき る. このことを確かめるために, yitと fitについて, 次のような変換を行う. ˜ yit = yit− B(i)ft (i≤ k), fjt (i > k; j = i− k), i = 1, . . . , q. ただし, B(i)は, 1× p のベクトルで, 行列 B の i 番目の行を表す. Λt, Ψt, W , R がそれぞ
れ対角行列であることから, 全ての i = 1, . . . , q について, 次式のように書ける. ˜ yit = {βi(zit− ci) +√zitit} exp(hit/2), hi,t+1 = µi+ φi(hit− µi) + ηit, it ηit ∼ N(0, Ωi), Ωi = 1 ρiσi ρiσi σi2 . 式 (2.4)–(2.6) から分かるように, B も ftも推定すべきパラメータや潜在変数であるた め, モデルを識別するためには何かしらの制約を置く必要がある. 文献の中では, いくつか の方法が提案されているが, ここでは, Aguilar and West (2000), Lopes and West (2004) にしたがって, B の一部のパラメータに制約を置く. 具体的には, B の最初の p 行の正方行 列部分が対角成分を 1 とする下三角行列になるようにする. すなわち, i = 1, . . . , p につい て, bii = 1 とし, また, i < p かつ i < j の行列の上三角の部分について, bij = 0 と仮定す る. この方法は B に制約を置くことで, それ以外の部分に制約を置かずに推定できるため,
簡単な方法であるが, 推定結果が ytの最初の p 個の変数の順序に依存するという留意点が
ある (この点について, 例えば, Fr¨uwirth-Schnatter and Lopes (2018) 等を参照).
2.3 Skew パラメータに対するスパース事前分布 非対称 t 分布付きの MSV モデルには, 非対称性のパラメータ βiが, コレスキー MSV モ デルでは k 個, ファクター MSV モデルでは k + 2 個, それぞれ含まれることになる. ytの 変数の個数, つまり k が非常に大きい場合, 非対称性のパラメータの数もそれに応じて多く なる. もっとも, 変数間の相関を勘案すれば, 全ての非対称性のパラメータが必要でないと いう可能性もある. この点について, 以下, シミュレーションにより確認する. 簡易的に, k = 3 の場合を考え, SV モデルの先行研究における実証分析の結果から, 次の ようなパラメータの値を用いる. SV モデルのパラメータについて, φi = 0.995, σi= 0.05, ρi =−0.5, µi = −11, νi = 8 とし, A の各要素は一様分布 U [−2, 0] から, B の各要素は U [1, 2] からその都度, 発生させる. βiについては, 後述のように設定する. コレスキー MSV モデルについては, ファクターの数を p = 1 とし, µ4=−10 と仮定する. µ4を (µ1, µ2, µ3) よりも高く設定しているのは, 実証分析で, ファクターのボラティリティの平均的な大きさ が, 各変数固有のショックのそれよりも高くなることが多いためである. T をサンプルサイ ズとして, t = 1, . . . , T の順に htおよび ytを発生させ, yi,1:T = (yi1, . . . , yiT) の歪度を計 算する. ここでは, T = 1000 とする. この計算を 1000 回繰り返す. 図 1(i) は, コレスキー MSV モデルについて, yitの歪度を 1000 回のシミュレーション ごとに計算し, 分布で表したとき, 平均 (丸印), 25%点から 75%点までの 50%範囲 (太い縦
(i)コレスキーMSVモデル (ii)ファクターMSVモデル 図 1 シミュレーションで得られた ytの歪度. 線) と, 10%点から 90%点までの 80%範囲 (細い縦線) をグラフにしたものである. まず, (a) β = (−1, −1, −1) と設定した場合, すなわち, 全ての非対称性パラメータが負である場合 は, 全ての yitの歪度が相応に負となる. 次に, (b) β = (0,−1, −1) とした場合, y2t, y3tの 式の誤差項の非対称性パラメータは負であるが, y1tの非対称性パラメータはゼロである ため, y2tと y3tの歪度が負となり, y1tの歪度はほぼゼロを中心に分布している. 一方, (c) β = (−1, 0, 0) とした場合, 最初の変数である y1tの歪度が負となり, それが変数間の相関を 通じて y2tと y3tに影響するため, y2tと y3tの誤差項の非対称性パラメータがゼロであって も, 歪度は負となることがわかる. 図 1(ii) は, ファクター MSV モデルについて, 同様のシミュレーションを行った結果である. (a) β = (−1, −1, −1, −1) の場合は, 全ての yitの歪度が負となる. (b) β = (−1, −1, −1, 0) のとき, すなわち, 個別の yitの誤差項の非対称性パラメータは負だが, ファクターの非対称性 パラメータがゼロの場合, yitの歪度は (a) の場合ほど負ではなくなる. これは, ファクターの ボラティリティが個別の yitの誤差項のそれよりも大きく, 左右対称な分布をもつファクター
の影響が左右非対称な個別の誤差項の影響に勝るためである. 一方, (c) β = (0, 0, 0,−1), すなわち, ファクターの非対称性パラメータのみ負にすると, 全ての yitの歪度は負となる. これらのシミュレーション結果から, 株価収益率にみられるような負の歪度を表現する ために, 全ての非対称性のパラメータが負であることは必要でない可能性がある. そこで, Nakajima (2017, 2020) は, ベイズ推定方法における縮小推定の手法の１つである, スパー ス事前分布 (sparsity prior) を非対称性パラメータに適用することを提案した. 代表的なス パース事前分布は, spike-and-slab 事前分布とも呼ばれ, 次式で定義される. βi ∼ κN(βi|0, τ02) + (1− κ)δ0(βi), 0 < κ < 1. ただし, δ0はディラック (Dirac) のデルタであり, βi= 0 で密度 1 をもつ. この事前分布は, 確率 κ で正規分布 (“slab”) にしたがい, 確率 1− κ でゼロの値 (“spike”) をとる. この縮小 推定の方法は, 様々なベイズ推定の研究で用いられている (例えば, George and McCulloch (1993, 1997), Clyde and George (2004) を参照).
3. ベイズ推定とマルコフ連鎖モンテカルロ法
SV モデルは潜在変数の個数が非常に多いため, 最尤法によるパラメータの推定が困難と なる場合が多い. このため, 先行研究では, SV モデルをベイズ推定の枠組みにおけるマル コフ連鎖モンテカルロ (Markov chain Monte Carlo; MCMC) 法を用いて, 推定すること が提案されている. 具体的には, モデルの事後分布からパラメータの確率標本を乱数発生 させ, 得られた標本をもとに統計的推測を行う. (MCMC について詳しくは, Koop (2003), Gamerman and Lopes (2006), 大森 (2001), 中妻 (2003), 伊庭他 (2001), 和合 (2005), 小 西他 (2008) 等を参照). 以下では, Nakajima (2017, 2020) で提案されている, MSV モデル の MCMC アルゴリズムを概説する. なお, Nakajima (2017) では, コレスキー MSV モデ ルの行列 A の各要素が AR(1) にしたがっていると仮定しているが, 本稿では簡単化のた め A の各要素が時間に依存しない形で定式化している. このため, MCMC アルゴリズムも Nakajima (2017) とはわずかに異なる. データとして, T 時点までの観測値, y1:T ={y1, . . . , yT} が得られた場合, 次の MCMC アルゴリズムを計算することにより, 全てのパラメータと潜在変数の同時事後確率密度関 数からの標本を得ることができる. コレスキー MSV モデル A1. 全てのパラメータと潜在変数の初期値を設定する.
A2. ボラティリティ hi,1:T と潜在変数 zi,1:T をサンプリングする (i = 1, . . . , k). A3. SV モデルのパラメータ θi={µi, φi, σi, ρi, νi} をサンプリングする (i = 1, . . . , k). A4. 非対称性パラメータ β とスパース事前分布のパラメータ κ をサンプリングする. A5. 変数間の相関を表すパラメータ A をサンプリングする. A6. A2 に戻る. ファクター MSV モデル B1. 全てのパラメータと潜在変数の初期値を設定する. B2. ボラティリティ hi,1:T と潜在変数 zi,1:T をサンプリングする (i = 1, . . . , q). B3. SV モデルのパラメータ θiをサンプリングする (i = 1, . . . , q). B4. 非対称性パラメータ β とスパース事前分布のパラメータ κ をサンプリングする. B5. ファクター ftと係数行列 B をサンプリングする. B6. B2 に戻る. 各ステップは, 各パラメータや潜在変数の条件付事後分布からのサンプリングである. 4 番目までは両方のモデルとも同様の手順で, 5 番目のサンプリングのみ異なる. 各ステップ の詳細は補論にまとめた. SV モデルのボラティリティ hi,1:T のサンプリングについて付言 すると, 他のパラメータと潜在変数を条件付ければ, 前述のとおり, i ごとに hi,1:T が独立 となり, i ごとにサンプリングを行うことができる. SV モデルのボラティリティの条件付 事後分布からのサンプリング法としては, 主に single-move sampler, multi-move sampler, mixture sampler の 3 つがしばしば使われている. multi-move sampler は Shephard and Pitt (1997) および Watanabe and Omori (2004), Omori and Watanabe (2008) によって, mixture sampler は Kim et al. (1998) および Omori et al. (2007) によってそれぞれ提案さ れた方法で, single-move sampler よりも効率的であることが知られている (詳しくは, 渡部 (2000), 大森・渡部 (2008) 等を参照). 非対称 t 分布付きの SV モデルについては, Nakajima and Omori (2012) が用いた multi-move sampler が有用あり, 効率的にボラティリティをサ ンプリングできることが知られている.
表 1 S&P500 業種別日次収益率の要約統計量 (2014/1/2 – 2017/12/29, T = 1007). 業種 英語名 平均 標準偏差 歪度 尖度 1 素材 Materials 0.0003 0.0098 −0.3296 4.7151 2 情報技術 Information Technology 0.0006 0.0095 −0.3699 5.7842 3 エネルギー Energy −0.0002 0.0129 −0.1023 5.1005 4 生活必需品 Consumer Staples 0.0003 0.0069 −0.2700 4.9068 4. S&P500 業種別インデックスの日次収益率の実証分析 本節では, 米国の株価指数 S&P500(スタンダード・プアーズ 500 種) における業種別イ ンデックスの日次収益率を用いた実証分析を行う. この日次収益率のデータを用いて, ファ クター MSV モデルを例に推定を行い, 非対称 t 分布がない場合とある場合で予測精度の比 較を行う. 4.1 データと事前分布 推定に使用した業種は表 1 の 4 つである (k = 4; i = 1, . . . , 4). 収益率 yitは, 終値を Pit として対数階差を yit= log Pit− log Pi,t−1と計算した. 推定に使用した期間は 2014 年 1 月 2 日から 2017 年 12 月 29 日までの 1,007 営業日である. 時系列データは図 2 に描かれて おり, 要約統計量は表 1 のとおりである. 全ての系列で歪度は−0.4∼−0.1 と負の値をとっ ており, 尖度は 4∼6 と, 正規分布の 3 と比べて大きい. なお, 推計には収益率 yitから各系 列のサンプル平均を引いたものを用いた. 推定するモデルの種類は, ファクター MSV モデルの次の 3 つである. • F-0: 左右対称な t 分布 (βi= 0) • F-S: 非対称 t 分布 (スパース事前分布はなし, κ = 1) • F-SS: 非対称 t 分布 (スパース事前分布あり) ファクターの数については, p = 1 とした. これは, データに単純な主成分分析を行うと, 1 つ目のファクターで系列の変動の 86%を説明することができる一方, 2 つ目のファクター は追加的に 6%の変動しか説明しないため, 1 つのファクターだけで十分であると判断した. 前述のとおり, ytの中の変数の順番は推定結果に影響を与えるため, その選択には注意 が必要である. 確固たる方法はないが, Nakajima (2017, 2020) では, 変数間の相関係数行 列を計算し, 各列の合計値が高い順, つまり, 他の変数との相関が平均的に高い変数から 並べる方法を用いている. 本稿でも同様の方法で順番を決め, 表 1 に掲載されている順を (y1t, . . . , y4t) とした.
図 2 S&P500 業種別日次収益率の時系列データ (2014/1/2 – 2017/12/29). 推計にあたって, パラメータの事前分布を次のように設定した. (φi+ 1)/2∼ B(20, 1.5), µi ∼ N(−11, 1), σ−2i ∼ G(20, 0.01), (ρi + 1)/2 ∼ B(1, 1), νi ∼ G(16, 0.8)I[νi > 4], βi∼ κN(βi|0, 10) + (1 − κ)δ0(βi), κ∼ B(2, 2), aj ∼ N(0, 10), bj ∼ N(0, 10). ただし, 分 布の B と G はそれぞれベータ分布, ガンマ分布を表す. MCMC 法による推定においては, 稼動検査期間 (burn-in period) として最初の 2,000 個 を捨てた後, 20,000 個のサンプルを発生させた. 本稿における数値計算は全て行列言語 Ox 6.0 (Doornik (2006)) を用いている. サンプリングの結果を確認すると, 標本自己相関関 数はいずれのパラメータについても十分に減衰していて, MCMC アルゴリズムの定常分布 への収束は早く, また標本経路も状態空間を万遍なく十分に訪れている. なお, サンプリン グのいくつかのステップで, Metropolis-Hasting (MH) アルゴリズムを用いているが, その MH アルゴリズムの採択確率は 80∼99%と高く, 効率的なサンプリングとなっていること がわかる. 4.2 モデルの推定結果 図 3 は, モデル F-SS のそれぞれについて, 推定されたパラメータの事後分布を, 事後平 均 (丸印) と 50%(太い縦線), 80%信用区間 (細い縦線) で示したものである. 横軸はファク ター (F1) と各変数 (Y1,. . .,Y4) を表す. φiの事後平均は全体的に 0.995 以上と 1 に近く, ボラティリティの持続性が観察される. ρiの事後平均はファクターが−0.75 と推定されて
図 3 各パラメータの事後分布 (F-SS モデル). 図 4 パラメータ βiの事後分布. おり, レバレッジ効果が生じていることがわかる. νiの事後平均は 5∼20 と低めであり, 誤 差項が裾の厚い分布となっていることがわかる. 図 4 は, 各モデルの βi の事後分布を示したものである. モデル F-S の推定値をみると, ファクター (F1) がゼロから離れた負の値をとっており, MSV モデルの誤差分布の非対称 性が現れている. 一方, それ以外の βiについては, ほとんどゼロのまわりを分布していると いえる. モデル F-SS の βiの事後分布をみると, モデル F-S でゼロまわりに分布していた
図 5 ボラティリティの事後平均 (F-SS モデル). βiはスパース事前分布の影響により, ゼロに縮小推定されていることがわかる. 図 5 は, ボラティリティの事後平均を図示したものである. ファクター (F1) のボラティ リティが比較的大きく, 変動も大きいことがわかる. 1, 2, 4 番目の変数のボラティリティは ファクターに比べて小さめである一方, 3 番目の変数のボラティリティは相対的に大きく, サンプルの後半ではファクターのボラティリティより大きくなっているのが特徴的である. 3 番目の変数は 4 つの業種別収益率の共通成分とは異なった動きをしていることが多く, 特 にサンプルの後半は固有の変動の方が共通成分より大きくなっていたことを示唆している. 4.3 予測精度の比較 モデルの予測精度を計測するために, 1 日先から 5 日先までの out-of-sample 予測の計算 を行う. 具体的には, まず, 2014 年 1 月 2 日から 2017 年 5 月 26 日までのデータ (T = 857) を用いて各モデルを推計し, 1 日先から 5 日先までの日次収益率を予測する. 次に, 5 日間 のデータを加えて, モデルを再推計し, 1∼5 日先までの日次収益率を予測する. これを 30 回繰り返し, 150 日分の予測が得られるまで続ける. 予測精度としては, 次式の対数予測密 度比 (Log predictive density ratio) を用いる.
LPDRt(d) = log { pM1(yt+d|y1:t) pM0(yt+d|y1:t) } . ただし, pM(yt+d|y1:t) はモデル M の d 日先の予測分布の対数密度関数について実現した 値 yt+dで評価したものである.
表 2 予測精度の推計結果. モデル 1日先 2日先 3日先 4日先 5日先 合計 対数予測分布密度 F-0 −404.7 −464.6 −380.9 −421.0 −413.9 −2085.1 F-S −529.1 −625.3 −495.6 −513.0 −491.3 −2654.3 F-SS −363.8 −425.3 −360.4 −388.8 −359.6 −1897.9 対数予測分布密度比(LPDR) F-S −124.4 −160.7 −114.7 −92.0 −77.4 −569.2 F-SS 40.9 39.3 20.5 32.2 54.3 187.2 表 2 は, 対数予測分布密度と, モデル F-0 を基準として計算した LPDR の結果である. 興 味深いことに, モデル F-S の LPDR は負となっており, 誤差項の分布が対称なモデル F-0 よりも, 分布が非対称なモデル F-S の方が, 予測精度が低いことがわかる. 一方, スパース 事前分布を用いたモデル F-SS は, LPDR が正となっており, モデル F-0 よりも予測精度が 高いことがわかる. このことから, 単純に分布の非対称性のパラメータを入れるだけでは, ファクターが表す分布の非対称性と, 各変数の誤差項の非対称性が適切に識別できない可 能性があり, そのために予測精度が低下していると考えられる. そこでスパース事前分布 を用いて, より parsimonious な構造を組み込むと, 予測精度が向上し, 分布が対称なモデル (F-0) よりも優れたモデルとなることがわかる. 5. 結論と今後の課題 本稿では, 非対称 t 分布付き MSV モデルを紹介し, その推定方法について解説を行った. 分布の非対称性の程度を表すパラメータについて縮小推定する方法についても議論を行っ た. 米国株価市場 S&P500 の日次収益率データを用いた実証分析の結果, 通常の左右対称 な t 分布よりも, 非対称 t 分布とスパース事前分布を用いた方が, MSV モデルの予測精度 が向上することが示された. 今後の課題としては, データや分析の問題意識に合わせて, MSV モデルの拡張を行うこと が考えられる. 例えば, 変数間の相関を表すパラメータである, A や B の各要素が, AR(1) 過程やマルコフスイッチングモデルなど時間を通じて変化する定式化が考えられる. また, ファクターが AR(1) などダイナミックな過程にしたがう定式化も可能である. フレキシブ ルすぎる定式化は推定が不安定となり得るが, 分布の非対称性を考慮しつつ, 多変量間のダ イナミクスを適切に捉えるために様々な工夫を行うことが今後の課題である.
補論. 事後分布からのサンプリング法 この補論では, 3 節で述べた MCMC アルゴリズムの詳細を各ステップごとに解説する. パラメータや潜在変数の条件付事後分布を π(θ| · ) と表記する. これは, θ について, それ以 外の全てのパラメータと潜在変数を条件としたときの事後確率密度関数を表している. な お, ボラティリティの初期値については, hi1∼ N(µi, σi2/(1− φ2i)) とする. (ステップ 2) ボラティリティを表す潜在変数 hi,1:T を効率的にサンプリングするために, multi-move
sampler を用いる (Shephard and Pitt (1997), Watanabe and Omori (2004), Omori and Watanabe (2008)). まず, ¯hit= hit− µi, γi = exp(µi/2), ˆyit= ˜yit/γi, ¯zit= zit− ciと定 義し, SV モデルの部分を次の状態空間モデルに書き換える. ˆ yit = {βiz¯it+ √ zitεit} exp(¯hit/2), t = 1, . . . , T, ¯ hi,t+1 = φi¯hit+ ηit, t = 0, . . . , T− 1. これは, ¯hitを潜在変数とする非線形ガウス型状態空間モデルの形をしており, Omori and
Watanabe (2008) が提案した multi-move sampler を用いて効率的に hi,1:Tをその条件付事後 分布からサンプリングすることができる. 具体的なサンプリング方法については, Nakajima and Omori (2012) の Appendix や, 中島・大森 (2011) の補論 B を参照されたい.
次に, zi,1:T のサンプリングは, 次の条件付事後分布を考える. π(zit| · ) ∝ g(zit)× z−((νi +1)/2+1) it exp ( − νi 2zit ) . ただし, g(zit) = exp { −(˜yit− βiz¯itehit/2)2 2zitehit −(¯hi,t+1− φ¯hit− ¯yit)2 2σ2 i(1− ρ2i) I(t < T ) } , および, ¯yit= ρiσi(˜yite−hit/2− βiz¯it)/√zitである. 上の条件付事後確率密度関数は, g(zit) 以外の部分が逆ガンマ分布の密度関数の核の形をしているため, 提案分布を z∗it∼ IG((νi+ 1)/2, νi/2) として, 候補点 zit∗ をサンプリングし, MH アルゴリズムを適用する. zit∗ の受容 確率は, 現在の値を zitとして, min{g(zit∗)/g(zit), 1} である. コレスキー MSV モデルやファクター MSV モデルの利点として, i ごとに並列して hi,1:T や zi,1:T のサンプリングを行うことができる. (ステップ 3) φiの事前分布を π(φi) と書く. φiの条件付事後分布は, 次のように書くこと
ができる. π(φi| · ) ∝ π(φi) √ 1− φ2 iexp { −(1− φ2i)¯h 2 i1 2σ2 i − T−1 ∑ t=1 (¯hi,t+1− φ¯hit− ¯yit)2 2σ2 i(1− ρ2i) } ∝ π(φi) √ 1− φ2 iexp { −(φi− µφi)2 2σ2 φi } . ただし, µφi = ∑T−1 t=1 (¯hi,t+1− ¯yit)¯hit ρ2 i¯h2i1+ ∑T−1 t=2 ¯h 2 it , σφi2 = σ 2 i(1− ρ2i) ρ2 i¯h2i1+ ∑T−1 t=2 ¯h 2 it . この事後分布の一部が正規分布の密度関数の核の形をしていることから, 提案分布を T N(−1,1)
(µφi, σ2φi) とおいて候補点 φ∗iを発生させ, MH アルゴリズムを適用する. ただし, T N(a,b)(µ, σ2) は平均 µ, 分散 σ2の定義域 (a, b) における切断正規分布を表す. この候補点 φ∗ i の受容確率 は, 現在の値を φiとして, min{(π(φ∗i) √ 1− φ∗2i )/(π(φi) √ 1− φ2 i), 1} である. σiと ρiは同時条件付事後分布からサンプリングを行う. ϑi= (σi, ρi)0と定義して, 事前 分布を π(ϑi) と設定すると, ϑiの条件付事後確率密度関数は次式で表される. π(ϑi| · ) ∝ π(ϑi)× σiT(1− ρ 2 i) T−1 2 exp { −(1− φ2i)¯h 2 i1 2σ2 i − T−1 ∑ t=1 (¯hi,t+1− φ¯hit− ¯yit)2 2σ2 i(1− ρ2i) } . この条件付事後確率密度関数から確率標本を直接サンプリングすることはできないため, MH アルゴリズムを用いる. 条件付事後分布から適当な提案分布が見つからないため, 条 件付事後分布のモード (またはモードの近似値) を数値的に求め, その値の周りで対数密度 関数をテイラー展開して得られる多変量正規分布を提案分布とする. 詳しくは, 中島・大森 (2011) の補論 B を参照されたい. なお, Ishihara et al. (2016) は, 共分散行列のパラメータ をサンプリングする際に, ここで述べた多変量正規近似による条件付事後分布の近似を用 いるのではなく, ウィシャート分布とそれを条件付けた多変量正規分布で提案分布を作っ ている. 本稿の場合でも, ガンマ分布と条件付きの正規分布で ϑiの提案分布を作成するこ とは可能である. µiの事前分布を N (µi0, v2i0) とおく. µiの条件付事後確率密度関数は次の式で表される. π(µi| · ) ∝ exp { −(µi− µi0)2 2v2 i0 −(1− φ2i)¯h 2 i1 2σ2 i − T∑−1 t=1 {¯hi,t+1− φ¯hit− ¯yit}2 2σ2 i(1− ρ2i) } ∝ exp { −(µi− ˆµi)2 2σ2 µi } .
ただし, σµi2 = { 1 v2 i0 +(1− ρ 2 i)(1− φ 2 i) + (T− 1)(1 − φi)2 σ2 i(1− ρ2i) }−1 , ˆ µi = σ2µi { µi0 v2 i0 +(1− ρ 2 i)(1− φ2i)hi1+ (1− φi) ∑T−1 t=1 (hi,t+1− φihit− ¯yit) σ2 i(1− ρ 2 i) } . このサンプリングは, µi| · ∼ N(ˆµi, σ2µi) となる. νiの事前分布を π(νi) と設定する. νiの条件付事後確率密度関数は次の式で表される. π(νi| · ) ∝ π(νi)× T ∏ t=1 (νi/2)νi/2 Γ(νi/2) z−νi/2 it exp ( − νi 2zit ) × exp { − T ∑ t=1 (˜yit− βiz¯itehit/2)2 2zitehit − T−1 ∑ t=1 (¯hi,t+1− φ¯hit− ¯yit)2 2σ2 i(1− ρ2i) } , νi > 4. この条件付事後確率密度関数は, 直接サンプリングを行うことができないため, MH アルゴ リズムを用いる. (σi, ρi) と同様に, 条件付事後確率密度関数のモードを用いた正規分布に よる提案分布でサンプリングを行う. (ステップ 4) まず, βiの事前分布を βi∼ N(0, τ02) と設定する. このとき, βiの条件付事後 分布は, βi| · ∼ N( ˆβi, ˆτi2) となる. ただし, ˆ τi2 = ( 1 τ2 0 + T ∑ t=1 ¯ zit2 zit(1− ρ2i) )−1 , ˆ βi = ˆτi2 T ∑ t=1 ¯ zit { yite−hit/2− √ zitρi(¯hi,t+1− φi¯hit)/σi } zit(1− ρ2i) . 次に, sparsity prior として, βi∼ κN(βi|0, τ02) + (1− κ)δ0(βi) と設定する. この場合, βiの 条件付事後確率密度関数は次式で表される. βi| · ∼ ˆκiN (βi| ˆβi, ˆτi2) + (1− ˆκi)δ0(βi).
ただし, ˆκi = κbi/(κbi + 1− κ), bi = exp( ˆβi2/2ˆτi2)ˆτi/τ0 である (詳しくは, George and McCulloch (1997) を参照).
κ の事前分布をベータ分布 B(aκ0, bκ0) とおく. βiがゼロでない i の数を nκとおくと, κ の事後分布は, κ|· ∼ B(aκ0+ nκ, bκ0+ k− nκ) となる. なお, ファクター MSV の場合は, 2 つ目の引数が bκ0+ q− nκとなる.
を a(i) = (a
i1, . . . , ai,i−1)0と書く (i = 2, . . . , k). a(i)の事前分布を, a(i) ∼ N(ai0, Wi0) と
仮定する. コレスキー MSV モデルは, a(i)について次式のように書き換えることができる. ˘ yit = a (i)0 t xit+ ˆεit, εˆit∼ N ( 0, zit(1− ρ2i) ) . ただし, xit= (y1te−h1t/2, . . . , yi−1,te−hi−1,t/2)0, ˘yit= yite−hit/2− βiz¯it− √zitρi(¯hi,t+1− φih¯it)/σiである. これを利用して, a(i)の条件付事後分布は, a(i)| · ∼ N(ˆa(i), ˆWi) となり, この正規分布からサンプリングすることができる. ただし, ˆ Wi = ( W−1i0 + T ∑ t=1 1 zit(1− ρ2i) xitx0it )−1 , ˆ a(i) = Wˆ i ( W−1i0 ai0+ T ∑ t=1 ˘ yit zit(1− ρ2i) xit ) . である. このサンプリングを, i = 2, . . . , k に対して行う. ファクター MSV モデルのファクター ftとファクター係数 B の条件付事後分布を考え るために, まず, ˆσ2 it= zit(1− ρ2i)ehitとおく. 他のパラメータと潜在変数を条件付けた場合 に, ファクター MSV モデルは次式に書き換えられる. yit = B(i)ft+ αit+ ˆuit, i = 1, . . . , k, fjt = αk+j,t+ ˆuk+j,t, j = 1, . . . , p. ただし, α`t = {β`(z`t− c`) +√z`tρ`(¯h`,t+1− φ`¯h`t)/σi}eh`t/2, ˆ u`t ∼ N(0, ˆσ2`t), である (` = 1, . . . , q). ファクター係数 B の i 行目で推定すべきパラメータの部分を bi = (Bi1, . . . , Bir)0と書 く. ただし, i = 2, . . . , k について, r = min(i− 1, p) である. 上式をさらに書き換えると次 式が得られる. ˇ yit = dit+ g0itbi+ ˆuit, uˆit∼ N(0, ˆσit2).
ただし, ˇyit= yit− αit, git= (f1t, . . . , frt)0, および, dit = fit (i = 2, . . . , p), 0 (i = p + 1, . . . , k), である. biの事前分布を, bi ∼ N(bi0, Zi0) と仮定すると, 条件付事後分布は正規分布とな り, bi| · ∼ N(ˆbi, ˆZi) で表される. ただし, ˆ Zi = ( Z−1i0 + T ∑ t=1 1 ˆ σ2 it gitg0it )−1 , ˆbi = ˆZi ( Z−1i0 bi0+ T ∑ t=1 (ˇyit− dit) ˆ σ2 it git ) , である. 次に, ファクター ftのサンプリングを考えるために, ˇyt= (ˇy1t, . . . , ˇykt)0, ˆut= (ˆu1t, . . . , ˆukt)0 と定義し, 次式を得る. ˇ yt = Bft+ ˆut, ft ∼ N(αt, ˆΨt). ただし, αt = (αk+1,t, . . . , αqt)0, ˆΨt= diag(ˆσ2k+1,t, . . . , ˆσ 2 qt) である. ftの条件付事後分布 は, ft| · ∼ N(ˆat, ˆVt) となる. ただし, ˆVt = ( ˆΨ −1 t + B0Σˆ −1 t B)−1, ˆat = ˆVt( ˆΨ −1 t αt+ B0Σˆ−1t yˇt), ˆΣt= diag(ˆσ21t, . . . , ˆσkt2 ) である. この正規分布からのサンプリングを t = 1, . . . , T に対して行う. 謝辞 本稿の作成に当たっては, 匿名の査読者 2 名と吉羽要直編集委員長, 渡部敏明教授 (一橋大 学) より大変有益なコメントを頂いた. ここに記して感謝したい. 参 考 文 献
Aas, K. and Haff, I. H. (2006). The generalized hyperbolic skew Student’s t-distribution, Journal of Financial Econometrics, 4, 275–309.
Aguilar, O. and West, M. (2000). Bayesian dynamic factor models and portfolio allocation, Journal of Business and Economic Statistics, 18, 338–357.
Asai, M., McAleer, M. and Yu, J. (2006). Multivariate stochastic volatility: A review, Econometric Reviews, 25, 145–175.
Barndorff-Nielsen, O. E. (1977). Exponentially decreasing distributions for the logarithm of particle size, Pro-ceedings of the Royal Society of London, Series A, 353, 401–419.
Black, F. (1976). Studies of stock market volatility changes, in Proceedings of the American Statistical Associ-ation, Business and Economic Statistics Section, 177–181.
Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity, Journal of Econometrics, 31, 307–327.
Chib, S., Nardari, F. and Shephard, N. (2006). Analysis of high dimensional multivariate stochastic volatility models, Journal of Econometrics, 134, 341–371.
Chib, S., Omori, Y. and Asai, M. (2009). Multivariate stochastic volatility, in Andersen, T. G., Davis, R. A., Kreiss, J.-P. and Mikosch, T. eds. Handbook of Financial Time Series, Springer-Verlag Berlin Heidelberg, 365–400.
Clyde, M. and George, E. I. (2004). Model uncertainty, Statistical Science, 19, 81–94.
Cogley, T. and Sargent, T. J. (2005). Drifts and volatilities: Monetary policies and outcomes in the post WWII U.S., Review of Economic Dynamics, 8, 262–302.
Doornik, J. A. (2006). Ox: Object Oriented Matrix Programming, London, Timberlake Consultants Press. Engle, R. F. (1982). Autoregressive conditional heteroskedasticity with estimates of the variance of U.K.
infla-tion, Econometrica, 50, 987–1008.
Fr¨uwirth-Schnatter, S. and Lopes, H. F. (2018). Parsimonious Bayesian factor analysis when the number of factors is unknown. Working Paper.
Gamerman, D. and Lopes, H. F. (2006). Markov Chain Monte Carlo. Stochastic Simulation for Bayesian Inference, Boca Raton, FL, Chapman & Hall/CRC, 2nd edition.
George, E. I. and McCulloch, R. E. (1993). Variable selection via Gibbs sampling, Journal of the American Statistical Association, 88, 881–889.
George, E. I. and McCulloch, R. E. (1997). Approaches for Bayesian variable selection, Statistica Sinica, 7, 339–373.
George, E. I., Sun, D. and Ni, S. (2008). Bayesian stochastic search for VAR model restrictions, Journal of Econometrics, 142, 553–580.
Ghysels, E., Harvey, A. C. and Renault, E. (2002). Stochastic volatility, in Rao, C. R. and Maddala G. S. eds. Statistical Methods in Finance, Amsterdam: North-Holland, 119–191.
Gouri´eroux, C., Jasiak, J. and Sufana, R. (2009). The Wishart autoregressive process of multivariate stochastic volatility, Journal of Econometrics, 150, 167–181.
伊庭幸人, 種村正美, 大森裕浩, 和合肇, 佐藤整尚, 高橋明彦 (2001). 『計算統計 II マルコフ連鎖モンテカルロ法と その周辺』岩波書店．
Ishihara, T., Omori, Y. and Asai, M. (2016). Matrix exponential stochastic volatility with cross leverage, Computational Statistics and Data Analysis, 100, 331–350.
Kim, S., Shephard, N. and Chib, S. (1998). Stochastic volatility: likelihood inference and comparison with ARCH models, Review of Economic Studies, 65, 361–393.
小西貞則, 越智義道, 大森裕浩 (2008). 『計算統計学の方法−ブートストラップ, EM アルゴリズム, MCMC』朝 倉書店．
Koop, G. (2003). Bayesian Econometrics, Chichester, Wiley.
Lopes, H. F. and Carvalho, C. M. (2007). Factor stochastic volatility with time varying loadings and Markov switching regimes, Journal of Statistical Planning and Inference, 137, 3082–3091.
Lopes, H. F. and West, M. (2004). Bayesian model assessment in factor analysis, Statistica Sinica, 14, 41–67. Lopes, H. F., McCulloch, R. E. and Tsay, R. (2012). Cholesky stochastic volatility models for high-dimensional
time series, Technical report, University of Chicago, Booth Business School. 増田弘毅 (2002). 「GIG 分布と GH 分布に関する解析」『統計数理』50, 165–199．
Nakajima, J. (2017). Bayesian analysis of multivariate stochastic volatility with skew distribution, Econometric Reviews, 36, 546–562.
Nakajima, J. (2020). Skew selection for factor stochastic volatility models, Journal of Applied Statistics, 47, 582–601.
中島上智, 大森裕浩 (2011). 「一般化双曲型非対称 t 分布を用いた確率的ボラティリティ変動モデルの推定と株価 収益率データへの応用」『日本統計学会誌』40, 1–28．
Nakajima, J. and Omori, Y. (2012). Stochastic volatility model with leverage and asymmetrically heavy-tailed error using GH skew Student’s t-distribution. Computational Statistics and Data Analysis, 56, 3690–3704. Nakajima, J. and Watanabe T. (2011). Bayesian analysis of time-varying parameter vector autoregressive model
with the ordering of variables for the Japanese economy and monetary policy. Global COE Hi-Stat Discussion Paper Series 196, Hitotsubashi University.
中妻照雄 (2003). 『ファイナンスのための MCMC 法によるベイズ分析』三菱経済研究所． 大森裕浩 (2001). 「マルコフ連鎖モンテカルロ法の最近の展開」『日本統計学会誌』31, 308–334．
大森裕浩, 渡部敏明 (2008). 「MCMC とその確率的ボラティリティ変動モデルへの応用」国友直人・山本拓（編） 『21 世紀の統計科学 I 社会・経済と統計科学』東京大学出版会, 第 9 章, 223–266．
Omori, Y. and Watanabe, T. (2008). Block sampler and posterior mode estimation for asymmetric stochastic volatility models, Computational Statistics and Data Analysis, 52, 2892–2910.
Omori, Y., Chib, S., Shephard, N. and Nakajima, J. (2007). Stochastic volatility with leverage: fast likelihood inference, Journal of Econometrics, 140, 425–449.
Pinheiro, J. C. and Bates, D. M. (1996). Unconstrained parametrizations for variance-covariance matrices, Statistics and Computing, 6, 289–296.
Pourahmadi, M. (1999). Joint mean-covariance models with application to longitudinal data: Unconstrained parametrisation, Biometrika, 86, 677–690.
Prause, K. (1999). The Generalized Hyperbolic models: estimation, financial derivatives and risk measurement. PhD dissertation, University of Freiburg.
Primiceri, G. E. (2005). Time varying structural vector autoregressions and monetary policy, Review of Eco-nomic Studies, 72, 821–852.
Shephard, N. (2005). Stochastic Volatility: Selected Readings, Oxford, Oxford University Press.
Shephard, N. and Pitt, M. (1997). Likelihood analysis of non-Gaussian measurement time series, Biometrika, 84, 653–667.
Smith, M. and Kohn, R. (2002). Parsimonious covariance matrix estimation for longitudinal data, Journal of the American Statistical Association, 97, 1141–1153.
和合肇 (編) (2005). 『ベイズ計量経済分析』東洋経済新報社． 渡部敏明 (2000). 『ボラティリティ変動モデル』朝倉書店．
Watanabe, T. and Omori, Y. (2004). A multi-move sampler for estimating non-Gaussian time series models: Comments on Shephard & Pitt (1997), Biometrika, 91, 246–248.
