2020/05/28

R - 実践編 4 時系列のモデル推定

時系列データに当てはまるモデルの推定を R で行ってみましょう。実践編 3 で用いたものと同じ時系列 TS を対象として、モデルを推定した様子を紹介したいと思います。

基本統計量の確認

モデルの推定に先立って、基本統計量を確認しておく必要があります。平均、分散の値や単位根が無いこと、explosive であることは、実践編 3 ですでに見ている通りなので、ここでは、標本自己相関関数、標本偏自己相関関数を確認します。
R には、acf() という関数が用意されています。時系列データから標本自己相関関数を計算、返却しつつ、コレログラムをプロットしてくれる関数です。また、pacf() 関数もあり、こちらは標本偏自己相関関数の計算・プロットです。なお、acf() に type = "partial" を指定しても pacf() を呼んでくれます。

どの程度の次数まで見れば良いか分からないので、とりあえずデータの件数分全部見てみることにします。
> acf(TS, lag.max=nrow(TS))
> pacf(TS, lag.max=nrow(TS))
ちょっと見えづらいですが、青破線が「自己相関や偏自己相関が 0 である」という帰無仮説の有意水準5%での棄却点を示しています。標本自己相関関数は長期間に渡ってゆるやかに振動しながら減衰しており、標本偏自己相関関数は低い水準で激しく振動しています。少ない次数では偏自己相関があるようにも見えます。

forecast パッケージには、tsdisplay() 関数があり、これは原系列と acf, pacf を1つのプロットとして出力するものです。
> tsdisplay(TS)
原系列は諸事情によりボカしています。すみません。lag.max を指定していないので、適当なラグが選択されています。これで見ると標本偏自己相関関数については50次あたりまで減衰しながらの正の相関を持っているように見えています。

モデルの選択

基本統計量の概観を元に、モデルを選択してみます。自己相関に切断が見られないため、MA モデルは適しないと思われます。切断無く減衰傾向があることから、AR または ARMA モデルの適用ができる可能性があります。また、偏自己相関にも切断は見られないため、AR モデルよりは ARMA モデルの適用を試みる価値が高いのではないかと思われます。
ARMA モデルの場合、次数の目安はここからは分からないので、いくつかフィッティングを行った上で情報量規準による選択を行うことになります。

ARMA (ARIMA) モデルのフィッティング

stats::ar() は AR モデルの次数を AIC に基づいて選択する機能を持っているのですが、ARMA モデルについての同様の関数は標準にはありません。
ARMA モデルの最適次数の選択は、forecast パッケージにある auto.arima() を使います。この関数は自己回帰和分移動平均 (Autoregressive integrated moving average / ARIMA) モデルへのフィッティングを行うものです。ARIMA モデルは非定常過程のためのモデルで、階差に対して ARMA モデルを適用するものになっています。$ARIMA(p, d, q)$は$d$階の階差に$ARMA(p, q)$を適用するモデルを意味します。$ARMA(p, q)$は$ARIMA(p, 0, q)$と同等なので、auto.arima() が選ぶモデルの中には ARMA が選択対象に含まれています。ちなみに$p, q$についても$0$が選択され得るので、AR, MA も含まれています。つまり、ARIMA モデルは、AR, MA, ARMA の選択や階差の取り方を全てパラメータとして押し込んだモデルになっているので、ARIMA のパラメータ推定を情報量規準に基づいて行うことで、モデルの選択も含めてワンストップでできることになります。
auto.arima() を実行してみます。
> TSarima <- auto.arima(coredata(TS))
> TSarima
Series: coredata(TS)
ARIMA(3,1,2)

Coefficients:
ar1 ar2 ar3 ma1 ma2
0.7562 -0.0543 0.0405 -1.1877 0.1992
s.e. 0.0705 0.0390 0.0078 0.0706 0.0694

sigma^2 estimated as 692537: log likelihood=-696194.6
AIC=1392401 AICc=1392401 BIC=1392457

auto.arima() は配列を渡さないとうまく動作しないので、coredata(TS) を渡しています。次数$p, q$は初期値では$2$から$5$の値を試すようになっています。この範囲では、$ARIMA(3, 1, 2)$が最適ということになったようです。
もう少し次数の幅を広げて試してみましょう。
> TSarima <- auto.arima(coredata(TS), start.p = 0, max.p = 50, start.q = 0, max.q = 50)
> TSarima
Series: coredata(TS)
ARIMA(1,1,3)

Coefficients:
ar1 ma1 ma2 ma3
0.7719 -1.2032 0.1512 0.0619
s.e. 0.0056 0.0067 0.0057 0.0048

sigma^2 estimated as 692511: log likelihood=-696193
AIC=1392396 AICc=1392396 BIC=1392443

$ARIMA(1, 1, 3)$が、$ARIMA(3, 1, 2)$よりもわずかに AIC の値が小さく、最適とのことになりました。ic="bic" を付けると、AIC にかえて BIC で評価するようになるが、これでも同じく$ARIMA(1, 1, 3)$が選択されました。
従って、$ARIMA(1, 1, 3)$がモデル化の有力な候補ということになります。

モデルの診断

$ARIMA(1, 1, 3)$が一番フィッティングが良かったからといって、これが正解のモデルとは限りません。あくまでモデルの候補が見付かったに過ぎません。モデルの診断を行う必要があります。自己相関のモデル化をしているわけですから、残差に自己相関が無いことを確認することになります。こちらの記事でも取りあげましたが、かばん検定といいます。R では Box.test(data, lag=m, fitdf=k) という関数が用意されています。$Q(m)$値が$\chi^2(m-k)$で検定されます。
今回は$ARIMA(1, 1, 3)$なのでパラメータ数$k$は$5$です。$Q(6)$から$Q(50)$まで確認してみることにします。
> sapply(6:50, function(l) Box.test(TSarima$residuals, lag=l, type="Ljung-Box", fitdf=5)$p.value)
[1] 6.119297e-07 2.136468e-07 8.512888e-08 0.000000e+00 0.000000e+00
[6] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[31] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
[41] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00

p値が極めて小さく、Ljung-Box 検定の帰無仮説である「自己相関が無いこと」は棄却されてしまいました。つまり自己相関が残っているとみられます。

結果として、ARIMA モデルへのフィッティングは失敗だったということになります。acf, pacf から自己相関があることはほぼ間違いないので、ARIMA の範囲でモデル化できるような構造の自己相関ではなかったということです。より複雑な他のモデルを検討する必要があるでしょう。