#download packages
library("Quandl")
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library("vars")
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library("gdata")
## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
##
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
##
## Attaching package: 'gdata'
## The following objects are masked from 'package:xts':
##
## first, last
## The following object is masked from 'package:stats':
##
## nobs
## The following object is masked from 'package:utils':
##
## object.size
library("stargazer")
##
## Please cite as:
## Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2. http://CRAN.R-project.org/package=stargazer
library("tseries")
Quandl.auth(" 3HQPzReHC1WfzuSG41gy")
## Warning: 'Quandl.auth' is deprecated.
## Use 'Quandl.api_key' instead.
## See help("Deprecated")
#download data for gdp, gdpd, and sp
gdp<-Quandl("FRED/GDPC1", type="zoo")
gdpd<-Quandl("FRED/GDPDEF", type="zoo")
sp<-Quandl("YAHOO/INDEX_GSPC", collapse="quarterly", type="zoo")
gdp<-window(gdp,start="1950 Q1", end="2015 Q4")
gdpd<-window(gdpd,start="1950 Q1", end="2015 Q4")
sp<-window(sp,start="1950 Q1", end="2015 Q4")
#obtain log for gdp, gdpd, and sp, and first difference of log
#gdp, gdpd, and sp
lgdp<-log(gdp)
dlgdp<-diff(lgdp,1)*100
lgdpd<-log(gdpd)
dlgdpd<-diff(lgdpd,1)*100
lsp<-log(sp$Close)
dlsp<-diff(lsp,1)*100
par(mfrow=c(1,2))
plot(gdp,xlab="Time", ylab="gdp")
plot(dlgdp,xlab="Time", ylab="dlgdp")
par(mfrow=c(1,2))
plot(gdpd,xlab="Time", ylab="gdpd")
plot(dlgdpd,xlab="Time", ylab="dlgdpd")
par(mfrow=c(1,2))
plot(lsp,xlab="Time", ylab="lsp")
plot(dlsp,xlab="Time", ylab="dlsp")
#a)Estimate a bivariate reduced form VAR for yt = (y1,t,y2,t)′ for the period 1951Q1-2014Q4, use information criteria to select number of lags. How large is the correlation of residuals in the model, and what are the implications?
# form dataset for VAR model
y<-cbind(dlsp,dlgdpd,dlgdp)
y<-na.trim(y)
y2<-y$dlsp-y$dlgdpd
y1<-y$dlgdp
y3<-cbind(y1,y2)
y3<-sweep(y3,2,apply(y3,2,mean))
plot(y3)
# choose lag by using AIC
y3<-window(y3,start="1951 Q1", end="2014 Q4")
VARselect(y3, lag.max=10)
## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 1 2
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) 3.781582 3.725213 3.736946 3.763633 3.781136 3.802590
## HQ(n) 3.816008 3.782589 3.817271 3.866909 3.907362 3.951766
## SC(n) 3.867078 3.867706 3.936436 4.020121 4.094621 4.173073
## FPE(n) 43.885536 41.480544 41.970897 43.107576 43.871083 44.825965
## 7 8 9 10
## AIC(n) 3.811918 3.840491 3.860155 3.877568
## HQ(n) 3.984044 4.035568 4.078181 4.118544
## SC(n) 4.239397 4.324968 4.401629 4.476039
## FPE(n) 45.250817 46.568888 47.501981 48.346834
# reduced model VAR(1)
var1<-VAR(y3,p=2,type="const")
summary(var1)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1, y2
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: -1184.933
## Roots of the characteristic polynomial:
## 0.4454 0.4454 0.304 0.304
## Call:
## VAR(y = y3, p = 2, type = "const")
##
##
## Estimation results for equation y1:
## ===================================
## y1 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1.l1 0.245040 0.061477 3.986 8.84e-05 ***
## y2.l1 0.023097 0.006405 3.606 0.000375 ***
## y1.l2 0.079349 0.059543 1.333 0.183872
## y2.l2 0.026300 0.006572 4.002 8.29e-05 ***
## const -0.018032 0.050452 -0.357 0.721085
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.8038 on 249 degrees of freedom
## Multiple R-Squared: 0.2309, Adjusted R-squared: 0.2185
## F-statistic: 18.68 on 4 and 249 DF, p-value: 1.91e-13
##
##
## Estimation results for equation y2:
## ===================================
## y2 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1.l1 0.256054 0.609111 0.420 0.6746
## y2.l1 0.106501 0.063457 1.678 0.0945 .
## y1.l2 -0.819353 0.589953 -1.389 0.1661
## y2.l2 -0.040573 0.065115 -0.623 0.5338
## const -0.007916 0.499874 -0.016 0.9874
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 7.964 on 249 degrees of freedom
## Multiple R-Squared: 0.02087, Adjusted R-squared: 0.00514
## F-statistic: 1.327 on 4 and 249 DF, p-value: 0.2605
##
##
##
## Covariance matrix of residuals:
## y1 y2
## y1 0.646 0.87
## y2 0.870 63.42
##
## Correlation matrix of residuals:
## y1 y2
## y1 1.0000 0.1359
## y2 0.1359 1.0000
#B)(b) Run the Granger causality tests for both variables What do the results suggest about the predictive power of the two variables?
# Granger causality
causality(var1,cause="y1")
## $Granger
##
## Granger causality H0: y1 do not Granger-cause y2
##
## data: VAR object var1
## F-Test = 0.96654, df1 = 2, df2 = 498, p-value = 0.3811
##
##
## $Instant
##
## H0: No instantaneous causality between: y1 and y2
##
## data: VAR object var1
## Chi-squared = 4.6072, df = 1, p-value = 0.03184
causality(var1,cause="y2")
## $Granger
##
## Granger causality H0: y2 do not Granger-cause y1
##
## data: VAR object var1
## F-Test = 15.762, df1 = 2, df2 = 498, p-value = 2.305e-07
##
##
## $Instant
##
## H0: No instantaneous causality between: y2 and y1
##
## data: VAR object var1
## Chi-squared = 4.6072, df = 1, p-value = 0.03184
#there is Granger Causality from S&P to GDP. so, the predictive power of two variables is S&P.
#C)Estimate two restricted VAR models - in first remove lags based on the Granger causality test from (b), in the second one remove all variables which have t-statistics lower than 2.
#based on Granger causality test eliminate lags of gdp
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: y1, y2
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: -1185.934
## Roots of the characteristic polynomial:
## 0.4297 0.2121 0.2121 0.1847
## Call:
## VAR(y = y3, p = 2, type = "const")
##
##
## Estimation results for equation y1:
## ===================================
## y1 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1.l1 0.245040 0.061477 3.986 8.84e-05 ***
## y2.l1 0.023097 0.006405 3.606 0.000375 ***
## y1.l2 0.079349 0.059543 1.333 0.183872
## y2.l2 0.026300 0.006572 4.002 8.29e-05 ***
## const -0.018032 0.050452 -0.357 0.721085
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.8038 on 249 degrees of freedom
## Multiple R-Squared: 0.2315, Adjusted R-squared: 0.2161
## F-statistic: 15 on 5 and 249 DF, p-value: 7.118e-13
##
##
## Estimation results for equation y2:
## ===================================
## y2 = y2.l1 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y2.l1 0.11089 0.06304 1.759 0.0798 .
## y2.l2 -0.04500 0.06304 -0.714 0.4760
## const 0.00313 0.49962 0.006 0.9950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 7.963 on 251 degrees of freedom
## Multiple R-Squared: 0.01327, Adjusted R-squared: 0.001474
## F-statistic: 1.125 on 3 and 251 DF, p-value: 0.3395
##
##
##
## Covariance matrix of residuals:
## y1 y2
## y1 0.646 0.87
## y2 0.870 63.91
##
## Correlation matrix of residuals:
## y1 y2
## y1 1.0000 0.1354
## y2 0.1354 1.0000
mat.r2<- matrix(1,nrow=2, ncol=5)
mat.r2[1,c(3)]<-0
mat.r2[2,c(1,2,3,4)]<-0
mat.r2
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 1 0 1 1
## [2,] 0 0 0 0 1
var1.r <- restrict(var1, method="manual", resmat=mat.r)
var1.r
##
## VAR Estimation Results:
## =======================
##
## Estimated coefficients for equation y1:
## =======================================
## Call:
## y1 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## y1.l1 y2.l1 y1.l2 y2.l2 const
## 0.24504017 0.02309677 0.07934933 0.02630040 -0.01803217
##
##
## Estimated coefficients for equation y2:
## =======================================
## Call:
## y2 = y2.l1 + y2.l2 + const
##
## y2.l1 y2.l2 const
## 0.110892062 -0.045000768 0.003129995
summary(var1.r)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1, y2
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: -1185.934
## Roots of the characteristic polynomial:
## 0.4297 0.2121 0.2121 0.1847
## Call:
## VAR(y = y3, p = 2, type = "const")
##
##
## Estimation results for equation y1:
## ===================================
## y1 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1.l1 0.245040 0.061477 3.986 8.84e-05 ***
## y2.l1 0.023097 0.006405 3.606 0.000375 ***
## y1.l2 0.079349 0.059543 1.333 0.183872
## y2.l2 0.026300 0.006572 4.002 8.29e-05 ***
## const -0.018032 0.050452 -0.357 0.721085
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.8038 on 249 degrees of freedom
## Multiple R-Squared: 0.2315, Adjusted R-squared: 0.2161
## F-statistic: 15 on 5 and 249 DF, p-value: 7.118e-13
##
##
## Estimation results for equation y2:
## ===================================
## y2 = y2.l1 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y2.l1 0.11089 0.06304 1.759 0.0798 .
## y2.l2 -0.04500 0.06304 -0.714 0.4760
## const 0.00313 0.49962 0.006 0.9950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 7.963 on 251 degrees of freedom
## Multiple R-Squared: 0.01327, Adjusted R-squared: 0.001474
## F-statistic: 1.125 on 3 and 251 DF, p-value: 0.3395
##
##
##
## Covariance matrix of residuals:
## y1 y2
## y1 0.646 0.87
## y2 0.870 63.91
##
## Correlation matrix of residuals:
## y1 y2
## y1 1.0000 0.1354
## y2 0.1354 1.0000
# based on t-statistics lower than 2
mat.r2<- matrix(1,nrow=2, ncol=5)
mat.r2[1,c(3)]<-0
mat.r2[2,c(1,2,3,4)]<-0
mat.r2
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 1 0 1 1
## [2,] 0 0 0 0 1
var1.r2 <- restrict(var1, method="manual", resmat=mat.r2)
summary(var1.r2)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1, y2
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: -1188.827
## Roots of the characteristic polynomial:
## 0.2734 0 0 0
## Call:
## VAR(y = y3, p = 2, type = "const")
##
##
## Estimation results for equation y1:
## ===================================
## y1 = y1.l1 + y2.l1 + y2.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1.l1 0.273418 0.057761 4.734 3.70e-06 ***
## y2.l1 0.022632 0.006405 3.533 0.000488 ***
## y2.l2 0.026628 0.006578 4.048 6.88e-05 ***
## const -0.019029 0.050525 -0.377 0.706772
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.805 on 250 degrees of freedom
## Multiple R-Squared: 0.226, Adjusted R-squared: 0.2136
## F-statistic: 18.25 on 4 and 250 DF, p-value: 3.611e-13
##
##
## Estimation results for equation y2:
## ===================================
## y2 = const
##
## Estimate Std. Error t value Pr(>|t|)
## const 0.001059 0.500968 0.002 0.998
##
##
## Residual standard error: 7.984 on 253 degrees of freedom
## Multiple R-Squared: 1.767e-08, Adjusted R-squared: -0.003953
## F-statistic: 4.469e-06 on 1 and 253 DF, p-value: 0.9983
##
##
##
## Covariance matrix of residuals:
## y1 y2
## y1 0.6506 0.8224
## y2 0.8224 64.7701
##
## Correlation matrix of residuals:
## y1 y2
## y1 1.0000 0.1267
## y2 0.1267 1.0000
#D)(d) Use stargazer package to show estimation results from all three models in (a) and (c).
lm2 <- var1$varresult
lm.r <- var1.r$varresult
lm.r2 <- var1.r2$varresult
stargazer(lm2$y1, lm2$y2, lm.r$y1, lm.r$y2, lm.r2$y1, lm.r2$y2,
type="html", column.labels=rep(colnames(y3), 3),
dep.var.labels.include=FALSE)
##
## <table style="text-align:center"><tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="6"><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="6" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>y1</td><td>y2</td><td>y1</td><td>y2</td><td>y1</td><td>y2</td></tr>
## <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td><td>(5)</td><td>(6)</td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">y1.l1</td><td>0.245<sup>***</sup></td><td>0.256</td><td>0.245<sup>***</sup></td><td></td><td>0.273<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.061)</td><td>(0.609)</td><td>(0.061)</td><td></td><td>(0.058)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">y2.l1</td><td>0.023<sup>***</sup></td><td>0.107<sup>*</sup></td><td>0.023<sup>***</sup></td><td>0.111<sup>*</sup></td><td>0.023<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.006)</td><td>(0.063)</td><td>(0.006)</td><td>(0.063)</td><td>(0.006)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">y1.l2</td><td>0.079</td><td>-0.819</td><td>0.079</td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.060)</td><td>(0.590)</td><td>(0.060)</td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">y2.l2</td><td>0.026<sup>***</sup></td><td>-0.041</td><td>0.026<sup>***</sup></td><td>-0.045</td><td>0.027<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.007)</td><td>(0.065)</td><td>(0.007)</td><td>(0.063)</td><td>(0.007)</td><td></td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left">const</td><td>-0.018</td><td>-0.008</td><td>-0.018</td><td>0.003</td><td>-0.019</td><td>0.001</td></tr>
## <tr><td style="text-align:left"></td><td>(0.050)</td><td>(0.500)</td><td>(0.050)</td><td>(0.500)</td><td>(0.051)</td><td>(0.501)</td></tr>
## <tr><td style="text-align:left"></td><td></td><td></td><td></td><td></td><td></td><td></td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>254</td><td>254</td><td>254</td><td>254</td><td>254</td><td>254</td></tr>
## <tr><td style="text-align:left">R<sup>2</sup></td><td>0.231</td><td>0.021</td><td>0.231</td><td>0.013</td><td>0.226</td><td>0.00000</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.219</td><td>0.005</td><td>0.216</td><td>0.001</td><td>0.214</td><td>-0.004</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.804 (df = 249)</td><td>7.964 (df = 249)</td><td>0.804 (df = 249)</td><td>7.963 (df = 251)</td><td>0.805 (df = 250)</td><td>7.984 (df = 253)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>18.685<sup>***</sup> (df = 4; 249)</td><td>1.327 (df = 4; 249)</td><td>15.001<sup>***</sup> (df = 5; 249)</td><td>1.125 (df = 3; 251)</td><td>18.251<sup>***</sup> (df = 4; 250)</td><td>0.00000 (df = 1; 253)</td></tr>
## <tr><td colspan="7" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="6" style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
#(e) Plot IRFs and FEVD for the VAR model based on Choleski decomposition. Afterwards reverse the order of the variables in the VAR and plot IRFs and FEVD based on Choleski decomposition for the
#alternative order. Does the order matter much in this particular case?
#first IRFs
var1.irfs <- irf(var1, n.ahead=40)
par(mfcol=c(2,2), cex=0.6)
plot(var1.irfs, plot.type="single")
#second FEVD
var1.fevd <- fevd(var1, n.ahead=40)
var1.fevd[[1]][c(1,4,8,40),]
## y1 y2
## [1,] 1.0000000 0.0000000
## [2,] 0.8607464 0.1392536
## [3,] 0.8591922 0.1408078
## [4,] 0.8591882 0.1408118
var1.fevd[[2]][c(1,4,8,40),]
## y1 y2
## [1,] 0.01847382 0.9815262
## [2,] 0.02623100 0.9737690
## [3,] 0.02642786 0.9735721
## [4,] 0.02642813 0.9735719
plot(var1.fevd)
plot(var1.fevd, addbars=8)
#reverse the order of the variables in the VAR
y4<-cbind(y2,y1)
y4<-window(y4,start="1951 Q1", end="2014 Q4")
VARselect(y4, lag.max=10)
## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 1 2
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) 3.781582 3.725213 3.736946 3.763633 3.781136 3.802590
## HQ(n) 3.816008 3.782589 3.817271 3.866909 3.907362 3.951766
## SC(n) 3.867078 3.867706 3.936436 4.020121 4.094621 4.173073
## FPE(n) 43.885536 41.480544 41.970897 43.107576 43.871083 44.825965
## 7 8 9 10
## AIC(n) 3.811918 3.840491 3.860155 3.877568
## HQ(n) 3.984044 4.035568 4.078181 4.118544
## SC(n) 4.239397 4.324968 4.401629 4.476039
## FPE(n) 45.250817 46.568888 47.501981 48.346834
# we found p=2 is good , has lowest AIC
var2<-VAR(y4,p=2,type="const")
#IRFs and FEVD based on Choleski decomposition
var2.irfs <- irf(var2, n.ahead=40)
par(mfcol=c(2,2), cex=0.6)
plot(var2.irfs, plot.type="single")
var2.fevd <- fevd(var2, n.ahead=40)
var2.fevd[[1]][c(1,4,8,40),]
## y2 y1
## [1,] 1.0000000 0.000000000
## [2,] 0.9935239 0.006476147
## [3,] 0.9934082 0.006591832
## [4,] 0.9934079 0.006592058
var2.fevd[[2]][c(1,4,8,40),]
## y2 y1
## [1,] 0.01847382 0.9815262
## [2,] 0.18059627 0.8194037
## [3,] 0.18207471 0.8179253
## [4,] 0.18207906 0.8179209
plot(var2.fevd)
