library("Quandl")
library("zoo")
library("xts")
library("urca")
library("vars")
library("dygraphs")
library("stargazer")
library("forecast")
Quandl.api_key("sfE3eBMHtyGwpQTDwkWQ")

Import data

GDP<-Quandl("FRED/GDPC96",collapse="quarter", type="zoo")
D<-Quandl("FRED/GDPDEF",type="zoo")
SP<-Quandl("YAHOO/INDEX_GSPC/CLOSE",collapse="quarterly",type="zoo")

a-Estimate a bivariate reduced form VAR

dGDP<-diff(log(GDP))
dD<-diff(log(D))
dSP<-diff(log(SP))
Y<-dSP-dD
plot(dGDP, type= "l",xlab= "Years", ylab="DRGDP", main="", major.format="%Y Q%q")

plot(Y, type= "l",xlab= "Years", ylab="Y", main="", major.format="%Y Q%q")

y.Q <- cbind(dGDP, Y)
y.Q <- window(y.Q, start="1961 Q1", end="2016 Q4")

VARselect(y.Q, lag.max=8, type="const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2 
## 
## $criteria
##                    1             2             3             4
## AIC(n) -1.478987e+01 -1.487484e+01 -1.485274e+01 -1.483480e+01
## HQ(n)  -1.475199e+01 -1.481171e+01 -1.476436e+01 -1.472117e+01
## SC(n)  -1.469611e+01 -1.471858e+01 -1.463397e+01 -1.455353e+01
## FPE(n)  3.774357e-07  3.466931e-07  3.544517e-07  3.608860e-07
##                    5             6             7             8
## AIC(n) -1.481497e+01 -1.479206e+01 -1.478291e+01 -1.475799e+01
## HQ(n)  -1.467609e+01 -1.462792e+01 -1.459352e+01 -1.454334e+01
## SC(n)  -1.447120e+01 -1.438577e+01 -1.431412e+01 -1.422670e+01
## FPE(n)  3.681415e-07  3.767189e-07  3.802409e-07  3.899152e-07

AIC at lag 2 is minimum (-1.485738e+01 ), so it will be used to calculate VAR in reduced form.

var1 <- VAR(y.Q, p=2, type="const")
var1
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation dGDP: 
## ========================================= 
## Call:
## dGDP = dGDP.l1 + Y.l1 + dGDP.l2 + Y.l2 + const 
## 
##     dGDP.l1        Y.l1     dGDP.l2        Y.l2       const 
## 0.216477458 0.016139290 0.168937034 0.024020293 0.004297188 
## 
## 
## Estimated coefficients for equation Y: 
## ====================================== 
## Call:
## Y = dGDP.l1 + Y.l1 + dGDP.l2 + Y.l2 + const 
## 
##      dGDP.l1         Y.l1      dGDP.l2         Y.l2        const 
##  0.603703071  0.101999176 -0.884164506 -0.078228109  0.009605538
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dGDP, Y 
## Deterministic variables: const 
## Sample size: 222 
## Log Likelihood: 1030.366 
## Roots of the characteristic polynomial:
## 0.3746 0.3243 0.257 0.257
## Call:
## VAR(y = y.Q, p = 2, type = "const")
## 
## 
## Estimation results for equation dGDP: 
## ===================================== 
## dGDP = dGDP.l1 + Y.l1 + dGDP.l2 + Y.l2 + const 
## 
##         Estimate Std. Error t value Pr(>|t|)    
## dGDP.l1 0.216477   0.064678   3.347 0.000963 ***
## Y.l1    0.016139   0.006032   2.675 0.008032 ** 
## dGDP.l2 0.168937   0.063182   2.674 0.008070 ** 
## Y.l2    0.024020   0.006122   3.924 0.000117 ***
## const   0.004297   0.000733   5.863 1.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.007215 on 217 degrees of freedom
## Multiple R-Squared: 0.2373,  Adjusted R-squared: 0.2232 
## F-statistic: 16.88 on 4 and 217 DF,  p-value: 4.602e-12 
## 
## 
## Estimation results for equation Y: 
## ================================== 
## Y = dGDP.l1 + Y.l1 + dGDP.l2 + Y.l2 + const 
## 
##          Estimate Std. Error t value Pr(>|t|)
## dGDP.l1  0.603703   0.727743   0.830    0.408
## Y.l1     0.101999   0.067875   1.503    0.134
## dGDP.l2 -0.884165   0.710912  -1.244    0.215
## Y.l2    -0.078228   0.068882  -1.136    0.257
## const    0.009606   0.008247   1.165    0.245
## 
## 
## Residual standard error: 0.08118 on 217 degrees of freedom
## Multiple R-Squared: 0.02376, Adjusted R-squared: 0.005766 
## F-statistic:  1.32 on 4 and 217 DF,  p-value: 0.2633 
## 
## 
## 
## Covariance matrix of residuals:
##           dGDP         Y
## dGDP 5.205e-05 9.627e-05
## Y    9.627e-05 6.590e-03
## 
## Correlation matrix of residuals:
##        dGDP      Y
## dGDP 1.0000 0.1644
## Y    0.1644 1.0000

The correlation between residuals of dGDP and Y is 0.1644. This means that there is some sort of causal relationship between two variables.

b- Granger causality tests

causality(var1, cause="dGDP")
## $Granger
## 
##  Granger causality H0: dGDP do not Granger-cause Y
## 
## data:  VAR object var1
## F-Test = 0.87077, df1 = 2, df2 = 434, p-value = 0.4194
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: dGDP and Y
## 
## data:  VAR object var1
## Chi-squared = 5.8402, df = 1, p-value = 0.01566
causality(var1, cause="Y")
## $Granger
## 
##  Granger causality H0: Y do not Granger-cause dGDP
## 
## data:  VAR object var1
## F-Test = 12.1, df1 = 2, df2 = 434, p-value = 7.696e-06
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Y and dGDP
## 
## data:  VAR object var1
## Chi-squared = 5.8402, df = 1, p-value = 0.01566

There is no causality between dGDP and Y in both directions.In other words S&P500 do not Granger cause GDP.

c- A restricted VAR model

mat.r <- matrix(1, nrow=2, ncol=5)
mat.r[2, c(1,3)] <- 0
mat.r
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    1    1    1    1
## [2,]    0    1    0    1    1
var1.r <- restrict(var1, method="manual", resmat=mat.r)
summary(var1.r)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dGDP, Y 
## Deterministic variables: const 
## Sample size: 222 
## Log Likelihood: 1029.454 
## Roots of the characteristic polynomial:
## 0.5333 0.3168 0.2774 0.2774
## Call:
## VAR(y = y.Q, p = 2, type = "const")
## 
## 
## Estimation results for equation dGDP: 
## ===================================== 
## dGDP = dGDP.l1 + Y.l1 + dGDP.l2 + Y.l2 + const 
## 
##         Estimate Std. Error t value Pr(>|t|)    
## dGDP.l1 0.216477   0.064678   3.347 0.000963 ***
## Y.l1    0.016139   0.006032   2.675 0.008032 ** 
## dGDP.l2 0.168937   0.063182   2.674 0.008070 ** 
## Y.l2    0.024020   0.006122   3.924 0.000117 ***
## const   0.004297   0.000733   5.863 1.68e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.007215 on 217 degrees of freedom
## Multiple R-Squared: 0.5877,  Adjusted R-squared: 0.5782 
## F-statistic: 61.86 on 5 and 217 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation Y: 
## ================================== 
## Y = Y.l1 + Y.l2 + const 
## 
##        Estimate Std. Error t value Pr(>|t|)
## Y.l1   0.107702   0.067362   1.599    0.111
## Y.l2  -0.076951   0.067129  -1.146    0.253
## const  0.007432   0.005489   1.354    0.177
## 
## 
## Residual standard error: 0.08113 on 219 degrees of freedom
## Multiple R-Squared: 0.02454, Adjusted R-squared: 0.01118 
## F-statistic: 1.836 on 3 and 219 DF,  p-value: 0.1414 
## 
## 
## 
## Covariance matrix of residuals:
##           dGDP         Y
## dGDP 5.205e-05 9.627e-05
## Y    9.627e-05 6.643e-03
## 
## Correlation matrix of residuals:
##        dGDP      Y
## dGDP 1.0000 0.1637
## Y    0.1637 1.0000
var1.r$restrictions
##      dGDP.l1 Y.l1 dGDP.l2 Y.l2 const
## dGDP       1    1       1    1     1
## Y          0    1       0    1     1
Acoef(var1.r)
## [[1]]
##        dGDP.l1       Y.l1
## dGDP 0.2164775 0.01613929
## Y    0.0000000 0.10770248
## 
## [[2]]
##       dGDP.l2        Y.l2
## dGDP 0.168937  0.02402029
## Y    0.000000 -0.07695147

d- Choleski decomposition

var1.irfs <- irf(var1, n.ahead=10)
par(mfcol=c(2,2), cex=0.6)
plot(var1.irfs, plot.type="single")

var1.fevd <- fevd(var1, n.ahead=10)
var1.fevd[[1]][c(1,3,6,10),]
##           dGDP         Y
## [1,] 1.0000000 0.0000000
## [2,] 0.8927329 0.1072671
## [3,] 0.8813672 0.1186328
## [4,] 0.8812738 0.1187262
var1.fevd[[2]][c(1,3,6,10),]
##            dGDP         Y
## [1,] 0.02701807 0.9729819
## [2,] 0.03617662 0.9638234
## [3,] 0.03661349 0.9633865
## [4,] 0.03661909 0.9633809
par(mfcol=c(2,2), cex=0.6)
plot(var1.fevd)

par(mfcol=c(2,2), cex=0.6)
plot(var1.fevd, addbars=8)

y.Q.ord2 <- cbind(Y, dGDP)
y.Q.ord2 <- window(y.Q.ord2, start="1961 Q1", end="2016 Q4")
var1.ord2 <- VAR(y.Q.ord2, p=1, type="const")
var1.ord1 <- VAR(y.Q, p=1, type="const")
var1.irfs.ord1 <- irf(var1.ord1, n.ahead=10)
var1.irfs.ord2 <- irf(var1.ord2, n.ahead=10)
par(mfcol=c(2,2), cex=0.6)
plot(var1.irfs.ord1, plot.type="single")

par(mfcol=c(2,2), cex=0.6)
plot(var1.irfs.ord2, plot.type="single")

The order does not matter much in this case.

e- Index for the United States

US <- Quandl("FRED/USSLIND", collapse="quarter", type="zoo") 
y.Q3 <- cbind(dGDP, Y, US)
y.Q3 <- window(y.Q3, start="1982 Q1", end="2016 Q4")
VARselect(y.Q3, lag.max=8, type="const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      1      1      1      1 
## 
## $criteria
##                    1             2             3             4
## AIC(n) -1.773523e+01 -1.771439e+01 -1.765990e+01 -1.756538e+01
## HQ(n)  -1.762874e+01 -1.752802e+01 -1.739366e+01 -1.721928e+01
## SC(n)  -1.747316e+01 -1.725576e+01 -1.700472e+01 -1.671365e+01
## FPE(n)  1.984771e-08  2.027073e-08  2.141814e-08  2.356595e-08
##                    5             6             7             8
## AIC(n) -1.748128e+01 -1.743725e+01 -1.748607e+01 -1.748680e+01
## HQ(n)  -1.705530e+01 -1.693140e+01 -1.690035e+01 -1.682121e+01
## SC(n)  -1.643299e+01 -1.619240e+01 -1.604467e+01 -1.584884e+01
## FPE(n)  2.567661e-08  2.689815e-08  2.570303e-08  2.579947e-08
var3 <- VAR(y.Q3, p=1, type="const")
summary(var3)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dGDP, Y, US 
## Deterministic variables: const 
## Sample size: 139 
## Log Likelihood: 649.381 
## Roots of the characteristic polynomial:
## 0.8175 0.05397 0.0003733
## Call:
## VAR(y = y.Q3, p = 1, type = "const")
## 
## 
## Estimation results for equation dGDP: 
## ===================================== 
## dGDP = dGDP.l1 + Y.l1 + US.l1 + const 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## dGDP.l1 0.0735616  0.0979529   0.751   0.4540    
## Y.l1    0.0107875  0.0056310   1.916   0.0575 .  
## US.l1   0.0035767  0.0007624   4.691 6.57e-06 ***
## const   0.0015192  0.0008108   1.874   0.0631 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.005138 on 135 degrees of freedom
## Multiple R-Squared: 0.3478,  Adjusted R-squared: 0.3333 
## F-statistic:    24 on 3 and 135 DF,  p-value: 1.648e-12 
## 
## 
## Estimation results for equation Y: 
## ================================== 
## Y = dGDP.l1 + Y.l1 + US.l1 + const 
## 
##         Estimate Std. Error t value Pr(>|t|)
## dGDP.l1 0.230105   1.524378   0.151    0.880
## Y.l1    0.055120   0.087632   0.629    0.530
## US.l1   0.005331   0.011865   0.449    0.654
## const   0.006458   0.012618   0.512    0.610
## 
## 
## Residual standard error: 0.07995 on 135 degrees of freedom
## Multiple R-Squared: 0.01024, Adjusted R-squared: -0.01176 
## F-statistic: 0.4653 on 3 and 135 DF,  p-value: 0.7069 
## 
## 
## Estimation results for equation US: 
## =================================== 
## US = dGDP.l1 + Y.l1 + US.l1 + const 
## 
##         Estimate Std. Error t value Pr(>|t|)    
## dGDP.l1 12.81179    8.00938   1.600 0.112026    
## Y.l1     1.40438    0.46043   3.050 0.002754 ** 
## US.l1    0.74244    0.06234  11.909  < 2e-16 ***
## const    0.24343    0.06630   3.672 0.000346 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.4201 on 135 degrees of freedom
## Multiple R-Squared: 0.7551,  Adjusted R-squared: 0.7497 
## F-statistic: 138.8 on 3 and 135 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           dGDP         Y       US
## dGDP 2.639e-05 0.0000737 0.001187
## Y    7.370e-05 0.0063923 0.009232
## US   1.187e-03 0.0092322 0.176471
## 
## Correlation matrix of residuals:
##        dGDP      Y     US
## dGDP 1.0000 0.1794 0.5500
## Y    0.1794 1.0000 0.2749
## US   0.5500 0.2749 1.0000

f- VAR from (e)

par(mfcol=c(2,2), cex=0.6)
var3.f <- predict(var3, n.ahead=4, ci=0.9)
plot(var3.f)

par(mfcol=c(2,2), cex=0.6)
fanchart(var3.f)

var3.f
## $dGDP
##             fcst         lower      upper          CI
## [1,] 0.007731432 -0.0007190564 0.01618192 0.008450488
## [2,] 0.007633004 -0.0015240471 0.01679006 0.009157051
## [3,] 0.007555994 -0.0019894145 0.01710140 0.009545409
## [4,] 0.007493178 -0.0023010954 0.01728745 0.009794273
## 
## $Y
##            fcst      lower     upper        CI
## [1,] 0.01738547 -0.1141241 0.1488950 0.1315095
## [2,] 0.01718090 -0.1146989 0.1490607 0.1318798
## [3,] 0.01704627 -0.1149203 0.1490129 0.1319666
## [4,] 0.01693810 -0.1150827 0.1489589 0.1320208
## 
## $US
##          fcst     lower    upper        CI
## [1,] 1.497875 0.8068976 2.188852 0.6909771
## [2,] 1.478985 0.5267156 2.431255 0.9522696
## [3,] 1.463412 0.3669462 2.559879 1.0964663
## [4,] 1.450675 0.2674320 2.633918 1.1832429

In VAR model the forecast for first quarter of 2017 is 0.7% which is lower than the FRBNY Staff Nowcast (2.7%), Economic Forecasting Survey (1.4%) and is higher than Atlanta FED forecast(0.5%).