« 読了：Butler & Denham (2000) PLS回帰の偏回帰係数は奇妙な縮小特性を示す | メイン | 読了：「ペリリュー」「近所の最果て」「マイ・ブロークン・マリコ」「さめない町の喫茶店」「７人のシェイクスピア」 »
2020年2月14日 (金)
時系列データを扱っていると、たまーに単位根検定をやりたくなることがあるんだけど、Rのパッケージがいろいろあって困る。
仕方がないので、目についたやつについてメモを取った。なんというか、恥をさらしているような気がしますが...
tseriesパッケージ
時系列分析パッケージの古手だと思う。ADF検定, KPSS検定, PP検定の関数を持っている。
adf.test(x, alternative, k): H0「xは単位根を持つ」のADF検定。モデルは定数と線形トレンドを含む。引数は:
- x: 時系列
- alternative: 対立仮説の指定。値は"stationary", "explosive"。後者はnonstatonaryってこと？
- k: ラグ次数。デフォルトではtrunc((length(x)-1)^(1/3))となるそうだ。k=0とするとDickey-Fuller検定になる。
えーっと、いっつも混乱するんだけど、ここで「モデルは定数と線形トレンドを含む」っていうのはどういうことなの？
話を単純にするためにラグ次数を1として、ここで「定数と線形トレンドを含む」といっているのは、差分時系列について
$\Delta y_t = \beta_1 + \beta_2 t + (\phi-1) y_{t-1} + e_t$
というドリフト＋トレンドつきモデルを考えて$H_0: \phi = 1$の検定をいたします、ってこと？ それとも、元の時系列に切片と1次の項をいれて
$y_t = \mu + \beta_1 t + \phi y_{t-1} + e_t$
$\Delta y_t = \beta_1 + (\phi-1) y_{t-1} + e_t$
つまり差分時系列はドリフトつきモデルです、ってことなの？ うぐぐぐぐ。
[02/16 追記: コードをざっと眺めたところ、どうやら前者、つまり$\Delta y_t$のモデルに定数項と1次の項がはいるという話らしい... しらんけど]
kpss.test(x, null, lshort): H0「xはレベル定常ないしトレンド定常である」のKPSS検定。引数は
- x: 単変量時系列
- null: H0はレベル定常かトレンド定常か。値は"Level", "Trend". しっかし、すごい引数名でびびるわ... nullになにかを代入する日が来るとは...
- lshort: 値は論理値。TRUEのとき、truncation lagパラメータが trunc(4+(n/100)^0.25), FALSEのとき trunc(12+(n/100)^0.25) となるのだそうだ。なんだかよくわからん、KPSS検定についてちゃんと勉強しなきゃ。
pp.test(x, alternative, lshort): H0「xは単位根を持つ」のPP検定。定数項と線形トレンドを含む。引数は
- x: 単変量時系列
- alternative: 対立仮説の指定。adf.test()と同じ。
- type: 検定のタイプ。値は"Z(alpha)", "Z(t_alpha)". ううむ、PP検定について不勉強で、なにいってんだかさっぱりわからん。
- lshort: kpss.test()と同じ。
urcaパッケージ.
単位根検定のためのパッケージとして最も有名なのはこれだと思う。
ur.df(y, type, lags, selectlags): ADF検定。引数は
- y: 時系列。
- type: "none"だと切片もトレンドもなし, "drift"だと切片をいれる、"trend"だと切片とトレンドをいれる。これは差分時系列のモデルのことをいってんですよね。
- lags: ラグ次数.
- selectlags. "Fixed"だとなにもしない。"AIC", "BIC"だとlagsに指定した以下の次数をAICなりBICなりで勝手に選んでくれる。
ur.ers(y, type, model, lag.max): Ellot, Rothenberg, & Stock の単位根検定だそうだ。なにそれってかんじだけど、どうやらADF-GLS検定のことらしい。ADF検定より検定力が高いって前に読んだことがあるぞ、理屈はさっぱり理解できなかったけど。引数は
- y: 時系列.
- type: 値は"DF-GLS", "P-test". 後者は誤差項の系列相関を考慮するのだそうだ。へー。
- mode: 値は"constant", "trend". これって元の時系列の話？それとも差分時系列の話？
- lag.max: よくわからんけど、ラグ次数の最大値だそうな。
ur.kpss(y, type, lags, use.lag): KPSS検定. 引数は
- y: 時系列.
- type: 値は"mu", "tau". 前者は切片のみ、後者は切片と線形トレンド。
- lags: 値は"short", "long", "nil"。"short"にするとラグ次数はほにゃらら(メモ省略)、"long"にするとラグ次数はほにゃらら、"nil"にする誤差項の指定なし。
- use.lag: lagsを使わずに、ここで最大ラグ次数を指定してもよい。
ur.pp(x, type, model, lags, use.lag): PP検定。引数は
- x: 時系列.
- type: 値は"Z-alpha", "Z-tau". なんだかわからん。
- model: 値は"constant", "trend".
- lags: 値は"short", "long". マニュアルに説明が書いてない...
- use.lag: lagsを使わずにここで次数を指定してもよい。
ur.sp(y, type, pol.deg, signif): Schmidt & Phillips単位根検定。黒住(2008)ではADF検定ではない「その他の検定」というところで名前だけ出てくる。いわく、ADF検定ってのはWaldタイプの検定なんだけど(なるほど、最尤推定量の差の検定なわけね)、これはラグランジュ乗数検定なんだってさ。そういわれても困るけどな！
引数は
- y: 時系列
- type: 値は"tau", "rho". そういう種類があるんだってさ。しらんがな。
- pol.deg: 多項式の次数だそうだ。値は1,2,3,4。
- signif: 有意水準。値は0.01, 0.05, 0.1.
- y: 時系列
- model: 値は"intercept", "trend", "both". potential breakが切片で起きるか、線形トレンドで起きるか、両方で起きるか。
- lag: ラグ次数の最大値。指定しなくてもいいらしい(えっ、どういうこと？)
forecastパッケージ
泣く子も黙る(?) 有名パッケージ。中の人Hyndman先生は、ただいまこのパッケージのtidyverse対応版であるfableパッケージを鋭意ご製作中らしいのだが、単位根検定関係はまだ移植していない模様。
ndiffs(x, alpha, test, type, max.d, ...): 超お手軽な単位根検定の関数。実はこれ、昨年社外で時系列解析のセミナーやった時にはじめて存在を知り、あまりの簡単さにがっくり膝をついた次第である。だってあれですよ、検定統計量とか一切無視で、単に「何階差分をとるべし」という整数値だけをぽろっと返してくるんですよ。それって、カツカレーくださいって言ったらカツカレーが出てきたようなものじゃないですか...(ちょっとちがうか)
引数は
- x: 時系列
- alpha: 有意水準. 0.01から0.1までの値。
- test: 値は"kpss", "adf", "pp"。
- type: 値は"level", "trend". 恥ずかしながら、これがよく理解できてないんですが... 差分時系列に切片だけはいるか(つまりドリフトつきモデルか)、$t$の1次項もはいるか(つまりトレンドつきモデルか)、ということでしょうか。ってつぶやいてないで、今度ちゃんと調べよう。
- max.d: ラグ次数の最大値。AICかなんかで自動選択してるんだろうなあ。これも今度調べてみよう。
- ...: 単位根検定に渡す引数。たぶんurcaパッケージに渡るんだと思う。
[02/16 追記: 上記の想像通りで、ndiffs()のコードをざっと眺めたところ、ur.df(), ur.kpss(), ur.pp()のいずれかをコールしている模様。次のような対応になっているようだ：
ndiffs(x, test="adf",type="level") → ur.df(x, type="drift")をコール
ndiffs(x, test="adf",type="trend") → ur.df(x, type="trend")をコール
ndiffs(x, test="kpss", type="level") → ur.kpss(x, type = "mu")をコール
ndiffs(x, test="kpss", type="trend") → ur.kpss(x, type="tau")をコール
ndiffs(x, test="pp", type="level") → ur.pp(x, type="Z-tau", model="constant")をコール
ndiffs(x, test="pp", type="trend")→ ur.pp(x, type="Z-tau", model="trend")をコール]
nsdiffs(x, alpha, test, max.D, ...): ndiffs()のSARIMA版、つまり、季節について何階差分をとるべきかを返す。testの値は"seas", "ch", "hegy", "oscb"。おっと、hegyってのがある... uroot::hegy()と同じことなのかなあ...
そのほかのパッケージ
uroot::hegy(x, ...): Hylleberg, Engle, Granger & Yoo の季節単位根検定、だそうだ。勉強不足であれですけど、要するにモデルに季節ダミーが入りますってこと？ それとも、SARIMA(p,d,q,P,D,Q)モデルを当てはめるときにDをどうするかってこと？ いやまてよ、それって実は同じことなの？ うーん、よくわからん。引数はいっぱいあるので省略する。
CADFtest::CADFtest(model, x, type, data, max.lag, min.lag.y, min.lag.X, max.lag.X, dname, criterion, ...): Covariate ADF単位根検定、だそうだ。どうやら、共変量をいれたADF検定らしい。へええ！そんなのあるんだ！ ってことはあれですか、時系列変数間で回帰するとき、yが単位根を持っていてもこの検定に通るなら差分取らなくていいよってことっすか？ それ助かるわー... J.Stat.Softwareに論文が出てるらしい、今度読んでみよう。
MultipleBubbles::ADF_FL(y, adflag, mflag): ラグ次数を固定したADF検定。adflagがラグ次数、mflagは1のとき切片のみ、2のとき切片とトレンド、3のとき両方なし。これ、関数のヘルプを見る限りふつうのADF検定なんだけど、このパッケージは「バブルの存在をチェックするためのADF検定」のパッケージなので、きっとなんか特別なことをやっているんだろうなあ。
MultipleBubbles::ADF_IC(y, adflag, mflag, IC): 上と同じだが、mflagは最大ラグ次数となり、IC=1とするとAIC, 2とするとBICで次数選択するらしい。
fUnitRootsパッケージ: チューリッヒのThe Rmetrics Accoc.というところが本やRパッケージをたくさん出していて、これもそのひとつ。ADF検定の関数を自前で持っている(unitrootTest(), adfTest())。ほかにurcaパッケージへのラッパーもある。
FinTS::Unitroot(x, trend, medhod, lags): これはTsay(2005)という本のコンパニオン・パッケージ。この関数はfUnitRootsパッケージへのラッパーだそうだ。
aTSAパッケージ: SASとよく似た出力を返す時系列分析パッケージなのだそうだ(ははは)。adf.test(), pp.test(), kpss.test()を持っている。
Rの時系列分析パッケージとしては、ほかにTSAという強力なやつがあるんだけど、単位根検定の関数は持っていないようだ。
雑記：データ解析 - 覚え書き：単位根検定のためのRパッケージ