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
Quandl.api_key('4-KG5x_Vo7rXzmZNAHch')
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
gdp1 <- Quandl("FRED/GDPC1", type="zoo")
gdp1 <- window(gdp1, start="1950 Q1", end="2015 Q4")
lgdp1 <- log(gdp1)
ldgdp1 <- diff(lgdp1,1)*100
par (mfrow=c(1,2))
plot(gdp1, xlab= "Period", ylab="Real GDP")
plot(lgdp1, xlab="Periods", ylab="Log GDP")
plot(ldgdp1, xlab ="Periods", ylab="Differenced log GDP")
def <- Quandl("FRED/GDPDEF", type="zoo")
def <- window(def,start="1950 Q1", end="2015 Q4")
ldef<-log(def)
lddef <- diff(ldef,1)*100
par(mfrow=c(1,2))
plot(def,xlab="Periods", ylab="GDP Deflator")
plot(lddef, xlab="Time", ylab="Differenced Log GDP deflator")
sp <- Quandl("YAHOO/INDEX_GSPC", collapse="quarterly", type="zoo")
sp <- window(sp,start="1950 Q1", end="2015 Q4")
lsp <- log(sp$Close)
ldsp <- diff(lsp,1)*100
par(mfrow=c(1,2))
plot(lsp,xlab="Period", ylab="Log of S&P 500")
plot(ldsp,xlab="Periods", ylab="Differenced log S&P")
yt <- cbind(lddef,ldgdp1,ldsp)
yt <- na.trim(yt)
y2t <- yt$ldsp-yt$ldgdp1
y1t <- yt$ldgdp1
ytyt<- cbind(y1t,y2t)
ytyt <-sweep(ytyt,2,apply(ytyt,2,mean))
plot(ytyt)
A
ytyt<- window(ytyt,start="1951 Q1", end="2014 Q4")
VARselect(ytyt, lag.max=10)
## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 2 2
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) 3.762278 3.702907 3.714041 3.742261 3.760352 3.782276
## HQ(n) 3.796704 3.760282 3.794367 3.845537 3.886578 3.931452
## SC(n) 3.847774 3.845400 3.913532 3.998749 4.073837 4.152759
## FPE(n) 43.046497 40.565498 41.020503 42.196047 42.968693 43.924552
## 7 8 9 10
## AIC(n) 3.790697 3.818828 3.837216 3.855601
## HQ(n) 3.962823 4.013904 4.055242 4.096577
## SC(n) 4.218176 4.303304 4.378690 4.454072
## FPE(n) 44.300670 45.570875 46.424734 47.296374
var2 <-VAR(ytyt, p=2, type="const")
B
causality(var2, cause="y1t")
## $Granger
##
## Granger causality H0: y1t do not Granger-cause y2t
##
## data: VAR object var2
## F-Test = 1.6375, df1 = 2, df2 = 498, p-value = 0.1955
##
##
## $Instant
##
## H0: No instantaneous causality between: y1t and y2t
##
## data: VAR object var2
## Chi-squared = 0.28736, df = 1, p-value = 0.5919
causality(var2, cause="y2t")
## $Granger
##
## Granger causality H0: y2t do not Granger-cause y1t
##
## data: VAR object var2
## F-Test = 15.73, df1 = 2, df2 = 498, p-value = 2.375e-07
##
##
## $Instant
##
## H0: No instantaneous causality between: y2t and y1t
##
## data: VAR object var2
## Chi-squared = 0.28736, df = 1, p-value = 0.5919
C
b0 <- diag(2)
diag(b0) <- NA
b0[1, 2] <- NA
b0
## [,1] [,2]
## [1,] NA NA
## [2,] 0 NA
svar2 <- SVAR(var2, estmethod="direct", Amat=b0, hessian=TRUE, method="BFGS")
## Warning in SVAR(var2, estmethod = "direct", Amat = b0, hessian = TRUE,
## method = "BFGS"): The A-model is just identified. No test possible.
summary(svar2)
##
## SVAR Estimation Results:
## ========================
##
## Call:
## SVAR(x = var2, estmethod = "direct", Amat = b0, hessian = TRUE,
## method = "BFGS")
##
## Type: A-model
## Sample size: 254
## Log Likelihood: -1187.027
## Method: direct
## Number of iterations: 65
## Convergence code: 0
##
## Estimated A matrix:
## y1t y2t
## y1t 1.245 0.004316
## y2t 0.000 0.128174
##
## Estimated standard errors for A matrix:
## y1t y2t
## y1t 0.05523 0.008045
## y2t 0.00000 0.005687
##
## Estimated B matrix:
## y1t y2t
## y1t 1 0
## y2t 0 1
##
## Covariance matrix of reduced form residuals (*100):
## y1t y2t
## y1t 64.62 -21.1
## y2t -21.10 6086.9
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
var2.r <- restrict(var2, method="manual", resmat=mat.r)
summary(var2.r)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1t, y2t
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: -1183.64
## Roots of the characteristic polynomial:
## 0.4915 0.2835 0.2835 0.2215
## Call:
## VAR(y = ytyt, p = 2, type = "const")
##
##
## Estimation results for equation y1t:
## ====================================
## y1t = y1t.l1 + y2t.l1 + y1t.l2 + y2t.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1t.l1 0.270010 0.061101 4.419 1.48e-05 ***
## y2t.l1 0.023622 0.006470 3.651 0.000319 ***
## y1t.l2 0.108866 0.059725 1.823 0.069533 .
## y2t.l2 0.026592 0.006637 4.007 8.14e-05 ***
## const -0.017890 0.050458 -0.355 0.723223
## ---
## 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.2313, Adjusted R-squared: 0.2159
## F-statistic: 14.99 on 5 and 249 DF, p-value: 7.315e-13
##
##
## Estimation results for equation y2t:
## ====================================
## y2t = y2t.l1 + y2t.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y2t.l1 0.07466 0.06287 1.187 0.236
## y2t.l2 -0.08036 0.06286 -1.278 0.202
## const 0.01872 0.49078 0.038 0.970
##
##
## Residual standard error: 7.822 on 251 degrees of freedom
## Multiple R-Squared: 0.01123, Adjusted R-squared: -0.0005878
## F-statistic: 0.9503 on 3 and 251 DF, p-value: 0.4169
##
##
##
## Covariance matrix of residuals:
## y1t y2t
## y1t 0.6462 0.2111
## y2t 0.2111 61.6710
##
## Correlation matrix of residuals:
## y1t y2t
## y1t 1.00000 0.03344
## y2t 0.03344 1.00000
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
var2.r2 <- restrict(var2, method="manual", resmat=mat.r2)
summary(var2.r2)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1t, y2t
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: -1186.845
## Roots of the characteristic polynomial:
## 0.3111 0 0 0
## Call:
## VAR(y = ytyt, p = 2, type = "const")
##
##
## Estimation results for equation y1t:
## ====================================
## y1t = y1t.l1 + y2t.l1 + y2t.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1t.l1 0.311101 0.057054 5.453 1.19e-07 ***
## y2t.l1 0.023040 0.006492 3.549 0.000462 ***
## y2t.l2 0.025591 0.006645 3.851 0.000149 ***
## const -0.019210 0.050686 -0.379 0.705017
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.8076 on 250 degrees of freedom
## Multiple R-Squared: 0.2211, Adjusted R-squared: 0.2086
## F-statistic: 17.74 on 4 and 250 DF, p-value: 7.843e-13
##
##
## Estimation results for equation y2t:
## ====================================
## y2t = const
##
## Estimate Std. Error t value Pr(>|t|)
## const 0.0175 0.4916 0.036 0.972
##
##
## Residual standard error: 7.835 on 253 degrees of freedom
## Multiple R-Squared: 5.01e-06, Adjusted R-squared: -0.003948
## F-statistic: 0.001267 on 1 and 253 DF, p-value: 0.9716
##
##
##
## Covariance matrix of residuals:
## y1t y2t
## y1t 0.6548 0.1318
## y2t 0.1318 62.3712
##
## Correlation matrix of residuals:
## y1t y2t
## y1t 1.00000 0.02062
## y2t 0.02062 1.00000
D
l2 <- var2$varresult
l.r <- var2.r$varresult
l.r2 <- var2.r2$varresult
stargazer(l2$y1t, l2$y2, l.r$y2t, l.r$y2t, l.r2$y1t, l.r2$y2t, type="html", column.labels = rep(colnames(ytyt), 3), dep.var.labels.include = F)
##
## <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>y1t</td><td>y2t</td><td>y1t</td><td>y2t</td><td>y1t</td><td>y2t</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">y1t.l1</td><td>0.270<sup>***</sup></td><td>0.081</td><td></td><td></td><td>0.311<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.061)</td><td>(0.593)</td><td></td><td></td><td>(0.057)</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">y2t.l1</td><td>0.024<sup>***</sup></td><td>0.069</td><td>0.075</td><td>0.075</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.063)</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">y1t.l2</td><td>0.109<sup>*</sup></td><td>-1.001<sup>*</sup></td><td></td><td></td><td></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.060)</td><td>(0.580)</td><td></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">y2t.l2</td><td>0.027<sup>***</sup></td><td>-0.082</td><td>-0.080</td><td>-0.080</td><td>0.026<sup>***</sup></td><td></td></tr>
## <tr><td style="text-align:left"></td><td>(0.007)</td><td>(0.064)</td><td>(0.063)</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.0002</td><td>0.019</td><td>0.019</td><td>-0.019</td><td>0.018</td></tr>
## <tr><td style="text-align:left"></td><td>(0.050)</td><td>(0.490)</td><td>(0.491)</td><td>(0.491)</td><td>(0.051)</td><td>(0.492)</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.024</td><td>0.011</td><td>0.011</td><td>0.221</td><td>0.00001</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.218</td><td>0.008</td><td>-0.001</td><td>-0.001</td><td>0.209</td><td>-0.004</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.804 (df = 249)</td><td>7.802 (df = 249)</td><td>7.822 (df = 251)</td><td>7.822 (df = 251)</td><td>0.808 (df = 250)</td><td>7.835 (df = 253)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>18.666<sup>***</sup> (df = 4; 249)</td><td>1.535 (df = 4; 249)</td><td>0.950 (df = 3; 251)</td><td>0.950 (df = 3; 251)</td><td>17.738<sup>***</sup> (df = 4; 250)</td><td>0.001 (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
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),]
## y1t y2t
## [1,] 1.0000000 0.0000000
## [2,] 0.8623911 0.1376089
## [3,] 0.8612114 0.1387886
## [4,] 0.8612085 0.1387915
var2.fevd[[2]][c(1,4,8,40),]
## y1t y2t
## [1,] 0.001132627 0.9988674
## [2,] 0.012758827 0.9872412
## [3,] 0.012942936 0.9870571
## [4,] 0.012943297 0.9870567
plot(var2.fevd)
plot(var2.fevd, addbars=8)
yy2<- cbind(y2t,y1t)
yy2 <- window(yy2,start="1951 Q1", end="2014 Q4")
VARselect (yy2, lag.max=8)
## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 1 2
##
## $criteria
## 1 2 3 4 5 6
## AIC(n) 3.751070 3.694427 3.706871 3.734600 3.752420 3.774346
## HQ(n) 3.785289 3.751458 3.786715 3.837256 3.877889 3.922627
## SC(n) 3.836072 3.836098 3.905210 3.989607 4.064095 4.142689
## FPE(n) 42.566703 40.222961 40.727412 41.873939 42.629080 43.577389
## 7 8
## AIC(n) 3.782360 3.810193
## HQ(n) 3.953453 4.004099
## SC(n) 4.207371 4.291873
## FPE(n) 43.932553 45.178626
var21 <- VAR(yy2,p=2,type="const")
var21.irfs <-irf(var21, n.ahead=30)
par(mfcol=c(2,2), cex=0.75)
plot(var21.irfs, plot.type="single")
var21.fevd <- fevd(var21, n.ahead=30)
var21.fevd[[1]][c(1,4,8,30),]
## y2t y1t
## [1,] 1.0000000 0.00000000
## [2,] 0.9889798 0.01102018
## [3,] 0.9888223 0.01117773
## [4,] 0.9888219 0.01117806
var21.fevd[[2]][c(1,4,8,30)]
## [1] 0.001132627 0.145494015 0.146667564 0.146670559
plot(var21.fevd)
plot(var21.fevd, addbars=6)