"mmlr" <- function(z,y,constant=T){
# This program performs multivariate multiple linear regression analysis.
# z: design matrix 
# constant: switch for the constant term of the regression model
# y: dependent variables
z=as.matrix(z)
n=nrow(z)
nx=ncol(z)
zc=z
if (constant) zc=cbind(rep(1,n),z)
p=ncol(zc)
y=as.matrix(y)
ztz=t(zc)%*%zc
zty=t(zc)%*%y
ztzinv=solve(ztz)
beta=ztzinv%*%zty
Beta=beta
res=y-zc%*%beta
sig=t(res)%*%res/(n-p)
co=kronecker(sig,ztzinv)
print("LSE of parameters")
print("  est   s.d.   t-ratio   prob")
par=beta[,1]
deg=n-p
p1=ncol(y)
if (p1 > 1){
for (i in 2:p1){
par=c(par,beta[,i])
}
}
iend=nrow(beta)*ncol(beta)
tmp=matrix(0,iend,4)
for (i in 1:iend){
sd=sqrt(co[i,i])
tt=par[i]/sd
pr=2*(1-pt(abs(tt),deg))
tmp[i,]=c(par[i],sd,tt,pr)
}
print(tmp,digits=3)
# testing
if (nx > 1){
sig=sig*(n-p)/n
det1=det(sig)
ztmp=z
print("Test the regressor")
for (j in 1:nx){
zc=ztmp[,2:nx]
if(constant) zc=cbind(rep(1,n),zc)
ztz=t(zc)%*%zc
zty=t(zc)%*%y
ztzinv=solve(ztz)
beta=ztzinv%*%zty
res=y-zc%*%beta
sig1=t(res)%*%res/n
det2=det(sig1)
tst=log(det1)-log(det2)
tst=-(n-p-0.5*ncol(y))*tst
deg=ncol(y)
pr=1-pchisq(tst,deg)
tmp=c(j,tst,pr)
print("regr    test-stat   p-value")
print(tmp,digits=4)
ztmp=cbind(ztmp[,2:nx],ztmp[,1])
}

}
mmlr <- list(beta=Beta,residuals=res,sigma=sig)
}