1. Tải dữ liệu
suppressMessages(library(vars))
suppressMessages(library(dplyr))
library(thongke.club)
suppressMessages(library(aTSA))
dulieu <-price_tom
dulieu <-as.data.frame(dulieu)
# Xét thời gian, khi không xét có thể dùng dulieu$bien
ao <-ts(dulieu$price_ao,start=(c(1993,1)),frequency=4)
bs <-ts(dulieu$price_bs,start=(c(1993,1)),frequency=4)
bl <-ts(dulieu$price_bl,start=(c(1993,1)),frequency=4)
xk <-ts(dulieu$price_xk,start=(c(1993,1)),frequency=4)
dulieu <-cbind(ao,bs,bl,xk)
plot(dulieu)
2. Unit root test
adf.test(ao)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -0.950 0.338
## [2,] 1 -0.584 0.470
## [3,] 2 -0.390 0.531
## [4,] 3 -0.379 0.535
## [5,] 4 -0.442 0.517
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -9.97 0.01
## [2,] 1 -8.05 0.01
## [3,] 2 -6.39 0.01
## [4,] 3 -5.79 0.01
## [5,] 4 -5.69 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -9.99 0.01
## [2,] 1 -8.09 0.01
## [3,] 2 -6.41 0.01
## [4,] 3 -5.78 0.01
## [5,] 4 -5.65 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
adf.test(diff(ao))
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -16.61 0.01
## [2,] 1 -13.17 0.01
## [3,] 2 -10.43 0.01
## [4,] 3 -8.58 0.01
## [5,] 4 -7.05 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -16.54 0.01
## [2,] 1 -13.11 0.01
## [3,] 2 -10.38 0.01
## [4,] 3 -8.54 0.01
## [5,] 4 -7.02 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -16.47 0.01
## [2,] 1 -13.06 0.01
## [3,] 2 -10.35 0.01
## [4,] 3 -8.53 0.01
## [5,] 4 -7.01 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
kpss.test(bs)
## KPSS Unit Root Test
## alternative: nonstationary
##
## Type 1: no drift no trend
## lag stat p.value
## 2 1.05 0.1
## -----
## Type 2: with drift no trend
## lag stat p.value
## 2 0.0402 0.1
## -----
## Type 1: with drift and trend
## lag stat p.value
## 2 0.039 0.1
## -----------
## Note: p.value = 0.01 means p.value <= 0.01
## : p.value = 0.10 means p.value >= 0.10
kpss.test(diff(bs))
## KPSS Unit Root Test
## alternative: nonstationary
##
## Type 1: no drift no trend
## lag stat p.value
## 2 0.0236 0.1
## -----
## Type 2: with drift no trend
## lag stat p.value
## 2 0.021 0.1
## -----
## Type 1: with drift and trend
## lag stat p.value
## 2 0.0201 0.1
## -----------
## Note: p.value = 0.01 means p.value <= 0.01
## : p.value = 0.10 means p.value >= 0.10
ur.ers(bl, type="P-test", model="trend")
##
## ###############################################################
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: 0.3958
ur.ers(diff(bl), type="P-test", model="trend")
##
## ###############################################################
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test #
## ###############################################################
##
## The value of the test statistic is: 11.2616
pp.test(xk)
## Phillips-Perron Unit Root Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag Z_rho p.value
## 4 -0.121 0.662
## -----
## Type 2: with drift no trend
## lag Z_rho p.value
## 4 -112 0.01
## -----
## Type 3: with drift and trend
## lag Z_rho p.value
## 4 -112 0.01
## ---------------
## Note: p-value = 0.01 means p.value <= 0.01
pp.test(diff(xk))
## Phillips-Perron Unit Root Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag Z_rho p.value
## 4 -140 0.01
## -----
## Type 2: with drift no trend
## lag Z_rho p.value
## 4 -140 0.01
## -----
## Type 3: with drift and trend
## lag Z_rho p.value
## 4 -140 0.01
## ---------------
## Note: p-value = 0.01 means p.value <= 0.01
3. Test xu hướng và mùa
stationary.test(ao)
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -0.950 0.338
## [2,] 1 -0.584 0.470
## [3,] 2 -0.390 0.531
## [4,] 3 -0.379 0.535
## [5,] 4 -0.442 0.517
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -9.97 0.01
## [2,] 1 -8.05 0.01
## [3,] 2 -6.39 0.01
## [4,] 3 -5.79 0.01
## [5,] 4 -5.69 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -9.99 0.01
## [2,] 1 -8.09 0.01
## [3,] 2 -6.41 0.01
## [4,] 3 -5.78 0.01
## [5,] 4 -5.65 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
stationary.test(diff(ao))
## Augmented Dickey-Fuller Test
## alternative: stationary
##
## Type 1: no drift no trend
## lag ADF p.value
## [1,] 0 -16.61 0.01
## [2,] 1 -13.17 0.01
## [3,] 2 -10.43 0.01
## [4,] 3 -8.58 0.01
## [5,] 4 -7.05 0.01
## Type 2: with drift no trend
## lag ADF p.value
## [1,] 0 -16.54 0.01
## [2,] 1 -13.11 0.01
## [3,] 2 -10.38 0.01
## [4,] 3 -8.54 0.01
## [5,] 4 -7.02 0.01
## Type 3: with drift and trend
## lag ADF p.value
## [1,] 0 -16.47 0.01
## [2,] 1 -13.06 0.01
## [3,] 2 -10.35 0.01
## [4,] 3 -8.53 0.01
## [5,] 4 -7.01 0.01
## ----
## Note: in fact, p.value = 0.01 means p.value <= 0.01
trend.test(ao)
##
## Approximate Cox-Stuart trend test
##
## data: ao
## D+ = 34, p-value = 0.1508
## alternative hypothesis: data have a increasing trend
trend.test(bs)
##
## Approximate Cox-Stuart trend test
##
## data: bs
## D+ = 30, p-value = 0.3964
## alternative hypothesis: data have a increasing trend
trend.test(bl)
##
## Approximate Cox-Stuart trend test
##
## data: bl
## D+ = 31, p-value = 0.3481
## alternative hypothesis: data have a increasing trend
trend.test(xk)
##
## Approximate Cox-Stuart trend test
##
## data: xk
## D- = 30, p-value = 0.5
## alternative hypothesis: data have a decreasing trend
4. Tìm lag phù hợp
#type = c("const", "trend", "both", "none")
chonlag <- VARselect(dulieu, lag.max = 12, type = "const")
chonlag$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
chonlag <- VARselect(dulieu, lag.max = 12, type = "trend")
chonlag$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 1 1 3
chonlag <- VARselect(dulieu, lag.max = 12, type = "both")
chonlag$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
chonlag <- VARselect(dulieu, lag.max = 12, type = "none")
chonlag$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 1 1 3
5. Chạy VAR với lag phù hợp
chonvar <- VAR(dulieu, p = 1, type = "none", season = NULL, exog = NULL)
summary(chonvar)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: ao, bs, bl, xk
## Deterministic variables: none
## Sample size: 119
## Log Likelihood: -1997.765
## Roots of the characteristic polynomial:
## 0.9947 0.202 0.06619 0.06619
## Call:
## VAR(y = dulieu, p = 1, type = "none", exogen = NULL)
##
##
## Estimation results for equation ao:
## ===================================
## ao = ao.l1 + bs.l1 + bl.l1 + xk.l1
##
## Estimate Std. Error t value Pr(>|t|)
## ao.l1 -0.09326 0.13436 -0.694 0.48899
## bs.l1 0.40073 0.12036 3.330 0.00117 **
## bl.l1 0.33857 0.12062 2.807 0.00587 **
## xk.l1 0.20490 0.09338 2.194 0.03023 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 19.19 on 115 degrees of freedom
## Multiple R-Squared: 0.976, Adjusted R-squared: 0.9751
## F-statistic: 1167 on 4 and 115 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation bs:
## ===================================
## bs = ao.l1 + bs.l1 + bl.l1 + xk.l1
##
## Estimate Std. Error t value Pr(>|t|)
## ao.l1 -0.04539 0.12711 -0.357 0.72168
## bs.l1 0.21060 0.11386 1.850 0.06694 .
## bl.l1 0.34897 0.11411 3.058 0.00277 **
## xk.l1 0.29763 0.08834 3.369 0.00103 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 18.16 on 115 degrees of freedom
## Multiple R-Squared: 0.9784, Adjusted R-squared: 0.9776
## F-statistic: 1301 on 4 and 115 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation bl:
## ===================================
## bl = ao.l1 + bs.l1 + bl.l1 + xk.l1
##
## Estimate Std. Error t value Pr(>|t|)
## ao.l1 -0.11942 0.14024 -0.852 0.396220
## bs.l1 0.49793 0.12563 3.964 0.000129 ***
## bl.l1 0.28626 0.12590 2.274 0.024835 *
## xk.l1 0.33072 0.09747 3.393 0.000949 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 20.03 on 115 degrees of freedom
## Multiple R-Squared: 0.9813, Adjusted R-squared: 0.9806
## F-statistic: 1506 on 4 and 115 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation xk:
## ===================================
## xk = ao.l1 + bs.l1 + bl.l1 + xk.l1
##
## Estimate Std. Error t value Pr(>|t|)
## ao.l1 -0.07994 0.15062 -0.531 0.596632
## bs.l1 0.56863 0.13493 4.214 5.01e-05 ***
## bl.l1 0.30586 0.13522 2.262 0.025582 *
## xk.l1 0.38502 0.10469 3.678 0.000359 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 21.52 on 115 degrees of freedom
## Multiple R-Squared: 0.9844, Adjusted R-squared: 0.9838
## F-statistic: 1809 on 4 and 115 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## ao bs bl xk
## ao 368.0 187.2 254.1 227.5
## bs 187.2 328.8 197.6 192.4
## bl 254.1 197.6 400.6 212.8
## xk 227.5 192.4 212.8 462.1
##
## Correlation matrix of residuals:
## ao bs bl xk
## ao 1.0000 0.5381 0.6619 0.5518
## bs 0.5381 1.0000 0.5445 0.4937
## bl 0.6619 0.5445 1.0000 0.4946
## xk 0.5518 0.4937 0.4946 1.0000
6. Các kiểm định phù hợp
# serial correlation test
serial.test(chonvar, lags.pt = 12, type = "PT.asymptotic") -> testchuoi
testchuoi
##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object chonvar
## Chi-squared = 196.11, df = 176, p-value = 0.1426
plot(testchuoi, names = "price_ao")
# Test phương sai thay đôi
#arch.test(chonvar, lags.multi = 12, multivariate.only = TRUE)
# Test phân phối chuẩn
normality.test(chonvar, multivariate.only = TRUE)
## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object chonvar
## Chi-squared = 7.1548, df = 8, p-value = 0.52
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object chonvar
## Chi-squared = 5.208, df = 4, p-value = 0.2666
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object chonvar
## Chi-squared = 1.9469, df = 4, p-value = 0.7455
# Tính ổn định
stability(chonvar, type = "OLS-CUSUM") -> bv.cusum
plot(bv.cusum)
7. Chạy Granger causality
causality(chonvar, cause="ao")
## $Granger
##
## Granger causality H0: ao do not Granger-cause bs bl xk
##
## data: VAR object chonvar
## F-Test = 0.2572, df1 = 3, df2 = 460, p-value = 0.8562
##
##
## $Instant
##
## H0: No instantaneous causality between: ao and bs bl xk
##
## data: VAR object chonvar
## Chi-squared = 40.942, df = 3, p-value = 6.726e-09
causality(chonvar, cause = "bs")
## $Granger
##
## Granger causality H0: bs do not Granger-cause ao bl xk
##
## data: VAR object chonvar
## F-Test = 7.4747, df1 = 3, df2 = 460, p-value = 6.783e-05
##
##
## $Instant
##
## H0: No instantaneous causality between: bs and ao bl xk
##
## data: VAR object chonvar
## Chi-squared = 33.334, df = 3, p-value = 2.738e-07
causality(chonvar, cause="bl")
## $Granger
##
## Granger causality H0: bl do not Granger-cause ao bs xk
##
## data: VAR object chonvar
## F-Test = 3.7904, df1 = 3, df2 = 460, p-value = 0.01047
##
##
## $Instant
##
## H0: No instantaneous causality between: bl and ao bs xk
##
## data: VAR object chonvar
## Chi-squared = 39.566, df = 3, p-value = 1.317e-08
causality(chonvar, cause = "xk")
## $Granger
##
## Granger causality H0: xk do not Granger-cause ao bs bl
##
## data: VAR object chonvar
## F-Test = 5.0544, df1 = 3, df2 = 460, p-value = 0.001873
##
##
## $Instant
##
## H0: No instantaneous causality between: xk and ao bs bl
##
## data: VAR object chonvar
## Chi-squared = 32.238, df = 3, p-value = 4.663e-07
8. Chạy impulse response
xkao <- irf(chonvar, impulse = "xk", response = "ao", n.ahead = 20, boot = TRUE)
plot(xkao, ylab = "Cú shocks", main = "Biến động giá tại ao khi giá xuất khẩu biến động")
9. Phân rã phương sai
phanra <- fevd(chonvar, n.ahead = 10)
phanra
## $ao
## ao bs bl xk
## [1,] 1.0000000 0.0000000 0.00000000 0.00000000
## [2,] 0.7952632 0.1330761 0.05009221 0.02156853
## [3,] 0.7061878 0.1800110 0.06480137 0.04899988
## [4,] 0.6526387 0.2101373 0.07248792 0.06473609
## [5,] 0.6172104 0.2297801 0.07778509 0.07522434
## [6,] 0.5919124 0.2438575 0.08153818 0.08269191
## [7,] 0.5729681 0.2543901 0.08435378 0.08828801
## [8,] 0.5582474 0.2625761 0.08654069 0.09263576
## [9,] 0.5464809 0.2691191 0.08828890 0.09611111
## [10,] 0.5368608 0.2744686 0.08971819 0.09895250
##
## $bs
## ao bs bl xk
## [1,] 0.2902056 0.7097944 0.00000000 0.00000000
## [2,] 0.3558199 0.5312280 0.06289827 0.05005381
## [3,] 0.3739703 0.4846302 0.07149094 0.06990855
## [4,] 0.3870915 0.4516933 0.07881188 0.08240337
## [5,] 0.3952158 0.4312909 0.08318067 0.09031258
## [6,] 0.4009914 0.4167871 0.08631809 0.09590339
## [7,] 0.4052575 0.4060737 0.08863002 0.10003875
## [8,] 0.4085468 0.3978134 0.09041359 0.10322622
## [9,] 0.4111584 0.3912549 0.09182950 0.10575716
## [10,] 0.4132824 0.3859210 0.09298109 0.10781553
##
## $bl
## ao bs bl xk
## [1,] 0.4386069 0.05016423 0.5112289 0.00000000
## [2,] 0.4183064 0.19532542 0.3392070 0.04716122
## [3,] 0.4244590 0.22922799 0.2751183 0.07119475
## [4,] 0.4265470 0.25194452 0.2378971 0.08361132
## [5,] 0.4281062 0.26581607 0.2144270 0.09165070
## [6,] 0.4291459 0.27554040 0.1981081 0.09720559
## [7,] 0.4299165 0.28266051 0.1861360 0.10128706
## [8,] 0.4305048 0.28811259 0.1769727 0.10440982
## [9,] 0.4309698 0.29241851 0.1697351 0.10687656
## [10,] 0.4313463 0.29590563 0.1638739 0.10887414
##
## $xk
## ao bs bl xk
## [1,] 0.3050816 0.0547969 0.01168467 0.6284369
## [2,] 0.3596398 0.2022823 0.03965625 0.3984217
## [3,] 0.3839973 0.2373541 0.06109845 0.3175502
## [4,] 0.3955339 0.2595431 0.07096356 0.2739594
## [5,] 0.4029221 0.2727152 0.07728507 0.2470776
## [6,] 0.4079149 0.2817807 0.08155318 0.2287513
## [7,] 0.4115395 0.2883333 0.08465255 0.2154747
## [8,] 0.4142857 0.2933029 0.08700061 0.2054109
## [9,] 0.4164390 0.2971987 0.08884182 0.1975205
## [10,] 0.4181725 0.3003352 0.09032405 0.1911682
plot(phanra)
10. Dự báo
uocluong <- predict(chonvar, n.ahead = 10, ci = 0.95)
plot(uocluong, names = "bl")
plot(uocluong, names = "xk")
fanchart(uocluong, names = "ao")
fanchart(uocluong, names = "bs")
