【R】時系列データに関する基本的な検定

時系列データに関する基本的な検定について調べてみました。 Granger の因果性検定は 【R】金利と日経平均株価の関係【Granger因果】 で扱ったので今回は省略しました。

環境は macOS 10.12, R 3.3.1 です。

はじめに

架空の日次PVデータを例に以下の検定を行ってみる。

名称 検定内容
ADF検定 単位根検定
PP検定 単位根検定
Engle-Granger検定 共和分検定
Phillips-Ouliaris検定 共和分検定
連検定 非独立性の検定
BDS検定 非独立性の検定
Jarque–Bera検定 正規性の検定
Ljung-Box検定 自己相関の検定

架空の日次PVデータ。

ADF検定

多くの時系列モデルでは定常過程を前提としているので, 時系列データに対してまず最初に単位根を確認することが多い。

ADF検定 (Augmented Dickey-Fuller test) は単位根検定の手法で DF検定 (Dickey–Fuller test) を拡張したものである。単位根検定は帰無仮説を単位根過程, 対立仮説を定常過程とした検定である。単位根過程は原系列が非定常過程で差分系列が定常過程の場合である。単位根過程は平均回帰性を持たないため t検定 が使えない。従って, ブラウン運動を用いた漸近分布 (DF分布) が検定統計量として使われる。
DF検定では AR(1) を真のモデルと仮定しているがこの仮定を緩めて自己相関を取り除く AR(p) まで許したのが ADF検定となる。

ADF検定は tseries::adf.test() で行う。 default の AR次数 は trunc((length(x)-1)^(1/3)) となっている。

> require(tseries)
> df <- read.csv("data/pv.csv", header = T)
> adf.test(df$pv)

	Augmented Dickey-Fuller Test

data:  df$pv
Dickey-Fuller = -3.9709, Lag order = 7, p-value = 0.01073
alternative hypothesis: stationary

> adf.test(seq(1,100))

	Augmented Dickey-Fuller Test

data:  seq(1, 100)
Dickey-Fuller = 1.7321, Lag order = 4, p-value = 0.99
alternative hypothesis: stationary

pvデータに対して有意水準5%で単位根過程であるという帰無仮説は棄却された。

非定常過程を定常過程に変換する方法として, 差分系列・対数差分系列・移動平均を取る方法などがある。 また, 経済・ファイナンスデータでは収益率が使われることが多い。

PP検定

PP検定 (Phillips–Perron test) は単位根検定でADF検定が AR(p) であるのに対してPP検定は自己相関や分散不均一性を持つことまで許す。

PP検定は tseries::pp.test() で行う。

> pp.test(df$pv)

	Phillips-Perron Unit Root Test

data:  df$pv
Dickey-Fuller Z(alpha) = -125.84, Truncation lag parameter = 5,
p-value = 0.01
alternative hypothesis: stationary

 警告メッセージ: 
 pp.test(df$pv) で:  p-value smaller than printed p-value

Engle-Granger検定

x_t と y_t が単位根過程で ax_t + by_t が定常過程となるような a と b が存在する場合, x_t と y_t の間には共和分 (cointegration) の関係があるという。つまり, 2つの単位根過程の線形和が定常過程になることを指す。

Engle-Grangerの共和分検定は, 線形回帰モデルの残差にADF検定を用いて行う。

> model <- lm(pv ~ ad + vacation, data = df)

> summary(model)

Call:
lm(formula = pv ~ ad + vacation, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-401.23 -184.76   44.24  120.24  612.24 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  446.755      9.188  48.622  < 2e-16 ***
ad            14.837      1.840   8.065 1.09e-14 ***
vacation    -199.612     38.044  -5.247 2.64e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 169.2 on 362 degrees of freedom
Multiple R-squared:  0.2081,	Adjusted R-squared:  0.2037 
F-statistic: 47.56 on 2 and 362 DF,  p-value: < 2.2e-16

> model.resid <- resid(model)
> plot(model.resid, type = "l")
> adf.test(model.resid)

	Augmented Dickey-Fuller Test

data:  model.resid
Dickey-Fuller = -3.9109, Lag order = 7, p-value = 0.01373
alternative hypothesis: stationary

有意水準5%で共和分がないという帰無仮説は棄却された。

Phillips-Ouliaris検定

Phillips-Ouliarisの共和分検定は urca::ca.po() で行う。

> require(urca)
> x <- data.frame(df$pv, df$ad, df$vacation)
> 
> pu <- ca.po(x, demean = "const", type = "Pu")
> 
> summary(pu)

######################################## 
# Phillips and Ouliaris Unit Root Test # 
######################################## 

Test of type Pu 
detrending of series with constant only 


Call:
lm(formula = z[, 1] ~ z[, -1])

Residuals:
    Min      1Q  Median      3Q     Max 
-401.23 -184.76   44.24  120.24  612.24 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)         446.755      9.188  48.622  < 2e-16 ***
z[, -1]df.ad         14.837      1.840   8.065 1.09e-14 ***
z[, -1]df.vacation -199.612     38.044  -5.247 2.64e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 169.2 on 362 degrees of freedom
Multiple R-squared:  0.2081,	Adjusted R-squared:  0.2037 
F-statistic: 47.56 on 2 and 362 DF,  p-value: < 2.2e-16


Value of test-statistic is: 137.5592 

Critical values of Pu are:
                  10pct    5pct    1pct
critical values 33.6955 40.5252 53.8731

連検定

連検定 (Run test) は2値データ系列のランダム性を検定する。連は同じ値の連続性のことで, ランダム(=独立)なら連の個数はサイズNに応じた平均と分散を持つ正規分布に従う。データの平均値が理論値から外れるほどランダムではないとする。

連検定は tseries::runs.test() で行う。

今回は原系列の差分を取り各値の符号から 1, -1 に2値化し連検定を行う。

> d <- diff(df$pv)

> x <- factor(sign(d[-which(d %in% 0)])) # binarize

> runs.test(x)

	Runs Test

data:  x
Standard Normal = -1.2316, p-value = 0.2181
alternative hypothesis: two.sided
> set.seed(123)

> rnd <- factor(sign(rnorm(1000, 0, 1)))  # randomness

> runs.test(rnd)

	Runs Test

data:  rnd
Standard Normal = 0.88913, p-value = 0.3739
alternative hypothesis: two.sided

pvデータの上下の変動がランダムであるという帰無仮説は棄却されなかった。

BDS検定

BDS検定 (W. A. Brock, W. Dechert and J. Scheinkman test) は時系列データの内部構造の偏り (ランダム性の有無) を検定する。偏りがある場合, 座標変換によって有形化するためデータは空間内の一部に局在し, ない場合(帰無仮説) はデータは空間全体に散らばる。

BDS検定は tseries::bds.test() で行う。

> bds.test(df$pv, m = 10)

	 BDS Test 

data:  df$pv 

Embedding dimension =   2  3  4  5  6  7  8  9 10 

Epsilon for close points =   94.7925 189.5851 284.3776 379.1701 

Standard Normal = 
       [ 94.7925 ] [ 189.5851 ] [ 284.3776 ] [ 379.1701 ]
[ 2 ]      26.2488      16.9095      12.0554      13.3231
[ 3 ]      34.7476      19.2542      11.0756      11.6600
[ 4 ]      46.0679      21.2769      11.3600      10.4822
[ 5 ]      70.4751      27.1976      12.1520       9.5538
[ 6 ]     124.7066      39.9101      12.2764       8.7043
[ 7 ]     218.3399      56.0526      14.2528       8.8265
[ 8 ]     408.5722      83.9699      19.0836       9.7265
[ 9 ]     773.5772     123.9684      24.1734      10.4397
[ 10 ]   1465.8982     183.4653      30.1376      10.9997

p-value = 
       [ 94.7925 ] [ 189.5851 ] [ 284.3776 ] [ 379.1701 ]
[ 2 ]            0            0            0            0
[ 3 ]            0            0            0            0
[ 4 ]            0            0            0            0
[ 5 ]            0            0            0            0
[ 6 ]            0            0            0            0
[ 7 ]            0            0            0            0
[ 8 ]            0            0            0            0
[ 9 ]            0            0            0            0
[ 10 ]           0            0            0            0

> rnd <- rnorm(1000, 0, 1)

> bds.test(rnd, m = 10)

	 BDS Test 

data:  rnd 

Embedding dimension =   2  3  4  5  6  7  8  9 10 

Epsilon for close points =  0.4958 0.9917 1.4875 1.9834 

Standard Normal = 
       [ 0.4958 ] [ 0.9917 ] [ 1.4875 ] [ 1.9834 ]
[ 2 ]     -0.2093    -0.8773    -0.5987    -0.4051
[ 3 ]     -0.4777    -1.2593    -0.8352    -0.7121
[ 4 ]     -0.6752    -1.0543    -0.5397    -0.5165
[ 5 ]     -1.1127    -1.3509    -0.7835    -0.8297
[ 6 ]      0.4928    -1.5359    -1.0823    -1.1382
[ 7 ]      0.6879    -1.6424    -1.3213    -1.3235
[ 8 ]     -0.3658    -1.8280    -1.5801    -1.5340
[ 9 ]      3.1748    -2.0609    -1.7665    -1.7532
[ 10 ]     2.8868    -2.2599    -1.9084    -1.9465

p-value = 
       [ 0.4958 ] [ 0.9917 ] [ 1.4875 ] [ 1.9834 ]
[ 2 ]      0.8342     0.3803     0.5494     0.6854
[ 3 ]      0.6328     0.2079     0.4036     0.4764
[ 4 ]      0.4996     0.2917     0.5894     0.6055
[ 5 ]      0.2658     0.1767     0.4333     0.4067
[ 6 ]      0.6222     0.1246     0.2791     0.2550
[ 7 ]      0.4915     0.1005     0.1864     0.1857
[ 8 ]      0.7146     0.0676     0.1141     0.1250
[ 9 ]      0.0015     0.0393     0.0773     0.0796
[ 10 ]     0.0039     0.0238     0.0563     0.0516

pvデータに対して内部構造の偏りがないという帰無仮説は棄却された。

Jarque–Bera検定

Jarque–Bera検定は正規分布に従う尖度と歪度を有しているかどうかを検定する。QQプロットによる可視化は見た目に分かりやすい利点があるが程度は数値化できない。

Jarque–Bera検定は tseries::jarque.bera.test() で行う。

今回は ARIMA(p,d,q) モデルの残差に対して Jarque–Bera検定 を行う。

> require(forecast)

> frame_files <- lapply(sys.frames(), function(x) x$ofile)

> frame_files <- Filter(Negate(is.null), frame_files)

> setwd(dirname(frame_files[[length(frame_files)]]))

> df <- read.csv("data/pv.csv", header = T)

> best <- auto.arima(
+     ts(df$pv),
+     max.p = 7, # AR degree
+     max.q = 7, # MA degree
+     ic = "aic",# model selection.
+     trace = T,  .... [TRUNCATED] 

 Fitting models using approximations to speed things up...

 ARIMA(2,0,2) with non-zero mean : 4693.874
 ARIMA(0,0,0) with non-zero mean : 4867.555
 ARIMA(1,0,0) with non-zero mean : 4750.954
 ARIMA(0,0,1) with non-zero mean : 4701.177
 ARIMA(0,0,0) with zero mean     : 5548.315
 ARIMA(1,0,2) with non-zero mean : 4693.331
 ARIMA(1,0,1) with non-zero mean : 4696.207
 ARIMA(1,0,3) with non-zero mean : 4692.652
 ARIMA(0,0,2) with non-zero mean : 4691.878
 ARIMA(0,0,2) with zero mean     : 5002.88
 ARIMA(0,0,3) with non-zero mean : 4690.069
 ARIMA(1,0,4) with non-zero mean : 4694.631
 ARIMA(0,0,3) with zero mean     : 4959.851
 ARIMA(0,0,4) with non-zero mean : 4691.273

 Now re-fitting the best model(s) without approximations...

 ARIMA(0,0,3) with non-zero mean : 4690.77

 Best model: ARIMA(0,0,3) with non-zero mean 


> summary(best)
Series: ts(df$pv) 
ARIMA(0,0,3) with non-zero mean 

Coefficients:
         ma1     ma2      ma3      mean
      0.7789  0.1217  -0.1181  443.7071
s.e.  0.0529  0.0683   0.0595   13.7238

sigma^2 estimated as 21921:  log likelihood=-2340.39
AIC=4690.77   AICc=4690.94   BIC=4710.27

Training set error measures:
                     ME     RMSE      MAE       MPE     MAPE      MASE
Training set 0.03101651 147.2451 111.2352 -15.06934 33.71827 0.8765312
                    ACF1
Training set 0.008160298

> plot(model)

> res <- model$res[-1]

> # qq-plot
> qqnorm(res)

> qqline(res)

> jarque.bera.test(res)

	Jarque Bera Test

data:  res
X-squared = 105.03, df = 2, p-value < 2.2e-16

残差の分布が歪度ゼロ, 尖度ゼロ (正規分布) であるという帰無仮説は棄却された。

Ljung-Box検定

Ljung-Box検定は時系列の自己相関の有無を検定する。帰無仮説は前後のデータが独立していること。

Ljung-Box検定は stats::Box.test() で行える。 Jarque–Bera検定同様, ARIMA(p,d,q) モデルの残差に対して検定を行う。

> k <- trunc((length(res)-1)^(1/3))

> Box.test(res, lag = k, type = "Ljung-Box", fitdf = 0)

	Box-Ljung test

data:  res
X-squared = 105.57, df = 7, p-value < 2.2e-16


> set.seed(123)

> x <- rnorm (1000)

> k <- trunc((length(x)-1)^(1/3))

> Box.test (x, lag = k, type = "Ljung-Box")

	Box-Ljung test

data:  x
X-squared = 4.7818, df = 9, p-value = 0.8529

残差に自己相関が無いという帰無仮説は棄却された。 従ってモデリングに改善の余地がありそうだ。

残差分析では正規性や自己相関以外にも, 平均ゼロ・分散一定であることを確認する。

Code は GitHub に置いた。

おわりに

単位根検定, 共和分検定については以下の書籍が詳しいです。


[1] 定常過程: 時系列の平均・分散が時刻tによらず一定のため全体で見ればその時系列の内部構造は変わらない性質. 平均は時間平均でなくアンサンブル平均.
[2] R による共和分分析
[3] ノンパラメトリック検定
[4] Rで計量時系列分析:単位根過程、見せかけの回帰、共和分、ベクトル誤差修正モデル
[5] Thoughts on the Ljung-Box test