"coint" <- function(x,lags=0,trend=1){
# Perform Johansen's co-integration test.
# trend = 0 means no constant
#       = 1 means with a constant
#       = 2 means with time trend.
#
T=nrow(x)
k=ncol(x)
y=as.matrix(x)
# setup the differenced data
dy=cbind(rep(0,k),t(diff(y)))
dy=t(dy)
#for (i in 1:k){
#dy[,i]=c(0,diff(y[,i]))
#}
ist = lags+1
ut=dy[(ist+1):T,]
vt=y[ist:(T-1),]
# Run regression to remove the impact of stationary part or trend.
nobe=T-ist
idx=0
if( (trend > 0) || (lags > 0)){
idim=trend + lags*k
##print(c(idim,nobe))
bx=matrix(0,nobe,idim)
if(trend >0){
bx[,1]=rep(1,nobe)
idx=1
if(trend > 1){
idx=2
bx[,2]=c(1:nobe)
}
}
if(lags > 0){
for (j in 1:lags){
bx[,(idx+1):(idx+k)]=dy[(ist-j+1):(T-j),]
idx=idx+k
}
}
xpx=t(bx)%*%bx
##print(xpx)
xpxinv=solve(xpx)
y1=dy[(ist+1):T,]
xpy=t(bx)%*%y1
b1=xpxinv%*%xpy
ut=y1-bx%*%b1
zm1=y[ist:(T-1),]
xpy=t(bx)%*%zm1
b2=xpxinv%*%xpy
vt=zm1-bx%*%b2
}
# Perform eigenvalues &eigenvector analysis
S00=t(ut)%*%ut/(T-lags)
S01=t(ut)%*%vt/(T-lags)
S11=t(vt)%*%vt/(T-lags)
print("S00-mtx")
print(S00)
print("S01-mtx")
print(S01)
print("S11-mtx")
print(S11)
S00inv=solve(S00)
S11inv=solve(S11)
A=S11inv%*%t(S01)
A1=S00inv%*%S01
AA=A%*%A1
print("The U-mtx")
print(AA)
mm=eigen(AA)
ev=mm$values
#print(ev)
evto=mm$vectors
lktr=rep(0,k)
lkmx=rep(0,k)
for (j in 1:k){
lkmx[j]=-nobe*log(1-ev[j])
}
#print(lkmx)
lktr[k]=lkmx[k]
for (j in 1:(k-1)){
lktr[(k-j)]=lktr[(k-j+1)]+lkmx[(k-j)]
}
#print(lktr)
rnk=c(1:k)-1
resu=cbind(rnk,ev,lktr,lkmx)
print("Rank, Eigenvalue, trace, & max-eig")
print(resu,digits=3)
#print(paste("selected order: aic =",aicor))
coint<-list(trace=lktr,maxeig=lkmx,eigvalue=ev,eigvector=evto)

}