Homework 5
Problem 2
The first step is to import Industrial Production Index, FRED/INDPRO, and Consumer Price Index, FRED/CPIAUCSL, using Quandl packages. We will denote FRED/INDPRO as \(IPI_{t}\) and FRED/CPIAUCSL* as as \(CPI_{t}\).
IPI <- Quandl("FRED/INDPRO", type="ts")
CPI <- Quandl("FRED/CPIAUCSL", type="ts")Next step is to construct the log of Industrial Production Index, \(y_{1,t}=log{IPI_{t}}\), and Consumer Price Index, \(y_{2,t}=log{CPI_{t}}\) as follow:
y1t <- log(IPI)
y2t <- log(CPI)The following figures shows the original and log time series of \(IND_{t}\) and \(CPI_{t}\).
par(mfrow=c(2,2), cex=0.5, mar=c(2,2,3,1))
plot(IPI, xlab="", ylab="", main="Industrial Production Index")
plot(y1t, xlab="", ylab="", main="Log Industrial Production Index")
plot(CPI, xlab="", ylab="", main="Consumer Price Index")
plot(y2t, xlab="", ylab="", main="Log Consumer Price Index")
(a) Unit Root Test
The next step is to test for the stationary for both \(y_{1,t}\) and \(y_{2,t}\).
y1t.adf <- adf.test(y1t)
y1t.adf
Augmented Dickey-Fuller Test
data: y1t
Dickey-Fuller = -2.7713, Lag order = 10, p-value = 0.2518
alternative hypothesis: stationaryy2t.adf <- adf.test(y2t)
y2t.adf
Augmented Dickey-Fuller Test
data: y2t
Dickey-Fuller = -1.2868, Lag order = 9, p-value = 0.8802
alternative hypothesis: stationaryFrom the ADF test we can see that the p-value of the test for both \(y_{1,t}\) and \(y_{2,t}\) are greater than 0.05, hence, we can conclude that both \(y_{1,t}\) and \(y_{2,t}\) are non-stationary.
The next step that we have to take is to take first difference of the series and retest the series for non-stationarity. We construct new series \(\Delta y_{1,t}\) and \(\Delta y_{2,t}\) using the following code
dy1t <- diff(y1t)
dy2t <- diff(y2t)Next, we redo the ADF test for both \(\Delta y_{1,t}\) and \(\Delta y_{2,t}\):
dy1t.adf <- adf.test(dy1t)
dy1t.adf
Augmented Dickey-Fuller Test
data: dy1t
Dickey-Fuller = -8.9515, Lag order = 10, p-value = 0.01
alternative hypothesis: stationarydy2t.adf <- adf.test(dy2t)
dy1t.adf
Augmented Dickey-Fuller Test
data: dy1t
Dickey-Fuller = -8.9515, Lag order = 10, p-value = 0.01
alternative hypothesis: stationaryFrom the ADF test above for \(\Delta y_{1,t}\) and \(\Delta y_{2,t}\), we can see that both of the series were stationary, indicated by p-values that are less than 0.05. We can now proceed to the next step to estimate the reduced form VAR model.
(b) Bivariate Reduced Form VAR
Before we can estimate the reduced form CAR for \(y_{t} = ( \Delta y_{1,t} , \Delta y_{2,t})'\), we need to select the number of lags. In practice we have several info criterion that we can use to determine the number of lags, but for now, we will use the AIC.
y <- cbind(dy1t, dy2t)
y <- na.trim(y)
var1 <- VAR(y, lag.max=8, ic="AIC", type="const")
var1
VAR Estimation Results:
=======================
Estimated coefficients for equation dy1t:
=========================================
Call:
dy1t = dy1t.l1 + dy2t.l1 + dy1t.l2 + dy2t.l2 + dy1t.l3 + dy2t.l3 + dy1t.l4 + dy2t.l4 + dy1t.l5 + dy2t.l5 + dy1t.l6 + dy2t.l6 + dy1t.l7 + dy2t.l7 + const
dy1t.l1 dy2t.l1 dy1t.l2 dy2t.l2 dy1t.l3
0.3119228703 0.2181774249 0.0952428204 -0.0930330560 0.0809581628
dy2t.l3 dy1t.l4 dy2t.l4 dy1t.l5 dy2t.l5
-0.0006986961 0.0404475882 -0.1110367505 -0.0773211659 -0.0249262474
dy1t.l6 dy2t.l6 dy1t.l7 dy2t.l7 const
0.0133328809 -0.2069329954 0.0297572787 -0.1086330658 0.0021726752
Estimated coefficients for equation dy2t:
=========================================
Call:
dy2t = dy1t.l1 + dy2t.l1 + dy1t.l2 + dy2t.l2 + dy1t.l3 + dy2t.l3 + dy1t.l4 + dy2t.l4 + dy1t.l5 + dy2t.l5 + dy1t.l6 + dy2t.l6 + dy1t.l7 + dy2t.l7 + const
dy1t.l1 dy2t.l1 dy1t.l2 dy2t.l2 dy1t.l3
-0.0058262613 0.4433806900 0.0290656163 0.0439567059 0.0067641397
dy2t.l3 dy1t.l4 dy2t.l4 dy1t.l5 dy2t.l5
0.0681148385 -0.0021679389 0.0404206993 -0.0199103348 0.0548126410
dy1t.l6 dy2t.l6 dy1t.l7 dy2t.l7 const
0.0241330113 0.0598067743 0.0038612591 0.0932928562 0.0004596373 summary(var1)
VAR Estimation Results:
=========================
Endogenous variables: dy1t, dy2t
Deterministic variables: const
Sample size: 822
Log Likelihood: 6452.474
Roots of the characteristic polynomial:
0.9028 0.7865 0.7371 0.7371 0.7223 0.7223 0.6597 0.6597 0.6512 0.6512 0.6193 0.6193 0.4736 0.4736
Call:
VAR(y = y, type = "const", lag.max = 8, ic = "AIC")
Estimation results for equation dy1t:
=====================================
dy1t = dy1t.l1 + dy2t.l1 + dy1t.l2 + dy2t.l2 + dy1t.l3 + dy2t.l3 + dy1t.l4 + dy2t.l4 + dy1t.l5 + dy2t.l5 + dy1t.l6 + dy2t.l6 + dy1t.l7 + dy2t.l7 + const
Estimate Std. Error t value Pr(>|t|)
dy1t.l1 0.3119229 0.0351873 8.865 < 2e-16 ***
dy2t.l1 0.2181774 0.1141490 1.911 0.05632 .
dy1t.l2 0.0952428 0.0367932 2.589 0.00981 **
dy2t.l2 -0.0930331 0.1245860 -0.747 0.45544
dy1t.l3 0.0809582 0.0369912 2.189 0.02891 *
dy2t.l3 -0.0006987 0.1242586 -0.006 0.99551
dy1t.l4 0.0404476 0.0370572 1.091 0.27538
dy2t.l4 -0.1110368 0.1239133 -0.896 0.37048
dy1t.l5 -0.0773212 0.0368772 -2.097 0.03633 *
dy2t.l5 -0.0249262 0.1232141 -0.202 0.83973
dy1t.l6 0.0133329 0.0368413 0.362 0.71752
dy2t.l6 -0.2069330 0.1208461 -1.712 0.08721 .
dy1t.l7 0.0297573 0.0350855 0.848 0.39661
dy2t.l7 -0.1086331 0.1116065 -0.973 0.33067
const 0.0021727 0.0004912 4.423 1.11e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.008703 on 807 degrees of freedom
Multiple R-Squared: 0.1932, Adjusted R-squared: 0.1792
F-statistic: 13.8 on 14 and 807 DF, p-value: < 2.2e-16
Estimation results for equation dy2t:
=====================================
dy2t = dy1t.l1 + dy2t.l1 + dy1t.l2 + dy2t.l2 + dy1t.l3 + dy2t.l3 + dy1t.l4 + dy2t.l4 + dy1t.l5 + dy2t.l5 + dy1t.l6 + dy2t.l6 + dy1t.l7 + dy2t.l7 + const
Estimate Std. Error t value Pr(>|t|)
dy1t.l1 -0.0058263 0.0108062 -0.539 0.58993
dy2t.l1 0.4433807 0.0350559 12.648 < 2e-16 ***
dy1t.l2 0.0290656 0.0112994 2.572 0.01028 *
dy2t.l2 0.0439567 0.0382611 1.149 0.25095
dy1t.l3 0.0067641 0.0113602 0.595 0.55173
dy2t.l3 0.0681148 0.0381606 1.785 0.07464 .
dy1t.l4 -0.0021679 0.0113805 -0.190 0.84897
dy2t.l4 0.0404207 0.0380545 1.062 0.28847
dy1t.l5 -0.0199103 0.0113252 -1.758 0.07912 .
dy2t.l5 0.0548126 0.0378398 1.449 0.14785
dy1t.l6 0.0241330 0.0113142 2.133 0.03323 *
dy2t.l6 0.0598068 0.0371126 1.611 0.10746
dy1t.l7 0.0038613 0.0107750 0.358 0.72017
dy2t.l7 0.0932929 0.0342750 2.722 0.00663 **
const 0.0004596 0.0001509 3.047 0.00239 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.002673 on 807 degrees of freedom
Multiple R-Squared: 0.4078, Adjusted R-squared: 0.3975
F-statistic: 39.69 on 14 and 807 DF, p-value: < 2.2e-16
Covariance matrix of residuals:
dy1t dy2t
dy1t 7.575e-05 7.511e-07
dy2t 7.511e-07 7.144e-06
Correlation matrix of residuals:
dy1t dy2t
dy1t 1.00000 0.03229
dy2t 0.03229 1.00000After we estimate the bivariate reduced VAR model above, next we can plot the impulse response function using the following code:
var1.irfs <- irf(var1, n.ahead=12)
par(mfcol=c(2,2), cex=0.6)
plot(var1.irfs, plot.type="single")
(c) Restricted SVAR
To imposed the restriction that demand shock do not affect industrial production in the long run, we can estimate SVAR model with Blanchard and Quah approach using the following code:
var2 <- BQ(var1)
summary(var2)
SVAR Estimation Results:
========================
Call:
BQ(x = var1)
Type: Blanchard-Quah
Sample size: 822
Log Likelihood: 6437.336
Estimated contemporaneous impact matrix:
dy1t dy2t
dy1t 0.007721 0.004017
dy2t -0.001156 0.002410
Estimated identified long run impact matrix:
dy1t dy2t
dy1t 0.01706 0.00000
dy2t -0.00277 0.01228
Covariance matrix of reduced form residuals (*100):
dy1t dy2t
dy1t 7.575e-03 7.511e-05
dy2t 7.511e-05 7.144e-04(d) Contemporaneous and Long-Run Impact
The long run cumulative effect of any demand shock on industrial production, \(\Delta y_{1,t}\), is 0. Meanwhile, the long run cumulative effect of a single positive one standard deviation productivity shocks on \(\Delta y_{1,t}\) is to increase it by 0.018%
(e) IRF for SVAR
var2.irfs <- irf(var2, n.ahead=40)
var2.irfs.c <- irf(var2, n.ahead=40, ci=.9, cumulative=TRUE)
var2.irfs$irf[[2]] <- -var2.irfs$irf[[2]]
var2.irfs$Lower[[2]] <- -var2.irfs$Lower[[2]]
var2.irfs$Upper[[2]] <- -var2.irfs$Upper[[2]]
var2.irfs.c$irf[[2]] <- -var2.irfs.c$irf[[2]]
var2.irfs.c$Lower[[2]] <- -var2.irfs.c$Lower[[2]]
var2.irfs.c$Upper[[2]] <- -var2.irfs.c$Upper[[2]]
setwd("C:/Users/Rullan Rinaldi/Google Drive/ECO - 5316/HW 5")
source("plotIRF.R")
plotIRF(var2.irfs.c, vnames="dy1t", vlabels="Industrial Production", slabels=c("Productivity Shock","Demand shock"))
plotIRF(var2.irfs.c, vnames="dy2t", vlabels="CPI", slabels=c("Productivity Shock","Demand shock"))
From the cumulative IRF function for industrial production, we can see that productivity shock has a positive cumulative effect. While the effect of demand shock has the peak at approximately 4 quarter and for each one standard deviation in demand shock, at the peak industrial production fall by -0.01%.
(f) FEVD for SVAR
plot( fevd(var2, n.ahead=40) ,addbars=10 )
Overall fluctuation in \(\Delta y_{1,t}\) is approximately 60% of the time can be explained by its own lagging value and the rest by \(\Delta y_{2,t}\). In the long-run the pattern does not change by much.
Overall fluctuation in \(\Delta y_{2,t}\) is approximately 20% of the time can be explained by its own lagging value and the rest by \(\Delta y_{1,t}\). In the long-run the pattern does not change by much.