library("Quandl")
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library("vars")
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library("gdata")
## gdata: read.xls support for 'XLS' (Excel 97-2004) files ENABLED.
## 
## gdata: read.xls support for 'XLSX' (Excel 2007+) files ENABLED.
## 
## Attaching package: 'gdata'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following object is masked from 'package:stats':
## 
##     nobs
## The following object is masked from 'package:utils':
## 
##     object.size
library("stargazer")
## 
## Please cite as:
##  Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2. http://CRAN.R-project.org/package=stargazer
library("tseries")
Quandl.auth(" 3HQPzReHC1WfzuSG41gy")
## Warning: 'Quandl.auth' is deprecated.
## Use 'Quandl.api_key' instead.
## See help("Deprecated")
#A)Test the log transformed time series, log of industrial production y1,t = logIPIt and log of consumer price index y2,t = logCPIt 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 ∆y1,t (which approximates the month-over-moth growth rate of the industrial production) and and ∆y2,t (which approximates the month-over-moth inflation rate).
ip<-Quandl("FRED/INDPRO", type="zoo")
cp<-Quandl("FRED/CPIAUCSL", type="zoo")
ip<-window(ip,start="JAN 1947", end="Feb 2016")
cp<-window(cp,start="JAN 1947", end="Feb 2016")
lip<-log(ip)
lcp<-log(cp)
dlip<-diff(lip,1)
dlcp<-diff(lcp,1)
par(mfrow=c(2,2))
plot(lip,xlab="Time", ylab="lip")
plot(lcp,xlab="Time", ylab="Log CP")
plot(dlip,xlab="Time", ylab="dlip")
plot(dlcp,xlab="Time", ylab="dlcp")

# first difference of log ip and log cp do not have unit root 
adf.test(lip)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lip
## Dickey-Fuller = -2.0144, Lag order = 9, p-value = 0.5722
## alternative hypothesis: stationary
adf.test(lcp)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  lcp
## Dickey-Fuller = -1.2868, Lag order = 9, p-value = 0.8802
## alternative hypothesis: stationary
adf.test(dlip)
## Warning in adf.test(dlip): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlip
## Dickey-Fuller = -8.0898, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
adf.test(dlcp)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dlcp
## Dickey-Fuller = -3.9607, Lag order = 9, p-value = 0.01096
## alternative hypothesis: stationary
#B)Estimate a bivariate reduced form VAR for yt = (∆y1,t,∆y2,t)′, use information criteria to select number of lags.
y<-cbind(dlip,dlcp)
y<-na.trim(y)
y<-sweep(y,2,apply(y,2,mean))
#get rid of missing observations and mean
y<-na.trim(y)
y<-sweep(y,2,apply(y,2,mean))
#choose number of lages based on AIC, we found p=1 is good (lowest AIC)
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
# estimate bivariate reduced form VAR
var1<-VAR(y,p=1, type="const")
summary(var1)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: dlip, dlcp 
## Deterministic variables: const 
## Sample size: 828 
## Log Likelihood: 6419.952 
## Roots of the characteristic polynomial:
## 0.5766 0.3856
## Call:
## VAR(y = y, p = 1, type = "const")
## 
## 
## Estimation results for equation dlip: 
## ===================================== 
## dlip = dlip.l1 + dlcp.l1 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## dlip.l1  3.859e-01  3.213e-02  12.012   <2e-16 ***
## dlcp.l1 -1.002e-02  8.844e-02  -0.113     0.91    
## const   -7.684e-06  3.077e-04  -0.025     0.98    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.008853 on 825 degrees of freedom
## Multiple R-Squared: 0.1489,  Adjusted R-squared: 0.1468 
## F-statistic: 72.15 on 2 and 825 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation dlcp: 
## ===================================== 
## dlcp = dlip.l1 + dlcp.l1 + const 
## 
##           Estimate Std. Error t value Pr(>|t|)    
## dlip.l1 -7.066e-03  1.034e-02  -0.683    0.495    
## dlcp.l1  5.762e-01  2.846e-02  20.244   <2e-16 ***
## const   -7.467e-06  9.902e-05  -0.075    0.940    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.002849 on 825 degrees of freedom
## Multiple R-Squared: 0.3321,  Adjusted R-squared: 0.3305 
## F-statistic: 205.1 on 2 and 825 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           dlip      dlcp
## dlip 7.838e-05 3.747e-07
## dlcp 3.747e-07 8.118e-06
## 
## Correlation matrix of residuals:
##         dlip    dlcp
## dlip 1.00000 0.01485
## dlcp 0.01485 1.00000
#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.