VAR_V1
03/08/2020
Framing the matrix
library(MASS)
sd.inf <-diff(inf, differences = 2)
sd.gdp <-diff(GDP, differences = 2)
sd.M1 <- diff(M1, differences = 2)
sd.rep <- diff(rep, differences = 2)
#bc.rep <- boxcox(rep ~ inf+GDP+M1)
#bc.rep
#lambda <- bc.rep$x # lambda values
#lik <- bc.rep$y # log likelihood values for SSE
#bc <- cbind(lambda, lik) # combine lambda and lik
#sorted_bc <- bc[order(-lik),] # values are sorted to identify the lambda value for the maximum log likelihood for obtaining minimum SSE
#head(sorted_bc, n = 10)
VAR_org <- window(ts.union(inf,GDP,M1,rep), start = c(2004,1), end = c(2019,4))
VAR_stat <- window(ts.union(sd.inf,sd.gdp,sd.M1,sd.rep), start = c(2004,1), end = c(2019,4))## Warning in window.default(x, ...): 'start' value not changed#plot.ts(VAR_org)
#plot.ts(VAR_stat)
plot(VAR_stat, nc = 2, xlab = "")
plot(VAR_org, nc = 2, xlab = "")
Making the series stationary, check ACF, PACF
library(urca)
library(tseries)## Warning: package 'tseries' was built under R version 3.6.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo#pp.test(GDP)
#pp.test(l.gdp)
summary(ur.ers(window(VAR_stat, start = c(2004, 1), end = c(2019, 4)),model = "trend", lag.max = 2))## Warning in window.default(x, ...): 'start' value not changed##
## ###############################################
## # Elliot, Rothenberg and Stock Unit Root Test #
## ###############################################
##
## Test of type DF-GLS
## detrending of series with intercept and trend
##
##
## Call:
## lm(formula = dfgls.form, data = data.dfgls)
##
## Residuals:
## Min 1Q Median 3Q Max
## -943607 -4509 -2060 3038 490430
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## yd.lag -3.62346 0.13379 -27.08 <2e-16 ***
## yd.diff.lag1 1.55640 0.09996 15.57 <2e-16 ***
## yd.diff.lag2 0.65761 0.04843 13.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 92820 on 242 degrees of freedom
## Multiple R-squared: 0.9058, Adjusted R-squared: 0.9046
## F-statistic: 775.7 on 3 and 242 DF, p-value: < 2.2e-16
##
##
## Value of test-statistic is: -27.0835
##
## Critical values of DF-GLS are:
## 1pct 5pct 10pct
## critical values -3.48 -2.89 -2.57#DF & ADF tests for stationarity in Y
#k is the number of lags, Augmented Dickey Fuller test for stationarity
adf.test(sd.inf,"stationary")## Warning in adf.test(sd.inf, "stationary"): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: sd.inf
## Dickey-Fuller = -8.4799, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryadf.test(sd.gdp,"stationary")## Warning in adf.test(sd.gdp, "stationary"): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: sd.gdp
## Dickey-Fuller = -7.5265, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryadf.test(sd.M1,"stationary")## Warning in adf.test(sd.M1, "stationary"): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: sd.M1
## Dickey-Fuller = -7.6893, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryadf.test(sd.rep,"stationary")## Warning in adf.test(sd.rep, "stationary"): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: sd.rep
## Dickey-Fuller = -5.1491, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationaryLag Length Selection for unrestricted VAR
VARselect(VAR_stat, lag.max = 8, type = "both")[["selection"]]## AIC(n) HQ(n) SC(n) FPE(n)
## 8 8 3 7#VARselect(VAR_org, lag.max = 8, type = "both")Lag lenth output = 3. Checking majorly AIC and SC (other name for BIC) We construct multivariate order 3 VAR model, VAR(3).
VAR <- VAR(VAR_stat, p = 3, type = "both")
#VAR <- VAR(VAR_org, p = 1, type = "both")
summary(VAR)##
## VAR Estimation Results:
## =========================
## Endogenous variables: sd.inf, sd.gdp, sd.M1, sd.rep
## Deterministic variables: both
## Sample size: 59
## Log Likelihood: -1594.224
## Roots of the characteristic polynomial:
## 0.9918 0.9918 0.9512 0.7693 0.7693 0.7615 0.7615 0.7078 0.7078 0.6905 0.6905 0.2366
## Call:
## VAR(y = VAR_stat, p = 3, type = "both")
##
##
## Estimation results for equation sd.inf:
## =======================================
## sd.inf = sd.inf.l1 + sd.gdp.l1 + sd.M1.l1 + sd.rep.l1 + sd.inf.l2 + sd.gdp.l2 + sd.M1.l2 + sd.rep.l2 + sd.inf.l3 + sd.gdp.l3 + sd.M1.l3 + sd.rep.l3 + const + trend
##
## Estimate Std. Error t value Pr(>|t|)
## sd.inf.l1 -8.177e-01 1.472e-01 -5.554 1.43e-06 ***
## sd.gdp.l1 9.730e-07 1.597e-06 0.609 0.54542
## sd.M1.l1 6.344e-07 9.499e-07 0.668 0.50763
## sd.rep.l1 -2.449e-02 3.420e-01 -0.072 0.94323
## sd.inf.l2 -5.931e-01 1.880e-01 -3.155 0.00286 **
## sd.gdp.l2 4.540e-06 1.828e-06 2.484 0.01679 *
## sd.M1.l2 -4.276e-08 1.333e-06 -0.032 0.97454
## sd.rep.l2 2.308e-01 3.333e-01 0.693 0.49213
## sd.inf.l3 -2.198e-01 1.649e-01 -1.333 0.18929
## sd.gdp.l3 3.563e-06 1.575e-06 2.262 0.02855 *
## sd.M1.l3 -1.719e-07 9.556e-07 -0.180 0.85802
## sd.rep.l3 -6.215e-01 3.577e-01 -1.738 0.08913 .
## const 1.232e-01 2.910e-01 0.423 0.67400
## trend -1.997e-03 7.842e-03 -0.255 0.80018
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.022 on 45 degrees of freedom
## Multiple R-Squared: 0.6903, Adjusted R-squared: 0.6008
## F-statistic: 7.715 on 13 and 45 DF, p-value: 1.01e-07
##
##
## Estimation results for equation sd.gdp:
## =======================================
## sd.gdp = sd.inf.l1 + sd.gdp.l1 + sd.M1.l1 + sd.rep.l1 + sd.inf.l2 + sd.gdp.l2 + sd.M1.l2 + sd.rep.l2 + sd.inf.l3 + sd.gdp.l3 + sd.M1.l3 + sd.rep.l3 + const + trend
##
## Estimate Std. Error t value Pr(>|t|)
## sd.inf.l1 -2.274e+03 7.533e+03 -0.302 0.7641
## sd.gdp.l1 -8.860e-01 8.172e-02 -10.842 3.89e-14 ***
## sd.M1.l1 -3.025e-02 4.860e-02 -0.622 0.5369
## sd.rep.l1 3.699e+04 1.750e+04 2.114 0.0401 *
## sd.inf.l2 9.275e+03 9.618e+03 0.964 0.3400
## sd.gdp.l2 -8.960e-01 9.352e-02 -9.581 1.95e-12 ***
## sd.M1.l2 1.059e-02 6.818e-02 0.155 0.8772
## sd.rep.l2 1.375e+04 1.705e+04 0.806 0.4244
## sd.inf.l3 -4.651e+03 8.436e+03 -0.551 0.5842
## sd.gdp.l3 -9.351e-01 8.059e-02 -11.604 4.03e-15 ***
## sd.M1.l3 3.316e-02 4.890e-02 0.678 0.5011
## sd.rep.l3 1.384e+04 1.830e+04 0.756 0.4533
## const 7.076e+03 1.489e+04 0.475 0.6370
## trend -8.556e+01 4.012e+02 -0.213 0.8321
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 52290 on 45 degrees of freedom
## Multiple R-Squared: 0.9245, Adjusted R-squared: 0.9027
## F-statistic: 42.41 on 13 and 45 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation sd.M1:
## ======================================
## sd.M1 = sd.inf.l1 + sd.gdp.l1 + sd.M1.l1 + sd.rep.l1 + sd.inf.l2 + sd.gdp.l2 + sd.M1.l2 + sd.rep.l2 + sd.inf.l3 + sd.gdp.l3 + sd.M1.l3 + sd.rep.l3 + const + trend
##
## Estimate Std. Error t value Pr(>|t|)
## sd.inf.l1 4.615e+04 2.132e+04 2.165 0.03572 *
## sd.gdp.l1 3.398e-01 2.312e-01 1.469 0.14867
## sd.M1.l1 -1.252e+00 1.375e-01 -9.102 9.11e-12 ***
## sd.rep.l1 -2.038e+04 4.952e+04 -0.412 0.68263
## sd.inf.l2 2.284e+04 2.722e+04 0.839 0.40572
## sd.gdp.l2 3.116e-01 2.646e-01 1.177 0.24521
## sd.M1.l2 -8.624e-01 1.929e-01 -4.470 5.25e-05 ***
## sd.rep.l2 -7.346e+02 4.826e+04 -0.015 0.98792
## sd.inf.l3 -2.425e+04 2.387e+04 -1.016 0.31517
## sd.gdp.l3 -4.276e-01 2.280e-01 -1.875 0.06728 .
## sd.M1.l3 -3.724e-01 1.384e-01 -2.691 0.00995 **
## sd.rep.l3 -5.978e+03 5.179e+04 -0.115 0.90862
## const 3.847e+03 4.214e+04 0.091 0.92767
## trend -7.483e+01 1.135e+03 -0.066 0.94774
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 148000 on 45 degrees of freedom
## Multiple R-Squared: 0.8096, Adjusted R-squared: 0.7546
## F-statistic: 14.72 on 13 and 45 DF, p-value: 4.041e-12
##
##
## Estimation results for equation sd.rep:
## =======================================
## sd.rep = sd.inf.l1 + sd.gdp.l1 + sd.M1.l1 + sd.rep.l1 + sd.inf.l2 + sd.gdp.l2 + sd.M1.l2 + sd.rep.l2 + sd.inf.l3 + sd.gdp.l3 + sd.M1.l3 + sd.rep.l3 + const + trend
##
## Estimate Std. Error t value Pr(>|t|)
## sd.inf.l1 -3.851e-02 6.274e-02 -0.614 0.5424
## sd.gdp.l1 6.105e-07 6.806e-07 0.897 0.3744
## sd.M1.l1 -3.686e-07 4.048e-07 -0.911 0.3673
## sd.rep.l1 9.295e-03 1.457e-01 0.064 0.9494
## sd.inf.l2 5.299e-02 8.010e-02 0.662 0.5116
## sd.gdp.l2 1.034e-06 7.788e-07 1.328 0.1909
## sd.M1.l2 -3.288e-07 5.678e-07 -0.579 0.5654
## sd.rep.l2 -3.791e-01 1.420e-01 -2.669 0.0105 *
## sd.inf.l3 1.609e-03 7.026e-02 0.023 0.9818
## sd.gdp.l3 7.938e-07 6.712e-07 1.183 0.2431
## sd.M1.l3 7.105e-09 4.072e-07 0.017 0.9862
## sd.rep.l3 -2.343e-01 1.524e-01 -1.537 0.1312
## const 9.539e-03 1.240e-01 0.077 0.9390
## trend -5.751e-04 3.341e-03 -0.172 0.8641
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.4355 on 45 degrees of freedom
## Multiple R-Squared: 0.2517, Adjusted R-squared: 0.03558
## F-statistic: 1.165 on 13 and 45 DF, p-value: 0.3351
##
##
##
## Covariance matrix of residuals:
## sd.inf sd.gdp sd.M1 sd.rep
## sd.inf 1.045e+00 -9460 4.062e+04 -4.973e-04
## sd.gdp -9.460e+03 2734611804 -1.010e+08 1.295e+03
## sd.M1 4.062e+04 -100963026 2.190e+10 1.162e+04
## sd.rep -4.973e-04 1295 1.162e+04 1.897e-01
##
## Correlation matrix of residuals:
## sd.inf sd.gdp sd.M1 sd.rep
## sd.inf 1.000000 -0.17700 0.26858 -0.001117
## sd.gdp -0.177002 1.00000 -0.01305 0.056857
## sd.M1 0.268578 -0.01305 1.00000 0.180313
## sd.rep -0.001117 0.05686 0.18031 1.000000serial.test(VAR, lags.pt=10, type="PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object VAR
## Chi-squared = 151.74, df = 112, p-value = 0.007412Forecast
#library(tsDyn)
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:MASS':
##
## select## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forecast)## Warning: package 'forecast' was built under R version 3.6.3library(ggplot2)
VAR.f <- predict(VAR,n.ahead=4)
VAR.f## $sd.inf
## fcst lower upper CI
## [1,] 2.004818 0.001679178 4.007956 2.003138
## [2,] 1.005095 -1.570513838 3.580704 2.575609
## [3,] -1.294545 -3.921476692 1.332386 2.626932
## [4,] -1.207838 -3.951789114 1.536114 2.743951
##
## $sd.gdp
## fcst lower upper CI
## [1,] -95535.72 -198029.13 6957.687 102493.4
## [2,] -194600.74 -333432.73 -55768.743 138832.0
## [3,] 57014.84 -86327.37 200357.048 143342.2
## [4,] 216812.36 68986.25 364638.469 147826.1
##
## $sd.M1
## fcst lower upper CI
## [1,] 122292.45 -167736.5 412321.4 290029.0
## [2,] -154529.60 -612082.4 303023.2 457552.8
## [3,] 63283.70 -448391.9 574959.3 511675.6
## [4,] -78327.44 -596936.7 440281.9 518609.3
##
## $sd.rep
## fcst lower upper CI
## [1,] 0.08561739 -0.7679576 0.9391924 0.8535750
## [2,] -0.11195840 -0.9814517 0.7575349 0.8694933
## [3,] -0.02429291 -0.9515375 0.9029517 0.9272446
## [4,] -0.11386996 -1.0653604 0.8376205 0.9514905plot(VAR.f, nc = 2, xlab = "")
fanchart(VAR.f, nc = 2, xlab = "")
forecast(VAR) %>% autoplot() + xlab("Year")
Arranging the forecasts (from stationary series) to be able to arrive at original series
#View(VAR.f)
C1 <- (VAR.f$fcst$sd.inf[,1])
C2 <- (VAR.f$fcst$sd.gdp[,1])
C3 <- (VAR.f$fcst$sd.M1[,1])
C4 <- (VAR.f$fcst$sd.rep[,1])
fd.inf <-diff(inf, differences = 1)
fd.gdp <-diff(GDP, differences = 1)
fd.M1 <- diff(M1, differences = 1)
fd.rep <- diff(rep, differences = 1)
VAR_fd <- cbind(fd.inf, fd.gdp, fd.M1, fd.rep)
VAR_ORG_STAT <- cbind (inf, GDP, M1, rep, fd.inf, fd.gdp, fd.M1, fd.rep, sd.inf, sd.gdp, sd.M1, sd.rep)
For_Stat <- cbind(for.inf = C1, for.gdp = C2, for.M1 = C3, for.rep = C4)
View(For_Stat)
STAT <- rbind(VAR_stat, For_Stat)
View(VAR_ORG_STAT)
dediff1 <- For_Stat[1,]+ VAR_ORG_STAT[64,5:8]
dediff2 <- dediff1 + VAR_ORG_STAT[64,1:4]
dediffrow2 <- matrix(0, nrow = 1, ncol = 4)
dediffMatrix <- dediff2
for (i in 2:4) {
dediffrow1 <- For_Stat[i,]+ dediff1
dediffrow2 <- dediffrow1 + dediff2
dediffMatrix <- rbind(dediffMatrix, dediffrow2)
}
dediffMatrix## for.inf for.gdp for.M1 for.rep
## dediffMatrix 127.6021 5330238 3852582 4.985617
## dediffrow2 130.7394 5292642 3862286 4.709276
## dediffrow2 128.4397 5544257 4080100 4.796942
## dediffrow2 128.5264 5704055 3938488 4.707365org_forecast <- rbind(VAR_org, dediffMatrix)
matrix(org_forecast, nrow = 68, ncol = 4)## [,1] [,2] [,3] [,4]
## [1,] 42.12162 762777.5 578693.9 6.000000
## [2,] 43.17190 726524.3 577137.7 6.000000
## [3,] 43.44828 750800.9 578205.9 6.000000
## [4,] 43.55884 847248.6 610005.9 6.000000
## [5,] 43.80759 878664.4 649765.6 6.000000
## [6,] 44.77495 826929.7 678003.6 6.000000
## [7,] 45.63175 845056.2 709833.2 6.250000
## [8,] 45.52120 964492.2 737321.9 6.500000
## [9,] 46.41377 1012499.2 826389.4 6.750000
## [10,] 47.56136 940813.3 806531.4 7.000000
## [11,] 48.58145 992085.5 847960.1 7.250000
## [12,] 48.70896 1125619.0 866989.5 7.500000
## [13,] 49.34651 1184566.8 967925.5 7.500000
## [14,] 50.74912 1112943.2 941001.9 7.750000
## [15,] 51.25916 1136971.8 983510.6 7.750000
## [16,] 51.76921 1305763.4 1023955.5 7.750000
## [17,] 53.17182 1371470.0 1155809.7 8.500000
## [18,] 55.33950 1317436.4 1108140.9 9.000000
## [19,] 56.48709 1349763.1 1143391.0 7.500000
## [20,] 56.61460 1447672.9 1128674.9 5.500000
## [21,] 57.88970 1447519.8 1259670.5 4.750000
## [22,] 61.84252 1415428.2 1247720.8 4.750000
## [23,] 64.01020 1479644.0 1311294.1 4.750000
## [24,] 65.28530 1680047.9 1330955.5 5.000000
## [25,] 65.79535 1824847.3 1489268.2 5.250000
## [26,] 68.21804 1723210.8 1486650.3 6.000000
## [27,] 69.87568 1765504.3 1519831.2 6.250000
## [28,] 71.15078 2020575.1 1591773.2 6.750000
## [29,] 71.66082 2181264.2 1638345.4 6.750000
## [30,] 74.46605 2043496.3 1592367.3 7.500000
## [31,] 75.74115 2029469.2 1582811.6 8.250000
## [32,] 76.25120 2244852.2 1697814.9 8.500000
## [33,] 78.92891 2418511.1 1737394.0 8.500000
## [34,] 81.73414 2313203.1 1807960.6 8.000000
## [35,] 83.39178 2356953.3 1746870.6 8.000000
## [36,] 85.17692 2548377.8 1812874.4 8.000000
## [37,] 87.34460 2725478.9 1897526.1 7.500000
## [38,] 90.53236 2565374.4 1971440.2 7.250000
## [39,] 92.18999 2686203.2 1910508.2 7.500000
## [40,] 91.04240 2921070.0 1989685.8 7.750000
## [41,] 93.33758 3060874.1 2059762.0 8.000000
## [42,] 96.65285 2914393.2 2147673.6 8.000000
## [43,] 96.78036 3053982.0 2114241.2 8.000000
## [44,] 97.03538 3184917.1 2188575.9 8.000000
## [45,] 98.82053 3314666.8 2292403.8 7.500000
## [46,] 101.11572 3243543.0 2344883.4 7.250000
## [47,] 103.02837 3355884.0 2349066.3 6.750000
## [48,] 102.51833 3476632.0 2458153.6 6.750000
## [49,] 104.94103 3695814.0 2602538.3 6.750000
## [50,] 106.47115 3639774.6 2665296.9 6.500000
## [51,] 105.83360 3754117.0 2828550.4 6.500000
## [52,] 104.94103 3875006.7 1994572.7 6.250000
## [53,] 106.47115 4122770.7 2681956.8 6.250000
## [54,] 109.02136 4007453.0 2701004.6 6.250000
## [55,] 109.78642 4172122.0 2873033.6 6.000000
## [56,] 109.91393 4326105.0 2886324.9 6.000000
## [57,] 110.67899 4592624.1 3267331.3 6.000000
## [58,] 115.14186 4551224.0 3189125.3 6.250000
## [59,] 115.39688 4647831.0 3291706.1 6.500000
## [60,] 117.69206 4817413.0 3278667.9 6.500000
## [61,] 120.11476 4954769.9 3710463.6 6.250000
## [62,] 122.53746 4918228.0 3549635.7 5.750000
## [63,] 125.34268 4920694.0 3646405.8 5.400000
## [64,] 125.47000 5173234.0 3688347.5 5.150000
## [65,] 127.60213 5330238.3 3852581.7 4.985617
## [66,] 130.73936 5292641.8 3862286.3 4.709276
## [67,] 128.43972 5544257.4 4080099.6 4.796942
## [68,] 128.52643 5704054.9 3938488.5 4.707365View(org_forecast)
F.INF <- ts(org_forecast[,1],start = c(2004,1), end = c(2020,4), frequency = 4)
F.GDP <- ts(org_forecast[,2],start = c(2004,1), end = c(2020,4), frequency = 4)
F.M1 <- ts(org_forecast[,3],start = c(2004,1), end = c(2020,4), frequency = 4)
F.REP <- ts(org_forecast[,4],start = c(2004,1), end = c(2020,4), frequency = 4)
org_forecast <- window(ts.union(F.INF, F.GDP, F.M1, F.REP), start = c(2004,1), end = c(2020,4))
plot.ts(org_forecast, nc = 2, xlab = "", col = "blue", main="Forecasted Values along with historical time series")
#lines(VAR_org, nc = 2, xlab = "", col = "blue")
zoom <- ts(org_forecast[-(1:64),], start = c(2020,1), end =c(2020,4), frequency = 4)
#ggplot(zoom, aes(x= (start = c(2018,1), end =c(2020,4), frequency = 4), y = zoom, colour = "red"))
plot.ts(zoom, nc = 2, xlab = "", type = "o",col = "red", main="Only Forecasted Values") 
######author: “Ramya Emandi”