Rで計量時系列分析：状態変化を伴うモデル（閾値モデル、平滑推移モデル、マルコフ転換モデル）

前回の記事までは多変量時系列モデルとしてのVARモデルを扱ってきました。今回は一旦このシリーズの最終回ということで、元の単変量時系列モデルに戻って「状態変化を伴うモデル」を扱ってみようと思います。

ということでもはや毎回恒例になってますが、使用テキストはいつもの沖本本です。

経済・ファイナンスデータの計量時系列分析 (統計ライブラリー)

作者: 沖本竜義
出版社/メーカー: 朝倉書店
発売日: 2010/02/01
メディア: 単行本
購入: 4人クリック: 101回
この商品を含むブログ (6件) を見る

ただし今回の最後に出てくるマルコフ転換モデルは沖本本の説明では不足と思われるので、Hamilton本もあった方が良いです。

Time Series Analysis

作者: James D. Hamilton
出版社/メーカー: Princeton Univ Pr
発売日: 1994/01/11
メディア: ハードカバー
購入: 1人クリック: 5回
この商品を含むブログ (8件) を見る

以下タイトルにのっとってRで各モデルの挙動を見ながらやっていきます。

必要なRパッケージ

{tsDyn}, {MSwM}をインストールしましょう。それから前回同様、演習用データを取ってきたいので{RFinanceYJ}もインストールしておいて下さい。ちなみに、{tsDyn}, {MSwM}ともにvignetteが内容が非常に充実していて極めて分かりやすいのでお薦めです。

演習用データは、とりあえず以下のような感じで{RFinanceYJ}を使って取ってきておくことにします。

> st1<-quoteStockTsData('4751.t',since='2012-01-01',date.end='2013-07-31')
> st2<-quoteStockTsData('2432.t',since='2012-01-01',date.end='2013-07-31')
> st3<-quoteStockTsData('3632.t',since='2012-01-01',date.end='2013-07-31')
> st4<-quoteStockTsData('4755.t',since='2012-01-01',date.end='2013-07-31')
> st5<-quoteStockTsData('3793.t',since='2012-01-01',date.end='2013-07-31')
> stock1<-ts(st1$adj_close)
> stock2<-ts(st2$adj_close)
> stock3<-ts(st3$adj_close)
> stock4<-ts(st4$adj_close)
> stock5<-ts(st5$adj_close)

毎回この5銘柄である理由は特にありません（笑）。強いて言えば「現職場は非上場なので」とか。。。

状態変化を伴うモデルとは何か

読んで字の如く、「途中で状態（レジーム：regime）が変わる時系列モデル」です。何を当たり前のこと言ってんだコラと思うかもしれませんが、ここで想定しているのは「ARIMA次数・係数までもが変わってしまうような時系列モデル」です。

つまり、途中で過程自体が変化してしまうようなケースを検出し、それぞれをモデリングするような方法論というわけです。実際に過程自体が切り替わるようなケースというのは特にファイナンスデータの世界では全く珍しくなくて*1、これをモデリングすると同時にそのタイミングもきちんと突き止められれば、むしろ異常値検出と同じ発想で対策を立てるのにも役立つことでしょう。

ということで、それらのモデルを見ていきましょう。沖本本では以下の全てのモデルについてAR(1)過程を例として説明していますが、もちろん普通にARIMA(p,1,q)過程に広く一般化できるので念のため。

閾値モデル

これは、ある状態変数 $s_t$ がある閾値を超えているかどうかによって、 $y_t$ が従うモデルが変化するモデルです。ここではTARモデル(threshold AR model)とそのファミリーを見ていきます。

まず、2状態TAR(1)モデルは沖本本pp.173-174に従えば以下のように書けます。

$\left \{ \begin{array} ~ y_t = \phi_{01} + \phi_{11} y_{t-1} + \sigma_1 \epsilon_t, s_t < c \\ ~ y_t = \phi_{02} + \phi_{12} y_{t-1} + \sigma_2 \epsilon_t, s_t \geq c \end{array} \right.$

なお $\epsilon_t \sim iid(0,1)$ ですが、これまた分布を状態に依存させることも可能です。この式のポイントは、 $\phi_{0k}, \phi_{1k}, \sigma_k$ のどれが状態変化するかによって意味合いが変わるというところを表せるところです。例えば $\phi_{0k}$ が状態変化するなら「期待値」が、 $\phi_{1k}$ なら自己相関構造が、 $\sigma_k$ なら分散が、それぞれ変動するモデルになります。

なおTARモデルは3状態以上にも拡張可能です。

$y_t = \phi_{0k} + \phi_{1k} y_{t-1} + \sigma_k \epsilon_t, c_{k-1} \leq s_t < c_k, k = 1,\cdots,M$

$-\infty = c_0 < c_1 < \cdots < c_{M-1} = \infty$

ところで $s_t$ をどう決めるかですが、沖本本では典型的なやり方として $s_t = y_{t-d}$ や $s_t = |y_{t-d}|$ を挙げています（線形トレンドとして $s_t=t/T$ もしくは $s_t=t$ とすることもある：ただしTは標本数）。この時、状態の推移は過去の $y_t$ 自身の値によって決まるので、自己励起型閾値モデル(self-exciting threshold model)とも呼ばれます。これを先ほどのTARモデルに適用したものを自己励起型TARモデル(self-exciting threshold AR model: SETAR)と呼ぶようです。

なお、閾値モデルの推定は条件付き最尤法によります。沖本本pp.174-175にざっくりとした説明がありますが、ここでは割愛してRでの実行例に行きましょう。setar(){tsDyn}関数で簡単にSETARモデルを推定できます。

> stock1.setar1<-setar(stock1,mL=2,mH=2,nthresh=1,include="const")

 1 T: Trim not respected:  0.8505155 0.1494845 from th: 223300

# AR次数は全てのレジームで2、閾値は1段階のみ、定数項あり

> stock1.setar2<-setar(stock1,mL=2,mM=2,mH=2,nthresh=2,include="const")

 1 T: Trim not respected:  0.8505155 0.1494845 from th: 223300

# 同じくAR次数は全てのレジームで2、閾値は2段階、定数項あり


> plot(stock1.setar1)
> plot(stock1.setar2)

f:id:TJO:20130821171856p:plain

f:id:TJO:20130821171907p:plain

自己励起型閾値モデルに基づいて、この株価時系列のレジームをlow / (middle) / highと切り分けた結果が見て取れるかと思います。また汎用のregime(){tsDyn}関数で細かいレジーム推移の時系列をベクトルとして得ることもできます。明確な閾値で判定できそうor実際に何かしらの閾値をまたいだ変化があると（ドメイン知識から）予想されるケースでは、有用な手法でしょう。

平滑推移モデル

閾値モデルの1つの問題点として、状態変化が離散的というか「カックン」と不連続に起きてしまうように仮定されている点があります。そこで、もっとこれを滑らかにしたらいいんじゃね？ということで提案されたのが平滑推移(ST)モデル(smooth-transition model)です。

モデルの概要として沖本本p.177の式を引用すると、

$y_t = (\phi_{01} + \phi_{11} y_{t-1})(1-G(s_t)) + (\phi_{02} + \phi_{12} y_{t-1})G(s_t) + \sigma \epsilon_t, \epsilon_t \sim iid(0,1)$

ここで $G(\cdot)$ は推移関数(transition function)と呼ばれる関数で、 $s_t$ は推移変数(transition variable)と呼ばれる変数です。式の形を見れば何となく分かるかと思いますが、このモデルは両極に2つの状態を持っていて、その間を推移関数＆変数に応じて滑らかに推移するモデルとして表現されています。

既に閾値モデルのところでも見ていますが、3状態以上にも容易に拡張でき、どの項まで状態変化を反映させるかも完全に自由に決められます。またモデル推定も閾値モデルの時と同様に条件付き最尤法で行います。

推移変数 $s_t$ は閾値モデル同様に選ぶことができます。そして推移関数ですが、これはその関数形を考慮してロジスティック型もしくは指数型が用いられます。ロジスティック型はまさにロジスティック分布のあの関数形と同じもので、極限として閾値モデルを取ることができるので、SETARモデルを内包しています。ちなみにこのロジスティック型の平滑推移モデルをLSTARモデルと呼びます。

Rではlstar(){tsDyn}関数で簡単に実践することができます。

> stock1.lstar<-lstar(stock1,m=2) # AR次数の最大値は2
Using maximum autoregressive order for low regime: mL = 2 
Using maximum autoregressive order for high regime: mH = 2 
Using default threshold variable: thDelay=0
Performing grid search for starting values...
Starting values fixed: gamma =  100 , th =  178026 ; SSE =  16511207226 
Grid search selected lower/upper bound gamma (was:  1 100 ]). 
					  Might try to widen bound with arg: 'starting.control=list(gammaInt=c(1,200))'
Optimization algorithm converged
Problem: singular hessian
Optimized values fixed for regime 2  : gamma =  100 , th =  178026 ; SSE =  16511207226 
> par(mfrow=c(2,1))
> plot(stock1)
> plot(regime(stock1.lstar))

f:id:TJO:20130821175926p:plain

実は結果を比べるとSETARの場合と同じなんですねｗ　あと、lstar(){tsDyn}関数は3状態以上のやり方が何故か見当たらないので、平滑推移でないとどうしても具合が悪いというケース以外ではsetar(){tsDyn}関数の方が使いやすいかもです。

マルコフ転換モデル

さぁ来ました、ラスボスのマルコフ転換モデルです。既に一度ブログ記事にまとめたことがありますが、改めて説明してみようと思います。

また、マルコフ転換モデルに限っては実は沖本本では記述が不足しているので、Hamilton本pp.690-696を参照するか、もしくは{MSwM}パッケージの元になったMatlabパッケージ{MS_Regress}のドキュメントを適宜参照して下さい。

何でこんなことを書くのかというと、沖本本には外部変数 $\mathbf{x}_t$ を導入した立式過程のコメントが一つもないからです。後で出てきますが、{MSwM}パッケージでは外部変数を説明変数として事後確率を求めるプロセスがあるので、ここは覚えておいた方が良いでしょう。

また{MSwM}パッケージのvignetteには原理の部分の説明が殆どないので、上記Matlab{MS_Regress}パッケージのドキュメントを参照して下さい。

基本的なアイデア

前回{MSwM}パッケージ紹介をした時にざっくり書きましたが、改めて再掲するとこんな感じです。

観測時系列が $n$ 個の状態に分かれ得ると仮定して、パラメータフリーなAR過程を $n$ 個置く
観測時系列は $n \times n$ の遷移確率行列で表されるマルコフ過程に従って状態変化すると仮定する
EMアルゴリズムを用いて、観測時系列が $n$ 通りある状態のいずれかであったかの事後確率列を求める

マルコフ過程についてはここでは特に取り上げません。適宜マルコフ過程に関する資料を参照するようにして下さい。またEMアルゴリズム本体についての説明もここでは割愛します。資料としては色々ありますが、例えばPRML邦訳版の下巻pp.146-171あたりを参照すると分かりやすいかと思います。

マルコフ転換モデルの原理

2状態MSAR(1)モデルは、沖本本p.182に従えば以下のように書けます。

$\left \{ \begin{array} ~ y_t = \phi_{01} + \phi_{11} y_{t-1} + \sigma_1 \epsilon_t, s_t = 1 \\ ~ y_t = \phi_{02} + \phi_{12} y_{t-1} + \sigma_2 \epsilon_t, s_t = 2 \end{array} \right.$

これをM状態に拡張させるのは非常に単純で、

$y_t = \phi_{0}(s_t) + \phi_{1}(s_t) y_{t-1} + \sigma (s_t) \epsilon_t$

となります。要は $s_t$ を何とかして求めてやろうという話で、そこで出力結果から遡って過去の状態列を同定していく隠れマルコフモデルのアイデアを取り入れてやろうというのがマルコフ転換モデルの基本的な原理です。

そこでHamiltonが考えたのがEMアルゴリズムを導入してやろうというアイデアで、基本的な原理がHamilton本pp.690-696に書いてあります。・・・が、これを全部説明するのは非常につらいのでここでは割愛します。問題はHamilton本pp.692-693で、ここに観測されている時系列データ $y_t$ 以外にさらに $s_t$ に関する何の情報もなく $y_t$ とも独立な説明変数 $x_t$ を入れることが要求されています。

実は僕もここだけはよく理解できてなくて*2憶測になってしまうんですが、EMアルゴリズムの初期値をここで与えているのかなぁと。。。Hamilton本の導出過程を見ていると、当たり前ながら必ず $x_t$ を含む条件付き確率から入っているので、そういう制約があってのことなんだろうとも思いますが。

MS_Regressのドキュメントの説明でも、p.9で唐突に以下の回帰式が出てくるので分かりにくいんですよね。

$y_t = {\displaystyle \sum_{i=1}^{N_{nS}}} \beta_i x_{i,t}^{nS} + {\displaystyle \sum_{j=1}^{N_S}} \phi_{i,S_t} x_{j,t}^{S} + \epsilon_t$

$\epsilon_t \sim P(\Phi_{S_t})$

ともあれ、気になる人は是非Hamilton本pp.690-696を読んで勉強して、是非正しい理解の仕方を僕に教えてくださいm(_ _)m　よろしくお願いします。。。

Rでの実践例

ということで、何はともあれ以前の記事を踏まえて{MSwM}パッケージの使い方を見て行こうかと思います。改めて手順を再掲すると、

事後確率列を求めたい時系列 $y$ を持ってくる
何でも良いので説明変数になり得る時系列 $x_1, x_2, \cdots, x_n$ を持ってくる
$y$ と $x_n$ とでデータフレームを作り、lm(){stats}関数で回帰モデルmodelを推定する
msmFit(){MSwM}関数にmodel、状態数k、推定パラメータリストswを与えて、マルコフ状態転換モデルを推定する

この $x_n$ ですが、何も良いものが思いつかなければ単調増加関数でも構いません。ただし、Hamilton本でのコメントとは裏腹に、ある程度隠れた状態変数に相関していそうな何かしらの時系列を持ってきた方が良いようです。僕は様々な周期の三角関数を $x_n$ に充てる、という怠惰なやり方をよくしていますｗ

modelは上述のMS_Regressのドキュメントp.9にある通りで、その回帰関係に基づいて $s_t$ の事後確率を推定するというものです。普通にlm(){stats}関数で求めるだけで大丈夫です。

推定パラメータリストswの与え方は、MS_regressのマニュアルpp.9-12に書いてあります。特に何も理由がなければ、最大値目一杯Trueを並べてしまって良いようです（実際にはどう選んでもあまり推定モデルに差が出ない）。

ともあれ動かしてみましょう。推定自体はmsmFit(){MSwM}関数でできます。

> x1<-ts(sin(seq(from=0,to=0.5*pi,length.out=390)))
> x2<-ts(sin(seq(from=0,to=pi,length.out=390)))
# 説明変数を三角関数で与えるようにした

f:id:TJO:20130823152302p:plain

こんな感じで回帰させて、マルコフ転換モデル推定にぶち込みます。

> model1<-lm(stock1~x1+x2)
> msmModel1<-msmFit(model1,k=2,sw=rep(T,4))
> plotProb(msmModel1,1)
> plotProb(msmModel1,2)
> plotProb(msmModel1,3)

f:id:TJO:20130821184818p:plain

f:id:TJO:20130821184830p:plain

f:id:TJO:20130821184839p:plain

という感じで、この株価時系列のアップレジームとダウンレジームが綺麗に分類できました。目視での印象とほぼ合っているのが見て取れるかと思います。

他の株価時系列でもやってみましょう。

> model5<-lm(stock5~x1+x2)
> msmModel5<-msmFit(model5,k=2,sw=rep(T,4))
> plotProb(msmModel5,1)
> plotProb(msmModel5,2)
> plotProb(msmModel5,3)

f:id:TJO:20130821185717p:plain

f:id:TJO:20130821185726p:plain

f:id:TJO:20130821185734p:plain

という感じで、多少（特に早い時期とかに）変な推定結果になっているところもありますが、概ねズドンと株価が跳ね上がったところも捉えられているように見受けられます。

発展：VARモデルのマルコフ転換モデル

TARをVARに拡張したTVARがあったということは、もちろんマルコフ転換モデルをVARに拡張することも可能です。ただ、これだけは沖本本にもHamilton本にも記述がないので、例えばこのドキュメントなどを参照して下さい。

ところで、VARモデルへの拡張は残念ながら{MSwM}パッケージでは対応していません。オリジナルのMatlab{MS_Regress}パッケージは対応していますが。。。Rで実践するには{MSBVAR}のインストールが必要です。

こちらは僕もまだ勉強中なので、参考までに{MSBVAR}のヘルプに掲載されている実行例を以下に示します。なお、ここで用いられているサンプルデータIsraelPalestineConflictとは、イスラエルとパレスチナとの長年にわたる紛争において、片方がもう片方に対して攻撃を仕掛けた回数をまとめたという物騒な代物です。。。

> data(IsraelPalestineConflict)  
> set.seed(123)
# MCMCがベースなので乱数の初期値を設定する

> xm <- msbvar(IsraelPalestineConflict, p=3, h=2,
+              lambda0=0.8, lambda1=0.15,
+              lambda3=1, lambda4=1, lambda5=0, mu5=0,
+              mu6=0, qm=12,
+              alpha.prior=matrix(c(10,5,5,9), 2, 2))
Initial Log Likelihood Value: -12573.35 
Starting blockwise optimization over 40 block optimizations...
Completed iteration 1 of 40. Log Likelihood Value: -12389.76
Completed iteration 2 of 40. Log Likelihood Value: -12360.08
# 中略
Completed iteration 40 of 40. Log Likelihood Value: -12358.54

N1 <- 1000
N2 <- 2000

x2 <- gibbs.msbvar(xm, N1=N1, N2=N2, permute=FALSE, Beta.idx=c(7,1))
Burn-in iteration :  1000 
Final iteration :  1000 
Final iteration :  2000 
# ギブスサンプラーで推定をかける。N1とN2はburn-inとかその辺

plot.SS(x2)