I first obtain quarterly data for Disposable Personal Income FRED/DPI and forPersonal Consumption Expenditures FED/PCEC. These two time series should be co-integrated.

First we load the “Quandl” packages and etc.

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

A) Plot Time Series and Conduct Unit Root Test

We plot both series in the same plot, and also plot their first differences.

DPI<-Quandl("FRED/DPI", type="zoo")
PCE<-Quandl("FRED/PCEC", type="zoo")
DPIl<-log(DPI)
PCEl<-log(PCE)
DPIdiff<-diff(DPIl,1)
PCEdiff<-diff(PCEl,1)
par(mfrow=c(2,2), cex=0.5, mar=c(2,2,3,1))
plot(DPIl, xlab="Year", ylab="", main="Log(PCE)")
plot(PCEl, xlab="Year", ylab="", main="Log(DPI")
plot(DPIdiff, xlab="Year", ylab="", main="dLog(PCE)")
plot(PCEdiff, xlab="Year", ylab="", main="dLog(DPI)")

Now we do unit root tests to test for stationarity.

summary( ur.ers(DPIl, type="P-test", lag.max=8, model="trend") )
## 
## ############################################### 
## # Elliot, Rothenberg and Stock Unit Root Test # 
## ############################################### 
## 
## Test of type P-test 
## detrending of series with intercept and trend 
## 
## Value of test-statistic is: 83.7733 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89
summary( ur.ers(PCEl, type="P-test", lag.max=8, model="trend") )
## 
## ############################################### 
## # Elliot, Rothenberg and Stock Unit Root Test # 
## ############################################### 
## 
## Test of type P-test 
## detrending of series with intercept and trend 
## 
## Value of test-statistic is: 63.1974 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89
summary( ur.ers(DPIdiff, type="P-test", lag.max=8, model="trend") )
## 
## ############################################### 
## # Elliot, Rothenberg and Stock Unit Root Test # 
## ############################################### 
## 
## Test of type P-test 
## detrending of series with intercept and trend 
## 
## Value of test-statistic is: 2.2919 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89
summary( ur.ers(PCEdiff, type="P-test", lag.max=8, model="trend") )
## 
## ############################################### 
## # Elliot, Rothenberg and Stock Unit Root Test # 
## ############################################### 
## 
## Test of type P-test 
## detrending of series with intercept and trend 
## 
## Value of test-statistic is: 1.7594 
## 
## Critical values of P-test are:
##                 1pct 5pct 10pct
## critical values 3.96 5.62  6.89

Now we can say that the first differences are stationary.

B) Perform Trace and Max Eigenvalue Tests for Cointegration

First we do the Trace test and then the Eigenvalue Test

y<-cbind(DPIdiff,PCEdiff)
y.CA <- ca.jo(y, ecdet="const", type="trace", K=2, spec="transitory", season=4)
summary(y.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 4.202666e-01 1.205877e-01 1.665335e-16
## 
## Values of teststatistic and critical values of test:
## 
##            test 10pct  5pct  1pct
## r <= 1 |  35.21  7.52  9.24 12.97
## r = 0  | 184.59 17.85 19.96 24.60
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               DPIdiff.l1  PCEdiff.l1   constant
## DPIdiff.l1  1.0000000000  1.00000000  1.0000000
## PCEdiff.l1 -0.9871040761  2.50569420  2.1485644
## constant   -0.0002429824 -0.05493107 -0.7093921
## 
## Weights W:
## (This is the loading matrix)
## 
##           DPIdiff.l1 PCEdiff.l1     constant
## DPIdiff.d -0.8848011 -0.1066346 8.042714e-18
## PCEdiff.d  0.4776245 -0.1085864 2.652189e-19
y.CA <- ca.jo(y, ecdet="const", type="eigen", K=2, spec="transitory", season=4)
summary(y.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , without linear trend and constant in cointegration 
## 
## Eigenvalues (lambda):
## [1] 4.202666e-01 1.205877e-01 1.665335e-16
## 
## Values of teststatistic and critical values of test:
## 
##            test 10pct  5pct  1pct
## r <= 1 |  35.21  7.52  9.24 12.97
## r = 0  | 149.38 13.75 15.67 20.20
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               DPIdiff.l1  PCEdiff.l1   constant
## DPIdiff.l1  1.0000000000  1.00000000  1.0000000
## PCEdiff.l1 -0.9871040761  2.50569420  2.1485644
## constant   -0.0002429824 -0.05493107 -0.7093921
## 
## Weights W:
## (This is the loading matrix)
## 
##           DPIdiff.l1 PCEdiff.l1     constant
## DPIdiff.d -0.8848011 -0.1066346 8.042714e-18
## PCEdiff.d  0.4776245 -0.1085864 2.652189e-19

The test statistic is greater than the critical value, so we reject the null. There is cointegration. It seems that there is just one cointegration relationship.

C. Test of Restricted Constant

lttest(y.CA, r=1)
## LR-test for no linear trend
## 
## H0: H*2(r<=1)
## H1: H2(r<=1)
## 
## Test statistic is distributed as chi-square
## with 1 degress of freedom
##         test statistic p-value
## LR test           0.04    0.84

Now we redo the trace test and the eigenvalue test

y.CA <- ca.jo(y, ecdet="trend", type="trace", K=2, spec="transitory", season=4)
summary(y.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: trace statistic , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1] 4.234132e-01 1.306315e-01 1.665335e-16
## 
## Values of teststatistic and critical values of test:
## 
##            test 10pct  5pct  1pct
## r <= 1 |  38.36 10.49 12.25 16.26
## r = 0  | 189.23 22.76 25.32 30.45
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               DPIdiff.l1   PCEdiff.l1   trend.l1
## DPIdiff.l1  1.000000e+00 1.000000e+00  1.0000000
## PCEdiff.l1 -9.586421e-01 2.633716e+00 -1.7927550
## trend.l1    6.897249e-06 9.428044e-05  0.0135177
## 
## Weights W:
## (This is the loading matrix)
## 
##           DPIdiff.l1 PCEdiff.l1      trend.l1
## DPIdiff.d -0.9155089 -0.1091323 -3.336542e-19
## PCEdiff.d  0.4684666 -0.1145687 -1.292973e-20
y.CA <- ca.jo(y, ecdet="trend", type="eigen", K=2, spec="transitory", season=4)
summary(y.CA)
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Test type: maximal eigenvalue statistic (lambda max) , with linear trend in cointegration 
## 
## Eigenvalues (lambda):
## [1] 4.234132e-01 1.306315e-01 1.665335e-16
## 
## Values of teststatistic and critical values of test:
## 
##            test 10pct  5pct  1pct
## r <= 1 |  38.36 10.49 12.25 16.26
## r = 0  | 150.87 16.85 18.96 23.65
## 
## Eigenvectors, normalised to first column:
## (These are the cointegration relations)
## 
##               DPIdiff.l1   PCEdiff.l1   trend.l1
## DPIdiff.l1  1.000000e+00 1.000000e+00  1.0000000
## PCEdiff.l1 -9.586421e-01 2.633716e+00 -1.7927550
## trend.l1    6.897249e-06 9.428044e-05  0.0135177
## 
## Weights W:
## (This is the loading matrix)
## 
##           DPIdiff.l1 PCEdiff.l1      trend.l1
## DPIdiff.d -0.9155089 -0.1091323 -3.336542e-19
## PCEdiff.d  0.4684666 -0.1145687 -1.292973e-20

So we can conclude that there is only 1 cointegrating relationship.

D. Estimate the Unrestricted VEC Model

unres.vec <- cajorls(y.CA, r=1)
unres.vec
## $rlm
## 
## Call:
## lm(formula = substitute(form1), data = data.mat)
## 
## Coefficients:
##              DPIdiff.d   PCEdiff.d 
## ect1         -9.155e-01   4.685e-01
## constant      1.415e-03  -8.719e-04
## sd1          -4.026e-05  -5.155e-05
## sd2          -3.293e-04   1.184e-03
## sd3          -1.626e-03  -1.235e-03
## DPIdiff.dl1  -1.094e-01  -2.273e-01
## PCEdiff.dl1  -3.196e-01  -4.122e-01
## 
## 
## $beta
##                     ect1
## DPIdiff.l1  1.000000e+00
## PCEdiff.l1 -9.586421e-01
## trend.l1    6.897249e-06

E. Hypothesis Testing

Perform the following three tests: 1) \(\beta_1=-1\) 2) \(\alpha_2=0\) 3) Joint hypothesis of 1 and 2

We first test the first one

rest.beta <- matrix(c(1,0,0,0,-1,0), c(3,2))
summary( blrtest(y.CA, H=rest.beta, r=1) )
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Estimation and testing under linear restrictions on beta 
## 
## The VECM has been estimated subject to: 
## beta=H*phi and/or alpha=A*psi
## 
##      [,1] [,2]
## [1,]    1    0
## [2,]    0   -1
## [3,]    0    0
## 
## Eigenvalues of restricted VAR (lambda):
## [1] 0.4203 0.1205
## 
## The value of the likelihood ratio test statistic:
## 1.49 distributed as chi square with 1 df.
## The p-value of the test statistic is: 0.22 
## 
## Eigenvectors, normalised to first column
## of the restricted VAR:
## 
##         [,1]   [,2]
## [1,]  1.0000 1.0000
## [2,] -0.9871 2.5051
## [3,]  0.0000 0.0000
## 
## Weights W of the restricted VAR:
## 
##              [,1]    [,2]
## DPIdiff.d -0.8848 -0.1067
## PCEdiff.d  0.4776 -0.1086

Now the second one

rest.mat <- matrix(c(1,0), c(2,1))
summary(alrtest(y.CA, A=rest.mat, r=1) )
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Estimation and testing under linear restrictions on beta 
## 
## The VECM has been estimated subject to: 
## beta=H*phi and/or alpha=A*psi
## 
##      [,1]
## [1,]    1
## [2,]    0
## 
## Eigenvalues of restricted VAR (lambda):
## [1] 0.378 0.000 0.000
## 
## The value of the likelihood ratio test statistic:
## 20.76 distributed as chi square with 1 df.
## The p-value of the test statistic is: 0 
## 
## Eigenvectors, normalised to first column
## of the restricted VAR:
## 
##                  [,1]
## RK.DPIdiff.l1  1.0000
## RK.PCEdiff.l1 -0.8044
## RK.trend.l1    0.0000
## 
## Weights W of the restricted VAR:
## 
##         [,1]
## [1,] -1.2023
## [2,]  0.0000
summary(ablrtest(y.CA, A=rest.mat, H=rest.beta, r=1) )
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Estimation and testing under linear restrictions on alpha and beta 
## 
## The VECM has been estimated subject to: 
## beta=H*phi and/or alpha=A*psi
## 
##      [,1] [,2]
## [1,]    1    0
## [2,]    0   -1
## [3,]    0    0
## 
## 
##      [,1]
## [1,]    1
## [2,]    0
## 
## Eigenvalues of restricted VAR (lambda):
## [1] 0.3703 0.0000
## 
## The value of the likelihood ratio test statistic:
## 24.13 distributed as chi square with 1 df.
## The p-value of the test statistic is: 0 
## 
## Eigenvectors, normalised to first column
## of the restricted VAR:
## 
##         [,1]
## [1,]  1.0000
## [2,] -0.8354
## [3,]  0.0000
## 
## Weights W of the restricted VAR:
## 
##         [,1]
## [1,] -1.1797
## [2,]  0.0000

Now the third joint one

summary( ablrtest(y.CA, A=rest.mat, H=rest.beta, r=1) )
## 
## ###################### 
## # Johansen-Procedure # 
## ###################### 
## 
## Estimation and testing under linear restrictions on alpha and beta 
## 
## The VECM has been estimated subject to: 
## beta=H*phi and/or alpha=A*psi
## 
##      [,1] [,2]
## [1,]    1    0
## [2,]    0   -1
## [3,]    0    0
## 
## 
##      [,1]
## [1,]    1
## [2,]    0
## 
## Eigenvalues of restricted VAR (lambda):
## [1] 0.3703 0.0000
## 
## The value of the likelihood ratio test statistic:
## 24.13 distributed as chi square with 1 df.
## The p-value of the test statistic is: 0 
## 
## Eigenvectors, normalised to first column
## of the restricted VAR:
## 
##         [,1]
## [1,]  1.0000
## [2,] -0.8354
## [3,]  0.0000
## 
## Weights W of the restricted VAR:
## 
##         [,1]
## [1,] -1.1797
## [2,]  0.0000

We can say that the second test is true

F) Forecasting

We convert the VEC model into a VAR model and forecast for 8 quarters.

y.VAR <- vec2var(y.CA, r=1)
y.VAR.for <- predict(y.VAR, n.ahead=8)
plot(y.VAR.for)