時系列データに関する基本的な検定について調べてみました。 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