########## Library + Data
##########

library(urca)
library(vars)


data<- read.table("C:/Users/Alexandre Plaquet/Desktop/Mémoire/Analysis/R/Master Thesis/data.txt", header = TRUE, sep=",")
attach(data)


# Shapiro Test for Normality

shapiro.test(lnoil)
shapiro.test(lngas)
shapiro.test(lncoal)

########## 5.2. Stationarity test ##########
##########

# 5.2.2. ADF


VARselect(lnoil, type="trend")
VARselect(lngas, type="trend")
VARselect(lncoal, type="const")

lnoil.adf=ur.df(lnoil,type='trend',selectlags="AIC")
summary(lnoil.adf)
lngas.adf=ur.df(lngas,type='trend',selectlags="AIC")
summary(lngas.adf)
lncoal.adf=ur.df(lncoal,type='drift',selectlags="AIC")
summary(lncoal.adf)

# 5.2.3. Phillips-Perron test

lnoil.pp=ur.pp(lnoil,type="Z-tau",model="trend",use.lag=4)
summary(lnoil.pp)
lngas.pp=ur.pp(lngas,type="Z-tau",model="trend",use.lag=4)
summary(lngas.pp)
lncoal.pp=ur.pp(lncoal,type="Z-tau",model="trend",use.lag=4)
summary(lngas.pp)

# 5.2.4. Zivot-Andrews test

lnoil.za=ur.za(lnoil,model="both",lag=1)
summary(lnoil.za)
AIC(eval(attributes(lnoil.za)$testreg))
BIC(eval(attributes(lnoil.za)$testreg))

lngas.za=ur.za(lngas,model="intercept",lag=3)
summary(lngas.za)
AIC(eval(attributes(lngas.za)$testreg))
BIC(eval(attributes(lngas.za)$testreg))

lncoal.za=ur.za(lncoal,model="intercept",lag=3)
summary(lncoal.za)
AIC(eval(attributes(lncoal.za)$testreg))
BIC(eval(attributes(lncoal.za)$testreg))

# 5.2.5. KPSS test

lnoil.kpss=ur.kpss(lnoil,type ="tau",use.lag=1)
summary(lnoil.kpss)
lngas.kpss=ur.kpss(lngas,type ="tau",use.lag=1)
summary(lngas.kpss)
lncoal.kpss=ur.kpss(lncoal,type ="mu",use.lag=1)
summary(lncoal.kpss)


########## 5.3. Cointegration tests ##########
##########

# 5.3.1. Order of integration (ADF)

dlnoil=diff(lnoil)
dlngas=diff(lngas)
dlncoal=diff(lncoal)

dlnoil.ct=ur.df(dlnoil,type='none',selectlags="AIC")
summary(dlnoil.ct)

dlngas.ct=ur.df(dlngas,type='none',selectlags="AIC")
summary(dlngas.ct)

dlncoal.ct=ur.df(dlncoal,type='none',selectlags="AIC")
summary(dlncoal.ct)


########## 5.3.2. Engle-Granger method
##########

### 5.3.2.1. Granger Causality

# Causality tests

# l=1

grangertest(dlnoil~dlngas, order=1)
grangertest(dlnoil~dlncoal, order=1)
grangertest(dlngas~dlnoil, order=1)
grangertest(dlngas~dlncoal, order=1)
grangertest(dlncoal~dlnoil, order=1)
grangertest(dlncoal~dlngas, order=1)

# l=2

grangertest(dlnoil~dlngas, order=2)
grangertest(dlnoil~dlncoal, order=2)
grangertest(dlngas~dlnoil, order=2)
grangertest(dlngas~dlncoal, order=2)
grangertest(dlncoal~dlnoil, order=2)
grangertest(dlncoal~dlngas, order=2)

# l=3

grangertest(dlnoil~dlngas, order=3)
grangertest(dlnoil~dlncoal, order=3)
grangertest(dlngas~dlnoil, order=3)
grangertest(dlngas~dlncoal, order=3)
grangertest(dlncoal~dlnoil, order=3)
grangertest(dlncoal~dlngas, order=3)

# l=4

grangertest(dlnoil~dlngas, order=4)
grangertest(dlnoil~dlncoal, order=4)
grangertest(dlngas~dlnoil, order=4)
grangertest(dlngas~dlncoal, order=4)
grangertest(dlncoal~dlnoil, order=4)
grangertest(dlncoal~dlngas, order=4)

# l=5

grangertest(dlnoil~dlngas, order=5)
grangertest(dlnoil~dlncoal, order=5)
grangertest(dlngas~dlnoil, order=5)
grangertest(dlngas~dlncoal, order=5)
grangertest(dlncoal~dlnoil, order=5)
grangertest(dlncoal~dlngas, order=5)

### 5.3.2.2. Cointegration tests

lnoil=ts(lnoil,start=c(1980,1), end=c(2014,12), frequency = 12)
lngas=ts(lngas,start=c(1980,1), end=c(2014,12), frequency = 12)
lncoal=ts(lncoal,start=c(1980,1), end=c(2014,12), frequency = 12)

oc.eq=summary(lm(lnoil~lncoal))
oc.eq

go.eq=summary(lm(lngas~lnoil))
go.eq

gc.eq=summary(lm(lngas~lncoal))
gc.eq

### Extraction and storage of residuals

error.oc.eq=ts(resid(oc.eq),start=c(1980,2), end=c(2014,12), frequency = 12)
error.go.eq=ts(resid(go.eq),start=c(1980,2), end=c(2014,12), frequency = 12)
error.gc.eq=ts(resid(gc.eq),start=c(1980,2), end=c(2014,12), frequency = 12)

### ADF test on residuals

ci.oc=ur.df(error.oc.eq, type='none',selectlags="AIC")
summary(ci.oc)

ci.go=ur.df(error.go.eq, type='none',selectlags="AIC")
summary(ci.go)

ci.gc=ur.df(error.gc.eq, type='none',selectlags="AIC")
summary(ci.gc)

### 5.3.2.3. ECM

dlnoil=diff(lnoil)
dlngas=diff(lngas)
dlncoal=diff(lncoal)

lerror.oc.eq=lag(error.oc.eq)
lerror.go.eq=lag(error.go.eq)
lerror.gc.eq=lag(error.gc.eq)

ecm.oc.eq=summary(lm(dlnoil~dlncoal+lerror.oc.eq))
ecm.oc.eq
ecm.go.eq=summary(lm(dlngas~dlnoil+lerror.go.eq))
ecm.go.eq
ecm.gc.eq=summary(lm(dlngas~dlncoal+lerror.gc.eq))
ecm.gc.eq


### 5.3.2.4. Diagnostic error tests

library(tseries)
library(FinTS)

jarque.bera.test(lag(error.oc.eq))
ArchTest(lag(error.oc.eq),lags=3)
Box.test(lag(error.oc.eq),lag=1,type="Ljung-Box")

jarque.bera.test(lag(lerror.go.eq))
ArchTest(lag(lerror.go.eq),lags=3)
Box.test(lag(lerror.go.eq),lag=1,type="Ljung-Box")

jarque.bera.test(lag(error.gc.eq))
ArchTest(lag(error.gc.eq),lags=3)
Box.test(lag(error.gc.eq),lag=1,type="Ljung-Box")

########## 5.3.3. Johasen's procedure

VARselect(data.frame(lnoil,lngas,lncoal))

### 5.3.3.1. Trace Test

summary(ca.jo(data.frame(lnoil,lngas,lncoal), type="trace", ecdet="const",K=2))

### 5.3.3.2. Maximum eigenvalue test

summary(ca.jo(data.frame(lnoil,lngas,lncoal), type="eigen", ecdet="const",K=2))

### 5.3.3.3. VECM

library(vars)

vecm.eigen<-ca.jo(data.frame(lnoil,lngas,lncoal), type="eigen",ecdet="const", K=2)

vecm.model<-cajorls(vecm.eigen,r=1)

vecm.model

## Calculation of t-values for alpha and beta
##
alpha <- coef(vecm.model$rlm)[1, ]
names(alpha) <- c("lnoil", "lngas", "lncoal")
alpha
beta <- vecm.model$beta
beta
resids <- resid(vecm.model$rlm)
N <- nrow(resids)
sigma <- crossprod(resids) / N
## t-stats for alpha (calculated by hand)
alpha.se <- sqrt(solve(crossprod(cbind(vecm.eigen@ZK %*% beta, vecm.eigen@Z1)))[1, 1] * diag(sigma))
names(alpha.se) <-  c("lnoil", "lngas", "lncoal")
alpha.t <- alpha / alpha.se
alpha.t
## Differ slightly from coef(summary(vecm.r1$rlm))
## due to degrees of freedom adjustment 
coef(summary(vecm.model$rlm))
## t-stats for beta
beta.se <- sqrt(diag(kronecker(solve(crossprod(vecm.eigen@RK[, -1])),
                               solve(t(alpha) %*% solve(sigma) %*% alpha))))
beta.t <- c(NA, beta[-1] / beta.se)
names(beta.t) <- rownames(vecm.model$beta)
beta.t


### 5.3.3.4 Structural Shift

summary(cajolst(data.frame(lnoil,lngas,lncoal),trend=FALSE, K=2))

vecm.break<-ca.jo(data.frame(lnoil,lngas,lncoal), type="eigen",ecdet="const", K=2)

vecm.model2<-cajorls(vecm.break,r=2)

vecm.model2

## Calculation of t-values for alpha and beta
##
alpha <- coef(vecm.model2$rlm)[1, ]
names(alpha) <- c("lnoil", "lngas", "lncoal")
alpha
beta <- vecm.model2$beta
beta
resids <- resid(vecm.model2$rlm)
N <- nrow(resids)
sigma <- crossprod(resids) / N
## t-stats for alpha (calculated by hand)
alpha.se <- sqrt(solve(crossprod(cbind(vecm.break@ZK %*% beta, vecm.break@Z1)))[1, 1] * diag(sigma))
names(alpha.se) <-  c("lnoil", "lngas", "lncoal")
alpha.t <- alpha / alpha.se
alpha.t
## Differ slightly from coef(summary(vecm.r1$rlm))
## due to degrees of freedom adjustment 
coef(summary(vecm.model2$rlm))
## t-stats for beta
beta.se <- sqrt(diag(kronecker(solve(crossprod(vecm.break@RK[, -1])),
                               solve(t(alpha) %*% solve(sigma) %*% alpha))))
beta.t <- c(NA, beta[-1] / beta.se)
names(beta.t) <- rownames(vecm.model2$beta)
beta.t

### 5.3.3.5. Diagnostic tests

vecm.level=vec2var(vecm.eigen,r=1)

normality.test(vecm.level)
arch.test(vecm.level)