Problem 1: Estimate VAR Model

(a)

library(tidyquant)
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
## Loading required package: PerformanceAnalytics
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
## Loading required package: quantmod
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
## Loading required package: tidyverse
## Warning: package 'tidyverse' was built under R version 3.4.4
## -- Attaching packages ----------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 2.2.1     v purrr   0.2.4
## v tibble  1.4.2     v dplyr   0.7.4
## v tidyr   0.8.0     v stringr 1.3.0
## v readr   1.1.1     v forcats 0.3.0
## -- Conflicts -------------------------------------------------------------- tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date()        masks base::date()
## x dplyr::filter()          masks stats::filter()
## x dplyr::first()           masks xts::first()
## x lubridate::intersect()   masks base::intersect()
## x dplyr::lag()             masks stats::lag()
## x dplyr::last()            masks xts::last()
## x lubridate::setdiff()     masks base::setdiff()
## x lubridate::union()       masks base::union()
library(ggplot2)
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 3.4.4
theme_set(theme_bw())

# obtain data on US Real GDP and GDP deflator and S&P500 Index


data1 <-
  tq_get("GDPC1", get = "economic.data",
         from = "1950-01-01", to = "2017-12-31")  %>%
  rename(GDP = price) %>%
  mutate(lGDP = log(GDP),
         dlGDP = lGDP - lag(lGDP),
         dlrGDP = 400*(dlGDP))


data2 <-
  tq_get("GDPDEF", get = "economic.data",
         from = "1950-01-01", to = "2017-12-31")  %>%
  rename(GDPDEF = price) %>%
  mutate(lGDPDEF = log(GDPDEF),
         dlGDPDEF = lGDPDEF - lag(lGDPDEF))


data3 <-
  tq_get("^GSPC", get = "stock.prices",
         from = "1950-01-01", to = "2017-12-31")  %>%
  select(date, adjusted) %>%
  mutate(qtryear = as.yearqtr(date)) %>%
  group_by(qtryear) %>%
  summarise(SP500 = mean(adjusted)) %>%
  ungroup()

(b)

dlrGDP <- ts(data1$dlrGDP)
plot(dlrGDP)

lSP500 <- log(data3$SP500)
dlSP500 <- ts(lSP500 - lag(lSP500))
plot(dlSP500)

dlGDPDEF <- ts(data2$dlGDPDEF)
plot(dlGDPDEF)

dlrSP500 <- 100*(dlSP500)-dlGDPDEF
plot(dlrSP500)

(c)

 # Estimate a bivariate reduced form VAR for the period 1990Q1-2017Q4

# Series for the period 1990Q1-2017Q4

dlrGDP1 <- ts(data1$dlrGDP[data1$date > "1990-01-01"], start=c(1990,1), frequency = 4)
plot(dlrGDP1)

dlGDPDEF1 <- ts(data2$dlGDPDEF[data2$date > "1990-01-01"], start=c(1990,1), frequency = 4)
plot(dlGDPDEF1)

dlSP5001 <- ts((lSP500 - lag(lSP500))[data3$qtryear > "1990 Q1"], start=c(1990,1), frequency = 4)
plot(dlSP5001)

dlrSP5001 <- 100*(dlSP5001)- dlGDPDEF1
plot(dlrSP5001)

library(vars)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: strucchange
## Loading required package: sandwich
## 
## Attaching package: 'strucchange'
## The following object is masked from 'package:stringr':
## 
##     boundary
## Loading required package: urca
## Loading required package: lmtest
library(timetk)

# VAR Model

v1 <- cbind(dlrSP5001, dlrGDP1)
VARselect(v1, lag.max = 8, type = "const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2 
## 
## $criteria
##                 1          2          3          4          5          6
## AIC(n)   4.813560   4.769170   4.831048   4.888774   4.958325   5.000042
## HQ(n)    4.875725   4.872777   4.976099   5.075267   5.186261   5.269421
## SC(n)    4.967040   5.024968   5.189167   5.349212   5.521083   5.665119
## FPE(n) 123.173438 117.839349 125.394905 132.909324 142.588142 148.821625
##                 7          8
## AIC(n)   5.059363   5.107884
## HQ(n)    5.370185   5.460149
## SC(n)    5.826760   5.977600
## FPE(n) 158.148462 166.330815
# VAR estimation with lag 2

var1 <- VAR(v1, p = 2, type = "const")
var1
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation dlrSP5001: 
## ============================================== 
## Call:
## dlrSP5001 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
## dlrSP5001.l1   dlrGDP1.l1 dlrSP5001.l2   dlrGDP1.l2        const 
##    0.4121764    0.2075374   -0.1081056   -0.1768264    1.2558179 
## 
## 
## Estimated coefficients for equation dlrGDP1: 
## ============================================ 
## Call:
## dlrGDP1 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
## dlrSP5001.l1   dlrGDP1.l1 dlrSP5001.l2   dlrGDP1.l2        const 
##   0.11705801   0.18302514   0.01586628   0.17558029   1.31433474
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dlrSP5001, dlrGDP1 
## Deterministic variables: const 
## Sample size: 109 
## Log Likelihood: -561.895 
## Roots of the characteristic polynomial:
## 0.4685 0.3568 0.3111 0.3111
## Call:
## VAR(y = v1, p = 2, type = "const")
## 
## 
## Estimation results for equation dlrSP5001: 
## ========================================== 
## dlrSP5001 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## dlrSP5001.l1   0.4122     0.1054   3.909 0.000165 ***
## dlrGDP1.l1     0.2075     0.2818   0.737 0.463039    
## dlrSP5001.l2  -0.1081     0.1099  -0.984 0.327546    
## dlrGDP1.l2    -0.1768     0.2684  -0.659 0.511416    
## const          1.2558     0.8619   1.457 0.148101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.648 on 104 degrees of freedom
## Multiple R-Squared: 0.1757,  Adjusted R-squared: 0.144 
## F-statistic: 5.541 on 4 and 104 DF,  p-value: 0.0004397 
## 
## 
## Estimation results for equation dlrGDP1: 
## ======================================== 
## dlrGDP1 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## dlrSP5001.l1  0.11706    0.03903   2.999  0.00339 ** 
## dlrGDP1.l1    0.18303    0.10430   1.755  0.08223 .  
## dlrSP5001.l2  0.01587    0.04068   0.390  0.69731    
## dlrGDP1.l2    0.17558    0.09934   1.768  0.08008 .  
## const         1.31433    0.31903   4.120 7.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 2.091 on 104 degrees of freedom
## Multiple R-Squared: 0.2806,  Adjusted R-squared: 0.253 
## F-statistic: 10.14 on 4 and 104 DF,  p-value: 5.676e-07 
## 
## 
## 
## Covariance matrix of residuals:
##           dlrSP5001 dlrGDP1
## dlrSP5001    31.895   5.128
## dlrGDP1       5.128   4.370
## 
## Correlation matrix of residuals:
##           dlrSP5001 dlrGDP1
## dlrSP5001    1.0000  0.4344
## dlrGDP1      0.4344  1.0000
plot(var1, nc = 4, lag.acf = 16, lag.pacf = 16)

# VAR estimation Using AIC information criteria 

varp1 <- VAR(v1, ic = "AIC", lag.max = 8, type = "const")
varp1
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation dlrSP5001: 
## ============================================== 
## Call:
## dlrSP5001 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
## dlrSP5001.l1   dlrGDP1.l1 dlrSP5001.l2   dlrGDP1.l2        const 
##    0.4121764    0.2075374   -0.1081056   -0.1768264    1.2558179 
## 
## 
## Estimated coefficients for equation dlrGDP1: 
## ============================================ 
## Call:
## dlrGDP1 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
## dlrSP5001.l1   dlrGDP1.l1 dlrSP5001.l2   dlrGDP1.l2        const 
##   0.11705801   0.18302514   0.01586628   0.17558029   1.31433474
summary(varp1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dlrSP5001, dlrGDP1 
## Deterministic variables: const 
## Sample size: 109 
## Log Likelihood: -561.895 
## Roots of the characteristic polynomial:
## 0.4685 0.3568 0.3111 0.3111
## Call:
## VAR(y = v1, type = "const", lag.max = 8, ic = "AIC")
## 
## 
## Estimation results for equation dlrSP5001: 
## ========================================== 
## dlrSP5001 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## dlrSP5001.l1   0.4122     0.1054   3.909 0.000165 ***
## dlrGDP1.l1     0.2075     0.2818   0.737 0.463039    
## dlrSP5001.l2  -0.1081     0.1099  -0.984 0.327546    
## dlrGDP1.l2    -0.1768     0.2684  -0.659 0.511416    
## const          1.2558     0.8619   1.457 0.148101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.648 on 104 degrees of freedom
## Multiple R-Squared: 0.1757,  Adjusted R-squared: 0.144 
## F-statistic: 5.541 on 4 and 104 DF,  p-value: 0.0004397 
## 
## 
## Estimation results for equation dlrGDP1: 
## ======================================== 
## dlrGDP1 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## dlrSP5001.l1  0.11706    0.03903   2.999  0.00339 ** 
## dlrGDP1.l1    0.18303    0.10430   1.755  0.08223 .  
## dlrSP5001.l2  0.01587    0.04068   0.390  0.69731    
## dlrGDP1.l2    0.17558    0.09934   1.768  0.08008 .  
## const         1.31433    0.31903   4.120 7.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 2.091 on 104 degrees of freedom
## Multiple R-Squared: 0.2806,  Adjusted R-squared: 0.253 
## F-statistic: 10.14 on 4 and 104 DF,  p-value: 5.676e-07 
## 
## 
## 
## Covariance matrix of residuals:
##           dlrSP5001 dlrGDP1
## dlrSP5001    31.895   5.128
## dlrGDP1       5.128   4.370
## 
## Correlation matrix of residuals:
##           dlrSP5001 dlrGDP1
## dlrSP5001    1.0000  0.4344
## dlrGDP1      0.4344  1.0000
plot(varp1, nc = 4, lag.acf = 16, lag.pacf = 16)

(d)

 ## Run the Granger causality tests for both variables (for var1)

causality(var1, cause = "dlrSP5001")
## $Granger
## 
##  Granger causality H0: dlrSP5001 do not Granger-cause dlrGDP1
## 
## data:  VAR object var1
## F-Test = 5.1594, df1 = 2, df2 = 208, p-value = 0.006503
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: dlrSP5001 and dlrGDP1
## 
## data:  VAR object var1
## Chi-squared = 17.302, df = 1, p-value = 3.189e-05
causality(var1, cause = "dlrGDP1")
## $Granger
## 
##  Granger causality H0: dlrGDP1 do not Granger-cause dlrSP5001
## 
## data:  VAR object var1
## F-Test = 0.38235, df1 = 2, df2 = 208, p-value = 0.6827
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: dlrGDP1 and dlrSP5001
## 
## data:  VAR object var1
## Chi-squared = 17.302, df = 1, p-value = 3.189e-05
# Granger causality H0: dlrSP5001 do not Granger-cause dlrGDP1: p-value = 0.006503, then reject Null Hypothesis
# Granger causality H0: dlrGDP1 do not Granger-cause dlrSP5001: p-value = 0.6827, then can not reject Null Hypothesis

## Run the Granger causality tests for both variables (for varp)  

causality(varp1, cause = "dlrSP5001")
## $Granger
## 
##  Granger causality H0: dlrSP5001 do not Granger-cause dlrGDP1
## 
## data:  VAR object varp1
## F-Test = 5.1594, df1 = 2, df2 = 208, p-value = 0.006503
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: dlrSP5001 and dlrGDP1
## 
## data:  VAR object varp1
## Chi-squared = 17.302, df = 1, p-value = 3.189e-05
causality(varp1, cause = "dlrGDP1")
## $Granger
## 
##  Granger causality H0: dlrGDP1 do not Granger-cause dlrSP5001
## 
## data:  VAR object varp1
## F-Test = 0.38235, df1 = 2, df2 = 208, p-value = 0.6827
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: dlrGDP1 and dlrSP5001
## 
## data:  VAR object varp1
## Chi-squared = 17.302, df = 1, p-value = 3.189e-05
# Granger causality H0: dlrSP5001 do not Granger-cause dlrGDP1: p-value = 0.006503, then reject Null Hypothesis
# Granger causality H0: dlrGDP1 do not Granger-cause dlrSP5001: p-value = 0.6827, then can not reject Null Hypothesis

(e)

# estimate a restricted VAR model in which you remove lags based on Granger causality test from (d)

mat.r <- matrix(1, nrow = 2, ncol = 5)
mat.r[1, c(2,4,5)] <- 0

# estimate a restricted VAR

varp1.r <- restrict(varp1, method = "manual", resmat = mat.r)
varp1.r
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation dlrSP5001: 
## ============================================== 
## Call:
## dlrSP5001 = dlrSP5001.l1 + dlrSP5001.l2 
## 
## dlrSP5001.l1 dlrSP5001.l2 
##   0.48221294  -0.07864958 
## 
## 
## Estimated coefficients for equation dlrGDP1: 
## ============================================ 
## Call:
## dlrGDP1 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
## dlrSP5001.l1   dlrGDP1.l1 dlrSP5001.l2   dlrGDP1.l2        const 
##   0.11705801   0.18302514   0.01586628   0.17558029   1.31433474
summary(varp1.r)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dlrSP5001, dlrGDP1 
## Deterministic variables: const 
## Sample size: 109 
## Log Likelihood: -565.563 
## Roots of the characteristic polynomial:
## 0.5204 0.3374 0.2804 0.2804
## Call:
## VAR(y = v1, type = "const", lag.max = 8, ic = "AIC")
## 
## 
## Estimation results for equation dlrSP5001: 
## ========================================== 
## dlrSP5001 = dlrSP5001.l1 + dlrSP5001.l2 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## dlrSP5001.l1  0.48221    0.09615   5.015 2.11e-06 ***
## dlrSP5001.l2 -0.07865    0.09607  -0.819    0.415    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.723 on 107 degrees of freedom
## Multiple R-Squared: 0.2047,  Adjusted R-squared: 0.1899 
## F-statistic: 13.77 on 2 and 107 DF,  p-value: 4.756e-06 
## 
## 
## Estimation results for equation dlrGDP1: 
## ======================================== 
## dlrGDP1 = dlrSP5001.l1 + dlrGDP1.l1 + dlrSP5001.l2 + dlrGDP1.l2 + const 
## 
##              Estimate Std. Error t value Pr(>|t|)    
## dlrSP5001.l1  0.11706    0.03903   2.999  0.00339 ** 
## dlrGDP1.l1    0.18303    0.10430   1.755  0.08223 .  
## dlrSP5001.l2  0.01587    0.04068   0.390  0.69731    
## dlrGDP1.l2    0.17558    0.09934   1.768  0.08008 .  
## const         1.31433    0.31903   4.120 7.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 2.091 on 104 degrees of freedom
## Multiple R-Squared: 0.6395,  Adjusted R-squared: 0.6221 
## F-statistic: 36.89 on 5 and 104 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           dlrSP5001 dlrGDP1
## dlrSP5001    32.303   5.128
## dlrGDP1       5.128   4.370
## 
## Correlation matrix of residuals:
##           dlrSP5001 dlrGDP1
## dlrSP5001    1.0000  0.4316
## dlrGDP1      0.4316  1.0000

(f)

# Use the VAR model to create a forecast for 2018Q1-2018Q4

# construct 1 to 4 quarter ahead forecast and its 90% confidence interval

library(forecast)
## Warning: package 'forecast' was built under R version 3.4.4
varp1.f <- predict(varp1, n.ahead = 4, ci = 0.9)
varp1.f
## $dlrSP5001
##          fcst     lower    upper        CI
## [1,] 3.201157 -6.088267 12.49058  9.289425
## [2,] 2.125751 -8.064282 12.31579 10.190034
## [3,] 1.833878 -8.382894 12.05065 10.216773
## [4,] 1.834796 -8.383515 12.05311 10.218311
## 
## $dlrGDP1
##          fcst      lower    upper       CI
## [1,] 3.053921 -0.3846471 6.492490 3.438568
## [2,] 2.832996 -0.9082330 6.574226 3.741229
## [3,] 2.668679 -1.3095256 6.646884 3.978205
## [4,] 2.548586 -1.4804880 6.577661 4.029074
summary(varp1.f)
##          Length Class  Mode   
## fcst       2    -none- list   
## endog    222    mts    numeric
## model     10    varest list   
## exo.fcst   0    -none- NULL
# fanchart - by default the step is 0.1

fanchart(varp1.f, lwd = 2)

autoplot(varp1.f) + geom_hline(yintercept = 0, linetype = "dashed") + theme_bw()

# Real GDP growth rate in 2018Q1 based on  
# (1) the Federal Bank of New York Nowcast = 2.91% 
# (2) the GDPNow Federal Bank of Atlanta forecast = 2.0%
# (3) the minimum, the average, and the maximum forecasts in the Wall Street Journal Economic Forecasting Survey 2.1%

# Here, real GDP growth rate in 2018Q1 = 3.05% which is higher than the other forecasts in (1), (2), and (3).

# Impulse-Response Functions (IRF) 

varp1.irfs <- irf(varp1, n.ahead = 40)
par(mfcol = c(2,2), cex = 0.8, mar = c(3,4,2,2))
plot(varp1.irfs, plot.type = "single", lwd = 2)

Problem 2: Estimate SVAR Model

(a)

# obtain data on labor productivity, measured as Nonfarm Business Sector: Real Output Per Hour of All Persons
# OPHNFB and for total hours worked, measured as Nonfarm Business Sector: Hours of All Persons HOANBS.

data4 <-
  tq_get("OPHNFB", get = "economic.data",
         from = "1947-01-01", to = "2017-12-31")  %>%
  rename(y1 = price) %>%
  mutate(ly1 = log(y1),
         dly1 = ly1 - lag(ly1)) 


data5 <-
  tq_get("HOANBS", get = "economic.data",
         from = "1947-01-01", to = "2017-12-31")  %>%
  rename(y2 = price) %>%
  mutate(ly2 = log(y2),
         dly2 = ly2 - lag(ly2)) 

(b)

# Test for the presence of unit root using ERS test.

# urersTest  Elliott-Rothenberg-Stock test for unit roots,  

library(urca)

ur.ers(data4$ly1)  # The value of the test statistic is: 3.5051 
## 
## ############################################################### 
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 3.5051
ur.ers(data4$dly1) # The value of the test statistic is: -1.8145
## 
## ############################################################### 
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -1.8145
ur.ers(data5$ly2) # The value of the test statistic is: 2.1617
## 
## ############################################################### 
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: 2.1617
ur.ers(data5$dly2) # The value of the test statistic is: -6.5027
## 
## ############################################################### 
## # Elliot, Rothenberg and Stock Unit Root / Cointegration Test # 
## ############################################################### 
## 
## The value of the test statistic is: -6.5027
# Series are stationary, then we can use them for estimating svar model.

(c)

# Estimate a bivariate reduced form VAR using AIC information criteria 

v2 <- ts(na.omit(cbind(data4$dly1, data5$dly2)))
varp2 <- VAR(v2, ic = "AIC", lag.max = 8, type = "const")
varp2
## 
## VAR Estimation Results:
## ======================= 
## 
## Estimated coefficients for equation Series.1: 
## ============================================= 
## Call:
## Series.1 = Series.1.l1 + Series.2.l1 + Series.1.l2 + Series.2.l2 + Series.1.l3 + Series.2.l3 + const 
## 
##  Series.1.l1  Series.2.l1  Series.1.l2  Series.2.l2  Series.1.l3 
## -0.061088229  0.063011021  0.051721762 -0.195943680  0.007129394 
##  Series.2.l3        const 
## -0.192711574  0.006194392 
## 
## 
## Estimated coefficients for equation Series.2: 
## ============================================= 
## Call:
## Series.2 = Series.1.l1 + Series.2.l1 + Series.1.l2 + Series.2.l2 + Series.1.l3 + Series.2.l3 + const 
## 
##   Series.1.l1   Series.2.l1   Series.1.l2   Series.2.l2   Series.1.l3 
##  0.1069386681  0.6265493767  0.1034210595 -0.0176821607  0.0874578733 
##   Series.2.l3         const 
## -0.0445960211 -0.0002338489
summary(varp2)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Series.1, Series.2 
## Deterministic variables: const 
## Sample size: 280 
## Log Likelihood: 1970.349 
## Roots of the characteristic polynomial:
## 0.7381 0.7381 0.4285 0.4285 0.4066 0.4066
## Call:
## VAR(y = v2, type = "const", lag.max = 8, ic = "AIC")
## 
## 
## Estimation results for equation Series.1: 
## ========================================= 
## Series.1 = Series.1.l1 + Series.2.l1 + Series.1.l2 + Series.2.l2 + Series.1.l3 + Series.2.l3 + const 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## Series.1.l1 -0.0610882  0.0588022  -1.039  0.29978    
## Series.2.l1  0.0630110  0.0709282   0.888  0.37512    
## Series.1.l2  0.0517218  0.0555934   0.930  0.35301    
## Series.2.l2 -0.1959437  0.0836382  -2.343  0.01986 *  
## Series.1.l3  0.0071294  0.0546710   0.130  0.89634    
## Series.2.l3 -0.1927116  0.0727477  -2.649  0.00854 ** 
## const        0.0061944  0.0007383   8.390 2.64e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.007891 on 273 degrees of freedom
## Multiple R-Squared: 0.1169,  Adjusted R-squared: 0.09753 
## F-statistic: 6.025 on 6 and 273 DF,  p-value: 6.166e-06 
## 
## 
## Estimation results for equation Series.2: 
## ========================================= 
## Series.2 = Series.1.l1 + Series.2.l1 + Series.1.l2 + Series.2.l2 + Series.1.l3 + Series.2.l3 + const 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## Series.1.l1  0.1069387  0.0502517   2.128   0.0342 *  
## Series.2.l1  0.6265494  0.0606145  10.337   <2e-16 ***
## Series.1.l2  0.1034211  0.0475095   2.177   0.0303 *  
## Series.2.l2 -0.0176822  0.0714764  -0.247   0.8048    
## Series.1.l3  0.0874579  0.0467213   1.872   0.0623 .  
## Series.2.l3 -0.0445960  0.0621694  -0.717   0.4738    
## const       -0.0002338  0.0006309  -0.371   0.7112    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.006743 on 273 degrees of freedom
## Multiple R-Squared: 0.4384,  Adjusted R-squared: 0.426 
## F-statistic: 35.52 on 6 and 273 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           Series.1  Series.2
## Series.1 6.226e-05 6.784e-06
## Series.2 6.784e-06 4.547e-05
## 
## Correlation matrix of residuals:
##          Series.1 Series.2
## Series.1   1.0000   0.1275
## Series.2   0.1275   1.0000

(d)

# Analyze effects of two types of shocks - technology shocks and demand shocks
# Use Blanchard and Quah approach to obtain an SVAR model

# Blanchard-Quah long run restriction: row 1 column 2 element of the cumulative effect matrix is 0

svar <- BQ(varp2)
summary(svar)
## 
## SVAR Estimation Results:
## ======================== 
## 
## Call:
## BQ(x = varp2)
## 
## Type: Blanchard-Quah 
## Sample size: 280 
## Log Likelihood: 1963.26 
## 
## Estimated contemporaneous impact matrix:
##           Series.1 Series.2
## Series.1  0.006496 0.004480
## Series.2 -0.003089 0.005994
## 
## Estimated identified long run impact matrix:
##           Series.1 Series.2
## Series.1  0.007188  0.00000
## Series.2 -0.002176  0.01376
## 
## Covariance matrix of reduced form residuals (*100):
##           Series.1  Series.2
## Series.1 0.0062260 0.0006784
## Series.2 0.0006784 0.0045470

(e)

# Report and interpret the contemporaneous impact and the long run impact matrices for the SVAR.

# In Estimated contemporaneous imapct matrix, the rows refer to two variables (dly1, dly2) and the columns to the  two shocks - technology shock and non-technology (Demand)  shock.  

# a. A positive one standard deviation technology shock increases and decreases dly1 and dly2 by 0.006% and 0.003% points, respectively. 

# b. A positive one standard deviation non-technology shock increases dly1 and dly2 by 0.004%, and 0.005% points, respectively. 

# Estimated Identified Long run Impact Matrix

# a. The long run effect of one standard deviation technology shock on dly1 is 0.007% points and single deviation will decrease dly2 by 0.002% points. 

# b. The long run effect of one standard deviation non-technological shock on dly1 will be 0.000% point and single deviation will increase dly2 by 0.014% point.

(f)

#  Plot the cumulative IRFs based on the SVAR model from (d) 

# standard non-cumulative IRFs

myIRF <- irf(svar, n.ahead=40, ci=.9)

# cumulative myIRFs

myIRF.c <- irf(svar, n.ahead=40, ci=.9, cumulative=TRUE)
plot( myIRF.c, plot.type="multiple")

# SVAR IRF from series 1 (cumulative) shows that initially series 1 increase is 0.006 due to one standard deviation shock in technology, decrease until the first quarter, and then starts increasing and eventually phases out completely. 

# SVAR IRF from series 1 (cumulative) shows that initially series 2 decrease is 0.003 due to one standard deviation shock in technology, decrease again until the second quarter, and then starts increasing and eventually phases out completely.  

# SVAR IRF from series 2 (cumulative) shows that initially series 1 increase is 0.0045 due to one standard deviation shock in demand, decrease until the seventh quarter, and then starts increasing and eventually phases out completely. 

# SVAR IRF from series 2 (cumulative) shows that initially series 2 increase is 0.005 due to one standard deviation shock in demand, increase again until the fourth quarter, and then starts decreasing and eventually phases out completely.  

(g)

# Compare your IRFs with Figure 2 from Gali (1999) AER.

library(imager)
## Warning: package 'imager' was built under R version 3.4.4
## Loading required package: plyr
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
## The following object is masked from 'package:lubridate':
## 
##     here
## Loading required package: magrittr
## Warning: package 'magrittr' was built under R version 3.4.4
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
## 
##     set_names
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Attaching package: 'imager'
## The following object is masked from 'package:magrittr':
## 
##     add
## The following object is masked from 'package:plyr':
## 
##     liply
## The following object is masked from 'package:vars':
## 
##     B
## The following object is masked from 'package:strucchange':
## 
##     boundary
## The following object is masked from 'package:stringr':
## 
##     boundary
## The following object is masked from 'package:tidyr':
## 
##     fill
## The following objects are masked from 'package:stats':
## 
##     convolve, spectrum
## The following object is masked from 'package:graphics':
## 
##     frame
## The following object is masked from 'package:base':
## 
##     save.image
library(png)

im <- load.image("C:/Users/Sima/Pictures/master_img-001.png")
plot(im)

# The results are quite different than in the Gali figure.

(h)

# Construct the FEVD for the SVAR model from (d).

library(forecast)

irfvc <- irf(svar, n.ahead = 40, cumulative=TRUE)
summary(irfvc)
##            Length Class  Mode     
## irf        2      -none- list     
## Lower      2      -none- list     
## Upper      2      -none- list     
## response   2      -none- character
## impulse    2      -none- character
## ortho      1      -none- logical  
## cumulative 1      -none- logical  
## runs       1      -none- numeric  
## ci         1      -none- numeric  
## boot       1      -none- logical  
## model      1      -none- character
vth <- fevd(svar, n.ahead=60)
summary(vth)
##          Length Class  Mode   
## Series.1 120    -none- numeric
## Series.2 120    -none- numeric
plot(vth)

# There are totally differences in the figures. By reversing the path of variables, the path of description forecast error variance of each variable by exogenous shocks to other variable is affected.