In this task we have to obtain monthly time series for Industrial Production Index FRED/INDPRO and for Consumer Price Index FRED/CPIAUCSL.
library(Quandl)
library(vars)
library(stargazer)
library(tseries)
Our first task is to test log transformed time series for the presence of unit root/stationarity using relevant tests.
ipi <- Quandl("FRED/INDPRO", type="zoo")
cpi <- Quandl("FRED/CPIAUCSL", type="zoo")
y1 <- log(ipi)
dy1 <- diff(log(ipi))
y2 <- log(cpi)
dy2 <- diff(log(cpi))
Now we can plot the graphs both graphs of the time-series
par(mfrow=c(1,2))
plot(y1, xlab="", ylab="", main="log(IPI)")
plot(dy1, xlab="", ylab="", main="diff.log(IPI)")
plot(y2, xlab="", ylab="", main="log(CPI)")
plot(dy2, xlab="", ylab="", main="diff.log(CPI)")
Now we can conduct unit root test for transformed series.
adf.test(y1)
##
## Augmented Dickey-Fuller Test
##
## data: y1
## Dickey-Fuller = -2.7713, Lag order = 10, p-value = 0.2518
## alternative hypothesis: stationary
adf.test(dy1)
## Warning in adf.test(dy1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dy1
## Dickey-Fuller = -8.9515, Lag order = 10, p-value = 0.01
## alternative hypothesis: stationary
adf.test(y2)
##
## Augmented Dickey-Fuller Test
##
## data: y2
## Dickey-Fuller = -1.2868, Lag order = 9, p-value = 0.8802
## alternative hypothesis: stationary
adf.test(dy2)
##
## Augmented Dickey-Fuller Test
##
## data: dy2
## Dickey-Fuller = -3.9607, Lag order = 9, p-value = 0.01096
## alternative hypothesis: stationary
Based on ADF test we can observe that both of the first differences are stationary. But we when it comes to the logs series is non stationary (have unit root).
In this task we have to Estimate a bivariate reduced form VAR for \(y_t=(\Delta y_{1,t}, \Delta y_{2,t})'\) and use information criteria to select number of lags.
y<-cbind(dy1,dy2)
y<-na.trim(y)
y<-sweep(y,2,apply(y,2,mean))
VARselect(y,lag.max=8)
## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 7 6 2 7
##
## $criteria
## 1 2 3 4
## AIC(n) -2.123881e+01 -2.128068e+01 -2.130123e+01 -2.131146e+01
## HQ(n) -2.122560e+01 -2.125867e+01 -2.127041e+01 -2.127183e+01
## SC(n) -2.120439e+01 -2.122331e+01 -2.122091e+01 -2.120818e+01
## FPE(n) 5.971746e-10 5.726866e-10 5.610375e-10 5.553301e-10
## 5 6 7 8
## AIC(n) -2.132000e+01 -2.132956e+01 -2.133296e+01 -2.133237e+01
## HQ(n) -2.127157e+01 -2.127232e+01 -2.126692e+01 -2.125752e+01
## SC(n) -2.119377e+01 -2.118038e+01 -2.116083e+01 -2.113729e+01
## FPE(n) 5.506083e-10 5.453727e-10 5.435220e-10 5.438459e-10
AIC shows that the number of the lags is equal to 7. And the model is VAR(7)
var <- VAR(y, p=7, type="const")
summary(var)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: dy1, dy2
## 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, p = 7, type = "const")
##
##
## Estimation results for equation dy1:
## ====================================
## dy1 = dy1.l1 + dy2.l1 + dy1.l2 + dy2.l2 + dy1.l3 + dy2.l3 + dy1.l4 + dy2.l4 + dy1.l5 + dy2.l5 + dy1.l6 + dy2.l6 + dy1.l7 + dy2.l7 + const
##
## Estimate Std. Error t value Pr(>|t|)
## dy1.l1 3.119e-01 3.519e-02 8.865 < 2e-16 ***
## dy2.l1 2.182e-01 1.141e-01 1.911 0.05632 .
## dy1.l2 9.524e-02 3.679e-02 2.589 0.00981 **
## dy2.l2 -9.303e-02 1.246e-01 -0.747 0.45544
## dy1.l3 8.096e-02 3.699e-02 2.189 0.02891 *
## dy2.l3 -6.987e-04 1.243e-01 -0.006 0.99551
## dy1.l4 4.045e-02 3.706e-02 1.091 0.27538
## dy2.l4 -1.110e-01 1.239e-01 -0.896 0.37048
## dy1.l5 -7.732e-02 3.688e-02 -2.097 0.03633 *
## dy2.l5 -2.493e-02 1.232e-01 -0.202 0.83973
## dy1.l6 1.333e-02 3.684e-02 0.362 0.71752
## dy2.l6 -2.069e-01 1.208e-01 -1.712 0.08721 .
## dy1.l7 2.976e-02 3.509e-02 0.848 0.39661
## dy2.l7 -1.086e-01 1.116e-01 -0.973 0.33067
## const 8.453e-06 3.036e-04 0.028 0.97779
## ---
## 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 dy2:
## ====================================
## dy2 = dy1.l1 + dy2.l1 + dy1.l2 + dy2.l2 + dy1.l3 + dy2.l3 + dy1.l4 + dy2.l4 + dy1.l5 + dy2.l5 + dy1.l6 + dy2.l6 + dy1.l7 + dy2.l7 + const
##
## Estimate Std. Error t value Pr(>|t|)
## dy1.l1 -5.826e-03 1.081e-02 -0.539 0.58993
## dy2.l1 4.434e-01 3.506e-02 12.648 < 2e-16 ***
## dy1.l2 2.907e-02 1.130e-02 2.572 0.01028 *
## dy2.l2 4.396e-02 3.826e-02 1.149 0.25095
## dy1.l3 6.764e-03 1.136e-02 0.595 0.55173
## dy2.l3 6.811e-02 3.816e-02 1.785 0.07464 .
## dy1.l4 -2.168e-03 1.138e-02 -0.190 0.84897
## dy2.l4 4.042e-02 3.805e-02 1.062 0.28847
## dy1.l5 -1.991e-02 1.133e-02 -1.758 0.07912 .
## dy2.l5 5.481e-02 3.784e-02 1.449 0.14785
## dy1.l6 2.413e-02 1.131e-02 2.133 0.03323 *
## dy2.l6 5.981e-02 3.711e-02 1.611 0.10746
## dy1.l7 3.861e-03 1.077e-02 0.358 0.72017
## dy2.l7 9.329e-02 3.428e-02 2.722 0.00663 **
## const -2.298e-05 9.323e-05 -0.246 0.80537
## ---
## 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:
## dy1 dy2
## dy1 7.575e-05 7.511e-07
## dy2 7.511e-07 7.144e-06
##
## Correlation matrix of residuals:
## dy1 dy2
## dy1 1.00000 0.03229
## dy2 0.03229 1.00000
In this part we will analyze effects of two shocks - productivity shocks and demand shocks.
myVAR<-VAR(y,ic="AIC", lag.max=8)
mySVAR<-BQ(myVAR)
summary(mySVAR)
##
## SVAR Estimation Results:
## ========================
##
## Call:
## BQ(x = myVAR)
##
## Type: Blanchard-Quah
## Sample size: 822
## Log Likelihood: 6437.336
##
## Estimated contemporaneous impact matrix:
## dy1 dy2
## dy1 0.007721 0.004017
## dy2 -0.001156 0.002410
##
## Estimated identified long run impact matrix:
## dy1 dy2
## dy1 0.01706 0.00000
## dy2 -0.00277 0.01228
##
## Covariance matrix of reduced form residuals (*100):
## dy1 dy2
## dy1 7.575e-03 7.511e-05
## dy2 7.511e-05 7.144e-04
Now we plot cumulative IRFs based on the SVAR model from part (c).
irfs <- irf(mySVAR, n.ahead=40, cumulative=TRUE)
plot(irfs, plot.type="single")
From the IRF we can see that a shock positively affect the IPI, it continues to grow for about a year untill it reaches a new higher level. And when CPI is shocked it positively affects the growth of GDP. Though, inflation will be higher all the proceeding years
Now we construct the FEVD for the SVAR and access how much of the overall fluctuations in \(y_{1,t}\) and \(y_{2,t}\) is explained in the short run by the two shocks? How about in the long run?
fevd <- fevd(mySVAR, n.ahead=40)
fevd[[1]][c(1,4,8,40),]
## dy1 dy2
## [1,] 0.7869807 0.2130193
## [2,] 0.7647388 0.2352612
## [3,] 0.7640897 0.2359103
## [4,] 0.7497064 0.2502936
fevd[[2]][c(1,4,8,40),]
## dy1 dy2
## [1,] 0.1871906 0.8128094
## [2,] 0.1771762 0.8228238
## [3,] 0.1639252 0.8360748
## [4,] 0.1455469 0.8544531
plot(fevd)
As it can be seen from the plots, more than 70% of fluctuations in \(\Delta y_{1,t}\) is explained by IPI and the rest by CPI. About 80% of fluctuations in \(\Delta y_{2,t}\) is explained by CPI and the rest by IPI.