library(urca)
library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: lmtest
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
library(readxl)
my_datae <- read_excel("Desktop/my_datae.xlsx")
summary(my_datae )
##      DATE              IndProd            RGDP      
##  Length:212         Min.   : 22.13   Min.   : 2800  
##  Class :character   1st Qu.: 41.55   1st Qu.: 4819  
##  Mode  :character   Median : 55.20   Median : 7079  
##                     Mean   : 60.75   Mean   : 7665  
##                     3rd Qu.: 86.42   3rd Qu.:10859  
##                     Max.   :100.51   Max.   :13665
indprod<-ts(my_datae$IndProd,start=c(1960,2012),frequency=4)
rgdp<-ts(my_datae$RGDP,start=c(1960,2012),frequency=4)

QUESTION 3(a.i)

The time series plots of growth rates for quarterly industrial production and real GDP below have clear trends indicating non-stationarity.The correlograms have almost all spikes sticking out the threshold band. Hence there exist strong significant serial autocorrelations in the series for both growth rates.

par(mfrow=c(2,2))
plot.ts(indprod)
plot.ts(rgdp)
acf(indprod)
acf(rgdp)

QUESTION 3(a.ii)

The time series plots for the growth rate by difference in the logs for quarterly industrial production and real GDP seem stationary since they do not show any trend nor fanning out patterns.The ACF for two series under consideration show significant autocorrelations.

indprod_grate <- diff(log(indprod))*400  #  industrial production  growth rate    
rgdp_grate <- diff(log(rgdp))*400  #  real GDP growth rate  
par(mfrow=c(2,2))
plot.ts(indprod_grate)
plot.ts(rgdp_grate)
acf(indprod_grate)
acf(rgdp_grate)

QUESTION 3(b.i)

By eyeballing the joint plot of the two series below they seem to have some co-movement behavior patttern.They are correlated in a way but we’re not sure which lag(s) show a significant correlation. In order that we ascertain the correlation coefficient of which lags are significant we perform a statistical analysis which is the unit-root test for stationarity purposes. The results from our unit root tests indicate that the test statistic for each series falls within their corresponding region.Thus,we reject the null hypothesis for each test and conclude that both series are stationary.

y_t <- data.frame(indprod_grate, rgdp_grate) 
colnames(y_t) <-c("Indprod_growth", "RGDP_growth") 
plot.ts(y_t)  

ur_indprod <- ur.df(y_t$Indprod_growth, type="drift", lags=6, selectlags="AIC")
summary(ur_indprod)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.8703  -2.4788   0.3411   2.5697  16.5072 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.21779    0.39952   3.048  0.00262 ** 
## z.lag.1     -0.46088    0.07333  -6.285 2.03e-09 ***
## z.diff.lag1  0.14524    0.07963   1.824  0.06969 .  
## z.diff.lag2 -0.05537    0.07242  -0.765  0.44543    
## z.diff.lag3  0.13817    0.06739   2.050  0.04164 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.884 on 199 degrees of freedom
## Multiple R-squared:  0.2398, Adjusted R-squared:  0.2245 
## F-statistic: 15.69 on 4 and 199 DF,  p-value: 3.528e-11
## 
## 
## Value of test-statistic is: -6.2846 19.7761 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81
ur_rgdp <- ur.df(y_t$RGDP_growth, type="drift", lags=6, selectlags="AIC")
summary(ur_rgdp)
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.434  -1.608  -0.169   1.778  13.482 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.59285    0.32599   4.886 2.10e-06 ***
## z.lag.1     -0.53967    0.07897  -6.834 9.67e-11 ***
## z.diff.lag  -0.18356    0.06902  -2.659  0.00846 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.156 on 201 degrees of freedom
## Multiple R-squared:  0.3544, Adjusted R-squared:  0.348 
## F-statistic: 55.17 on 2 and 201 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -6.8336 23.3701 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.46 -2.88 -2.57
## phi1  6.52  4.63  3.81

QUESTION 3(b.ii)

We examine cross-correlation matrices. At the contemporaneous level at lag zero we see that the two series are signficantly correlated.We also see from the cross correlation coefficient function plots that industrial production growth rate significantly affect lagged real GDP growth rate at lags 1,2 and 3. Similarly, it is clear that real GDP growth rate significantly affect lagged industrial production growth rate at lags 1 and 2.Thus there exists a significant dynamic relation between the the two series.

corr_coff<-MTS::ccm(y_t)
## [1] "Covariance matrix:"
##                Indprod_growth RGDP_growth
## Indprod_growth           40.9        16.8
## RGDP_growth              16.8        11.9
## CCM at lag:  0 
##       [,1]  [,2]
## [1,] 1.000 0.763
## [2,] 0.763 1.000
## Simplified matrix: 
## CCM at lag:  1 
## + + 
## + + 
## CCM at lag:  2 
## + + 
## + + 
## CCM at lag:  3 
## + + 
## . . 
## CCM at lag:  4 
## . . 
## . . 
## CCM at lag:  5 
## . . 
## . . 
## CCM at lag:  6 
## . . 
## . . 
## CCM at lag:  7 
## . . 
## . . 
## CCM at lag:  8 
## - - 
## . . 
## CCM at lag:  9 
## . . 
## . . 
## CCM at lag:  10 
## . . 
## . . 
## CCM at lag:  11 
## . . 
## . . 
## CCM at lag:  12 
## . . 
## - .

## Hit Enter for p-value plot of individual ccm:

corr_coff$ccm
##           [,1]      [,2]      [,3]      [,4]         [,5]        [,6]
## [1,] 1.0000000 0.5970837 0.2821397 0.1575315 -0.002976263 -0.10476709
## [2,] 0.7631361 0.4402353 0.1827462 0.0834439  0.027048272 -0.04788363
## [3,] 0.7631361 0.5064331 0.3294065 0.2155684  0.064855122 -0.02705649
## [4,] 1.0000000 0.3467701 0.2745193 0.1057918  0.097306456 -0.04817318
##               [,7]        [,8]        [,9]       [,10]      [,11]       [,12]
## [1,] -7.050802e-02 -0.08456491 -0.17984834 -0.11273667 0.01058445 -0.02454379
## [2,] -1.822087e-03  0.01004278 -0.05597467 -0.01501415 0.01278770  0.02076477
## [3,] -4.705929e-02 -0.11418336 -0.18798530 -0.08470278 0.02536258 -0.02728550
## [4,]  4.097336e-05 -0.05760423 -0.05556675  0.04463847 0.04302910  0.02026571
##            [,13]
## [1,] -0.10164357
## [2,] -0.14118511
## [3,] -0.09409534
## [4,] -0.11389876
corr_coff$pvalue
##  [1] 4.440892e-16 2.034113e-05 1.875974e-02 4.683336e-01 1.030981e-01
##  [6] 6.400030e-01 1.530101e-01 1.253148e-02 5.406957e-02 9.504526e-01
## [11] 8.850382e-01 3.218145e-01

We compute the multivariate Ljung-Box Q statistics (adj=0 by default)

MTS::mq(y_t,8) 
## Ljung-Box Statistics:  
##        m       Q(m)     df    p-value
## [1,]   1.0      78.7     4.0        0
## [2,]   2.0     105.9     8.0        0
## [3,]   3.0     117.8    12.0        0
## [4,]   4.0     121.4    16.0        0
## [5,]   5.0     129.2    20.0        0
## [6,]   6.0     131.7    24.0        0
## [7,]   7.0     138.5    28.0        0
## [8,]   8.0     151.4    32.0        0

QUESTION 3(b.iii)

We further perform a portmanteau test corroborate the existence of correlation between the two series.The p-values of the test statistics in the output for the multivariate Ljung-Box test are less than 0.05. We have no properties of white noise exhibitted here. Hence it is clear that there exists a significant correlation between the series

QUESTION 3(c.i)

We perform order selection by information criteria Our best choice of lag order using the AIC criterion is 2.

lag_sct<-VARselect(y_t,type="const")
lag_sct
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2 
## 
## $criteria
##                 1          2          3          4          5          6
## AIC(n)   4.867067   4.837764   4.856137   4.855807   4.880278   4.908180
## HQ(n)    4.906967   4.904265   4.949238   4.975508   5.026579   5.081081
## SC(n)    4.965673   5.002108   5.086218   5.151625   5.241833   5.335473
## FPE(n) 129.939787 126.189513 128.533965 128.499719 131.696047 135.441703
##                 7          8          9         10
## AIC(n)   4.920353   4.921189   4.920558   4.943244
## HQ(n)    5.119855   5.147290   5.173260   5.222546
## SC(n)    5.413384   5.479957   5.545063   5.633486
## FPE(n) 137.127217 137.276716 137.234350 140.439149
lag_sct$select
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2

QUESTION 3(c.ii)

VAR_fit <- VAR(ts(y_t), p=2,type="const") 
summary(VAR_fit)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Indprod_growth, RGDP_growth 
## Deterministic variables: const 
## Sample size: 209 
## Log Likelihood: -1092.434 
## Roots of the characteristic polynomial:
## 0.4777 0.4777 0.3456 0.05199
## Call:
## VAR(y = ts(y_t), p = 2, type = "const")
## 
## 
## Estimation results for equation Indprod_growth: 
## =============================================== 
## Indprod_growth = Indprod_growth.l1 + RGDP_growth.l1 + Indprod_growth.l2 + RGDP_growth.l2 + const 
## 
##                   Estimate Std. Error t value Pr(>|t|)    
## Indprod_growth.l1  0.54837    0.09497   5.774 2.85e-08 ***
## RGDP_growth.l1     0.24477    0.15548   1.574  0.11696    
## Indprod_growth.l2 -0.26049    0.09053  -2.878  0.00443 ** 
## RGDP_growth.l2     0.38044    0.15625   2.435  0.01576 *  
## const              0.10249    0.57175   0.179  0.85792    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 4.994 on 204 degrees of freedom
## Multiple R-Squared: 0.3921,  Adjusted R-squared: 0.3802 
## F-statistic:  32.9 on 4 and 204 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation RGDP_growth: 
## ============================================ 
## RGDP_growth = Indprod_growth.l1 + RGDP_growth.l1 + Indprod_growth.l2 + RGDP_growth.l2 + const 
## 
##                   Estimate Std. Error t value Pr(>|t|)    
## Indprod_growth.l1  0.25360    0.05782   4.386 1.85e-05 ***
## RGDP_growth.l1     0.03916    0.09467   0.414  0.67957    
## Indprod_growth.l2 -0.17297    0.05512  -3.138  0.00195 ** 
## RGDP_growth.l2     0.26835    0.09514   2.821  0.00527 ** 
## const              1.86204    0.34814   5.349 2.37e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 3.041 on 204 degrees of freedom
## Multiple R-Squared: 0.2366,  Adjusted R-squared: 0.2216 
## F-statistic:  15.8 on 4 and 204 DF,  p-value: 2.773e-11 
## 
## 
## 
## Covariance matrix of residuals:
##                Indprod_growth RGDP_growth
## Indprod_growth          24.94      10.284
## RGDP_growth             10.28       9.245
## 
## Correlation matrix of residuals:
##                Indprod_growth RGDP_growth
## Indprod_growth         1.0000      0.6773
## RGDP_growth            0.6773      1.0000
stargazer(VAR_fit[["varresult"]],type="text")
## 
## ===========================================================
##                                    Dependent variable:     
##                                ----------------------------
##                                             y              
##                                     (1)            (2)     
## -----------------------------------------------------------
## Indprod_growth.l1                 0.548***      0.254***   
##                                   (0.095)        (0.058)   
##                                                            
## RGDP_growth.l1                     0.245          0.039    
##                                   (0.155)        (0.095)   
##                                                            
## Indprod_growth.l2                -0.260***      -0.173***  
##                                   (0.091)        (0.055)   
##                                                            
## RGDP_growth.l2                    0.380**       0.268***   
##                                   (0.156)        (0.095)   
##                                                            
## const                              0.102        1.862***   
##                                   (0.572)        (0.348)   
##                                                            
## -----------------------------------------------------------
## Observations                        209            209     
## R2                                 0.392          0.237    
## Adjusted R2                        0.380          0.222    
## Residual Std. Error (df = 204)     4.994          3.041    
## F Statistic (df = 4; 204)        32.895***      15.803***  
## ===========================================================
## Note:                           *p<0.1; **p<0.05; ***p<0.01
roots(VAR_fit, modulus=FALSE) # eigenvalues of the companion coefficient matrix
## [1]  0.44054728+0.1846212i  0.44054728-0.1846212i -0.34556009+0.0000000i
## [4]  0.05199498+0.0000000i
roots(VAR_fit)
## [1] 0.47766818 0.47766818 0.34556009 0.05199498

We write down the models for each series

Let {\(y_{1t}\)} be the series for the industrial production growth rate and {\(y_{2t}\)} be the series for the GDP growth rate. Then we have the following equations

\[\begin{equation}\label{eq:AR} y_{1t} = 0.10249 + 0.54837y_{(1,t-1)} + 0.24477y_{(2,t-1)} -0.26049 y_{(1,t-2)}+ 0.38044y_{(2,t-2)}+\varepsilon_t\\ y_{2t} = 1.86204 + 0.25360y_{(1,t-1)} + 0.03916y_{(2,t-1)} -0.17297 y_{(1,t-2)}+ 0.26835y_{(2,t-2)}+\varepsilon_t \end{equation}\]

The question here is, are all the coefficient significant ? check the p-values We see from the out put that the constant term and the coefficient of \(y_{(2,t-1)}\) for the equation whose response variable is \(y_{1t}\) are insignificant hence we can consider taking them out of the equation. Also, we see that the coefficient of \(y_{(1,t-1)}\) for the equation whose response variable is \(y_{2t}\) is insignificant. Thus we take the insignificant coefficients out of the equations.

We also see that VAR model is stationary since the eigen values are less than one in modulus. Thus the eigen values are all within the unit circle .We conclude that model is a stable system and hence we can proceed to perform further analysis.

res <- serial.test(VAR_fit, lags.pt=8, type="PT.adjusted")
res
## 
##  Portmanteau Test (adjusted)
## 
## data:  Residuals of VAR object VAR_fit
## Chi-squared = 30.4, df = 24, p-value = 0.1718
par(mfrow=c(2,1))
plot(res, names="Indprod_growth") 

plot(res, names="RGDP_growth")

The residual plot for each model of both series behave like those of white noises. Hence it is okay to say our model is adequate.

QUESTION 3(d) Below are the the 4-step-ahead forecasts of the growth rates for industrial production and real GDP.

(fcast <- predict(VAR_fit, n.ahead=4))
## $Indprod_growth
##          fcst      lower    upper       CI
## [1,] 2.622130  -7.165161 12.40942  9.78729
## [2,] 1.819861  -9.898883 13.53861 11.71874
## [3,] 2.223234 -10.118929 14.56540 12.34216
## [4,] 2.422800 -10.069670 14.91527 12.49247
## 
## $RGDP_growth
##          fcst     lower    upper       CI
## [1,] 3.256079 -2.703342 9.215501 5.959422
## [2,] 2.316833 -4.203471 8.837136 6.520304
## [3,] 2.834511 -3.937070 9.606092 6.771581
## [4,] 2.843796 -3.960685 9.648277 6.804481
plot(fcast)

fanchart(fcast)

QUESTION 3(e)

grg_1<-causality(VAR_fit, cause="Indprod_growth")
grg_1$Granger
## 
##  Granger causality H0: Indprod_growth do not Granger-cause RGDP_growth
## 
## data:  VAR object VAR_fit
## F-Test = 11.739, df1 = 2, df2 = 408, p-value = 1.104e-05
grg_2<-causality(VAR_fit, cause="RGDP_growth")
grg_2$Granger
## 
##  Granger causality H0: RGDP_growth do not Granger-cause Indprod_growth
## 
## data:  VAR object VAR_fit
## F-Test = 3.9628, df1 = 2, df2 = 408, p-value = 0.01975

We see from the outputs that the series do Granger cause each other since the p-values are less 0.05 level of signoficance in each case.Thus we reject the Null hypothesis for the test in each case.

QUESTION 3(e)

irf_gdp <- irf(VAR_fit, impulse="RGDP_growth", response="Indprod_growth", n.ahead=20,boot=TRUE,run=1000,ci=0.95)
plot(irf_gdp, ylab="Indprod_growth",main="Indprod growth rate response to RGDP growth rate shock")

irf_spread <- irf(VAR_fit, impulse="Indprod_growth", response="RGDP_growth", n.ahead=20,boot=TRUE,run=1000,ci=0.95)
plot(irf_spread, ylab="RGDP_growth",main="RGDP growth rate response to Indprod growth rate shock ")

plot(irf(VAR_fit)) # examine all  IRFs

Since the confidence interval bands do not wrap around the zero line, We conclude that the industrial production growth rate response to RGDP growth rate’s shock is significant. Thus if there is a shock to industrial production growth rate then the real GDP growth rate will respond positively.