"growth" <- function(da,nv,tp,q){
# This program tests the q-order growth curve model.
# The da is the data matrix.
# nv is the vector of sample sizes of the groups.
# q is the polynomial order.
# tp: time vector
# Written by Ruey S. Tsay on April 24, 2008
growth=NULL
N = dim(da)[1]
p = dim(da)[2]
g=length(nv)
Sp=matrix(0,p,p)
idx=0
xbar=NULL
for (i in 1:g){
ni=nv[i]
x=da[(idx+1):(idx+ni),]
ave=apply(x,2,mean)
xbar=cbind(xbar,ave)
S=cov(x)
Sp=(ni-1)*S+Sp
idx=idx+ni
}
W=Sp
Sp=Sp/(N-g)
Si=solve(Sp)
beta=matrix(0,(q+1),g)
semtx=matrix(0,(q+1),g)
B=matrix(1,p,1)
for (j in 1:q){
B=cbind(B,tp^j)
}
idx=0
Wq=matrix(0,p,p)
k=(N-g)*(N-g-1)/((N-g-p+q)*(N-g-p+q+1))
for (i in 1:g){
ni=nv[i]
x=da[(idx+1):(idx+ni),]
S=t(B)%*%Si%*%B
Sinv=solve(S)
d1=xbar[,i]
xy=t(B)%*%Si%*%d1
bhat=Sinv%*%xy
beta[,i]=bhat
semtx[,i]=sqrt((k/ni)*diag(Sinv))
fit=t(B%*%bhat)
err=x-kronecker(fit,matrix(1,ni,1))
err=as.matrix(err)
Wq=t(err)%*%err+Wq
idx=idx+ni
}
Lambda=det(W)/det(Wq)
Tst=-(N-(p-q+g)/2)*log(Lambda)
pv=1-pchisq(Tst,(p-q-1)*g)
growth=cbind(growth,c(Tst,pv))
row.names(growth)=c("LR-stat","p.value")
print("Growth curve model")
print("Order: ")
print(q)
print("Beta-hat: ")
print(beta,digits=4)
print("Standard errors: ")
print(semtx,digits=4)
print("W")
print(W,digits=4)
print("Wq")
print(Wq,digits=4)
print("Lambda:")
print(Lambda)
print("Test result:")
print(growth,digits=4)
}