date()


#########################################################################
#Generating data                                                         
#The k samples are stored in one vector y=(y_11,...,y_1n1,...,y_k1,...,  
#y_knk)' and the covariates are stored in one vector x similarly.        
#The vector SIZES contains the sizes for each sample. k is the number of 
#samples.                                                              
#########################################################################
r1<-matrix(integer(15000),5000,3)
r2<-matrix(integer(15000),5000,3)
r3<-matrix(integer(15000),5000,3)

sizes<-c(16,24,12,20,18,15)
n <- sum(sizes)
k <- length(sizes)

ja<-matrix(c(1,1),2,1)%*%matrix(c(1,1),1,2)/2
jb<-matrix(c(1,1,1),3,1)%*%matrix(c(1,1,1),1,3)/3
pa<-diag(c(1,1))-ja
pb<-diag(c(1,1,1))-jb
h1.mat<-kronecker(pa,jb)
h1<-h1.mat[1,1]
h2.mat<-kronecker(ja,pb)
h2<-h2.mat[1,1]
h3.mat<-kronecker(pa,pb)
h3<-h3.mat[1,1]

d<-diag(1/(sizes-1))
m<-solve(diag(rep(1,6))+2*d)
set.seed(368)

group <- rep(1:k,sizes)
#group

allsizes <- rep(sizes,sizes)
#allsizes


an <- 9/allsizes^0.26
#an


for(i in 1:5000)
{
index<-sample(1:length(tbbm),n)
x<-htvel[index]
y<-tbbm[index]



#########################################################################
#Calculattion of the test statistic;
#########################################################################


 midrank <- rank(y)  #The midranks of the responses
 
 x.mat <- t(matrix(x,n,n))
 k.mat <- (x.mat-x)/an
 k.mat <- replace(k.mat, abs(k.mat)>1,1)
 k.mat <- (15/16)*(1-k.mat^2)^2
 #k.mat

 temp1 <- rowsum(k.mat,group)
 temp <- replace(temp1,temp1==0,1)

 tmp <- apply(temp,2,rep,times=sizes)
 w.mat <- k.mat/tmp
 #w.mat

 wr.mat<- (1/n)*w.mat*midrank

 int.pre <- apply(wr.mat,1,mean)
 int <- rowsum(int.pre,group)
# int

########################################################################
#calculate \hat \delta_{1,i}^2 and \hat \delta_{1,i_1,i_2},i=1,...,k
########################################################################

 tp1.mat <- rowsum(wr.mat,group)
# v1.mat<-(n-1)/n*var(t(tp1.mat))

#########################################################################
#calculate \hat \delta_{2,i}^2, i=1,...,k
#########################################################################

 wrr.mat <- wr.mat * midrank
 tp4.mat <- rowsum(wrr.mat,group)
 uij.vec <- (apply(w.mat,1,mean))^2
 tp5.mat <- t(tp4.mat)*uij.vec
 tp6.mat <- rowsum(tp5.mat,group)
 v21.vec <- diag(tp6.mat)
 wrsq.mat <- (rowsum(wr.mat,group))^2
 tp7.mat <- n*rowsum(t(wrsq.mat)*uij.vec,group)
 v22.vec <- diag(tp7.mat)
 v2.mat <- diag(v21.vec-v22.vec)
# v2.mat

#########################################################################
#calculate the estimate of the covariance matrix
#########################################################################

# v.mat <- v1.mat + v2.mat
# v.mat

#########################################################################
#calculate the statistic
#########################################################################
v.mat<-sizes/n*v2.mat

qh<-t(int)%*%h1.mat%*%int
 f1<-h1^2*((sum(diag(v2.mat)))^2-2*sum(diag(d%*%m%*%v2.mat%*%v2.mat)))/ sum(
     diag(h1.mat%*%v2.mat%*%h1.mat%*%v2.mat))
 f2<-((sum(diag(solve(d)%*%v.mat)))^2-2*sum(diag(solve(d)%*%m%*%v.mat%*%
     v.mat)))/sum(diag(solve(d)%*%v.mat%*%v.mat))
 testst<- n/(h1*sum(diag(v2.mat)))*qh
 qt<-qf(c(0.9,0.95,0.99),f1,f2)
 r1[i,]<-(testst>=qt)

qh<-t(int)%*%h2.mat%*%int
 f1<-h2^2*((sum(diag(v2.mat)))^2-2*sum(diag(d%*%m%*%v2.mat%*%v2.mat)))/ sum(
     diag(h2.mat%*%v2.mat%*%h2.mat%*%v2.mat))
 f2<-((sum(diag(solve(d)%*%v.mat)))^2-2*sum(diag(solve(d)%*%m%*%v.mat%*%
     v.mat)))/sum(diag(solve(d)%*%v.mat%*%v.mat))
 testst<- n/(h2*sum(diag(v2.mat)))*qh
 qt<-qf(c(0.9,0.95,0.99),f1,f2)
 r2[i,]<-(testst>=qt)

qh<-t(int)%*%h3.mat%*%int
 f1<-h3^2*((sum(diag(v2.mat)))^2-2*sum(diag(d%*%m%*%v2.mat%*%v2.mat)))/ sum(
     diag(h3.mat%*%v2.mat%*%h3.mat%*%v2.mat))
 f2<-((sum(diag(solve(d)%*%v.mat)))^2-2*sum(diag(solve(d)%*%m%*%v.mat%*%
     v.mat)))/sum(diag(solve(d)%*%v.mat%*%v.mat))
 testst<- n/(h3*sum(diag(v2.mat)))*qh
 qt<-qf(c(0.9,0.95,0.99),f1,f2)
 r3[i,]<-(testst>=qt)
}
print("Treatment effect:")
apply(r1,2,mean)
print("Age effect:")
apply(r2,2,mean)
print("Interaction:")
apply(r3,2,mean)