"mq_res" <- function(x,lag,p){
# Compute multivariate Ljung-Box test statistics 
# for residuals of a VAR(p)
#
if(!is.matrix(x))x=as.matrix(x)
nr=nrow(x)
nc=ncol(x)
g0=var(x)
ginv=solve(g0)
qm=0.0
print("m,       Q(m) and  p-value:")
df = 0
for (i in 1:lag){
  x1=x[(i+1):nr,]
  x2=x[1:(nr-i),]
  g = cov(x1,x2)
  g = g*(nr-i-1)/(nr-1)
  h=t(g)%*%ginv%*%g%*%ginv
  qm=qm+nr*nr*sum(diag(h))/(nr-i)
  df=df+nc*nc
  dfp = df-p*nc*nc
  if (i > p) {
     pv=1-pchisq(qm,dfp)
     print(c(i,qm,pv))
               }
}
}