Vector Autoregression (VAR) Model
Ibikunle Gabriel
2026-03-21
What is a VAR Model?
A Vector Autoregression (VAR) model is a time series model in which multiple variables depend on their own past values each other’s past values. So instead of modeling one variable, you model a system of variables together.
Loading required packages
Set and confirm the working directory
getwd()## [1] "C:/Users/Hp/Documents/ARDL"load the data into the R environment
data <- read.csv("C:\\Users\\Hp\\Documents\\ARDL\\ECON DEVELOPMENT DATA.csv")Select the Variables needed and convert t time series object
data_var <- data[, c("ECNDEV","INF", "EDU", "K", "HEA", "UNEMP")]data_var <- ts(data_var)Determine Optimal Lag Length
lag_selection <- VARselect(data_var, lag.max = 10, type = "const")
lag_selection## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 4 4 4 4
##
## $criteria
## 1 2 3 4 5 6 7 8 9 10
## AIC(n) -1.9636674 -5.40392011 -2.064357e+02 -Inf -Inf -Inf -Inf -Inf -Inf -Inf
## HQ(n) -1.5554748 -4.64584824 -2.053277e+02 -Inf -Inf -Inf -Inf -Inf -Inf -Inf
## SC(n) 0.1273704 -1.52056424 -2.007600e+02 -Inf -Inf -Inf -Inf -Inf -Inf -Inf
## FPE(n) 0.1689179 0.02023512 8.740723e-86 0 0 0 0 0 0 0lag_selection$selection## AIC(n) HQ(n) SC(n) FPE(n)
## 4 4 4 4Estimate VAR Model
var_model <- VAR(data_var, p = 2, type = "const")
summary(var_model)##
## VAR Estimation Results:
## =========================
## Endogenous variables: ECNDEV, INF, EDU, K, HEA, UNEMP
## Deterministic variables: const
## Sample size: 28
## Log Likelihood: -174.992
## Roots of the characteristic polynomial:
## 1.311 1.164 1.164 0.9641 0.9641 0.5181 0.5181 0.5041 0.5041 0.3905 0.1493 0.1493
## Call:
## VAR(y = data_var, p = 2, type = "const")
##
##
## Estimation results for equation ECNDEV:
## =======================================
## ECNDEV = ECNDEV.l1 + INF.l1 + EDU.l1 + K.l1 + HEA.l1 + UNEMP.l1 + ECNDEV.l2 + INF.l2 + EDU.l2 + K.l2 + HEA.l2 + UNEMP.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## ECNDEV.l1 4.496e-01 2.992e-01 1.503 0.154
## INF.l1 1.224e-04 2.610e-04 0.469 0.646
## EDU.l1 -3.384e-02 2.388e-02 -1.417 0.177
## K.l1 1.674e-07 1.350e-06 0.124 0.903
## HEA.l1 1.394e-03 7.440e-03 0.187 0.854
## UNEMP.l1 -2.856e-03 1.238e-02 -0.231 0.821
## ECNDEV.l2 1.547e-01 2.487e-01 0.622 0.543
## INF.l2 -2.342e-04 2.691e-04 -0.870 0.398
## EDU.l2 3.147e-02 2.354e-02 1.337 0.201
## K.l2 3.866e-07 1.236e-06 0.313 0.759
## HEA.l2 7.201e-03 7.448e-03 0.967 0.349
## UNEMP.l2 5.235e-03 2.916e-02 0.180 0.860
## const 2.374e-01 1.584e-01 1.499 0.155
##
##
## Residual standard error: 0.01049 on 15 degrees of freedom
## Multiple R-Squared: 0.9236, Adjusted R-squared: 0.8624
## F-statistic: 15.1 on 12 and 15 DF, p-value: 3.022e-06
##
##
## Estimation results for equation INF:
## ====================================
## INF = ECNDEV.l1 + INF.l1 + EDU.l1 + K.l1 + HEA.l1 + UNEMP.l1 + ECNDEV.l2 + INF.l2 + EDU.l2 + K.l2 + HEA.l2 + UNEMP.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## ECNDEV.l1 3.040e+01 1.059e+02 0.287 0.778046
## INF.l1 3.892e-01 9.244e-02 4.210 0.000758 ***
## EDU.l1 1.039e+01 8.457e+00 1.229 0.238073
## K.l1 -1.146e-04 4.781e-04 -0.240 0.813876
## HEA.l1 1.297e+00 2.634e+00 0.492 0.629607
## UNEMP.l1 2.074e+00 4.383e+00 0.473 0.642916
## ECNDEV.l2 1.784e+01 8.806e+01 0.203 0.842200
## INF.l2 -7.891e-02 9.529e-02 -0.828 0.420625
## EDU.l2 -1.047e+01 8.337e+00 -1.256 0.228333
## K.l2 3.471e-04 4.377e-04 0.793 0.440072
## HEA.l2 -2.041e+00 2.638e+00 -0.774 0.451096
## UNEMP.l2 -2.987e+00 1.033e+01 -0.289 0.776318
## const -6.642e+00 5.608e+01 -0.118 0.907293
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 3.713 on 15 degrees of freedom
## Multiple R-Squared: 0.7497, Adjusted R-squared: 0.5494
## F-statistic: 3.743 on 12 and 15 DF, p-value: 0.009093
##
##
## Estimation results for equation EDU:
## ====================================
## EDU = ECNDEV.l1 + INF.l1 + EDU.l1 + K.l1 + HEA.l1 + UNEMP.l1 + ECNDEV.l2 + INF.l2 + EDU.l2 + K.l2 + HEA.l2 + UNEMP.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## ECNDEV.l1 -2.851e+00 3.982e+00 -0.716 0.485110
## INF.l1 -4.170e-03 3.475e-03 -1.200 0.248780
## EDU.l1 1.397e+00 3.179e-01 4.395 0.000522 ***
## K.l1 3.960e-08 1.797e-05 0.002 0.998271
## HEA.l1 -3.901e-03 9.904e-02 -0.039 0.969101
## UNEMP.l1 1.425e-01 1.648e-01 0.865 0.400629
## ECNDEV.l2 5.939e-01 3.310e+00 0.179 0.860010
## INF.l2 -2.324e-03 3.582e-03 -0.649 0.526321
## EDU.l2 -4.331e-01 3.134e-01 -1.382 0.187298
## K.l2 7.732e-06 1.645e-05 0.470 0.645146
## HEA.l2 -1.859e-02 9.915e-02 -0.188 0.853760
## UNEMP.l2 -3.792e-02 3.882e-01 -0.098 0.923470
## const 2.886e+00 2.108e+00 1.369 0.191203
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.1396 on 15 degrees of freedom
## Multiple R-Squared: 0.9984, Adjusted R-squared: 0.9971
## F-statistic: 768.1 on 12 and 15 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation K:
## ==================================
## K = ECNDEV.l1 + INF.l1 + EDU.l1 + K.l1 + HEA.l1 + UNEMP.l1 + ECNDEV.l2 + INF.l2 + EDU.l2 + K.l2 + HEA.l2 + UNEMP.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## ECNDEV.l1 -8.605e+04 6.169e+04 -1.395 0.183401
## INF.l1 -2.939e+01 5.383e+01 -0.546 0.593146
## EDU.l1 -7.624e+02 4.925e+03 -0.155 0.879041
## K.l1 -4.120e-02 2.784e-01 -0.148 0.884328
## HEA.l1 1.292e+03 1.534e+03 0.842 0.412819
## UNEMP.l1 6.569e+02 2.552e+03 0.257 0.800376
## ECNDEV.l2 -2.075e+03 5.128e+04 -0.040 0.968251
## INF.l2 5.168e+01 5.549e+01 0.931 0.366444
## EDU.l2 1.016e+03 4.855e+03 0.209 0.837121
## K.l2 1.143e+00 2.549e-01 4.484 0.000437 ***
## HEA.l2 2.005e+02 1.536e+03 0.131 0.897889
## UNEMP.l2 1.220e+04 6.013e+03 2.028 0.060700 .
## const -3.094e+04 3.266e+04 -0.947 0.358448
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 2162 on 15 degrees of freedom
## Multiple R-Squared: 0.9941, Adjusted R-squared: 0.9894
## F-statistic: 210.9 on 12 and 15 DF, p-value: 1.871e-14
##
##
## Estimation results for equation HEA:
## ====================================
## HEA = ECNDEV.l1 + INF.l1 + EDU.l1 + K.l1 + HEA.l1 + UNEMP.l1 + ECNDEV.l2 + INF.l2 + EDU.l2 + K.l2 + HEA.l2 + UNEMP.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## ECNDEV.l1 -9.409e-02 9.768e+00 -0.010 0.9924
## INF.l1 2.358e-03 8.524e-03 0.277 0.7859
## EDU.l1 5.852e-01 7.798e-01 0.750 0.4646
## K.l1 2.350e-05 4.408e-05 0.533 0.6017
## HEA.l1 5.129e-01 2.429e-01 2.111 0.0519 .
## UNEMP.l1 3.892e-01 4.041e-01 0.963 0.3508
## ECNDEV.l2 7.197e+00 8.120e+00 0.886 0.3894
## INF.l2 -3.013e-03 8.787e-03 -0.343 0.7364
## EDU.l2 -3.216e-01 7.688e-01 -0.418 0.6816
## K.l2 -7.152e-06 4.036e-05 -0.177 0.8617
## HEA.l2 -3.826e-01 2.432e-01 -1.573 0.1365
## UNEMP.l2 -1.073e+00 9.521e-01 -1.127 0.2776
## const -6.383e+00 5.171e+00 -1.234 0.2361
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.3424 on 15 degrees of freedom
## Multiple R-Squared: 0.9317, Adjusted R-squared: 0.877
## F-statistic: 17.04 on 12 and 15 DF, p-value: 1.355e-06
##
##
## Estimation results for equation UNEMP:
## ======================================
## UNEMP = ECNDEV.l1 + INF.l1 + EDU.l1 + K.l1 + HEA.l1 + UNEMP.l1 + ECNDEV.l2 + INF.l2 + EDU.l2 + K.l2 + HEA.l2 + UNEMP.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## ECNDEV.l1 5.822e+00 8.930e+00 0.652 0.5243
## INF.l1 2.316e-03 7.792e-03 0.297 0.7703
## EDU.l1 1.069e-01 7.128e-01 0.150 0.8827
## K.l1 -5.977e-05 4.030e-05 -1.483 0.1587
## HEA.l1 1.358e-01 2.221e-01 0.611 0.5501
## UNEMP.l1 9.375e-01 3.694e-01 2.538 0.0227 *
## ECNDEV.l2 7.906e-01 7.423e+00 0.107 0.9166
## INF.l2 3.809e-03 8.032e-03 0.474 0.6422
## EDU.l2 -1.166e-01 7.028e-01 -0.166 0.8704
## K.l2 2.400e-05 3.689e-05 0.651 0.5252
## HEA.l2 1.559e-01 2.223e-01 0.701 0.4939
## UNEMP.l2 2.974e-01 8.704e-01 0.342 0.7374
## const -5.983e+00 4.727e+00 -1.266 0.2250
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.313 on 15 degrees of freedom
## Multiple R-Squared: 0.858, Adjusted R-squared: 0.7445
## F-statistic: 7.555 on 12 and 15 DF, p-value: 0.0002285
##
##
##
## Covariance matrix of residuals:
## ECNDEV INF EDU K HEA UNEMP
## ECNDEV 1.100e-04 -3.304e-03 -0.0006976 -8.493e-01 -1.192e-05 3.549e-04
## INF -3.304e-03 1.379e+01 -0.0035223 -2.653e+03 -4.324e-01 -9.179e-02
## EDU -6.976e-04 -3.522e-03 0.0194867 1.653e+01 1.077e-02 -2.528e-02
## K -8.493e-01 -2.653e+03 16.5345369 4.676e+06 1.487e+02 1.479e+02
## HEA -1.192e-05 -4.324e-01 0.0107678 1.487e+02 1.172e-01 1.077e-02
## UNEMP 3.549e-04 -9.179e-02 -0.0252769 1.479e+02 1.077e-02 9.797e-02
##
## Correlation matrix of residuals:
## ECNDEV INF EDU K HEA UNEMP
## ECNDEV 1.00000 -0.084847 -0.476571 -0.03745 -0.00332 0.10814
## INF -0.08485 1.000000 -0.006795 -0.33036 -0.34008 -0.07897
## EDU -0.47657 -0.006795 1.000000 0.05477 0.22528 -0.57851
## K -0.03745 -0.330358 0.054773 1.00000 0.20087 0.21846
## HEA -0.00332 -0.340076 0.225284 0.20087 1.00000 0.10047
## UNEMP 0.10814 -0.078970 -0.578506 0.21846 0.10047 1.00000Check Stability of VAR
stability <- roots(var_model)
stability## [1] 1.3107490 1.1639995 1.1639995 0.9641346 0.9641346 0.5180960 0.5180960
## [8] 0.5040750 0.5040750 0.3905212 0.1493322 0.1493322plot(stability)
Diagnostic Tests Serial Correlation
serial.test(var_model, lags.pt = 16, type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var_model
## Chi-squared = 444.65, df = 504, p-value = 0.973Heteroskedasticity
arch.test(var_model, lags.multi = 5)##
## ARCH (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 483, df = 2205, p-value = 1Normality
normality.test(var_model)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 273.02, df = 12, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 71.264, df = 6, p-value = 2.25e-13
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 201.76, df = 6, p-value < 2.2e-16Impulse Response Function (IRF)
irf_result <- irf(var_model,
impulse = "ECNDEV",
response = "HEA",
n.ahead = 10,
boot = TRUE)
plot(irf_result)
irf_result1 <- irf(var_model,
impulse = "ECNDEV",
response = "HEA",
n.ahead = 10,
boot = TRUE)
plot(irf_result1)
irf_result2 <- irf(var_model,
impulse = "ECNDEV",
response = "INF",
n.ahead = 10,
boot = TRUE)
plot(irf_result2)
irf_result3 <- irf(var_model,
impulse = "ECNDEV",
response = "EDU",
n.ahead = 10,
boot = TRUE)
plot(irf_result3)
irf_result4 <- irf(var_model,
impulse = "ECNDEV",
response = "K",
n.ahead = 10,
boot = TRUE)
plot(irf_result4)
irf_result5 <- irf(var_model,
impulse = "ECNDEV",
response = "UNEMP",
n.ahead = 10,
boot = TRUE)
plot(irf_result5)
Get all responses automatically
irf_all <- irf(var_model, n.ahead = 10, boot = TRUE)
plot(irf_all)
Forecast Error Variance Decomposition (FEVD)
fevd_result <- fevd(var_model, n.ahead = 10)
fevd_result## $ECNDEV
## ECNDEV INF EDU K HEA UNEMP
## [1,] 1.0000000 0.000000000 0.00000000 0.0000000000 0.000000000 0.000000000
## [2,] 0.9240525 0.001224893 0.07176226 0.0001887892 0.000324996 0.002446521
## [3,] 0.8730834 0.007579594 0.07913266 0.0085836480 0.028190350 0.003430379
## [4,] 0.8409085 0.008411385 0.08536294 0.0125830828 0.049453646 0.003280435
## [5,] 0.8249933 0.009117972 0.08291506 0.0237581799 0.056139050 0.003076391
## [6,] 0.8207409 0.009275980 0.08118352 0.0277322315 0.057761633 0.003305701
## [7,] 0.7915177 0.012032488 0.07972605 0.0503929856 0.059157056 0.007173686
## [8,] 0.7611158 0.014011624 0.08082358 0.0559750982 0.066059202 0.022014733
## [9,] 0.6345697 0.018117531 0.09101642 0.0826559507 0.078614292 0.095026057
## [10,] 0.5285712 0.016653161 0.10151658 0.0702967077 0.105690572 0.177271767
##
## $INF
## ECNDEV INF EDU K HEA UNEMP
## [1,] 0.00719897 0.9928010 0.00000000 0.0000000000 0.00000000 0.00000000
## [2,] 0.01605428 0.8970241 0.05647169 0.0001343425 0.01749320 0.01282239
## [3,] 0.01645516 0.8533097 0.07718949 0.0162770362 0.01742150 0.01934713
## [4,] 0.01616027 0.8400272 0.08294921 0.0162069642 0.02245323 0.02220310
## [5,] 0.01573922 0.7616886 0.09626122 0.0706546777 0.02122847 0.03442778
## [6,] 0.01466856 0.7103527 0.11881093 0.0661328655 0.03529558 0.05473936
## [7,] 0.01867821 0.4804761 0.16162543 0.1096816143 0.05835396 0.17118467
## [8,] 0.02672120 0.4007953 0.16569930 0.1095389836 0.08872974 0.20851550
## [9,] 0.06645039 0.2612462 0.17388032 0.1165843893 0.10402452 0.27781421
## [10,] 0.09649070 0.2212483 0.14978832 0.1676201177 0.11984153 0.24501101
##
## $EDU
## ECNDEV INF EDU K HEA UNEMP
## [1,] 0.2271201 0.002246888 0.7706330 0.000000000 0.000000000 0.00000000
## [2,] 0.3234770 0.014213818 0.6406598 0.001823896 0.001342623 0.01848284
## [3,] 0.3594336 0.033584121 0.5456631 0.004487685 0.003083817 0.05374763
## [4,] 0.3686043 0.048061121 0.4522017 0.006487850 0.008480082 0.11616490
## [5,] 0.3396989 0.055911855 0.3582512 0.007283218 0.025049952 0.21380486
## [6,] 0.2806432 0.055831616 0.2757579 0.006027481 0.055680436 0.32605934
## [7,] 0.2196508 0.051071874 0.2147083 0.004877613 0.093427665 0.41626367
## [8,] 0.1841135 0.045515489 0.1755019 0.007599013 0.128180967 0.45908914
## [9,] 0.1818487 0.041400189 0.1545760 0.016787236 0.151573391 0.45381451
## [10,] 0.1995295 0.038194994 0.1563249 0.035519159 0.155190242 0.41524123
##
## $K
## ECNDEV INF EDU K HEA UNEMP
## [1,] 0.001402711 0.11205272 0.0005753641 0.8859692 0.00000000 0.000000000
## [2,] 0.119925696 0.10311805 0.0028302583 0.7349281 0.03531231 0.003885586
## [3,] 0.022725602 0.05448733 0.1248998844 0.4641275 0.04821965 0.285539986
## [4,] 0.018073360 0.03832436 0.1406183321 0.3431576 0.09542850 0.364397811
## [5,] 0.025734068 0.02443298 0.1729743140 0.2413898 0.10066674 0.434802076
## [6,] 0.038719909 0.02006230 0.1542581799 0.2332823 0.14285735 0.410820011
## [7,] 0.093967338 0.01760554 0.1569295154 0.1913172 0.13871084 0.401469610
## [8,] 0.120796973 0.02042137 0.1218689058 0.2793461 0.14225042 0.315316246
## [9,] 0.191896776 0.02233985 0.1097926892 0.2761499 0.12705703 0.272763782
## [10,] 0.151570406 0.02677570 0.1029251138 0.3897253 0.08586752 0.243135922
##
## $HEA
## ECNDEV INF EDU K HEA UNEMP
## [1,] 1.102115e-05 0.11668348 0.05587648 0.007355668 0.8200734 0.00000000
## [2,] 5.089487e-03 0.11420251 0.04649747 0.051244742 0.7394720 0.04349381
## [3,] 4.267307e-03 0.09110582 0.17658308 0.062447885 0.5938353 0.07176056
## [4,] 3.852598e-03 0.07741748 0.20605072 0.079154821 0.5307375 0.10278686
## [5,] 3.331666e-02 0.06086019 0.22247892 0.073595988 0.4406822 0.16906601
## [6,] 6.100222e-02 0.07069721 0.18209071 0.182141921 0.3636971 0.14037081
## [7,] 1.367408e-01 0.06360354 0.16660518 0.179010692 0.3273188 0.12672103
## [8,] 9.326719e-02 0.05290282 0.14402589 0.265987046 0.2050579 0.23875914
## [9,] 8.967680e-02 0.04601682 0.12790794 0.274583364 0.1933804 0.26843470
## [10,] 5.350613e-02 0.03285330 0.14550089 0.247284221 0.1389434 0.38191200
##
## $UNEMP
## ECNDEV INF EDU K HEA UNEMP
## [1,] 0.01169396 0.004906612 0.3649054 0.05147368 0.04889305 0.5181273
## [2,] 0.04848687 0.003680035 0.3178366 0.04042969 0.08575721 0.5038096
## [3,] 0.11286904 0.002081372 0.2697037 0.02578600 0.12605460 0.4635053
## [4,] 0.18014006 0.004800206 0.2146572 0.06845749 0.15359037 0.3783547
## [5,] 0.29045192 0.003888556 0.1797790 0.05495739 0.15371669 0.3172064
## [6,] 0.28679183 0.009318426 0.1558274 0.14024699 0.11263396 0.2951814
## [7,] 0.28617728 0.008074213 0.1546987 0.11808846 0.09992394 0.3330374
## [8,] 0.14813801 0.008617612 0.1828993 0.12967457 0.09008513 0.4405854
## [9,] 0.11319952 0.008916827 0.1761997 0.13024841 0.11971314 0.4517224
## [10,] 0.12098774 0.007667358 0.1812266 0.10758187 0.12932528 0.4532112# Plot FEVD
plot(fevd_result)
Forecast (Optional)
———— VAR with retail dataset—————–
df1 <- read.csv("C:\\Users\\Hp\\Documents\\ARDL\\retail_sales_dataset.csv")# Select the Variables
data_var <- df1[, c("Age", "Price", "TotalAmount","Quantity")]# Convert to time series object
data_var <- ts(data_var)# Determine Optimal Lag Length
lag_selection <- VARselect(data_var, lag.max = 10, type = "const")
lag_selection## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1
##
## $criteria
## 1 2 3 4 5
## AIC(n) 2.671711e+01 2.673857e+01 2.675771e+01 2.677408e+01 2.679010e+01
## HQ(n) 2.675473e+01 2.680630e+01 2.685553e+01 2.690200e+01 2.694812e+01
## SC(n) 2.681605e+01 2.691667e+01 2.701496e+01 2.711048e+01 2.720566e+01
## FPE(n) 4.009530e+11 4.096529e+11 4.175686e+11 4.244628e+11 4.313240e+11
## 6 7 8 9 10
## AIC(n) 2.680735e+01 2.682702e+01 2.685120e+01 2.686279e+01 2.687021e+01
## HQ(n) 2.699547e+01 2.704524e+01 2.709951e+01 2.714120e+01 2.717873e+01
## SC(n) 2.730207e+01 2.740089e+01 2.750422e+01 2.759497e+01 2.768155e+01
## FPE(n) 4.388384e+11 4.475658e+11 4.585319e+11 4.638969e+11 4.673779e+11lag_selection$selection## AIC(n) HQ(n) SC(n) FPE(n)
## 1 1 1 1# Estimate VAR Model
var_model <- VAR(data_var, p = 2, type = "const")
summary(var_model)##
## VAR Estimation Results:
## =========================
## Endogenous variables: Age, Price, TotalAmount, Quantity
## Deterministic variables: const
## Sample size: 998
## Log Likelihood: -18970.316
## Roots of the characteristic polynomial:
## 0.2103 0.2103 0.2016 0.2016 0.1689 0.159 0.1496 0.1496
## Call:
## VAR(y = data_var, p = 2, type = "const")
##
##
## Estimation results for equation Age:
## ====================================
## Age = Age.l1 + Price.l1 + TotalAmount.l1 + Quantity.l1 + Age.l2 + Price.l2 + TotalAmount.l2 + Quantity.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## Age.l1 0.0233961 0.0318247 0.735 0.462
## Price.l1 -0.0061575 0.0055454 -1.110 0.267
## TotalAmount.l1 0.0021984 0.0020271 1.084 0.278
## Quantity.l1 0.8032704 0.5251267 1.530 0.126
## Age.l2 0.0189613 0.0316846 0.598 0.550
## Price.l2 -0.0020545 0.0055497 -0.370 0.711
## TotalAmount.l2 -0.0001209 0.0020295 -0.060 0.952
## Quantity.l2 0.1217929 0.5260773 0.232 0.817
## const 37.8686558 2.6540425 14.268 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 13.65 on 989 degrees of freedom
## Multiple R-Squared: 0.01291, Adjusted R-squared: 0.004925
## F-statistic: 1.617 on 8 and 989 DF, p-value: 0.1157
##
##
## Estimation results for equation Price:
## ======================================
## Price = Age.l1 + Price.l1 + TotalAmount.l1 + Quantity.l1 + Age.l2 + Price.l2 + TotalAmount.l2 + Quantity.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## Age.l1 0.453050 0.442326 1.024 0.306
## Price.l1 -0.011991 0.077074 -0.156 0.876
## TotalAmount.l1 0.004325 0.028175 0.153 0.878
## Quantity.l1 -1.271694 7.298659 -0.174 0.862
## Age.l2 -0.279887 0.440380 -0.636 0.525
## Price.l2 0.113265 0.077134 1.468 0.142
## TotalAmount.l2 -0.040796 0.028207 -1.446 0.148
## Quantity.l2 -2.047672 7.311871 -0.280 0.779
## const 179.277244 36.888145 4.860 1.36e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 189.7 on 989 degrees of freedom
## Multiple R-Squared: 0.006211, Adjusted R-squared: -0.001828
## F-statistic: 0.7726 on 8 and 989 DF, p-value: 0.6271
##
##
## Estimation results for equation TotalAmount:
## ============================================
## TotalAmount = Age.l1 + Price.l1 + TotalAmount.l1 + Quantity.l1 + Age.l2 + Price.l2 + TotalAmount.l2 + Quantity.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## Age.l1 1.26686 1.30855 0.968 0.333
## Price.l1 -0.15286 0.22801 -0.670 0.503
## TotalAmount.l1 0.05815 0.08335 0.698 0.486
## Quantity.l1 -23.17184 21.59182 -1.073 0.283
## Age.l2 -0.76089 1.30279 -0.584 0.559
## Price.l2 0.25126 0.22819 1.101 0.271
## TotalAmount.l2 -0.07776 0.08345 -0.932 0.352
## Quantity.l2 -0.28847 21.63090 -0.013 0.989
## const 484.97068 109.12718 4.444 9.82e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 561.3 on 989 degrees of freedom
## Multiple R-Squared: 0.00423, Adjusted R-squared: -0.003825
## F-statistic: 0.5251 on 8 and 989 DF, p-value: 0.8382
##
##
## Estimation results for equation Quantity:
## =========================================
## Quantity = Age.l1 + Price.l1 + TotalAmount.l1 + Quantity.l1 + Age.l2 + Price.l2 + TotalAmount.l2 + Quantity.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## Age.l1 2.469e-03 2.647e-03 0.933 0.3511
## Price.l1 -4.996e-04 4.613e-04 -1.083 0.2791
## TotalAmount.l1 2.319e-04 1.686e-04 1.375 0.1694
## Quantity.l1 -7.958e-02 4.368e-02 -1.822 0.0688 .
## Age.l2 1.110e-03 2.635e-03 0.421 0.6737
## Price.l2 6.144e-05 4.616e-04 0.133 0.8941
## TotalAmount.l2 2.695e-05 1.688e-04 0.160 0.8732
## Quantity.l2 -2.337e-03 4.376e-02 -0.053 0.9574
## const 2.532e+00 2.208e-01 11.471 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.135 on 989 degrees of freedom
## Multiple R-Squared: 0.004891, Adjusted R-squared: -0.003158
## F-statistic: 0.6076 on 8 and 989 DF, p-value: 0.772
##
##
##
## Covariance matrix of residuals:
## Age Price TotalAmount Quantity
## Age 186.3417 -95.023 -444.0 -0.3423
## Price -95.0227 35997.142 90732.4 3.8719
## TotalAmount -444.0476 90732.380 315036.1 238.1905
## Quantity -0.3423 3.872 238.2 1.2892
##
## Correlation matrix of residuals:
## Age Price TotalAmount Quantity
## Age 1.00000 -0.03669 -0.05796 -0.02209
## Price -0.03669 1.00000 0.85202 0.01797
## TotalAmount -0.05796 0.85202 1.00000 0.37375
## Quantity -0.02209 0.01797 0.37375 1.00000# Check Stability of VAR
stability <- roots(var_model)
stability## [1] 0.2103055 0.2103055 0.2015807 0.2015807 0.1689167 0.1590061 0.1495568
## [8] 0.1495568roots_var <- roots(var_model)# Plot empty coordinate system
plot(Re(roots_var), Im(roots_var),
xlim = c(-1.5, 1.5),
ylim = c(-1.5, 1.5),
xlab = "Real Part",
ylab = "Imaginary Part",
main = "Stability of VAR Model",
pch = 19,
col = "blue")
# Add unit circle
theta <- seq(0, 2*pi, length = 500)
lines(cos(theta), sin(theta), col = "red", lwd = 2)
# Add horizontal and vertical axis
abline(h = 0, v = 0, col = "gray")
# Add grid
grid(col = "lightgray")
# Save the plot automatically to your working directory
ggsave("Stability_plot.png", width = 12, height = 8, dpi = 300)Diagnostic Tests
## Serial Correlation
serial.test(var_model, lags.pt = 16, type = "PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var_model
## Chi-squared = 230.91, df = 224, p-value = 0.3615## Heteroskedasticity
arch.test(var_model, lags.multi = 5)##
## ARCH (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 515.05, df = 500, p-value = 0.3112# Normality
normality.test(var_model)## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 372.87, df = 8, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 91.201, df = 4, p-value < 2.2e-16
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object var_model
## Chi-squared = 281.67, df = 4, p-value < 2.2e-16Impulse Response Function (IRF)
# Plot of just two variables
irf_result <- irf(var_model,
impulse = "Price",
response = "TotalAmount",
n.ahead = 10,
boot = TRUE)
plot(irf_result)
# Plot of just two variables
irf_result1 <- irf(var_model,
impulse = "Price",
response = "TotalAmount",
n.ahead = 10,
boot = TRUE)
plot(irf_result1)
# Plot of just two variables
irf_result2 <- irf(var_model,
impulse = "Price",
response = "Age",
n.ahead = 10,
boot = TRUE)
plot(irf_result2)
# Plot of just two variables
irf_result3 <- irf(var_model,
impulse = "Price",
response = "Quantity",
n.ahead = 10,
boot = TRUE)
plot(irf_result3)
# all responses automatically
irf_all <- irf(var_model, n.ahead = 10, boot = TRUE)
plot(irf_all)
# all responses automatically
irf_all <- irf(var_model, n.ahead = 10, boot = TRUE)
plot(irf_all)
# A more publishable plot
irf_all <- irf(var_model,
n.ahead = 10,
boot = TRUE,
ci = 0.95)
irf_df <- do.call(rbind, lapply(names(irf_all$irf), function(imp) {
responses <- irf_all$irf[[imp]]
lower <- irf_all$Lower[[imp]]
upper <- irf_all$Upper[[imp]]
df <- as.data.frame(responses)
df$Horizon <- 0:(nrow(df)-1)
df$Impulse <- imp
df_long <- pivot_longer(df,
cols = -c(Horizon, Impulse),
names_to = "Response",
values_to = "IRF")
df_long$Lower <- as.vector(lower)
df_long$Upper <- as.vector(upper)
df_long
}))
ggplot(irf_df, aes(x = Horizon, y = IRF)) +
geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "skyblue", alpha = 0.3) +
geom_line(color = "darkblue", size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
facet_grid(Response ~ Impulse) +
labs(title = "Impulse Response Functions",
x = "Forecast Horizon",
y = "Response") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
strip.text = element_text(face = "bold"),
panel.grid.minor = element_blank()
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Save the plot automatically to your working directory
ggsave("IRF_plot.png", width = 12, height = 8, dpi = 300)par(mar = c(2, 2, 2, 2)) # default margins
# Forecast Error Variance Decomposition (FEVD)
fevd_result <- fevd(var_model, n.ahead = 10)
fevd_result## $Age
## Age Price TotalAmount Quantity
## [1,] 1.0000000 0.000000e+00 0.000000000 0.000000000
## [2,] 0.9890525 5.565483e-05 0.008549489 0.002342345
## [3,] 0.9880066 1.071511e-03 0.008579137 0.002342780
## [4,] 0.9879797 1.073369e-03 0.008601706 0.002345189
## [5,] 0.9879720 1.074118e-03 0.008608651 0.002345220
## [6,] 0.9879720 1.074125e-03 0.008608672 0.002345223
## [7,] 0.9879720 1.074124e-03 0.008608699 0.002345225
## [8,] 0.9879720 1.074124e-03 0.008608699 0.002345225
## [9,] 0.9879720 1.074124e-03 0.008608699 0.002345225
## [10,] 0.9879720 1.074124e-03 0.008608699 0.002345225
##
## $Price
## Age Price TotalAmount Quantity
## [1,] 0.001346097 0.9986539 0.000000e+00 0.000000e+00
## [2,] 0.002397323 0.9975698 2.187126e-06 3.070533e-05
## [3,] 0.002658790 0.9925876 4.670521e-03 8.307429e-05
## [4,] 0.002658767 0.9925751 4.671623e-03 9.454105e-05
## [5,] 0.002658831 0.9925638 4.681592e-03 9.574703e-05
## [6,] 0.002658831 0.9925638 4.681604e-03 9.576336e-05
## [7,] 0.002658831 0.9925638 4.681606e-03 9.576387e-05
## [8,] 0.002658831 0.9925638 4.681606e-03 9.576387e-05
## [9,] 0.002658831 0.9925638 4.681606e-03 9.576387e-05
## [10,] 0.002658831 0.9925638 4.681606e-03 9.576387e-05
##
## $TotalAmount
## Age Price TotalAmount Quantity
## [1,] 0.003358835 0.7232874 0.2733537 0.000000000
## [2,] 0.004272292 0.7217804 0.2727836 0.001163693
## [3,] 0.004575520 0.7205677 0.2736912 0.001165589
## [4,] 0.004575498 0.7205656 0.2736900 0.001168937
## [5,] 0.004575639 0.7205590 0.2736956 0.001169716
## [6,] 0.004575639 0.7205590 0.2736956 0.001169734
## [7,] 0.004575639 0.7205590 0.2736956 0.001169735
## [8,] 0.004575639 0.7205590 0.2736956 0.001169735
## [9,] 0.004575639 0.7205590 0.2736956 0.001169735
## [10,] 0.004575639 0.7205590 0.2736956 0.001169735
##
## $Quantity
## Age Price TotalAmount Quantity
## [1,] 0.0004877706 0.0002949575 0.4685011 0.5307162
## [2,] 0.0012586736 0.0004543710 0.4665118 0.5317752
## [3,] 0.0013870615 0.0008552249 0.4662929 0.5314648
## [4,] 0.0013870529 0.0008575298 0.4662936 0.5314618
## [5,] 0.0013870629 0.0008575297 0.4662943 0.5314611
## [6,] 0.0013870636 0.0008575311 0.4662943 0.5314611
## [7,] 0.0013870636 0.0008575313 0.4662943 0.5314611
## [8,] 0.0013870636 0.0008575313 0.4662943 0.5314611
## [9,] 0.0013870636 0.0008575313 0.4662943 0.5314611
## [10,] 0.0013870636 0.0008575313 0.4662943 0.5314611# Plot FEVD
plot(fevd_result)
# A more publishable plot
fevd_df <- do.call(rbind, lapply(names(fevd_result), function(var) {
df <- as.data.frame(fevd_result[[var]])
df$Horizon <- 1:nrow(df)
df$Variable <- var
df
}))
fevd_long <- melt(fevd_df,
id.vars = c("Horizon", "Variable"),
variable.name = "Shock",
value.name = "Contribution")
ggplot(fevd_long, aes(x = Horizon, y = Contribution, fill = Shock)) +
geom_area(position = "stack", alpha = 0.9) +
facet_wrap(~Variable, scales = "free_y") +
labs(title = "Forecast Error Variance Decomposition",
x = "Forecast Horizon",
y = "Variance Contribution",
fill = "Shock Variable") +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
legend.position = "bottom",
strip.text = element_text(face = "bold")
)+
scale_fill_brewer(palette = "Set2")
# Save the plot automatically to your working directory
ggsave("FEVD_plot.png", width = 10, height = 6, dpi = 300)# Forecast (Optional)
forecast_var <- predict(var_model, n.ahead = 2)
par(mar = c(2, 2, 2, 2))
plot(forecast_var)
# Forecast variance plot with ggplot
forecast_var <- predict(var_model, n.ahead = 10, ci = 0.95)
library(dplyr)
library(tidyr)
forecast_df <- do.call(rbind, lapply(names(forecast_var$fcst), function(var){
df <- as.data.frame(forecast_var$fcst[[var]])
df$Horizon <- 1:nrow(df)
df$Variable <- var
df
}))
library(ggplot2)
ggplot(forecast_df, aes(x = Horizon, y = fcst)) +
geom_ribbon(
aes(ymin = lower, ymax = upper),
fill = "grey70",
alpha = 0.4
) +
geom_line(
color = "green",
linewidth = 1
) +
facet_wrap(~Variable, scales = "free_y") +
labs(
title = "VAR Forecasts with 95% Confidence Intervals",
x = "Forecast Horizon",
y = "Forecast Value"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
strip.text = element_text(face = "bold")
)
forecast_var <- predict(var_model, n.ahead = 10, ci = 0.95)
forecast_df <- do.call(rbind, lapply(names(forecast_var$fcst), function(var){
df <- as.data.frame(forecast_var$fcst[[var]])
df$Horizon <- 1:nrow(df)
df$Variable <- var
df
}))
hist_df <- data.frame(data_var)
hist_df$Time <- 1:nrow(hist_df)
hist_long <- tidyr::pivot_longer(
hist_df,
cols = -Time,
names_to = "Variable",
values_to = "Value"
)
hist_df <- data.frame(data_var)
hist_df$Time <- 1:nrow(hist_df)
hist_long <- tidyr::pivot_longer(
hist_df,
cols = -Time,
names_to = "Variable",
values_to = "Value"
)
forecast_df$Time <- max(hist_long$Time) + forecast_df$Horizon
ggplot() +
# historical data
geom_line(
data = hist_long,
aes(x = Time, y = Value),
color = "black",
linewidth = 0.8
) +
# confidence interval
geom_ribbon(
data = forecast_df,
aes(x = Time, ymin = lower, ymax = upper),
fill = "grey70",
alpha = 0.4
) +
# forecast line
geom_line(
data = forecast_df,
aes(x = Time, y = fcst),
color = "blue",
linewidth = 1,
linetype = "dashed"
) +
facet_wrap(~Variable, scales = "free_y") +
labs(
title = "VAR Forecast with Historical Data",
x = "Time",
y = "Value"
) +
theme_classic(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
strip.text = element_text(face = "bold")
)
References
Pesaran, M. H., Shin, Y., & Smith, R. J. (2001). Bounds testing approaches to the analysis of level relationships. Journal of Applied Econometrics, 16(3), 289–326.https://doi.org/10.1002/jae.616
Pesaran, M. H., & Smith, R. (1995). Estimating long-run relationships from dynamic heterogeneous panels. Journal of Econometrics, 68(1), 79–113.https://doi.org/10.1016/0304-4076(94)01644-F