library("zoo")
library("forecast")
library("tseries")
library("urca")
library ("Quandl")
library("vars")
library("Quandl")
Quandl.api_key("DLk9RQrfTVkD4UTKc7op")
library("VARsignR")
library("stargazer")

Problem 2

Obtain monthly time series for Industrial Production Index FRED/INDPRO and for Consumer Price Index FRED/CPIAUCSL.

Introduction and Data :

In this problem, data is analysed about Industrial Production Index and Consumer Price Index from 1947 to 2015. The monthly data for this problem is collectd from Quandl.

IPI_d<-Quandl("FRED/INDPRO", type="zoo")
IPI <- window(IPI_d[,2], start=1947, end=2015)
str(IPI_d)
## 'zooreg' series from Jan 1919 to Feb 2016
##   Data: num [1:1166] 5.06 4.83 4.7 4.78 4.81 ...
##   Index: Class 'yearmon'  num [1:1166] 1919 1919 1919 1919 1919 ...
##   Frequency: 12
head(IPI)
## Jan 1947 Feb 1947 Mar 1947 Apr 1947 May 1947 Jun 1947 
##  14.1973  14.2811  14.3650  14.2532  14.3091  14.3091
tail(IPI)
## Aug 2014 Sep 2014 Oct 2014 Nov 2014 Dec 2014 Jan 2015 
## 105.1989 105.5750 105.6441 106.6868 106.5182 105.9906
tsdisplay(IPI, lag.max = 100, xlab=" Years 1947-2015", ylab= "Index 2007=100 Seasonally Adjusted ", main ="Industrial Production Index from 1947 to 2015,Monthly")

CPIUC_d<-Quandl("FRED/CPIAUCSL", type="zoo")
CPIUC <- window(CPIUC_d, start=1947, end=2015)

str(CPIUC)
## 'zooreg' series from Jan 1947 to Jan 2015
##   Data: num [1:817] 21.5 21.6 22 22 21.9 ...
##   Index: Class 'yearmon'  num [1:817] 1947 1947 1947 1947 1947 ...
##   Frequency: 12
head(CPIUC)
## Jan 1947 Feb 1947 Mar 1947 Apr 1947 May 1947 Jun 1947 
##    21.48    21.62    22.00    22.00    21.95    22.08
tail(CPIUC)
## Aug 2014 Sep 2014 Oct 2014 Nov 2014 Dec 2014 Jan 2015 
##  237.163  237.510  237.651  237.261  236.464  234.954
tsdisplay(CPIUC, lag.max = 100, xlab=" Years 1947-2015", ylab= "Index 1982-84=100 Seasonally Adjusted,  ", main ="Consumer Price Index for All Urban Consumers: All Items (USA Inflation) from 1947 to 2015,Monthly")

The monthly data for Industrial Production Index from 1947 to 2015 has increasing trend with some volatility. The PACF does not have peaks, while ACF peaks decays very slowly.

The monthly data for Consumer Price Index for All Urban Consumers from 1947 to 2015 has mostly increasing trend. The PACF does not have peaks, while ACF peaks decays very slowly.

(a) Test the log transformed time series, log of industrial production \(y_{1,t} = logIPI_t\) and log of consumer price index \(y_{2,t} = log CPI_t\) for the presence of unit root/stationarity using either ADF, ERS or KPSS tests. Afterwards apply the same unit root/stationarity test also to the first differences \(\Delta y_{1,t}\) (which approximates the month-over-moth growth rate of the industrial production) and and \(\Delta y_{2,t}\) (which approximates the month-over-moth inflation rate).

l.IPI= log(IPI)
diff.IPI=diff(IPI)
adf.l.IPI=adf.test(l.IPI)
adf.l.IPI
## 
##  Augmented Dickey-Fuller Test
## 
## data:  l.IPI
## Dickey-Fuller = -2.17, Lag order = 9, p-value = 0.5063
## alternative hypothesis: stationary
adf.diff.IPI=adf.test(diff.IPI)
adf.diff.IPI
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff.IPI
## Dickey-Fuller = -6.979, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
l.CPIUC= log(CPIUC)
diff.CPIUC=diff(CPIUC)
adf.l.CPIUC=adf.test(l.CPIUC)
adf.l.CPIUC
## 
##  Augmented Dickey-Fuller Test
## 
## data:  l.CPIUC
## Dickey-Fuller = -1.2692, Lag order = 9, p-value = 0.8876
## alternative hypothesis: stationary
adf.diff.CPIUC=adf.test(diff.CPIUC)
adf.diff.CPIUC
## 
##  Augmented Dickey-Fuller Test
## 
## data:  diff.CPIUC
## Dickey-Fuller = -5.6394, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
Industrial Production Index

The p-value for log transformed series is more than 0.05, while first difference is less than 0.05. The first difference is non-stationary.

Consumer Price Index for All Urban Consumers

The p-value for log transformed series is more than 0.05, while first difference is less than 0.05. The first difference is non-stationary.

(b) Estimate a bivariate reduced form VAR for \(y_t = (\Delta y_{1,t}, \Delta y_{2,t})'\), use information criteria to select number of lags.

y1= diff.IPI
y2= diff.CPIUC
y.Q <- cbind(y1, y2)
VARselect(y.Q, lag.max=10, type="const")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     10     10      3     10 
## 
## $criteria
##                  1           2           3           4           5
## AIC(n) -3.90360229 -3.98736075 -4.01608779 -4.03829624 -4.05073610
## HQ(n)  -3.89018901 -3.96500529 -3.98479014 -3.99805640 -4.00155407
## SC(n)  -3.86867363 -3.92914631 -3.93458758 -3.93351025 -3.92266433
## FPE(n)  0.02016913  0.01854861  0.01802335  0.01762751  0.01740961
##                  6           7           8           9          10
## AIC(n) -4.06089057 -4.06057378 -4.05781474 -4.07207329 -4.10123388
## HQ(n)  -4.00276636 -3.99350738 -3.98180615 -3.98712251 -4.00734092
## SC(n)  -3.90953303 -3.88593046 -3.85988565 -3.85085842 -3.85673324
## FPE(n)  0.01723376  0.01723927  0.01728697  0.01704232  0.01655263
var1 <- VAR(y.Q, p=10, type="const")
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: y1, y2 
## Deterministic variables: const 
## Sample size: 806 
## Log Likelihood: -592.532 
## Roots of the characteristic polynomial:
## 0.9246 0.8693 0.8693 0.8584 0.8584 0.8164 0.8164 0.8112 0.8112 0.8057 0.7836 0.7836 0.7666 0.7666 0.7611 0.7611 0.6972 0.6972 0.6363 0.2422
## Call:
## VAR(y = y.Q, p = 10, type = "const")
## 
## 
## Estimation results for equation y1: 
## =================================== 
## y1 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + y1.l3 + y2.l3 + y1.l4 + y2.l4 + y1.l5 + y2.l5 + y1.l6 + y2.l6 + y1.l7 + y2.l7 + y1.l8 + y2.l8 + y1.l9 + y2.l9 + y1.l10 + y2.l10 + const 
## 
##        Estimate Std. Error t value Pr(>|t|)    
## y1.l1   0.16551    0.03549   4.664 3.65e-06 ***
## y2.l1   0.18183    0.04810   3.780 0.000169 ***
## y1.l2   0.13752    0.03599   3.821 0.000143 ***
## y2.l2  -0.08089    0.05565  -1.453 0.146506    
## y1.l3   0.13885    0.03723   3.730 0.000205 ***
## y2.l3   0.01020    0.05577   0.183 0.854900    
## y1.l4   0.10816    0.03757   2.879 0.004102 ** 
## y2.l4   0.04892    0.05568   0.879 0.379882    
## y1.l5  -0.03529    0.03782  -0.933 0.351137    
## y2.l5  -0.09677    0.05595  -1.730 0.084094 .  
## y1.l6   0.03089    0.03797   0.814 0.416153    
## y2.l6  -0.05232    0.05598  -0.935 0.350297    
## y1.l7  -0.04842    0.03789  -1.278 0.201635    
## y2.l7   0.01009    0.05570   0.181 0.856291    
## y1.l8   0.03269    0.03747   0.873 0.383181    
## y2.l8  -0.04999    0.05572  -0.897 0.369944    
## y1.l9   0.05928    0.03679   1.611 0.107505    
## y2.l9   0.02884    0.05480   0.526 0.598809    
## y1.l10 -0.01865    0.03614  -0.516 0.605880    
## y2.l10 -0.18500    0.04765  -3.882 0.000112 ***
## const   0.09784    0.02359   4.147 3.74e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.4099 on 785 degrees of freedom
## Multiple R-Squared: 0.2136,  Adjusted R-squared: 0.1936 
## F-statistic: 10.66 on 20 and 785 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation y2: 
## =================================== 
## y2 = y1.l1 + y2.l1 + y1.l2 + y2.l2 + y1.l3 + y2.l3 + y1.l4 + y2.l4 + y1.l5 + y2.l5 + y1.l6 + y2.l6 + y1.l7 + y2.l7 + y1.l8 + y2.l8 + y1.l9 + y2.l9 + y1.l10 + y2.l10 + const 
## 
##         Estimate Std. Error t value Pr(>|t|)    
## y1.l1   0.026951   0.026509   1.017  0.30962    
## y2.l1   0.572865   0.035929  15.944  < 2e-16 ***
## y1.l2   0.145697   0.026880   5.420 7.92e-08 ***
## y2.l2  -0.115343   0.041570  -2.775  0.00566 ** 
## y1.l3   0.029186   0.027805   1.050  0.29420    
## y2.l3  -0.005237   0.041655  -0.126  0.89997    
## y1.l4  -0.041939   0.028064  -1.494  0.13547    
## y2.l4   0.119621   0.041592   2.876  0.00414 ** 
## y1.l5  -0.089712   0.028251  -3.176  0.00155 ** 
## y2.l5  -0.093677   0.041791  -2.242  0.02527 *  
## y1.l6  -0.005527   0.028365  -0.195  0.84555    
## y2.l6   0.082916   0.041815   1.983  0.04773 *  
## y1.l7   0.061482   0.028299   2.173  0.03011 *  
## y2.l7   0.071738   0.041607   1.724  0.08507 .  
## y1.l8   0.043253   0.027987   1.545  0.12263    
## y2.l8  -0.029670   0.041622  -0.713  0.47615    
## y1.l9  -0.082868   0.027480  -3.016  0.00265 ** 
## y2.l9   0.032423   0.040931   0.792  0.42851    
## y1.l10 -0.000556   0.026991  -0.021  0.98357    
## y2.l10  0.147139   0.035594   4.134 3.95e-05 ***
## const   0.045043   0.017622   2.556  0.01077 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3062 on 785 degrees of freedom
## Multiple R-Squared: 0.396,   Adjusted R-squared: 0.3806 
## F-statistic: 25.73 on 20 and 785 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           y1        y2
## y1  0.168029 -0.005433
## y2 -0.005433  0.093747
## 
## Correlation matrix of residuals:
##          y1       y2
## y1  1.00000 -0.04329
## y2 -0.04329  1.00000

The information criteria suggestes lag 10.

(c) Suppose that we want to analyze effects of two shocks - productivity shocks and demand shocks. Use Blanchard and Quah approach to obtain an SVAR model where we impose the condition that demand shocks do not affect industrial production y1,t in the long run.

# demean the data
y <- sweep(y.Q, 2, apply(y.Q, 2, mean))
# estimate reduced form VAR
myVAR <- VAR(y, ic="SC", lag.max=10)
# Blanchard-Quah long run restriction: row 1 column 2 element of the cumulative effect matrix is 0
mySVAR <- BQ(myVAR)
summary(mySVAR)
## 
## SVAR Estimation Results:
## ======================== 
## 
## Call:
## BQ(x = myVAR)
## 
## Type: Blanchard-Quah 
## Sample size: 813 
## Log Likelihood: -661.523 
## 
## Estimated contemporaneous impact matrix:
##          y1       y2
## y1  0.41534 -0.02961
## y2 -0.00482  0.31840
## 
## Estimated identified long run impact matrix:
##        y1     y2
## y1 0.9209 0.0000
## y2 0.2090 0.7086
## 
## Covariance matrix of reduced form residuals (*100):
##        y1     y2
## y1 17.339 -1.143
## y2 -1.143 10.140

(d) Report and interpret the contemporaneous impact and the long run impact matrices for the SVAR.

Interpretation

The contemporaneous impact

Here, on impact a positive one standard deviation technology shock increases Industrial Production Index by 0.415% and lowers Consumer Price Index for All Urban Consumers by 0.00482 percentage points.
A negative one standard deviation non-technology shock lowers Industrial Production Index on impact by 0.0296%, increases Consumer Price Index for All Urban Consumers by 0.3184 percentage points.

The long run impact matrices

The long run cumulative effect of any non-technology shock on Industrial Production Index is 0 (This is the long run constraint we imposed) In the long run cumulative effect of a single positive one standard deviation technology shocks on Industrial Production Index is to increase it by 0.9209%

(e) Plot the cumulative IRFs based on the SVAR model from (c) and interpret them/explain what say about the effects of the two types of shocks on industrial production and prices.

# standard non-cumulative IRFs
myIRF <- irf(mySVAR, n.ahead=20, ci=.9)
# cumulative IRFs
myIRF.c <- irf(mySVAR, n.ahead=20, ci=.9, cumulative=TRUE)
# change signs for non-technology shock to show effects of a positive shock, not a negative one
myIRF$irf[[2]] <- -myIRF$irf[[2]]
myIRF$Lower[[2]] <- -myIRF$Lower[[2]]
myIRF$Upper[[2]] <- -myIRF$Upper[[2]]
myIRF.c$irf[[2]] <- -myIRF.c$irf[[2]]
myIRF.c$Lower[[2]] <- -myIRF.c$Lower[[2]]
myIRF.c$Upper[[2]] <- -myIRF.c$Upper[[2]]
# graph with cumulative IRF for dlrGDP but a standard non-cumulative IRF for UR
par(mfrow=c(2,2), cex=.9 )
plot(myIRF.c, vnames="diff.IPI", vlabels="diff(industry index)",slabels=c("technology shock","non-technology shock"))

plot(myIRF, vnames="diff.CPIUC", vlabels="diff(CPI)",slabels=c("technology shock","non-technology shock"))

(f) Construct the FEVD for the SVAR. How much of the overall fluctuations in \(\Delta y_{1,t}\) and \(\Delta y_{2,t}\) is explained in the short run by the two shocks? How about in the long run?

var1.fevd <- fevd(mySVAR, n.ahead=20)
var1.fevd[[1]][c(1,3,6,10),]
##             y1         y2
## [1,] 0.9949421 0.00505786
## [2,] 0.9867129 0.01328715
## [3,] 0.9872113 0.01278871
## [4,] 0.9873121 0.01268790
var1.fevd[[2]][c(1,3,6,10),]
##                y1        y2
## [1,] 0.0002291574 0.9997708
## [2,] 0.0111226948 0.9888773
## [3,] 0.0286390823 0.9713609
## [4,] 0.0352540506 0.9647459
plot(var1.fevd)