対象データの基本統計量
今回は、多変量時系列データに対して、VAR モデルのパラメータを推定してみましょう。TS という多変量時系列データを用意しました。TS は、X と Y の2つの変量を持っています。基本統計量は以下の通り。
> class(TS)
[1] "xts" "zoo"
> nrow(TS)
[1] 85496
> frequency(TS)
[1] 1
> summary(TS)
Index X Y
Min. :1970-01-01 00:00:00.00 Min. :-61884.00 Min. :-3.423e-04
1st Qu.:1970-01-01 05:56:13.75 1st Qu.: 0.00 1st Qu.: 0.000e+00
Median :1970-01-01 11:52:27.50 Median : 0.00 Median : 0.000e+00
Mean :1970-01-01 11:52:27.50 Mean : -0.01 Mean :-9.850e-08
3rd Qu.:1970-01-01 17:48:41.25 3rd Qu.: 0.00 3rd Qu.: 0.000e+00
Max. :1970-01-01 23:44:55.00 Max. : 65967.00 Max. : 4.028e-04
frequency が 1 の xts データで、85496 行あります。秒刻みで1日分に少し足りないくらいのデータだと思ってください。X も Y も四分位点 Q1, Q2 (Median), Q3 がほぼ 0 なのに最大、最小の絶対値はずいぶん大きい、explosive な感じのするデータです。
単位根検定
X, Y それぞれに対し、Philllips-Perron 検定による単位根検定をしてみます。PP.test() を使います。
> PP.test(TS$X) Phillips-Perron Unit Root Test data: TS$X Dickey-Fuller = -637.5825, Truncation lag parameter = 21, p-value = 0.01 > PP.test(TS$Y) Phillips-Perron Unit Root Test data: TS$Y Dickey-Fuller = -263.3456, Truncation lag parameter = 21, p-value = 0.01Phillips-Perron 検定では、X, Y ともに p-value が小さく、単位根があるとは言い難いとの結論になっています。
※ 実は X も Y も単位根過程らしき時系列の1階の階差です。
標本自己相関の構造
自己相関を確認しておきましょう。> acf(TS, lag.max=10)
> pacf(TS, lag.max=10)
acf() を見ると、X, Y ともに、3～4次あたりまで自己相関を持っているように見えます。減衰傾向もあります。pacf() を見ると、X と Y の間の偏自己相関が多少ありそうにも見えます。
一応、VAR モデルを当ててみる価値はありそうです。
VAR モデルの選択
R で VAR モデルを扱うには、vars パッケージを使います。 VAR モデルの選択は VARselect() によって行うことができます。100次まで探索させてみましょう。> VARselect(TS, 100)
$selection
AIC(n) HQ(n) SC(n) FPE(n)
96 74 48 96
$criteria
1 2 3 4 5
AIC(n) -8.1011924695 -8.1383711520 -8.1530109267 -8.1631402110 -8.1740660487
HQ(n) -8.1009915707 -8.1380363206 -8.1525421628 -8.1625375145 -8.1733294197
SC(n) -8.1005351749 -8.1372756610 -8.1514772394 -8.1611683273 -8.1716559686
FPE(n) 0.0003031774 0.0002921126 0.0002878673 0.0002849661 0.0002818696
(略)
46 47 48 49 50
AIC(n) -8.2401103055 -8.2402657046 -8.2410426603 -8.241164891 -8.2413203039
HQ(n) -8.2338824420 -8.2339039086 -8.2345469317 -8.234535230 -8.2345567103
SC(n) -8.2197341738 -8.2194513765 -8.2197901358 -8.219474170 -8.2191913867
FPE(n) 0.0002638551 0.0002638141 0.0002636093 0.000263577 0.0002635361
(略)
71 72 73 74 75
AIC(n) -8.2458743893 -8.2463429275 -8.2464925460 -8.2466685854 -8.2466740306
HQ(n) -8.2362982121 -8.2366328178 -8.2366485037 -8.2366906106 -8.2365621233
SC(n) -8.2145433480 -8.2145736898 -8.2142851120 -8.2140229550 -8.2135902038
FPE(n) 0.0002623386 0.0002622157 0.0002621765 0.0002621304 0.0002621289
(略)
96 97 98 99 100
AIC(n) -8.2483858131 -8.2483204228 -8.2483074155 -8.2483632402 -8.2483824885
HQ(n) -8.2354613222 -8.2352619994 -8.2351150595 -8.2350369518 -8.2349222675
SC(n) -8.2060998623 -8.2055962756 -8.2051450719 -8.2047627003 -8.2043437522
FPE(n) 0.0002616806 0.0002616977 0.0002617011 0.0002616865 0.0002616815
$selection が各種情報量規準に照らしたときの最適次数を示しています。AIC では 96、SC (SIC と同じです) では 48 が最適次数とのことになりました。HQ はハンナン=クイン情報量規準という、また少しペナルティ関数の違う情報量規準です。FPE は赤池最終予測誤差のことで、AIC の原型的なものです。次数は 48 か 96 が良さそうとなったので、VAR(48), VAR(96) それぞれのフィッティングを行い、予測値を出させてみましょう。
> TSvar48 <- VAR(TS, 48)
> TSvar96 <- VAR(TS, 96)
> predict(TSvar48)
$X
fcst lower upper CI
[1,] 28.90594185 -1603.502 1661.314 1632.408
[2,] -80.28080566 -1852.209 1691.648 1771.929
[3,] 83.28752937 -1712.583 1879.158 1795.871
[4,] 20.07676960 -1779.406 1819.560 1799.483
[5,] -28.53638886 -1831.497 1774.424 1802.961
[6,] -5.67486663 -1811.140 1799.790 1805.465
[7,] -5.53840932 -1811.336 1800.259 1805.798
[8,] 20.36609880 -1785.637 1826.369 1806.003
[9,] -14.88094661 -1821.103 1791.341 1806.222
[10,] -0.09496698 -1807.907 1807.717 1807.812
$Y
fcst lower upper CI
[1,] -4.179815e-07 -3.857079e-05 3.773483e-05 3.815281e-05
[2,] 1.886498e-06 -3.654026e-05 4.031326e-05 3.842676e-05
[3,] -1.274244e-06 -3.972669e-05 3.717821e-05 3.845245e-05
[4,] -2.815922e-07 -3.875275e-05 3.818957e-05 3.847116e-05
[5,] -7.214299e-07 -3.919800e-05 3.775514e-05 3.847657e-05
[6,] -7.825950e-07 -3.925958e-05 3.769439e-05 3.847698e-05
[7,] 1.102569e-07 -3.837391e-05 3.859442e-05 3.848416e-05
[8,] -2.166975e-07 -3.870116e-05 3.826777e-05 3.848447e-05
[9,] -5.326749e-08 -3.854202e-05 3.843549e-05 3.848875e-05
[10,] -6.753158e-07 -3.918076e-05 3.783013e-05 3.850545e-05
> predict(TSvar96)
$X
fcst lower upper CI
[1,] 33.770859 -1592.898 1660.440 1626.669
[2,] -63.865710 -1833.106 1705.375 1769.241
[3,] 77.636350 -1716.435 1871.708 1794.072
[4,] 14.859046 -1783.052 1812.770 1797.911
[5,] -94.458953 -1895.953 1707.035 1801.494
[6,] 6.612096 -1797.559 1810.783 1804.171
[7,] 55.755082 -1748.780 1860.290 1804.535
[8,] 40.833596 -1763.947 1845.614 1804.780
[9,] -74.713039 -1879.741 1730.315 1805.028
[10,] -10.848328 -1817.510 1795.813 1806.661
$Y
fcst lower upper CI
[1,] 1.514097e-06 -3.660781e-05 3.963600e-05 3.812191e-05
[2,] 3.693606e-06 -3.470163e-05 4.208884e-05 3.839524e-05
[3,] -1.648351e-06 -4.006920e-05 3.677250e-05 3.842085e-05
[4,] 1.439977e-06 -3.699968e-05 3.987963e-05 3.843965e-05
[5,] -2.497274e-06 -4.094177e-05 3.594722e-05 3.844450e-05
[6,] -2.431337e-07 -3.868807e-05 3.820180e-05 3.844493e-05
[7,] -4.750051e-07 -3.892673e-05 3.797672e-05 3.845172e-05
[8,] -5.696669e-07 -3.902169e-05 3.788235e-05 3.845202e-05
[9,] 1.452356e-07 -3.831102e-05 3.860149e-05 3.845625e-05
[10,] -6.050617e-07 -3.907691e-05 3.786678e-05 3.847184e-05
X の予測値は似てるけど、Y は全然違いますね。係数行列を見てみても、高次のパラメータがあまり減衰していません。
残差の自己相関は無くなっていますが、モデル化としてはあまり良くない結果に思えます。状態空間モデルの選択
dse パッケージの estBlackBox() を用いると、状態空間モデルの選択をしてくれます。状態空間モデルでは、時系列を入力と出力に分けて考える必要があります。dse では、TSdata というオブジェクトで入出力のある時系列データを保持することができるので、まず TS から TSdata を作りましょう。今回の分析対象である TS は、実は Y が入力で X は応答であることを期待しています。以下のように input, output を指定して変換します。> DSETS <- TSdata(input=TS$Y, output=TS$X)
そして、以下のようにモデルの選択を実行します。
> DSEEST <- estBlackBox(DSETS)
> summary(DSEEST)
neg. log likelihood = 697537.7 sample length = 85496
X
RMSE 845.3734
innovations form state space: nested model a la Mittnik
inputs : Y
outputs: X
input dimension = 1 state dimension = 12 output dimension = 1
theoretical parameter space dimension = 36
180 actual parameters 0 non-zero constants
Initial values not specified.
上記のようなモデルが選択されました。しかし、残差の ACF, PACF を見てみると、かなり強い自己相関が残っていて、モデル化できたとは言えませんでした。残念。