Этот текст является вольным пересказом работы B. Pfaff (2004)
Подключаем необходимые пакеты:
library(car)
library(tseries)
library(zoo)
library(vars)
Данная работа повторяет регрессионный анализ рынка труда в Канаде за период 1980-2000 гг., проведенное в работе Breitung, et.al (2004). Для начала загрузим необходимые данные и посмотрим на них:
data("Canada")
summary(Canada)
## e prod rw U
## Min. :929 Min. :401 Min. :386 Min. : 6.70
## 1st Qu.:935 1st Qu.:405 1st Qu.:424 1st Qu.: 7.78
## Median :946 Median :406 Median :444 Median : 9.45
## Mean :944 Mean :408 Mean :441 Mean : 9.32
## 3rd Qu.:950 3rd Qu.:411 3rd Qu.:461 3rd Qu.:10.61
## Max. :962 Max. :418 Max. :470 Max. :12.77
Здесь
e — логарифм занятостиprod — производительность трудаrw — реальная заработная платаU — безработицаplot(Canada, nc = 2, col = "blue", xlab = "", mar = rep(2, 4), oma = rep(0.5,
4))
Для оценки теста Дикки-Фуллера в R используется следующая команда из пакета urca :
adf1 <- summary(ur.df(Canada[, "prod"], type = "trend", lag = 10, selectlags = "AIC"))
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.2253 -0.4163 0.0389 0.3785 1.6823
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.45958 15.32762 1.66 0.101
## z.lag.1 -0.06334 0.03826 -1.66 0.102
## tt 0.01161 0.00715 1.62 0.109
## z.diff.lag 0.30295 0.11931 2.54 0.013 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.673 on 69 degrees of freedom
## Multiple R-squared: 0.105, Adjusted R-squared: 0.0657
## F-statistic: 2.69 on 3 and 69 DF, p-value: 0.0531
##
##
## Value of test-statistic is: -1.656 2.119 1.454
##
## 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
Однако, читателю неинтересно такое нагромождение таблиц. Поэтому напишем функцию для компактного представления результатов теста ДФ в одной таблице:
adf.table <- function(adf.mult) {
i <- 0
name <- names(adf.mult)
Tab <- matrix(0, length(name), 7)
for (n in name) {
i <- i + 1
r.mult <- adf.mult[n]
names(r.mult) <- "a"
cval <- r.mult$a@cval[1, ]
test <- r.mult$a@teststat[[1]]
lags <- r.mult$a@lags
type <- r.mult$a@model
Tab[i, ] <- c(n, type, lags, round(test, 3), cval)
}
D <- data.frame(Tab)
names(D) <- c("Variable", "Det. term", "Lags", "Test Value", "1%_Cv", "5%_Cv",
"10%_Cv")
D
}
adf1.trend <- lapply(Canada, ur.df, type = "trend", lags = 12, selectlags = "BIC")
Tab <- adf.table(adf1.trend)
Tab
## Variable Det. term Lags Test Value 1%_Cv 5%_Cv 10%_Cv
## 1 e trend 12 -2.058 -4.04 -3.45 -3.15
## 2 prod trend 12 -1.594 -4.04 -3.45 -3.15
## 3 rw trend 12 -0.447 -4.04 -3.45 -3.15
## 4 U trend 12 -1.686 -4.04 -3.45 -3.15
Как видно из таблицы, все переменные имеют хотя бы один единичный корень. Протестируем первые разности:
adf2 <- lapply(diff(Canada), ur.df, type = "drift", lags = 12, selectlags = "BIC")
names(adf2) <- paste("d", names(adf2), sep = ".")
Tab2 <- adf.table(adf2)
Tab2
## Variable Det. term Lags Test Value 1%_Cv 5%_Cv 10%_Cv
## 1 d.e drift 12 -3.998 -3.51 -2.89 -2.58
## 2 d.prod drift 12 -4.524 -3.51 -2.89 -2.58
## 3 d.rw drift 12 -2.586 -3.51 -2.89 -2.58
## 4 d.U drift 12 -3.659 -3.51 -2.89 -2.58
\( VAR \) в своей самой базовой форме состоит из \( K \) энодогенных переменных \( \pmb{\mathbf{y_t}} = (y_{1t}, y_{2t}, \ldots , y_{Kt}) \). Тогда процесс \( VAR(p) \) выглядит так: \[ \pmb{\mathbf{y_t}} = A(1)\pmb{\mathbf{y_{t-1}}} + \ldots + A(p)\pmb{\mathbf{y_{t-1}}} + \varepsilon_{t} \] Итак, по результатам тестов на единичный корень, можно предположить интегрированность \( I(1) \) во всех переменных. Таким образом, будем оценивать \( VAR(p) \) в первых разностях. Прежде, чем оценивать сам \( VAR \), необходимо определить его порядок. Для этого используются критерии AIC, Ханна-Квина, SC и FPE:
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 6
## AIC(n) -6.27258 -6.636670 -6.771177 -6.634609 -6.398132 -6.307705
## HQ(n) -5.97843 -6.146420 -6.084828 -5.752160 -5.319584 -5.033057
## SC(n) -5.53656 -5.409968 -5.053794 -4.426546 -3.699388 -3.118280
## FPE(n) 0.00189 0.001319 0.001166 0.001363 0.001782 0.002044
## 7 8
## AIC(n) -6.070727 -6.06160
## HQ(n) -4.599979 -4.39475
## SC(n) -2.390622 -1.89081
## FPE(n) 0.002769 0.00306
Как видно из данной таблицы, в зависимости от выбора критерия, \( p \) может изменяться от 1 до 3. В работе были построены и протестированы все 3 варианта. Повторим их действия для \( p=1 \):
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.96314 0.01291 0.21109 -0.03909 16.24341 0.04613
##
##
## 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 trend
## 0.19465 1.23892 0.62301 -0.06776 -278.76121 -0.04066
##
##
## 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 trend
## -0.12319 -0.24844 0.39158 0.06581 259.98201 0.03452
##
##
## 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 trend
## -0.22309 -0.05104 -0.36864 0.94891 163.02453 0.07142
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.95 0.95 0.904 0.751
## 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.1947 0.0361 5.39 7.5e-07 ***
## e.l1 1.2389 0.0863 14.35 < 2e-16 ***
## U.l1 0.6230 0.1693 3.68 0.00043 ***
## rw.l1 -0.0678 0.0283 -2.40 0.01899 *
## const -278.7612 75.1830 -3.71 0.00039 ***
## trend -0.0407 0.0197 -2.06 0.04238 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.47 on 77 degrees of freedom
## Multiple R-Squared: 0.997, Adjusted R-squared: 0.997
## F-statistic: 6.09e+03 on 5 and 77 DF, p-value: <2e-16
##
##
##
## Covariance matrix of residuals:
## prod e U rw
## prod 0.46952 0.0677 -0.0413 0.00214
## e 0.06767 0.2210 -0.1320 -0.08279
## U -0.04128 -0.1320 0.1216 0.06374
## rw 0.00214 -0.0828 0.0637 0.59317
##
## Correlation matrix of residuals:
## prod e U rw
## prod 1.00000 0.210 -0.173 0.00406
## e 0.21008 1.000 -0.805 -0.22869
## U -0.17275 -0.805 1.000 0.23731
## rw 0.00406 -0.229 0.237 1.00000
plot(p1ct, nc = 1, names = "e")
e <- p1ct$varresult$e$residuals
qqnorm(e)
Данный тест использует следующую статистику: \[ LB test = T(T+2) \sum_{h=1}^{T/4} \frac{1}{(T-h)} trace\{ \hat{\Omega_h}' \hat{\Omega}^{-1} \hat{\Omega_h}' \hat{\Omega}^{-1} \} \sim \chi^{2}_{T/4-k+1} \] Нулевая гипотеза — отсутствие АК.
ser11 <- serial.test(p1ct, lags.pt = 16, type = "PT.asymptotic")
ser11$serial
##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 233.5, df = 240, p-value = 0.606
Для небольших выборок и для достаточно больших \( h \) следует использовать скорректированный тест LB: \[ LB test = T^{2} \sum_{h=1}^{T/4} \frac{1}{(T-h)} trace\{ \hat{\Omega_h}' \hat{\Omega}^{-1} \hat{\Omega_h}' \hat{\Omega}^{-1} \} \sim \chi^{2}_{T/4-k+1} \]
ser21 <- serial.test(p1ct, lags.pt = 16, type = "PT.adjusted")
ser21$serial
##
## Portmanteau Test (adjusted)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 256.9, df = 240, p-value = 0.2167
Для проверки АК используется также Breusch-Godfrey LM-тест, построенный на основе регрессии вида: \[ \hat{\varepsilon_t} = A_1 x_{t-1} + A_2 x_{t-2} + \ldots + A_k x_{t-k} + A_{\varepsilon} \hat{\varepsilon}_{t-h} + \tilde{\varepsilon_t} \] Нулевая гипотеза: \( A_{\varepsilon_t} = 0 \) Для каждого \( h \) лага остатков оцениваемая статистика выглядит следующим образом: \[ LM(h) = - (T-p(k+1) -\frac{1}{2})ln \left( \frac{| \hat{\Omega}(h) |}{| \hat{\Omega} |}\right) \sim \chi^{2}_{h K^{2}} \] В общем виде: \[ LM_h = T(K-trace( \tilde{\Omega}_{R}^{-1} \tilde{\Omega}_{UR})) \sim \chi^{2}_{h K^{2}} \] где \( \tilde{\Omega}_{R} \) — ковариационная матрица модели с ограничением, \tilde{\Omega}_{UR} — неограниченная модель:
ser31 <- serial.test(p1ct, lags.pt = 16, type = "BG")
ser31$serial
##
## Breusch-Godfrey LM test
##
## data: Residuals of VAR object p1ct
## Chi-squared = 118.3, df = 80, p-value = 0.003523
Edgerton and Shukur предложили коррекцию LM-теста для небольших выборок: \[ LMF_h = \frac{1-(1-R^{2})^{1/r}}{(1-R^{2})^{1/r}} \frac{Nr-q}{Km} \sim F(h K^{2}, int(Nr - q)) \]
ser41 <- serial.test(p1ct, lags.pt = 16, type = "ES")
ser41$serial
##
## Edgerton-Shukur F test
##
## data: Residuals of VAR object p1ct
## F statistic = 1.743, df1 = 80, df2 = 215, p-value = 0.000849
arch1 <- arch.test(p1ct, lags.multi = 5)
arch1$arch.mul
##
## ARCH (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 570.1, df = 500, p-value = 0.01606
plot(arch1, names = "e")
norm.t <- normality.test(p1ct)
norm.t
## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 9.919, df = 8, p-value = 0.2708
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 6.356, df = 4, p-value = 0.1741
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 3.563, df = 4, p-value = 0.4684
## $jb.mul
## $jb.mul$JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 9.919, df = 8, p-value = 0.2708
##
##
## $jb.mul$Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 6.356, df = 4, p-value = 0.1741
##
##
## $jb.mul$Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object p1ct
## Chi-squared = 3.563, df = 4, p-value = 0.4684
plot(norm.t)
plot(stability(p1ct), nc = 2)
Для того, чтобы оценить VECM модель, необходимо протестировать данные на коинтеграцию. Ниже представлен trace-test Йохансена: \[ \lambda_{trace} (r_0) = - T \sum_{j=r_0 + 1}^{k} log (1-\hat{\lambda}_{j}) \] где
summary(ca.jo(Canada, type = "trace", ecdet = "trend", K = 3, spec = "transitory"))
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 4.505e-01 1.963e-01 1.677e-01 4.647e-02 2.632e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.85 10.49 12.25 16.26
## r <= 2 | 18.72 22.76 25.32 30.45
## r <= 1 | 36.42 39.06 42.44 48.45
## r = 0 | 84.92 59.14 62.99 70.05
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.l1 1.00000 1.0000 1.0000 1.0000 1.0000
## e.l1 -0.02385 1.2947 -2.8832 4.2418 -8.2904
## U.l1 3.16875 3.4037 -7.4262 6.8414 -12.5578
## rw.l1 1.83528 -0.3331 1.3979 -0.1394 2.4466
## trend.l1 -1.30156 -0.2303 -0.5093 -1.5926 0.2831
##
## Weights W:
## (This is the loading matrix)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.d -0.006535 -0.02763 -0.070975 -0.014754 1.077e-11
## e.d -0.008503 0.11414 -0.008157 0.003988 7.400e-12
## U.d -0.004719 -0.06154 0.020719 -0.006557 -4.664e-12
## rw.d -0.046213 -0.14580 -0.016945 0.011896 6.952e-12
summary(ca.jo(Canada, type = "eigen", ecdet = "trend", K = 2, spec = "transitory"))
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 4.484e-01 2.324e-01 1.313e-01 4.878e-02 9.509e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 4.10 10.49 12.25 16.26
## r <= 2 | 11.54 16.85 18.96 23.65
## r <= 1 | 21.69 23.11 25.54 30.34
## r = 0 | 48.78 29.12 31.46 36.65
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.l1 1.0000 1.0000 1.00000 1.000 1.000
## e.l1 2.7132 -6.3190 0.49616 16.334 -10.369
## U.l1 8.8369 -15.2683 1.48063 25.774 -16.048
## rw.l1 -0.3716 3.1817 -0.04085 -2.546 4.927
## trend.l1 -0.4178 -0.9336 -0.26593 -3.414 -1.753
##
## Weights W:
## (This is the loading matrix)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.d 0.023156 -0.02833 -0.10915 -0.006296 -4.785e-13
## e.d 0.005602 -0.01739 0.08679 -0.001019 -4.386e-13
## U.d -0.019277 0.01382 -0.03696 -0.002277 4.920e-13
## rw.d -0.084619 -0.02739 -0.07798 0.003985 -1.032e-13
Так, гипотеза о равенстве \( r = 1 \) не отвергается на 10% уровне значимости.
vecm <- ca.jo(Canada[, c("prod", "e", "U", "rw")], type = "trace", ecdet = "trend",
K = 3, spec = "transitory")
summary(vecm)
##
## ######################
## # Johansen-Procedure #
## ######################
##
## Test type: trace statistic , with linear trend in cointegration
##
## Eigenvalues (lambda):
## [1] 4.505e-01 1.963e-01 1.677e-01 4.647e-02 2.632e-17
##
## Values of teststatistic and critical values of test:
##
## test 10pct 5pct 1pct
## r <= 3 | 3.85 10.49 12.25 16.26
## r <= 2 | 18.72 22.76 25.32 30.45
## r <= 1 | 36.42 39.06 42.44 48.45
## r = 0 | 84.92 59.14 62.99 70.05
##
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.l1 1.00000 1.0000 1.0000 1.0000 1.0000
## e.l1 -0.02385 1.2947 -2.8832 4.2418 -8.2904
## U.l1 3.16875 3.4037 -7.4262 6.8414 -12.5578
## rw.l1 1.83528 -0.3331 1.3979 -0.1394 2.4466
## trend.l1 -1.30156 -0.2303 -0.5093 -1.5926 0.2831
##
## Weights W:
## (This is the loading matrix)
##
## prod.l1 e.l1 U.l1 rw.l1 trend.l1
## prod.d -0.006535 -0.02763 -0.070975 -0.014754 1.077e-11
## e.d -0.008503 0.11414 -0.008157 0.003988 7.400e-12
## U.d -0.004719 -0.06154 0.020719 -0.006557 -4.664e-12
## rw.d -0.046213 -0.14580 -0.016945 0.011896 6.952e-12
vecm.r1 <- cajorls(vecm, r = 1)
vecm.r1
## $rlm
##
## Call:
## lm(formula = substitute(form1), data = data.mat)
##
## Coefficients:
## prod.d e.d U.d rw.d
## ect1 -0.00654 -0.00850 -0.00472 -0.04621
## constant 8.27481 10.33131 5.68783 55.46912
## prod.dl1 0.23444 0.20095 -0.13892 -0.07449
## e.dl1 -0.24654 0.82156 -0.64685 -0.63408
## U.dl1 -0.97987 0.00338 -0.19113 0.06314
## rw.dl1 0.00471 -0.07849 0.01726 -0.01208
## prod.dl2 -0.02952 0.04827 -0.00291 -0.25194
## e.dl2 -0.58047 -0.45969 -0.01974 0.08120
## U.dl2 -0.12810 -0.10341 -0.26269 -0.23001
## rw.dl2 -0.19026 -0.09583 0.08035 -0.15739
##
##
## $beta
## ect1
## prod.l1 1.00000
## e.l1 -0.02385
## U.l1 3.16875
## rw.l1 1.83528
## trend.l1 -1.30156
Для того, чтобы теперь перейти к структурной форме VAR, необходимо решить проблему идентификации. В только что построенной VECM необходимо \( \frac{1}{2} k(k-1) = 6 \) линейно независимых ограничений. Т.к. по результатам теста о коинтеграции можно предположить \( r=1 \), существует \( k^{*} = r(k-r)=3 \) LR соотношений и \( r=1 \) SR соотношений. Так как матрица неполного ранга, в ней присутствует \( r k^{*} = 3 \) линейно независимых соотношений. Таким образом, необходимо еще \( \frac{1}{2} k^{*}(k^{*}-1) = 3 \) ограничений.
Оценим SVEC с учетом данных ограничений:
vecm <- ca.jo(Canada[, c("prod", "e", "U", "rw")], type = "trace", ecdet = "trend",
K = 3, spec = "transitory")
SR <- matrix(NA, nrow = 4, ncol = 4) # Насильно создаем SR матрицу ограничений
SR[4, 2] <- 0 # Предположение 3
LR <- matrix(NA, nrow = 4, ncol = 4) # Насильно создаем LR матрицу ограничений
LR[1, 2:4] <- 0 # Предположение 2
LR[2:4, 4] <- 0 # Предположение 1
svec <- SVEC(vecm, LR = LR, SR = SR, r = 1, lrtest = FALSE, boot = TRUE, runs = 100)
summary(svec)
##
## SVEC Estimation Results:
## ========================
##
## Call:
## SVEC(x = vecm, LR = LR, SR = SR, r = 1, lrtest = FALSE, boot = TRUE,
## runs = 100)
##
## Type: B-model
## Sample size: 81
## Log Likelihood: -161.838
## Number of iterations: 10
##
## Estimated contemporaneous impact matrix:
## prod e U rw
## prod 0.5840 0.0743 -0.15258 0.0690
## e -0.1203 0.2614 -0.15510 0.0898
## U 0.0253 -0.2672 0.00549 0.0498
## rw 0.1117 0.0000 0.48377 0.4879
##
## Estimated standard errors for impact matrix:
## prod e U rw
## prod 0.0873 0.1259 0.236 0.0693
## e 0.0723 0.0602 0.190 0.0463
## U 0.0637 0.0502 0.056 0.0350
## rw 0.1465 0.0000 0.694 0.0807
##
## Estimated long run impact matrix:
## prod e U rw
## prod 0.791 0.000 0.000 0
## e 0.202 0.577 -0.492 0
## U -0.159 -0.341 0.141 0
## rw -0.153 0.596 -0.250 0
##
## Estimated standard errors for long-run matrix:
## prod e U rw
## prod 0.165 0.000 0.000 0
## e 0.264 0.161 0.619 0
## U 0.137 0.084 0.167 0
## rw 0.214 0.147 0.297 0
##
## Covariance matrix of reduced form residuals (*100):
## prod e U rw
## prod 37.464 -2.10 -0.251 2.51
## e -2.096 11.49 -6.927 -4.47
## U -0.251 -6.93 7.454 2.98
## rw 2.509 -4.47 2.978 48.46
В данной работе проблема идентификации решилась накладыванием априорных ограничений. Однако, что делать, если ограничений оказывается слишком много, и система становится переопределенной? В следующем разделе проводится тестирование ограничений с помощью LR-теста.
LR[3, 3] <- 0
svec.oi <- update(svec, LR = LR, lrtest = TRUE, boot = FALSE)
str(svec.oi$LRover)
## List of 5
## $ statistic: Named num 6.07
## ..- attr(*, "names")= chr "Chi^2"
## $ parameter: Named num 1
## ..- attr(*, "names")= chr "df"
## $ p.value : Named num 0.0137
## ..- attr(*, "names")= chr "Chi^2"
## $ method : chr "LR overidentification"
## $ data.name: chr "vecm"
## - attr(*, "class")= chr "htest"
Таким образом, видно, что на 5% уровне значимости нулевая гипотеза отвергается.
Строим impulse-response functions:
svec.irf <- irf(svec, n.ahead = 48, boot = TRUE, mar = rep(2, 4), oma = rep(0.5,
4))
plot(svec.irf, nc = 2)
vec2var(vecm)
##
## Coefficient matrix of lagged endogenous variables:
##
## A1:
## prod.l1 e.l1 U.l1 rw.l1
## prod 1.2279 -0.2464 -1.00058 -0.007287
## e 0.1924 1.8218 -0.02357 -0.094097
## U -0.1436 -0.6467 0.79392 0.008603
## rw -0.1207 -0.6330 -0.08330 0.903103
##
##
## A2:
## prod.l2 e.l2 U.l2 rw.l2
## prod -0.2640 -0.3339 0.85177 -0.19497
## e -0.1527 -1.2813 -0.10679 -0.01734
## U 0.1360 0.6271 -0.07156 0.06309
## rw -0.1774 0.7153 -0.29315 -0.14531
##
##
## A3:
## prod.l3 e.l3 U.l3 rw.l3
## prod 0.029520 0.58047 0.1281 0.19026
## e -0.048273 0.45969 0.1034 0.09583
## U 0.002909 0.01974 0.2627 -0.08035
## rw 0.251940 -0.08120 0.2300 0.15739
##
##
## Coefficient matrix of deterministic regressor(s).
##
## constant trend.l1
## prod 8.275 0.008506
## e 10.331 0.011068
## U 5.688 0.006142
## rw 55.469 0.060149
Для стандартного \( VAR(p) \) прогнозные значения создаются следующими коммандами:
pred <- predict(p1ct, n.ahead = 48)
fanchart(pred, nc = 2)
Для \( VECM \), \( SVEC \) и \( SVAR \) прогноз строится следующим способом:
fevd <- fevd(svec, n.ahead = 48)