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)