#***************************************************************
# Cointegration and Vector Error Correction Model
# install.packages("dynlm")
# install.packages("fUnitRoots")
# install.packages("dynlm")
#---------------------------------------------------------------
library("dynlm")
library("tseries")
library("urca")        #cointegration test
library("vars")
library("quantmod")    # To get data for symbols
library("fUnitRoots")  # Unit tes

# Data Generation for 2 variable cointegration Example

set.seed(12345678)
e1 <- rnorm(250, mean = 0, sd = 0.5)
e2 <- rnorm(250, mean = 0, sd = 0.5)
u.ar3 <- arima.sim(model =
         list(ar = c(0.6, -0.2, 0.1)), n = 250,
         innov = e1)
y2 <- cumsum(e2)
y1 <- u.ar3 + 0.75*y2
ymax <- max(c(y1, y2))
ymin <- min(c(y1, y2))
layout(matrix(1:2, nrow = 2, ncol = 1))
plot(y1, xlab = "", ylab = "", ylim =
        c(ymin, ymax), main =
        "Cointegrated System")
lines(y2, col = "green")
plot(u.ar3, ylab = "", xlab = "", main =
      "Cointegrating Residuals")
abline(h = 0, col = "red")

 y <- y1
 df20 <- adf.test(y, alternative = "stationary", k = 0)
 r0 <- c(df20$statistic,df20$p.value)
 
 df22 <- adf.test(y, alternative = "stationary", k = 2)
 r01 <- c(df22$statistic,df22$p.value)
 
 df24 <- adf.test(y, alternative = "stationary", k = 4)
 r1 <- c(df24$statistic,df24$p.value)
 
 pp2 <- pp.test(y,lshort=TRUE,type="Z(t_alpha)")
 r2 <- c(pp2$statistic,pp2$p.value)
 
 kp2 <- kpss.test(y,lshort=TRUE,null = c("Level", "Trend"))
 r3 <- c(kp2$statistic,kp2$p.value)
 
 table1 <- cbind(r0,r01,r1,r2,r3)
 table1
 
 # Cointegration Test
 # Engle and Granger Cointegreation Test 
 reg1 =lm (y1~y2)
 adfTest(reg1$residuals, type="nc")
         
 ## And now lets use adf test on spread (which is actually residuals of regression above step)
 
# Cointegration Test

 y.mat <- data.frame(y1, y2)

 #Phillips & Ouliaris Cointegration Test
 co.Pu <- ca.po(y.mat, demean = "none", type = "Pu")
 co.Pz <- ca.po(y.mat, demean = "none", type = "Pz")
 summary(co.Pu)
 summary(co.Pz)
 
 #Johansen Procedure for VAR
 vecm1 <- ca.jo(y.mat, type = "eigen", spec = "transitory")
 summary(vecm1)
 vecm2 <- ca.jo(y.mat, type = "trace", spec = "transitory")
 summary(vecm2)
 vecm.r2 <- cajorls(vecm1, r = 1) 
 vecm.r2
 
 vecm.level <- vec2var(vecm1, r = 1)
 vecm.level
 vecm.pred <- predict(vecm.level,n.ahead = 10)
 fanchart(vecm.pred)
 vecm.irf <- irf(vecm.level, impulse = 'y2',
                 response = 'y1', boot = FALSE)
 plot(vecm.irf)
 vecm.fevd <- fevd(vecm.level) #Forecast Error Variance Decomposition
 plot(vecm.fevd)
 
 vecm.norm <- normality.test(vecm.level)
 vecm.arch <- arch.test(vecm.level)
 vecm.serial <- serial.test(vecm.level)
 vecm.norm 
 vecm.arch
 vecm.serial 
 
 # ECM model Estimation  
 y.mat <- as.ts(y.mat)
 lr  <- lm(y1 ~ y2, data=y.mat)
summary(lr)
err <- (resid(lr))
n1 <- length(err)
err <- ts(err[1:n1-1])
dy1 <- diff(y.mat[,1])
dy2 <- diff(y.mat[,2])
ecmdat <- cbind(dy1, dy2, err)
ecmdat <- ts(ecmdat)

ecmdat <- na.omit(ecmdat)
colnames(ecmdat) <- c("dy1","dy2","ect")
ecm <- dynlm(dy1 ~ L(ect, 1) + L(dy1, 1) + L(dy2, 1) , data = ecmdat)
print(summary(ecm))