"t2chart" <- function(da){
# The data matrix is "da".
# Written by Ruey S. Tsay on April 10, 2008.
 nr = dim(da)[1]
 nc = dim(da)[2]
 cm=matrix(colMeans(da),1,nc)
 o=matrix(rep(1,nr),nr,1)
 dev=da-kronecker(o,cm)
 dev=as.matrix(dev)      
 s=t(dev)%*%dev/(nr-1)             
 si=solve(s)
 t2=diag(dev%*%si%*%t(dev)) 
 ul1=qchisq(0.95,nc)
 ul2=qchisq(0.99,nc)
 yl=max(t2,ul1,ul2)
 plot(t2,type='b',ylim=c(0,yl))
 abline(h=c(ul1),lty=2)
 abline(h=c(ul2),lty=2)
 title(main='The limits are 95% and 99% quantiles')
 list(Tsq=t2)
}