読了: Rodriguez-Rondon & Dufour (2024) MSTestパッケージをインストールしなさい。マルコフ・スイッチング・モデルのレジーム数についての本当の検定というものをお見せしますよ

　このたび仕事で時系列のレジーム・スイッチング・モデルの話に取り組んでいたんだけど、困ったのが、ある時系列についてわざわざレジーム・スイッチを考える必要があるのかないのか、どうやって判断すんの？ という問題であった。(あまりに初歩的な疑問で、なんだか恥をさらしているような気がする…)

　きっとなにかしら検定のような方法があるはずだ。でも、直観的には、単にレジーム数1のモデルと2のモデルのあいだで尤度比検定してはいけないような気がする。しかし、なぜそう思うのか自分でもうまく説明できない。
　RにはMSTestというパッケージがあって、レジーム数についての検定をご提供します、って書いてある。お、これじゃん。ところが、出力の見方がさっぱりわからない。思い余ってChatGPTくんに尋ねてみると、いいえ、RのMSTestパッケージとはマルチ・ステージ検定のパッケージです、マルコフ・スイッチング・モデルのレジーム数については、2つのモデルの間での尤度比検定が可能ですと断言し、Rのコードと出力例を示してくる。しかし、どうみてもシンプル過ぎるコードと雑な出力で、疲れている私の目から見ても疑わしい。ひとこと「ほんとですか？」と問い返すと、しばらく考えたのちに「申し訳ありません、検定はMSTestパッケージによって可能です」と反省のそぶりを見せるが、今度はありもしない関数の使い方を説明してくる。
　キレるよ？ 夜中にブチギレちゃうよ？？？

　結局時間切れとなり、課題そのものは検定ではなくてもっと定性的な方法で解決したけれど(そっちのほうが良かったと思う)、出力の読み方がわからんというのはいささか後味が悪かった。

Rodriguez-Rondon, G., Dufour, J.M. (2024) MSTest: An R-Package for Testing Markov Switching Models. arXiv:2411.08188.

　というわけで、用事はもう済んだんだけど、後学のために開発者謹製のプレプリントを読んでみた。
　名前の綴りをちゃんとみてなかったので、なんとなく「若手院生の恋人同士が仲良く書いてたりしてね…」なんて想像しながら読んでいたのだが、調べてみたところ第二著者はMarieさんではなく、Jean-Marie Dufour、高名な計量経済学者、御年76歳。騙された…(騙してません)

1. イントロダクション
　[マルコフ・スイッチング・モデルというのがありまして…広く使われてまして…]
　MSモデルではレジーム数を事前に決めないといけないので、レジーム数についての検定に関心が持たれる。従来提案されている方法はたいてい、帰無仮説が線形モデル、対立仮説が2レジームのMSモデルである。これまでにあった提案として… [列挙されている。略]。
　というわけで、MSTestというパッケージを作りましたよ。

2. マルコフ・スイッチング・モデル
　MSTestパッケージが想定するのは、平均か分散がマルコフ過程\(S_t\)に支配されているというモデルである。マルコフ・スイッチングモデルのうち自己回帰係数が含まれていないモデルのことを隠れマルコフモデルというけれど、これも想定内である。外生的説明変数があるのも想定内である。
　[マルコフ・スイッチング・モデルと隠れマルコフモデルってそういう違いがあるの？！ てっきり同じかと思ってた。それはともかく、えーと、外生的な説明変数の係数がスイッチする奴はカバーされないってことだよね。自己回帰係数がスイッチする奴はカバーされるのだろうか？]

　レジーム数を\(M\)とし、遷移行列を\(\mathbf{P}\)、ただし\(p_{ij} = P(S_t = j | S_{t-1} = i)\)とする。\(M=2\)ならエルゴード確率\(\pi = (\pi_1, \pi_2)^\top\)は$$ \pi_1 = \frac{1 – p_{22}}{2 – p_{11} – p_{22}} $$ となる。一般的には、$$ \mathbf{A} = \left[ \begin{array}{c} \mathbf{I}_M – \mathbf{P} \\ \mathbf{1}^\top \end{array} \right] $$ とし、\(\mathbf{I}_{M+1}\)の右端の列を\(\mathbf{e}_{M+1}\)として、$$ \mathbf{\pi} = (\mathbf{A}^\top \mathbf{A})^{-1} \mathbf{A}^\top \mathbf{e}_{M+1} $$ となる。[しれっと書いてあるけど、そうなんですか？ Hamilton(1994)がreferされている。あの分厚い本だよね… 手に届くところに置いてあるけどいまあれを持ち上げる気力がない]

　MSTestパッケージが考えるモデルは以下のとおり。

  • マルコフ・スイッチング自己回帰モデル。$$ y_t = \mu_{S_t} + \sum_{k=1}^p \phi_k (y_{t-1} – \mu_{S_t – 1}) + Z_t \beta_z + \sigma_{S_t} \epsilon_t $$ [あー、やっぱし自己回帰係数はスイッチしないんだ… そういうもんか…]
    　MSTestパッケージでは、このモデルをMSARXmdlと呼び、\(Z_t\)を抜いたのをMSARmdlという。現バージョンでは\(\epsilon_t\)は正規分布に従う。
    　尤度はどうなるかというと… [めんどくさいのでメモ省略]
  • マルコフ・スイッチングVARモデル。VARに拡張して、$$ \mathbf{y}_t = \mathbf{\mu}_{S_t} + \mathbf{\Phi}_1 (\mathbf{y}_{t-1} -\mathbf{\mu}_{S_{t-1}}) + \cdots + \mathbf{\Phi}_p (\mathbf{y}_{t-p} -\mathbf{\mu}_{S_{t-p}}) + Z_t \beta + \Sigma^{1/2}_{S_t} \epsilon_t $$ MSTestパッケージでゃこれをMSVARXmdl, MSVARmdlという。詳しい話はKrolzig(1997, 書籍)を勝手にみてちょうだい。
  • 隠れマルコフモデル。マルコフ・スイッチングモデルのラグ項を含まないやつで、時系列以外にも使う。$$ y_t = \mu_{S_t} + Z_t \beta + \sigma_{S_T} \epsilon_t $$ って感じね。HMmdlと呼ぶ。

　ふつうマルコフ・スイッチング・モデルはEMアルゴリズムとかベイジアン手法とかカルマン・フィルタとかで推定される。すごく単純なモデルだったら最尤推定できるけど。MSTestパッケージではEM, MLEを指定できる。

3. レジーム数についての仮説検定

3.1 尤度比検定
Hansen(1992 J.AppliedEcon.) いわく、ふつうの尤度比検定はうまくいかない。理由:

  • 尤度関数は帰無仮説と大域最適推定値がある領域で局所二次関数になっていないといけないが、帰無仮説の下でいくつかのパラメータが識別されないので、それらに関してはフラットになってしまう。
  • スコアは正でないといけないが、帰無仮説の最尤推定量ではゼロになってしまう。
  • 遷移確率のようなパラメータは0や1をとりうるので、パラメータ境界問題が生じる。
  • 尤度表面に複数の局所最適解があって、帰無仮説と大域最適解とは同じ領域にないかもしれない。[??? それって問題なの？ レジーム数1と2の検定だとして、レジーム数2の解が真の最尤推定値かどうかわからんとしても、レジーム数1よりましかどうかわかれば御の字じゃないの]

　そこでHansenはこういう方法を考えた。[説明を読んだんだけど、もうなにがなんだか… 尤度比をある関数と見立て、その漸近分布を考えて、標準化して上限をとるというような話だった]
　欠点: (1)標準化した尤度比ではなくてその上限を与えているだけだけなので、保守的かもしれない。(2)計算的に大変。
　でもこのやりかたもパッケージに載せました。ベンチマークとして。

　[以下、尤度比検定路線の既存案のうちパッケージに載せなかった方法についての解説が続く。もちろん読んでません…！]
　… という経緯を踏まえ、著者らは最大化モンテカルロ尤度比検定(MMC-LRT)と局所モンテカルロ尤度比検定(LMC-LRT)というのを開発した。ごくかんたんに説明すると… [どうみてもかんたんじゃない。ここからおよそ3頁スキップする]

3.2 積率ベースの検定
　レジーム数が1か2かという検定についていえば、最小二乗残差の平均・分散・尖度・歪度に焦点をあてて… [以下、話が私の理解の及ばない世界へと飛び去っているのを眺める私であった。1頁スキップ]

3.3 レジーム・スイッチングの最適性検定
　レジーム数の検定は、ランダム係数とマルコフ・スイッチング・モデルのあいだでのパラメータの一致性についての最適性検定と捉えることもできる。つまり、観察できない定常な変数\(\eta_t\)があって$$ H_0: \theta_t = \theta_0$$ $$ H_1: \theta_t = \theta_0 + \eta_t$$だと考えるわけだ。[あー、なるほど… レジームのあるモデルについて推定するのではなくて、レジーム数1の線形モデルのパラメータがほんとに固定か、それとも時変しちゃっているか、という検定をやるのかな？ レジーム数の検定としてみるとちょっとリベラルすぎるような気がするんだけど]

　[以下、わからんなりに細かくメモする]
　Carrasco, Hu, & Ploberger(2014)はこう考えた。\(\eta_t = chS_t\)とする。\(c\)は変化の大きさのスカラー、\(h\)はベクトルで対立仮説の方向を表す。\(S_t = \rho S_{t-1} + \epsilon_t\)とし、\(\epsilon_t\)はiidとする。\(-1 \lt \rho \lt 1\)と制約する(\(S_t\)は平均0で\(\pm 1/(|1-\rho|)\)の範囲を動く)。検定から見ると\(c^2, h, \rho\)が局外パラメータ。
　このとき、\(\mu_{2,t}\)は\(c, h, \rho, \theta, l_t\)の関数として書ける。[式は超めんどくさいしなぜそうなるのか理解できないので省略。要するにレジームがあるときの時変切片をなんらか式で書けるということだと思う。\(l_t\)ってなんなのかよくわかんない]
　ここからsupTSを求めることができる。[これが検定統計量なんだろうけど、式はまっっったく理解できないので省略。そもそもなんのsupなのかさえ理解できない]
　すでに述べたように[知らんがな]、この方法は局外パラメータの分布を通じたブートストラッピングを必要とする。従って、局外パラメータの事前分布を決める必要がある。MSTestパッケージでは一様分布を使っている。しかし、パラメータ\(c^2\)は必ずしも上によってboundされていないので一様分布は適切でないかもしれない。[さっぱりわからん]
　そこでCarrascoらはexpTSという検定統計量も考えた。[もちろん式はなんだかさっぱりわからない。結局のところ、supTSというのとexpTSというのと、どっちを使えばいいんだろうか]

4. MSTestパッケージ
4.1 データセット
[説明用のデータセットの紹介。略]

4.2 シミュレーション
[シミュレーション用の関数の紹介。simuXXX()という名前になっていて、XXXにはNorm, AR, ARX, VAR, VARX, MSAR, MSARX, MSVAR, MSVARX, HMMがはいる。それぞれの機能は名前から想像できる通りだが、NormとHMMには外生変数も入れられる]

4.3 モデル推定
[推定用の関数の紹介。XXXmdl()という名前になっていて、XXXにはN, AR, ARX, VAR, VARX, MSAR, MSARX, MSVAR, MSVARX, HMがはいる。機能は名前から想像できる通りだが、NとはNormのこと、HMとはHMMのこと]

4.4　仮説検定
[このパッケージの本丸、レジーム数の検定用の関数。XXXTest()という名前になっていて、XXXには以下がはいる。

  • LMCMR: 著者らが提案した局所モンテカルロ尤度比検定(LMC-LRT)
  • MMCLR: 著者らが提案した最大化モンテカルロ尤度比検定(MMC-LRT)
  • DLMC: Dufor & Luger (2017)の積率ベース(moment-based)検定。ちょっと待ってよ、Cってなんの頭文字なの…?
  • DLMMC: 積率ベース検定の最大化モンテカルロ(MMC)版。
  • CHP: Carrasco, Hu, & Ploberger(2014)の最適性検定。
  • HLR: Hansen(1992)の尤度比検定。

以下8ページにわたって出力例と説明が続く。愚かなパッケージユーザである自分向けにメモすると、

  • LMCMR, MMCLR: \(H_0, H_1\)のレジーム数を指定する(典型的には1,2)。出力の最後に\(H_0\)の下での\(p\)値がでてくる。つまりこれがちっちゃかったら\(H_0\)が棄却されるわけね。
  • DLMC, DLMMC : \(H_0, H_1\)はレジーム数1,2。DLMCは計算がめっちゃはやい。LMC_minとLMC_prod、ないしMMC-minとMMC-prodっていう行が出てきて、理屈のところを読み飛ばしたものでどうちがうのかわからんが、出力例をみる限り結果は似たようなもんらしい。その右端に\(H_0\)の下での\(p\)値がでてくる。
  • CHP: これも\(H_0, H_1\)はレジーム数1,2らしい。\(H_1\)では平均のみスイッチするのか、分散もスイッチするのかを指定できる。supTSとexpTSという行が出てきて、3章で説明があったけど、残念ながらどうちがうのかよくわからない。右端に\(H_0\)の下での\(p\)値がでてくる。
  • HLR。\(H_0, H_1\)はレジーム数1,2じゃないかと思う。平均のみスイッチするのか、分散もスイッチするのかを指定できる。Mが0,1,2,3,4のときのp値が出てくるんだけど、このMってのはどうやらレジーム数のことじゃなくて、Hansenが提案した標準化の方法のことらしい。ぶひー、不親切だ！！！ (いらだちのあまり机をひっくり返す)]

　
5. 計算例
[略。3つのデータに全手法をあてはめた結果が載っているんだけど、p値が結構違う…]

6. 結論
[略]
—————–
　というわけで、なにがなんだかぜんぜんわからんが、基本的には\(H_0, H_1\)はレジーム数1,2であり(手法によってはほかも指定できる)、\(H_1\)についてまじでパラメータ推定する方法と実はやらない方法があるけれど、まあいずれにせよ出力の最後のほうの行の右端がp値らしい。知らんけど。
　よっしゃ、今日はこのくらいにしといたるわ！ (ぼこぼこに殴られふらつきながら)