I first obtain quarterly time series data for Industrial Production Index using FRED/INDPRO and for the Consumer Price Index FRED/CPIAUCSL.

First we load the “Quandl” packages and etc.

library("Quandl")
library("vars")
library("gdata")
library("stargazer")
library("tseries")

A) Test for Unit Root on Log Transformed and the First Differences

First let us find the log transformed and first differences of each variable.

ipi<-Quandl("FRED/INDPRO", type="zoo")
cpi<-Quandl("FRED/CPIAUCSL", type="zoo")
ipi<-window(ipi,start="JAN 1947", end="Feb 2016")
cpi<-window(cpi,start="JAN 1947", end="Feb 2016")
lipi<-log(ipi)
lcpi<-log(cpi)
dlipi<-diff(lipi,1)
dlcpi<-diff(lcpi,1)
par(mfrow=c(2,2))
plot(lipi,xlab="Time", ylab="Log IPI")
plot(lcpi,xlab="Time", ylab="Log CPI")
plot(dlipi,xlab="Time", ylab="First Diff Log IPI")
plot(dlcpi,xlab="Time", ylab="First Diff Log CPI")

We can clearly see that log IPI and log CPI are not stationary. However, the first differences of both variables do look somewhat stationary. Lets conduct unit root tests on both.

adf.test(lipi)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lipi
## Dickey-Fuller = -2.0144, Lag order = 9, p-value = 0.5722
## alternative hypothesis: stationary
adf.test(lcpi)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lcpi
## Dickey-Fuller = -1.2868, Lag order = 9, p-value = 0.8802
## alternative hypothesis: stationary
adf.test(dlipi)
## Warning in adf.test(dlipi): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlipi
## Dickey-Fuller = -8.0898, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
adf.test(dlcpi)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlcpi
## Dickey-Fuller = -3.9607, Lag order = 9, p-value = 0.01096
## alternative hypothesis: stationary

We see that the log variables have a unit root, and we reject that there is a unit root for the first differences of the log of IPI and CPI.

B) Estimate a Bivariate Reduced Form VAR for \(Y_t\)

We first estimate the following:

y<-cbind(dlipi,dlcpi)
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
var1<-VAR(y,p=7, type="const")
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dlipi, dlcpi 
## 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 dlipi: 
## ====================================== 
## dlipi = dlipi.l1 + dlcpi.l1 + dlipi.l2 + dlcpi.l2 + dlipi.l3 + dlcpi.l3 + dlipi.l4 + dlcpi.l4 + dlipi.l5 + dlcpi.l5 + dlipi.l6 + dlcpi.l6 + dlipi.l7 + dlcpi.l7 + const 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## dlipi.l1  3.119e-01  3.519e-02   8.865  < 2e-16 ***
## dlcpi.l1  2.182e-01  1.141e-01   1.911  0.05632 .  
## dlipi.l2  9.524e-02  3.679e-02   2.589  0.00981 ** 
## dlcpi.l2 -9.303e-02  1.246e-01  -0.747  0.45544    
## dlipi.l3  8.096e-02  3.699e-02   2.189  0.02891 *  
## dlcpi.l3 -6.987e-04  1.243e-01  -0.006  0.99551    
## dlipi.l4  4.045e-02  3.706e-02   1.091  0.27538    
## dlcpi.l4 -1.110e-01  1.239e-01  -0.896  0.37048    
## dlipi.l5 -7.732e-02  3.688e-02  -2.097  0.03633 *  
## dlcpi.l5 -2.493e-02  1.232e-01  -0.202  0.83973    
## dlipi.l6  1.333e-02  3.684e-02   0.362  0.71752    
## dlcpi.l6 -2.069e-01  1.208e-01  -1.712  0.08721 .  
## dlipi.l7  2.976e-02  3.509e-02   0.848  0.39661    
## dlcpi.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 dlcpi: 
## ====================================== 
## dlcpi = dlipi.l1 + dlcpi.l1 + dlipi.l2 + dlcpi.l2 + dlipi.l3 + dlcpi.l3 + dlipi.l4 + dlcpi.l4 + dlipi.l5 + dlcpi.l5 + dlipi.l6 + dlcpi.l6 + dlipi.l7 + dlcpi.l7 + const 
## 
##            Estimate Std. Error t value Pr(>|t|)    
## dlipi.l1 -5.826e-03  1.081e-02  -0.539  0.58993    
## dlcpi.l1  4.434e-01  3.506e-02  12.648  < 2e-16 ***
## dlipi.l2  2.907e-02  1.130e-02   2.572  0.01028 *  
## dlcpi.l2  4.396e-02  3.826e-02   1.149  0.25095    
## dlipi.l3  6.764e-03  1.136e-02   0.595  0.55173    
## dlcpi.l3  6.811e-02  3.816e-02   1.785  0.07464 .  
## dlipi.l4 -2.168e-03  1.138e-02  -0.190  0.84897    
## dlcpi.l4  4.042e-02  3.805e-02   1.062  0.28847    
## dlipi.l5 -1.991e-02  1.133e-02  -1.758  0.07912 .  
## dlcpi.l5  5.481e-02  3.784e-02   1.449  0.14785    
## dlipi.l6  2.413e-02  1.131e-02   2.133  0.03323 *  
## dlcpi.l6  5.981e-02  3.711e-02   1.611  0.10746    
## dlipi.l7  3.861e-03  1.077e-02   0.358  0.72017    
## dlcpi.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:
##           dlipi     dlcpi
## dlipi 7.575e-05 7.511e-07
## dlcpi 7.511e-07 7.144e-06
## 
## Correlation matrix of residuals:
##         dlipi   dlcpi
## dlipi 1.00000 0.03229
## dlcpi 0.03229 1.00000

So while we use the BIC to select 7 lags, the output does not look pretty.

C) Obtain an SVAR Model Using Blanchard and Quah

I use the condition that demand shocks do not affect industrial production in the long run. I impose a B0 matrix with a restriction on the upper right half. I assume that CPI would be a demand shock.

svar1<-VAR(y,ic="AIC", lag.max=8)
mysvar1<-BQ(svar1)
summary(mysvar1)
## 
## SVAR Estimation Results:
## ======================== 
## 
## Call:
## BQ(x = svar1)
## 
## Type: Blanchard-Quah 
## Sample size: 822 
## Log Likelihood: 6437.336 
## 
## Estimated contemporaneous impact matrix:
##           dlipi    dlcpi
## dlipi  0.007721 0.004017
## dlcpi -0.001156 0.002410
## 
## Estimated identified long run impact matrix:
##          dlipi   dlcpi
## dlipi  0.01706 0.00000
## dlcpi -0.00277 0.01228
## 
## Covariance matrix of reduced form residuals (*100):
##           dlipi     dlcpi
## dlipi 7.575e-03 7.511e-05
## dlcpi 7.511e-05 7.144e-04

D) Interpret Results

When we look at the contemporaneous matrix, we see that IPI and CPI both have positive effects on the IPI. However, on the CPI, IPI has a negative affect, while CPI has a positive effect.

The long run impact matrix of has a similar effect, but stronger. IPI influences itself, while CPI does not. CPI is negatively influenced by IPI and positive influenced by CPI.

E) Plot the Cumulative IRFs and Interpret

So I plot the cumulative IRFs of the SVAR model.

var1.irfs <- irf(var1, n.ahead=40, cumulative=TRUE)
par(mfcol=c(2,2))
plot(var1.irfs, plot.type="single")

I see that a IPI shock is permanent on IPI and a CPI shock is temporary and in fact, turns into a negative shock on IPI. When we look at CPI we see that an IPI shock has a small but persistent effect on CPI. Also, a shock to CPI has a somewhat short term effect on CPI. This makes sense as demand for IPI would be hurt if prices rise. We would also expect a small increase in CPI from increased IPI, but it would not be as large.

F) Construct the FEVD

var1.fevd <- fevd(var1, n.ahead=40)
var1.fevd[[1]][c(1,4,8,40),]
##          dlipi       dlcpi
## [1,] 1.0000000 0.000000000
## [2,] 0.9955300 0.004469961
## [3,] 0.9836736 0.016326403
## [4,] 0.9612843 0.038715696
var1.fevd[[2]][c(1,4,8,40),]
##            dlipi     dlcpi
## [1,] 0.001042582 0.9989574
## [2,] 0.012803637 0.9871964
## [3,] 0.022576285 0.9774237
## [4,] 0.036513574 0.9634864
plot(var1.fevd)

plot(var1.fevd, addbars=8)

The overall fluctuations are mostly explained by internal shocks, and very little of it is explained by a shock from each other. This holds in both the short term and the long term.

This may be because CPI is made up of many components and the basket of goods only have a small percentage of items that can be directly counted from industrial production. For instance, many goods are imported. For IPI demand, it may also have small fluctuations from CPI because business demand may be a bigger factor on IPI than consumer demand.