#************************************************************ 
# Vector Error Correction Models GARCH  
# install.packages("fGarch")
# install.packages("AER")
#------------------------------------------------------------


library("dynlm")
library("tseries")
library("urca")  #cointegration test
library("vars")
library("quantmod")  # To get data for symbols
library("fUnitRoots") # Unit tes

getSymbols("MCOILWTICO", src="FRED")
oil <- MCOILWTICO
getSymbols("GASREGCOVM", src="FRED")
gas <- GASREGCOVM

dat1 <- merge(oil,gas)
dat1 <- window(dat1, start="2001-01-01",end="2016-12-31")
colnames(dat1) <- c("oil","gas")
par(mfrow=c(1,1))
plot(dat1, type='l', 
     main="Log of Crude Oil Prices vs. Regular Conventional Gas Price", 
     col=c("blue","red"), ylim=c(-1, 5))

# Unit Root Test
utest <-function(y) { 
  df20 <- adf.test(y, alternative = "stationary", k = 0)
  r0 <- c(df20$statistic,df20$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 <- rbind(r0,r1,r2,r3)
  colnames(table1) <- c("stat","p-value")
  rownames(table1) <- c("DF","ADF","PP","KPSS")
  table1
  }
utest(dat1$gas)
utest(dat1$oil)
n <- nrow(dat1$gas)
dg <- diff(dat1$gas) 
dg <- dg[2:n]
utest(dg)
do <- diff(dat1$oil) 
do <- do[2:n]
utest(do)

# Cointegration Test
# Engle and Granger Cointegreation Test 
reg1 =lm (gas~oil,data=dat1)
adfTest(reg1$residuals, type="nc")
utest(reg1$residuals)

## And now lets use adf test on spread (which is actually residuals of regression above step)

# Cointegration Test

y.mat <- data.frame(dat1$gas, dat1$oil)

#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 = 'oil',
                response = 'gas', 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 
lr  <- lm(gas ~ oil, data=dat1)
summary(lr)