[R]R言語で多変量時系列解析
Rで、多変量時系列解析であるベクトル自己回帰（VAR）をやってみる。
とりあえず、動かしてみる
> library(vars) > data("Canada") > library(vars) > data("Canada") > print(Canada[1:5,]) prod e U rw [1,] 405.3665 929.6105 7.53 386.1361 [2,] 404.6398 929.8040 7.70 388.1358 [3,] 403.8149 930.3184 7.47 390.5401 [4,] 404.2158 931.4277 7.27 393.9638 [5,] 405.0467 932.6620 7.37 396.7647 > summary(Canada) e prod rw U Min. :928.6 Min. :401.3 Min. :386.1 Min. : 6.700 1st Qu.:935.4 1st Qu.:404.8 1st Qu.:423.9 1st Qu.: 7.782 Median :946.0 Median :406.5 Median :444.4 Median : 9.450 Mean :944.3 Mean :407.8 Mean :440.8 Mean : 9.321 3rd Qu.:950.0 3rd Qu.:410.7 3rd Qu.:461.1 3rd Qu.:10.607 Max. :961.8 Max. :418.0 Max. :470.0 Max. :12.770
プロット
> plot(Canada, nc = 2, xlab = "")
ADF検定で単位根検定
系列データと差分系列データの単位根検定を行う。
- 系列データ
> adf1 <- summary(ur.df(Canada[, "prod"], type = "trend", lags = 2)) > adf1 ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression trend Call: lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.19924 -0.38994 0.04294 0.41914 1.71660 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.415228 15.309403 1.987 0.0506 . z.lag.1 -0.075791 0.038134 -1.988 0.0505 . tt 0.013896 0.006422 2.164 0.0336 * z.diff.lag1 0.284866 0.114359 2.491 0.0149 * z.diff.lag2 0.080019 0.116090 0.689 0.4927 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6851 on 76 degrees of freedom Multiple R-squared: 0.1354, Adjusted R-squared: 0.08993 F-statistic: 2.976 on 4 and 76 DF, p-value: 0.02438 Value of test-statistic is: -1.9875 2.3 2.3817 Critical values for test statistics: 1pct 5pct 10pct tau3 -4.04 -3.45 -3.15 phi2 6.50 4.88 4.16 phi3 8.73 6.49 5.47
- 差分系列データ
> adf2 <- summary(ur.df(diff(Canada[, "prod"]), type = "drift", lags=1)) > adf2 ############################################### # Augmented Dickey-Fuller Test Unit Root Test # ############################################### Test regression drift Call: lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag) Residuals: Min 1Q Median 3Q Max -2.05124 -0.39530 0.07819 0.41109 1.75129 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.11534 0.08029 1.437 0.155 z.lag.1 -0.68893 0.13350 -5.160 1.83e-06 *** z.diff.lag -0.04274 0.11275 -0.379 0.706 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6971 on 78 degrees of freedom Multiple R-squared: 0.3615, Adjusted R-squared: 0.3451 F-statistic: 22.08 on 2 and 78 DF, p-value: 2.526e-08 Value of test-statistic is: -5.1604 13.3184 Critical values for test statistics: 1pct 5pct 10pct tau2 -3.51 -2.89 -2.58 phi1 6.70 4.71 3.86
同様に他の系列データ、差分系列データについても行う
VARモデル選択
VARSelect関数で、それぞれの指標にあわせた最適なラグ次数を算出することが出来る。
> VARselect(Canada, lag.max = 8, type = "both") $selection AIC(n) HQ(n) SC(n) FPE(n) 3 2 1 3 $criteria 1 2 3 4 5 AIC(n) -6.272579064 -6.636669705 -6.771176872 -6.634609210 -6.398132246 HQ(n) -5.978429449 -6.146420347 -6.084827770 -5.752160366 -5.319583658 SC(n) -5.536558009 -5.409967947 -5.053794411 -4.426546046 -3.699388378 FPE(n) 0.001889842 0.001319462 0.001166019 0.001363175 0.001782055 6 7 8 AIC(n) -6.307704843 -6.070727259 -6.06159685 HQ(n) -5.033056512 -4.599979185 -4.39474903 SC(n) -3.118280272 -2.390621985 -1.89081087 FPE(n) 0.002044202 0.002768551 0.00306012
VARモデル構築
> Canada <- Canada[, c("prod", "e", "U", "rw")] > p1ct <- VAR(Canada, p = 1, type = "both") > p1ct VAR Estimation Results: ======================= Estimated coefficients for equation prod: ========================================= Call: prod = prod.l1 + e.l1 + U.l1 + rw.l1 + const + trend prod.l1 e.l1 U.l1 rw.l1 const trend 0.96313671 0.01291155 0.21108918 -0.03909399 16.24340747 0.04613085 Estimated coefficients for equation e: ====================================== Call: e = prod.l1 + e.l1 + U.l1 + rw.l1 + const + trend prod.l1 e.l1 U.l1 rw.l1 const 0.19465028 1.23892283 0.62301475 -0.06776277 -278.76121138 trend -0.04066045 Estimated coefficients for equation U: ====================================== Call: U = prod.l1 + e.l1 + U.l1 + rw.l1 + const + trend prod.l1 e.l1 U.l1 rw.l1 const -0.12319201 -0.24844234 0.39158002 0.06580819 259.98200967 trend 0.03451663 Estimated coefficients for equation rw: ======================================= Call: rw = prod.l1 + e.l1 + U.l1 + rw.l1 + const + trend prod.l1 e.l1 U.l1 rw.l1 const -0.22308744 -0.05104397 -0.36863956 0.94890946 163.02453066 trend 0.07142229
eを求める式について概要を出力
> summary(p1ct, equation = "e") VAR Estimation Results: ========================= Endogenous variables: prod, e, U, rw Deterministic variables: both Sample size: 83 Log Likelihood: -207.525 Roots of the characteristic polynomial: 0.9504 0.9504 0.9045 0.7513 Call: VAR(y = Canada, p = 1, type = "both") Estimation results for equation e: ================================== e = prod.l1 + e.l1 + U.l1 + rw.l1 + const + trend Estimate Std. Error t value Pr(>|t|) prod.l1 0.19465 0.03612 5.389 7.5e-07 *** e.l1 1.23892 0.08632 14.353 < 2e-16 *** U.l1 0.62301 0.16927 3.681 0.000430 *** rw.l1 -0.06776 0.02828 -2.396 0.018991 * const -278.76121 75.18295 -3.708 0.000392 *** trend -0.04066 0.01970 -2.064 0.042378 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.4701 on 77 degrees of freedom Multiple R-Squared: 0.9975, Adjusted R-squared: 0.9973 F-statistic: 6088 on 5 and 77 DF, p-value: < 2.2e-16 Covariance matrix of residuals: prod e U rw prod 0.469517 0.06767 -0.04128 0.002141 e 0.067667 0.22096 -0.13200 -0.082793 U -0.041280 -0.13200 0.12161 0.063738 rw 0.002141 -0.08279 0.06374 0.593174 Correlation matrix of residuals: prod e U rw prod 1.000000 0.2101 -0.1728 0.004057 e 0.210085 1.0000 -0.8052 -0.228688 U -0.172753 -0.8052 1.0000 0.237307 rw 0.004057 -0.2287 0.2373 1.000000
プロットしてみる
> plot(p1ct, names = "e")
参考
Rで学ぶ回帰分析と単位根検定
View more presentations from teramo nagi
