library(tseries) 
library(urca) 

data = read.table("672_2014_texas.txt", header = T) 
a = ts(data[,1], start = 1990, frequency = 12)
d = ts(data[,2], start = 1990, frequency = 12)
h = ts(data[,3], start = 1990, frequency = 12)
s = ts(data[,4], start = 1990, frequency = 12)


plot(d,type="l", ylab = "log price")
lines(h, lty = 2)

# bivariate VAR(2) for d(allas) and h(ouston)
ob = length(d)
co = rep(1, ob)
d.1 = c(NA, d[1:ob-1])
dd = d - d.1
dd.1 = c(NA, dd[1:ob-1])
h.1 = c(NA, h[1:ob-1])
dh = h - h.1
dh.1 = c(NA, dh[1:ob-1])

# obtain the residual
mo = lm(dd~dd.1+dh.1)
u1 = mo$res
mo = lm(dh~dd.1+dh.1)
u2 = mo$res
u = cbind(u1, u2)
uu = crossprod(u,u)

mo = lm(d.1~dd.1+dh.1)
w1 = mo$res
mo = lm(h.1~dd.1+dh.1)
w2 = mo$res
w = cbind(w1, w2)
ww = crossprod(w,w)

wu = crossprod(w,u)

# obtain the eigenvalues 
lam = eigen(wu%*%solve(uu)%*%t(wu)%*%solve(ww))
lam[1]


# obtain the LR tests (Trace Test)
LR1 = -(168-2)*(log(1-0.001074729)+log(1-0.244975013))
LR2 = -(168-2)*(log(1-0.001074729))
LR1
LR2


# built-in command ca.jo
joh = ca.jo(data[,c(2,3)],ecdet = "const", type="trace", K=2, spec="transitory")
summary(joh)