Homework 5 Problem 1
Eva Nurwita
Direction
Obtain quarterly time series for real GDP FRED/GDPC1, GDP deflator FRED/GDPDEF and quarterly closing value of S&P 500 Index YAHOO/INDEX_GSPC, estimate the VAR model and identified each models.
(a) Estimate a bivariate reduced form VAR for \(Y_{t}= (y_{1,t},y_{2,t})'\)
First, we need to define all the variables that we are going to use in this problem.
rgdp <- Quandl("FRED/GDPC1", type="zoo")
defgdp <- Quandl("FRED/GDPDEF", type="zoo")
sp500 <- Quandl("YAHOO/INDEX_GSPC", type="zoo")
qsp500 <- aggregate(sp500[,4],as.yearqtr,tail,1)
y1t <- diff(log(rgdp))
ddefgdp <- diff(log(defgdp))
dsp500 <- diff(log(qsp500[,4]))
y2t <- dsp500-ddefgdpBefore we estimate the VAR model, we can plot the datasets to see the movement of the data and to give initial description for our work
par(mfrow=c(1,2), cex=0.8)
plot(y1t, xlab="", ylab="", main="Growth Rate of GDP", col="black")
plot(y2t, xlab="", ylab="", main="Inflation Adjusted Return", col="red", ylim=c(-0.14,0.14))
Now we estimate the var model for the period 1951 Quarter 1 to 2014 Quarter 4 and using the information criteria to determine the lags, we can see that the most suitable lags to choose from the AIC will be equal to two, thus we will estimate the reduced form VAR model at lag 2. From the VAR estimation, we can see hat correlation of residuals is 0.1349 which means that there is an indication to instantaneous causality, for further analysis, we need to use Granger causality.
y.Q <- cbind(y1t, y2t)
y.Q <- window(y.Q, start="1951 Q1", end="2014 Q4")
VARselect(y.Q, lag.max=8, type="const")## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 1 2
##
## $criteria
## 1 2 3 4
## AIC(n) -1.465126e+01 -1.470486e+01 -1.469229e+01 -1.466603e+01
## HQ(n) -1.461704e+01 -1.464783e+01 -1.461245e+01 -1.456337e+01
## SC(n) -1.456626e+01 -1.456319e+01 -1.449395e+01 -1.441102e+01
## FPE(n) 4.335513e-07 4.109277e-07 4.161332e-07 4.272205e-07
## 5 6 7 8
## AIC(n) -1.464883e+01 -1.462725e+01 -1.461845e+01 -1.459013e+01
## HQ(n) -1.452337e+01 -1.447897e+01 -1.444736e+01 -1.439622e+01
## SC(n) -1.433716e+01 -1.425891e+01 -1.419344e+01 -1.410845e+01
## FPE(n) 4.346529e-07 4.441705e-07 4.481426e-07 4.610779e-07var1 <- VAR(y.Q, p=2, type="const")
summary(var1)##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1t, y2t
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: 1154.493
## Roots of the characteristic polynomial:
## 0.4454 0.4454 0.304 0.304
## Call:
## VAR(y = y.Q, 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.2450402 0.0614771 3.986 8.84e-05 ***
## y2t.l1 0.0230968 0.0064046 3.606 0.000375 ***
## y1t.l2 0.0793493 0.0595434 1.333 0.183872
## y2t.l2 0.0263004 0.0065720 4.002 8.29e-05 ***
## const 0.0046278 0.0007233 6.398 7.75e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.008038 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 y2t:
## ====================================
## y2t = y1t.l1 + y2t.l1 + y1t.l2 + y2t.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1t.l1 0.256054 0.609111 0.420 0.6746
## y2t.l1 0.106501 0.063457 1.678 0.0945 .
## y1t.l2 -0.819353 0.589953 -1.389 0.1661
## y2t.l2 -0.040573 0.065115 -0.623 0.5338
## const 0.013836 0.007167 1.931 0.0547 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.07964 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:
## y1t y2t
## y1t 6.46e-05 0.000087
## y2t 8.70e-05 0.006342
##
## Correlation matrix of residuals:
## y1t y2t
## y1t 1.0000 0.1359
## y2t 0.1359 1.0000(b) Granger Causality Test
To further test about the instantaneous causality existence, Granger Causality Test for both variable is conducted below:
causality(var1, cause="y1t")## $Granger
##
## Granger causality H0: y1t do not Granger-cause y2t
##
## data: VAR object var1
## F-Test = 0.96654, df1 = 2, df2 = 498, p-value = 0.3811
##
##
## $Instant
##
## H0: No instantaneous causality between: y1t and y2t
##
## data: VAR object var1
## Chi-squared = 4.6072, df = 1, p-value = 0.03184causality(var1, cause="y2t")## $Granger
##
## Granger causality H0: y2t do not Granger-cause y1t
##
## data: VAR object var1
## F-Test = 15.762, df1 = 2, df2 = 498, p-value = 2.305e-07
##
##
## $Instant
##
## H0: No instantaneous causality between: y2t and y1t
##
## data: VAR object var1
## Chi-squared = 4.6072, df = 1, p-value = 0.03184The granger causality test showing that there is no instantaneous causality, y1t does not improve the forecasting performance of y2t and so does otherwise, we need to separate/restrict the VAR estimation.
(c) Restricted VAR Model
First we estimate the restricted model to eliminate the lags
mat.r <- matrix(1, nrow=2, ncol=5)
mat.r[1, c(2,4)] <- 0
varp.r <- restrict(var1, method="manual", resmat=mat.r)
summary(varp.r)##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1t, y2t
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: 1139.086
## Roots of the characteristic polynomial:
## 0.4867 0.2014 0.2014 0.1565
## Call:
## VAR(y = y.Q, p = 2, type = "const")
##
##
## Estimation results for equation y1t:
## ====================================
## y1t = y1t.l1 + y1t.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1t.l1 0.3301391 0.0628607 5.252 3.21e-07 ***
## y1t.l2 0.0761755 0.0628199 1.213 0.226
## const 0.0044907 0.0007616 5.896 1.19e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.008497 on 251 degrees of freedom
## Multiple R-Squared: 0.4906, Adjusted R-squared: 0.4845
## F-statistic: 80.59 on 3 and 251 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation y2t:
## ====================================
## y2t = y1t.l1 + y2t.l1 + y1t.l2 + y2t.l2 + const
##
## Estimate Std. Error t value Pr(>|t|)
## y1t.l1 0.256054 0.609111 0.420 0.6746
## y2t.l1 0.106501 0.063457 1.678 0.0945 .
## y1t.l2 -0.819353 0.589953 -1.389 0.1661
## y2t.l2 -0.040573 0.065115 -0.623 0.5338
## const 0.013836 0.007167 1.931 0.0547 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.07964 on 249 degrees of freedom
## Multiple R-Squared: 0.03656, Adjusted R-squared: 0.01721
## F-statistic: 1.89 on 5 and 249 DF, p-value: 0.09665
##
##
##
## Covariance matrix of residuals:
## y1t y2t
## y1t 7.278e-05 0.000087
## y2t 8.700e-05 0.006342
##
## Correlation matrix of residuals:
## y1t y2t
## y1t 1.0000 0.1281
## y2t 0.1281 1.0000varp.r$restrictions## y1t.l1 y2t.l1 y1t.l2 y2t.l2 const
## y1t 1 0 1 0 1
## y2t 1 1 1 1 1Acoef(varp.r)## [[1]]
## y1t.l1 y2t.l1
## y1t 0.3301391 0.0000000
## y2t 0.2560542 0.1065012
##
## [[2]]
## y1t.l2 y2t.l2
## y1t 0.07617549 0.00000000
## y2t -0.81935297 -0.04057338And then we estimate the restricted VAR model removing all variables which have t-statistics lower than 2
varp.r.ser <- restrict(var1, method="ser", thresh=2.0)
summary(varp.r.ser)##
## VAR Estimation Results:
## =========================
## Endogenous variables: y1t, y2t
## Deterministic variables: const
## Sample size: 254
## Log Likelihood: 1150.599
## Roots of the characteristic polynomial:
## 0.2734 0 0 0
## Call:
## VAR(y = y.Q, 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.2734178 0.0577607 4.734 3.70e-06 ***
## y2t.l1 0.0226323 0.0064051 3.533 0.000488 ***
## y2t.l2 0.0266285 0.0065775 4.048 6.88e-05 ***
## const 0.0050198 0.0006618 7.585 6.55e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.00805 on 250 degrees of freedom
## Multiple R-Squared: 0.5446, Adjusted R-squared: 0.5374
## F-statistic: 74.76 on 4 and 250 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation y2t:
## ====================================
## y2t = const
##
## Estimate Std. Error t value Pr(>|t|)
## const 0.01017 0.00501 2.03 0.0434 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.07984 on 253 degrees of freedom
## Multiple R-Squared: 0.01602, Adjusted R-squared: 0.01213
## F-statistic: 4.119 on 1 and 253 DF, p-value: 0.04344
##
##
##
## Covariance matrix of residuals:
## y1t y2t
## y1t 6.506e-05 8.224e-05
## y2t 8.224e-05 6.477e-03
##
## Correlation matrix of residuals:
## y1t y2t
## y1t 1.0000 0.1267
## y2t 0.1267 1.0000varp.r.ser$restrictions## y1t.l1 y2t.l1 y1t.l2 y2t.l2 const
## y1t 1 1 0 1 1
## y2t 0 0 0 0 1Acoef(varp.r.ser)## [[1]]
## y1t.l1 y2t.l1
## y1t 0.2734178 0.02263225
## y2t 0.0000000 0.00000000
##
## [[2]]
## y1t.l2 y2t.l2
## y1t 0 0.02662845
## y2t 0 0.00000000(d) Using Stargazer Package to Show All Three Estimation Result
First, loading the source:
setwd("C:/Users/evanu/OneDrive/Time Series/Homework 5")
source("stargazerX.r")
source("stargazerX-internal.r")And then We compile our stargazer table using the following command:
lm1 <- var1$varresult
lmp <- varp.r$varresult
lms <- varp.r.ser$varresult
stargazerX(lm1$y1t, lm1$y2t, lmp$y1t, lmp$y2t, type="text", title="",
align=TRUE, header=FALSE, model.numbers=TRUE,
dep.var.caption="", result="asis", column.labels=c("VAR1-Y1t","VAR1-Y2t","VARP-Y1t", "VARP-Y2t"))##
## ============================================================================================================
## y
## VAR1-Y1t VAR1-Y2t VARP-Y1t VARP-Y2t
## (1) (2) (3) (4)
## ------------------------------------------------------------------------------------------------------------
## y1t.l1 0.245*** 0.256 0.330*** 0.256
## (0.061) (0.609) (0.063) (0.609)
##
## y2t.l1 0.023*** 0.107* 0.107*
## (0.006) (0.063) (0.063)
##
## y1t.l2 0.079 -0.819 0.076 -0.819
## (0.060) (0.590) (0.063) (0.590)
##
## y2t.l2 0.026*** -0.041 -0.041
## (0.007) (0.065) (0.065)
##
## const 0.005*** 0.014* 0.004*** 0.014*
## (0.001) (0.007) (0.001) (0.007)
##
## ------------------------------------------------------------------------------------------------------------
## Observations 254 254 254 254
## R2 0.231 0.021 0.491 0.037
## Adjusted R2 0.219 0.005 0.485 0.017
## Residual Std. Error 0.008 (df = 249) 0.080 (df = 249) 0.008 (df = 251) 0.080 (df = 249)
## F Statistic 18.685*** (df = 4; 249) 1.327 (df = 4; 249) 80.588*** (df = 3; 251) 1.890* (df = 5; 249)
## ============================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
##
## ====
## asis
## ----stargazerX(lms$y1t, lms$y2t, type="text", title="",
align=TRUE, header=FALSE, model.numbers=TRUE,
dep.var.caption="", result="asis", column.labels=c("VARP.SER-Y1t","VARP.SER-Y2t"))##
## =================================================================
## y
## VARP.SER-Y1t VARP.SER-Y2t
## (1) (2)
## -----------------------------------------------------------------
## y1t.l1 0.273***
## (0.058)
##
## y2t.l1 0.023***
## (0.006)
##
## y2t.l2 0.027***
## (0.007)
##
## const 0.005*** 0.010**
## (0.001) (0.005)
##
## -----------------------------------------------------------------
## Observations 254 254
## R2 0.545 0.016
## Adjusted R2 0.537 0.012
## Residual Std. Error 0.008 (df = 250) 0.080 (df = 253)
## F Statistic 74.756*** (df = 4; 250) 4.119** (df = 1; 253)
## =================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
##
## ====
## asis
## ----(e) Impulse Response Function
First we have to plot the IRF and FEVD of the VAR model based on Choleski Decomposition
varp.irfs <- irf(var1, n.ahead=10)
par(mfcol=c(2,2), cex=0.6)
plot(varp.irfs, plot.type="single")
And the FEVD
varp.fevd <- fevd(var1, n.ahead=40)
varp.fevd[[1]][c(1,4,8,40),]## y1t y2t
## [1,] 1.0000000 0.0000000
## [2,] 0.8607464 0.1392536
## [3,] 0.8591922 0.1408078
## [4,] 0.8591882 0.1408118varp.fevd[[2]][c(1,4,8,40),]## y1t y2t
## [1,] 0.01847382 0.9815262
## [2,] 0.02623100 0.9737690
## [3,] 0.02642786 0.9735721
## [4,] 0.02642813 0.9735719plot(varp.fevd)
Now if we reversed the order of the variable,
y.Qord2 <- cbind(y2t, y1t)
y.Qord2 <- window(y.Qord2, start="1951 Q1", end="2014 Q4")
Var1.ord2 <- VAR(y.Qord2, p=2, type="const")
varp.irfs <- irf(Var1.ord2, n.ahead=10)
par(mfcol=c(2,2), cex=0.6)
plot(varp.irfs, plot.type="single")
varp.fevd <- fevd(Var1.ord2, n.ahead=40)
varp.fevd[[1]][c(1,4,8,40),]## y2t y1t
## [1,] 1.0000000 0.000000000
## [2,] 0.9935239 0.006476147
## [3,] 0.9934082 0.006591832
## [4,] 0.9934079 0.006592058varp.fevd[[2]][c(1,4,8,40),]## y2t y1t
## [1,] 0.01847382 0.9815262
## [2,] 0.18059627 0.8194037
## [3,] 0.18207471 0.8179253
## [4,] 0.18207906 0.8179209plot(varp.fevd)
There is change in the impulse response function and also the FEVD, which tells us that order does matter in estimating VAR model