library(MASS)
 
 rm(list=ls(all=TRUE))
 
 p1<-2
 p2<-2
 m<-12
 
 NO<-c(50,250,50,100,200,150,50,150,100,150,100,50)
 N<-sum(NO)
 
 no<-c(1,5,1,2,4,3,1,3,2,3,2,1)
 n<-sum(no)
 
 to<-3*no
 sumt<-sum(to)
 
 I<-m+1
 s<-rep(0,I)
 S<-rep(0,I)
 st<-rep(0,I)
 st<-rep(0,I)
 s[1]<-1
 S[1]<-1
 st[1]<-1
 
 for(i in 2:I){
   S[i]<-S[i-1]+NO[i-1]
   s[i]<-s[i-1]+no[i-1]
   st[i]<-st[i-1]+to[i-1] 
 }
 s
 S
 st
 
 MC<-5000
 
 beta0<-100
 beta1<-rep(0.1,p1) 
 beta2<-rep(2,p2)
 vare0<-100
 varv0<-16
 vareta0<-diag(25, p2)
 mux0<-rep(194,p2)
 varx0<-diag(2737,p2) 
 
 runif(1)


 seed<-.Random.seed[1]
 set.seed(seed, kind ="default")
 
 mse.p.j<-matrix(0,MC,m)  
 mse.wp.j<-matrix(0,MC,m)
 mse1.p.j<-matrix(0,MC,m)  
 mse1.wp.j<-matrix(0,MC,m)
 mse2.p.j<-matrix(0,MC,m)  
 mse2.wp.j<-matrix(0,MC,m)
 
 
 MSE.p<-matrix(0,MC,m) 
 G1.p<-matrix(0,MC,m)
 MSE1.p<-matrix(0,MC,m)
 MSE2.p<-matrix(0,MC,m)
 MSEc.p<-matrix(0,MC,m)
 
 Bias.p<-matrix(0,MC,m)   
 Bias.mse.p<-matrix(0,MC,m)   
 Bias.mse.wp<-matrix(0,MC,m)   
 
 B0.p<-rep(0,MC)
 B1.p<-matrix(NA,MC,p1)
 B2.p<-matrix(NA,MC,p2)
 VARv.p<-rep(0,MC)
 VARe.p<-rep(0,MC)
 VAReta.p<-matrix(NA,MC,p2)
 VARx.p<-array(NA,c(p2,p2,MC))
 MUx.p<-matrix(NA,MC,p2)
 
 #VAReta.p.j<-matrix(0,MC,m)
 #VARe.p.j<-matrix(0,MC,m)
 #VARv.p.j<-matrix(0,MC,m)
 #VARx.p.j<--matrix(0,MC,m)
 #MUx.p.j<--matrix(0,MC,m)
 #B0new=matrix(0,MC,m)
 #B1new=matrix(0,MC,m)
 #B2new=matrix(0,MC,m)
 
 G.hat.p.j=matrix(0,MC,m)
 G1.hat.p.j=matrix(0,MC,m)
 #------------------------------------------------------
 #First generate once the free-error covarite w_ij: 
 W<-matrix(0,p1,N)
 W.bar<-matrix(0,p1,m)
 for(i in 1:m){
   A<-S[i]
   B<-S[i+1]-1
   for(j in 1:NO[i]){
     W[,A+j-1]<-rnorm(p1,1,1) 
   }
   W.bar[1,i]<-mean(W[1,A:B])
   W.bar[2,i]<-mean(W[2,A:B])
 }
 #------------------------------------------------------
 for(r in 1:MC){    #start loop of simulation
   x.i<-matrix(0,p2,m)
   for(i in 1:m){
     x.i[,i]<-mvrnorm(1, mux0, varx0) 
   }
   #print("A")
   
   Y<-rep(0,N)
   random2<-rnorm(m,0,1)
   for(i in 1:m){
     random3<-rnorm(NO[i],0,1)
     A<-S[i]
     for(j in 1:NO[i]){
       Y[A+j-1]<-beta0+beta1%*%(W[,A+j-1])+beta2%*%x.i[,i]+(sqrt(varv0))*random2[i]+(sqrt(vare0))*random3[j]
     }}
   
   X.ran<-matrix(0,p2,sumt)
   for(i in 1:m){
     a<-st[i]
     for(j in 1:to[i]){
       mu0<-rep(0,p2)
       X.ran[,a+j-1]<-x.i[,i]+mvrnorm(1, mu0, vareta0)
     }
   }
   #########
   #Now we need to sample from Y and W:  
   y.ran<-rep(0,n)
   w.ij.ran<-matrix(0,p1,n) 
   w.i.bar.ran<-matrix(0,p1,m)
   w.bar.ran<-rep(0,p1)
   for(i in 1:m){
     A<-S[i]
     B<-S[i+1]-1
     no.ran<-sample(A:B,no[i],replace=FALSE)
     a<-s[i]
     b<-s[i+1]-1
     for(j in 1:no[i]){
       y.ran[a+j-1]<-Y[no.ran[j]]
       w.ij.ran[,a+j-1]<-W[,no.ran[j]]
     }
     w.i.bar.ran[1,i]<-mean(w.ij.ran[1,a:b])
     w.i.bar.ran[2,i]<-mean(w.ij.ran[2,a:b])
     w.bar.ran<-w.bar.ran+(no[i]/n)*w.i.bar.ran[,i]
   }
   ################################  
   X.X<-matrix(0,p2+1,p2+1)
   weight<-rep(0,m)
   for(i in 1:m){
     a<-st[i]
     b<-st[i+1]-1
     X.X<-X.X+c(1,mean(X.ran[1,a:b]), mean(X.ran[2,a:b]))%*%t(c(1,mean(X.ran[1,a:b]),mean(X.ran[2,a:b])))} #when p2=2
   #print("E")
   for(i in 1:m){
     a<-st[i]
     b<-st[i+1]-1
     temp<-c(1,mean(X.ran[1,a:b]), mean(X.ran[2,a:b])) #when p2=2
     weight[i]<-1-t(temp)%*%solve(X.X)%*%temp}
   #print("F")
   ###############################  
   ybar.ran<-rep(0,m)
   Xbar.ran<-matrix(0,p2,m) 
   sum.y<-0
   sum.X<-rep(0,p2)
   for(i in 1:m){
     a<-s[i]
     b<-s[i+1]-1
     at<-st[i]
     bt<-st[i+1]-1
     ybar.ran[i]<-mean(y.ran[a:b])
     Xbar.ran[1,i]<-mean(X.ran[1,at:bt])  #when p2=2
     Xbar.ran[2,i]<-mean(X.ran[2,at:bt])
     sum.y<-sum.y+no[i]*ybar.ran[i]   
     sum.X[1]<-sum.X[1]+to[i]*Xbar.ran[1,i]   #when p2=2
     sum.X[2]<-sum.X[2]+to[i]*Xbar.ran[2,i]  
   } 
   ybarr<-sum.y/n
   XBARR<-sum.X/sumt #estimate of mu_x
   MUx.p[r,]<-XBARR
   #####################################  
   mu<-rep(0,m)
   for(i in 1:m){
     A<-S[i]
     B<-S[i+1]-1
     mu[i]<-mean(Y[A:B])}
   
   sum.XX<-rep(0,p2) 
   mXX<-rep(0,p2)
   for(i in 1:m){
     a<-st[i]
     for(k in 1:to[i]){
       sum.XX[1]<-sum.XX[1]+(X.ran[1,a+k-1]-Xbar.ran[1,i])^2
       sum.XX[2]<-sum.XX[2]+(X.ran[2,a+k-1]-Xbar.ran[2,i])^2
     }} 
   mXX[1]<-sum.XX[1]/(sumt-m) 
   mXX[2]<-sum.XX[2]/(sumt-m) 
   
   sum.yy<-0
   sum.ww<-matrix(0,p1,p1)
   msb.w<-matrix(0,p1,p1)
   for(i in 1:m){
     a<-s[i]
     for(k in 1:no[i]){
       sum.yy<-sum.yy+(y.ran[a+k-1]-ybar.ran[i])^2
       sum.ww<-sum.ww+(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])
     }
     msb.w<-msb.w+(no[i]/(m-1))*(w.i.bar.ran[,i]-w.bar.ran)%*%t(w.i.bar.ran[,i]-w.bar.ran)
   } 
   myy<-sum.yy/(n-m)
   mWW<-sum.ww/(n-m)
   ####################################
   vareta.p<-mXX   #estimate of var(eta) when p2=2
   VAReta.p[r,]<-vareta.p 
   #################################### 
   varx.p<-matrix(0,p2,p2)
   sum.X<-matrix(0,p2,p2)
   #sum.A<-matrix(0,p2,p2)
   tem.g1<-0
   for(i in 1:m){
     sum.X<-sum.X+to[i]*(Xbar.ran[,i]-XBARR)%*%t(Xbar.ran[,i]-XBARR)
     #sum.A<-sum.A+diag(mXX,p2)/to[i]
     tem.g1<-tem.g1+to[i]^2
   }
   g1m<-sumt-tem.g1/sumt
   varx.p<-((1/(m-1))*sum.X-mXX)*(m-1)/g1m #(1/m)*sum.A
 
 if(varx.p[1,1]<=0 | varx.p[2,2]<=0){
   varx.p.check<-varx.p
   varx.p.check[1,1]<-ifelse(varx.p[1,1]<=0,mXX[1]/m^{2/3},varx.p[1,1])
   varx.p.check[2,2]<-ifelse(varx.p[2,2]<=0,mXX[2]/m^{2/3},varx.p[2,2])
   sigma.check.av<-(1/2)*(varx.p.check[1,1]+varx.p.check[2,2])
   if(eigen(varx.p.check)$values[1]<0){
     a1<-(2*sigma.check.av)/((eigen(varx.p.check)$values[2])+sigma.check.av/(sqrt(m)))
     eig1<-a1*sigma.check.av/(sqrt(m))
     eig2<-a1*(eigen(varx.p.check)$values[2])
   }
   if(eigen(varx.p.check)$values[2]<0){
     a2<-(2*sigma.check.av)/((eigen(varx.p.check)$values[1])+sigma.check.av/(sqrt(m)))
     eig2<-a2*sigma.check.av/(sqrt(m))
     eig1<-a2*(eigen(varx.p.check)$values[1])
   }   
   varx.p<-eig1*(eigen(varx.p.check)$vectors[1,])%*%t(eigen(varx.p.check)$vectors[1,])+
     eig2*(eigen(varx.p.check)$vectors[2,])%*%t(eigen(varx.p.check)$vectors[2,])
 }
 
 VARx.p[,,r]<-varx.p #estimate of var(x)
 ################################### 
 tem.w<-matrix(0,p1,p1)
 tem.y<-rep(0,p1)
 for(i in 1:m){
   a<-s[i]
   for(j in 1:no[i]){
     tem.w<-tem.w+(w.ij.ran[,a+j-1]-w.bar.ran)%*%t(w.ij.ran[,a+j-1]) 
     tem.y<-tem.y+(w.ij.ran[,a+j-1]-w.bar.ran)*y.ran[a+j-1] 
   }  
 }
 beta1.p<-solve(tem.w)%*%tem.y
 beta1.p[1]<-ifelse(beta1.p[1]<0.01,0.01, beta1.p[1])
 beta1.p[2]<-ifelse(beta1.p[2]<0.01,0.01, beta1.p[2])
 beta1.p[1]<-ifelse(beta1.p[1]>1,1, beta1.p[1])
 beta1.p[2]<-ifelse(beta1.p[2]>1,1, beta1.p[2])
 
 B1.p[r,]<-beta1.p
 #beta1.p<rep(0,p1)   #for cehcking
 #################################
 sum.b2<-rep(0,p2)
 tem.g1<-0
 sum.b22<-0
 tem.g2<-0
 temm<-0
 
 for(i in 1:m){
   a<-s[i]
   for(j in 1:no[i]){
     temm<-temm+(to[i]*(1-to[i]/sumt)*t(w.i.bar.ran[,i])%*%solve(tem.w)%*%(w.ij.ran[,a+j-1]-w.bar.ran))
   }
   sum.b2<-sum.b2+to[i]*ybar.ran[i]*(Xbar.ran[,i]-XBARR)
   tem.g1<-tem.g1+to[i]^2
   tem.g2<-tem.g2+no[i]^2  
   sum.b22<-sum.b22+no[i]*(ybar.ran[i]-ybarr)^2
 }
 g1m<-sumt-tem.g1/sumt 
 g2m<-n-tem.g2/n
 
 beta2.p<-solve(sum.X/(m-1)-mXX)%*%(sum.b2/(m-1)) #g1m-temm)) 
 beta0.p<-ybarr-t(beta1.p)%*%w.bar.ran-t(beta2.p)%*%XBARR
 
 B2.p[r,]<-beta2.p
 B0.p[r]<-beta0.p
 ################################# 
 tem.w.bar<-0
 tem.w1<-0
 for(i in 1:m){
   a<-s[i]
   for(j in 1:no[i]){
     tem.w.bar<-tem.w.bar+(t(beta1.p)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]))^2  
     tem.w1<-tem.w1+no[i]*(t(beta1.p)%*%(w.i.bar.ran[,i]-w.bar.ran))^2
   }  
 }
 varv.p<-(((sum.b22/(m-1)) - myy- t(beta1.p)%*% (msb.w-mWW) %*%beta1.p )*(m-1)/g2m) - t(beta2.p)%*%(sum.b2/(g1m-0)) #-(1/g2m)*(tem.w1-((m-1)/(n-m))*tem.w.bar)
 #varv.p<-ifelse(varv.p<=0, t(beta2.p)%*%(sum.b2/(g1m-temm)), varv.p )   
 if(varv.p<=0)varv.p<-0.0001 #estimate of var(u)
 VARv.p[r]<-varv.p
 ################################## 
 tem.yw<-0
 tem.Aw1<-matrix(0,p1,p1)
 tem.Bw1<-matrix(0,p1,p1)
 A.w<-0
 B.w<-0
 for(i in 1:m){
   a<-s[i]
   for(j in 1:no[i]){
     tem.yw<-tem.yw+((y.ran[a+j-1]-ybar.ran[i])-t(beta1.p)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]))^2 
     tem.Aw1<-tem.Aw1+(w.ij.ran[,a+j-1]-w.bar.ran)%*%t(w.ij.ran[,a+j-1]) 
   }
   tem.Bw1<-tem.Bw1+(no[i]^2)*(w.i.bar.ran[,i]-w.bar.ran)%*%t(w.i.bar.ran[,i]-w.bar.ran)
 }
 
 tem.Aw2<-0
 for(i in 1:m){
   a<-s[i]
   for(j in 1:no[i]){
     tem.Aw2<-tem.Aw2+t(w.ij.ran[,a+j-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]) 
     B.w<-B.w+t(w.ij.ran[,a+j-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%tem.Bw1%*%solve(tem.Aw1)%*%(w.ij.ran[,a+j-1]-w.i.bar.ran[,i]) 
   }  
 }
 A.w<-(n-m)-tem.Aw2
 
 vare.p<-tem.yw/(n-m-p1)   #(1/A.w)*(tem.yw-(varv.p+t(beta2.p)%*%(sum.b2/(g1m-temm)))*B.w)
 #if(vare.p<=0)vare.p<-0 #estimate of var(e)
 VARe.p[r]<-vare.p
 ###############################################################################################################################
 muhat.p<-rep(0,m)
 mu.b.p<-rep(0,m)
 for(i in 1:m){
   fi<-(1-(no[i]/NO[i]))
   Ai<-vareta0/to[i]
   Hi<-diag(p2)+Ai%*%solve(varx0)
   Di.b<-vare0/(vare0+no[i]*(varv0+beta2%*%solve(Hi)%*%Ai%*%beta2))
   mu.b.p[i]<-(1-fi*Di.b)*(ybar.ran[i]+t(beta1)%*%(W.bar[,i]-w.i.bar.ran[,i]))+fi*Di.b*(beta0+t(beta1)%*%W.bar[,i]+t(beta2)%*%mux0)+fi*Di.b*(t(beta2)%*%varx0%*%solve(varx0+Ai)%*%(Xbar.ran[,i]-mux0) )
   
   Ai.p<-diag(mXX,p2)/to[i]
   varx.p<-diag(c(varx.p[1,1],varx.p[2,2]))
   Hi.p<-diag(p2)+Ai.p%*%solve(varx.p) 
   Di.p<-vare.p/(vare.p+no[i]*(varv.p+t(beta2.p)%*%solve(Hi.p)%*%Ai.p%*%beta2.p))
   muhat.p[i]<-(1-fi*Di.p)*(ybar.ran[i]+t(beta1.p)%*%(W.bar[,i]-w.i.bar.ran[,i]))+fi*Di.p*(beta0.p+t(beta1.p)%*%W.bar[,i]+t(beta2.p)%*%XBARR)+fi*Di.p*(t(beta2.p)%*%varx.p%*%solve(varx.p+Ai.p)%*%(Xbar.ran[,i]-XBARR) )
   
   MSE.p[r,i]<-(muhat.p[i]-mu[i])^2
   MSE1.p[r,i]<-(mu.b.p[i]-mu[i])^2
   MSE2.p[r,i]<-(muhat.p[i]-mu.b.p[i])^2
   MSEc.p[r,i]<-(muhat.p[i]-mu.b.p[i])*(mu.b.p[i]-mu[i])
   Bias.p[r,i]<-muhat.p[i]-mu[i]          
 }
 #-----------------------------------------------------------------
 g1.p<-rep(0,m)
 for(i in 1:m){
   Ai.p<-diag(mXX,p2)/to[i]
   Hi.p<-diag(p2)+Ai.p%*%solve(varx.p) 
   Di.p<-vare.p/no[i]
   Mi.p<-solve(Ai.p)+solve(varx.p)
   g1.p[i]<-(Di.p*(varv.p+t(beta2.p)%*%solve(Mi.p)%*%beta2.p))/(Di.p+varv.p+t(beta2.p)%*%solve(Mi.p)%*%beta2.p)
 }
 G1.p[r,]<-g1.p
 ##########################Jackknife method:#######################
 beta0.p.j<-rep(0,m)
 beta1.p.j<-matrix(0,p1,m)
 beta2.p.j<-matrix(0,p2,m)
 vareta.p.j<-matrix(0,p2,m)
 vare.p.j<-rep(0,m)
 varv.p.j<-rep(0,m)
 mux.p.j<-matrix(0,p2,m)
 varx.p.j<-array(0,c(p2,p2,m))
 
 XBARR.j<-matrix(0,p2,m)
 ybarr.j<-rep(0,m)
 w.bar.j<-matrix(0,p1,m)
 
 #print(r)
 for(j in 1:m){
   
   sum.y.j<-0
   sum.X.j<-rep(0,p2)
   sum.w.j<-rep(0,p1)
   for(i in 1:m){if (i !=j){
     sum.y.j<-sum.y.j+no[i]*ybar.ran[i]
     sum.X.j<-sum.X.j+to[i]*Xbar.ran[,i]
     sum.w.j<-sum.w.j+no[i]*w.i.bar.ran[,i]
   } }
   ybarr.j[j]<-sum.y.j/(n-no[j])
   XBARR.j[,j]<-sum.X.j/(sumt-to[j])
   w.bar.j[,j]<-sum.w.j/(n-no[j])
   
   sum.XX.j<-rep(0,p2)
   for(i in 1:m){if(i !=j){
     a<-st[i]
     for(k in 1:to[i]){
       sum.XX.j<-sum.XX.j+(X.ran[,a+k-1]-Xbar.ran[,i])^2}
   } }
   mXX.j<-sum.XX.j/((sumt-to[j])-(m-1))
   
   
   sum.yy.j<-0
   sum.ww.j<-matrix(0,p1,p1)
   msb.w.j<-matrix(0,p1,p1)
   for(i in 1:m){if(i !=j){
     a<-s[i]
     for(k in 1:no[i]){
       sum.yy.j<-sum.yy.j+(y.ran[a+k-1]-ybar.ran[i])^2
       sum.ww.j<-sum.ww.j+(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])
     }
     msb.w.j<-msb.w.j+(no[i]/(m-1-1))*(w.i.bar.ran[,i]-w.bar.j[,j])%*%t(w.i.bar.ran[,i]-w.bar.j[,j])
   } }
   myy.j<-sum.yy.j/((n-no[j])-(m-1))
   mWW.j<-sum.ww.j/((n-no[j])-(m-1))
   #print(r)
   
   mux.p.j[,j]<-XBARR.j[,j]
   vareta.p.j[,j]<-mXX.j             
   ################################### 
   sum.X<-matrix(0,p2,p2)
   sum.A<-matrix(0,p2,p2)
   tem.g1<-0
   for(i in 1:m){if(i !=j){
     sum.X<-sum.X+to[i]*(Xbar.ran[,i]-XBARR.j[,j])%*%t(Xbar.ran[,i]-XBARR.j[,j])
     sum.A<-sum.A+diag(mXX.j,p2)/to[i]
     tem.g1<-tem.g1+to[i]^2
   }}
   g1m.j<-(sumt-to[j])-tem.g1/(sumt-to[j]) 
   varx.p.j[,,j]<-((1/((m-1)-1))*sum.X-mXX.j)*(m-1-1)/g1m.j #(1/(m-1))*sum.A
 
 if(varx.p.j[1,1,j]<=0 | varx.p.j[2,2,j]<=0){
   varx.p.check<-varx.p.j[,,j]
   varx.p.check[1,1]<-ifelse(varx.p.j[1,1,j]<=0,mXX.j[1]/(m-1)^{2/3},varx.p.j[1,1,j])
   varx.p.check[2,2]<-ifelse(varx.p.j[2,2,j]<=0,mXX.j[2]/(m-1)^{2/3},varx.p.j[2,2,j])
   sigma.check.av<-(1/2)*(varx.p.check[1,1]+varx.p.check[2,2])
   if(eigen(varx.p.check)$values[1]<0){
     aa<-(2*sigma.check.av)/((eigen(varx.p.check)$values[2])+sigma.check.av/(sqrt(m-1)))
     eig1<-aa*sigma.check.av/(sqrt(m-1))
     eig2<-aa*(eigen(varx.p.check)$values[2])
   }
   if(eigen(varx.p.check)$values[2]<0){
     ab<-(2*sigma.check.av)/((eigen(varx.p.check)$values[1])+sigma.check.av/(sqrt(m-1)))
     eig2<-ab*sigma.check.av/(sqrt(m-1))
     eig1<-ab*(eigen(varx.p.check)$values[1])
   }   
   varx.p.j[,,j]<-eig1*(eigen(varx.p.check)$vectors[1,])%*%t(eigen(varx.p.check)$vectors[1,])+
     eig2*(eigen(varx.p.check)$vectors[2,])%*%t(eigen(varx.p.check)$vectors[2,])
 }  
 #################################### 
 tem.w<-matrix(0,p1,p1)
 tem.y<-rep(0,p1)
 for(i in 1:m){if(i !=j){
   a<-s[i]
   for(k in 1:no[i]){
     tem.w<-tem.w+(w.ij.ran[,a+k-1]-w.bar.j[,j])%*%t(w.ij.ran[,a+k-1]) 
     tem.y<-tem.y+(w.ij.ran[,a+k-1]-w.bar.j[,j])*y.ran[a+k-1] 
   }  
 }}
 beta1.p.j[,j]<-solve(tem.w)%*%tem.y
 beta1.p.j[1,j]<-ifelse(beta1.p.j[1,j]<0.01,0.01, beta1.p.j[1,j])
 beta1.p.j[2,j]<-ifelse(beta1.p.j[2,j]<0.01,0.01, beta1.p.j[2,j])
 beta1.p.j[1,j]<-ifelse(beta1.p.j[1,j]>1,1, beta1.p.j[1,j])
 beta1.p.j[2,j]<-ifelse(beta1.p.j[2,j]>1,1, beta1.p.j[2,j])
 
 #beta1.p.j[,j]<-rep(0,p1)  #for checking
 #################################
 sum.b2.j<-rep(0,p2)
 #   tem.g1<-0
 sum.b22.j<-0
 tem.g2<-0
 temm<-0
 
 for(i in 1:m){if(i !=j){
   a<-s[i]
   for(k in 1:no[i]){
     temm<-temm+(to[i]*(1-to[i]/(sumt-to[j]))*t(w.i.bar.ran[,i])%*%solve(tem.w)%*%(w.ij.ran[,a+k-1]-w.bar.j[,j]))
   }
   sum.b2.j<-sum.b2.j+to[i]*ybar.ran[i]*(Xbar.ran[,i]-XBARR.j[,j])
   #    tem.g1<-tem.g1+to[i]^2
   tem.g2<-tem.g2+no[i]^2  
   sum.b22.j<-sum.b22.j+no[i]*(ybar.ran[i]-ybarr.j[j])^2
 }}
 #   g1m.j<-(sumt-to[j])-tem.g1/(sumt-to[j]) 
 g2m.j<-(n-no[j])-tem.g2/(n-no[j])
 
 beta2.p.j[,j]<-solve(sum.X/(m-1-1)-mXX.j)%*%(sum.b2.j/(m-1-1)) #g1m.j-temm))  
 beta0.p.j[j]<-ybarr.j[j]-t(beta1.p.j[,j])%*%w.bar.j[,j]-t(beta2.p.j[,j])%*%XBARR.j[,j]
 ################################# 
 tem.w.bar<-0
 tem.w1<-0
 for(i in 1:m){if(i !=j){
   a<-s[i]
   for(k in 1:no[i]){
     tem.w.bar<-tem.w.bar+(t(beta1.p.j[,j])%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]))^2  
     tem.w1<-tem.w1+no[i]*(t(beta1.p.j[,j])%*%(w.i.bar.ran[,i]-w.bar.j[,j]))^2
   }  
 }}
 varv.p.j[j]<-(((sum.b22.j/((m-1)-1)) - myy.j- t(beta1.p.j[,j])%*% (msb.w.j-mWW.j) %*%beta1.p.j[,j]  )*((m-1)-1)/g2m.j) - t(beta2.p.j[,j])%*%(sum.b2.j/(g1m.j-0)) #- (1/g2m.j)*(tem.w1-(((m-1)-1)/((n-no[j])-(m-1)))*tem.w.bar)
 if(varv.p.j[j]<=0)varv.p.j[j]<-0.0001 #estimate of var(u)
 ################################## 
 tem.yw<-0
 tem.Aw1<-matrix(0,p1,p1)
 tem.Bw1<-matrix(0,p1,p1)
 A.w<-0
 B.w<-0
 for(i in 1:m){if(i !=j){
   a<-s[i]
   for(k in 1:no[i]){
     tem.yw<-tem.yw+((y.ran[a+k-1]-ybar.ran[i])-t(beta1.p.j[,j])%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]))^2 
     tem.Aw1<-tem.Aw1+(w.ij.ran[,a+k-1]-w.bar.j[,j])%*%t(w.ij.ran[,a+k-1]) 
   }
   tem.Bw1<-tem.Bw1+(no[i]^2)*(w.i.bar.ran[,i]-w.bar.j[,j])%*%t(w.i.bar.ran[,i]-w.bar.j[,j])
 }}
 
 tem.Aw2<-0
 for(i in 1:m){if(i !=j){
   a<-s[i]
   for(k in 1:no[i]){
     tem.Aw2<-tem.Aw2+t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]) 
     B.w<-B.w+t(w.ij.ran[,a+k-1]-w.i.bar.ran[,i])%*%solve(tem.Aw1)%*%tem.Bw1%*%solve(tem.Aw1)%*%(w.ij.ran[,a+k-1]-w.i.bar.ran[,i]) 
   }  
 }}
 A.w<-((n-no[j])-(m-1))-tem.Aw2
 
 vare.p.j[j]<-tem.yw/((n-no[j])-(m-1)-p1)   #(1/A.w)*(tem.yw-(varv.p.j[j]+t(beta2.p.j[,j]) %*%(sum.b2.j/(g1m.j-temm)) )*B.w)
 #  if(vare.p.j[j]<=0)vare.p.j[j]<-0 #estimate of var(e)
 ######
 }# loop of j
 #######################
 g.hat.p.j<-rep(0,m)
 g1.hat.p.j<-rep(0,m)
 g.hat.wp.j<-rep(0,m)
 g1.hat.wp.j<-rep(0,m)
 c1<-matrix(0,m,m)
 for(z in 1:m){
   fi<-((NO[z]-no[z])/NO[z])
   for(i in 1:m){
     #print("H1")
     Ai.p.j<-diag(mXX.j,p2)/to[z]
     varx.p.j[,,i]<-diag(c(varx.p.j[1,1,i], varx.p.j[2,2,i]))
     Hi.p.j<-diag(p2)+Ai.p.j%*%solve(varx.p.j[,,i]) 
     Di.p.j<-vare.p.j[i]/(vare.p.j[i]+no[z]*(varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Hi.p.j)%*%Ai.p.j%*%beta2.p.j[,i]))
     Di1.p.j<-vare.p.j[i]/no[z]
     Mi.p.j<-solve(Ai.p.j)+solve(varx.p.j[,,i])
     g1.hat.p.j[z]<-g1.hat.p.j[z]+(Di1.p.j*(varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]))/(Di1.p.j+varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]) - g1.p[z]
     g1.hat.wp.j[z]<-g1.hat.wp.j[z]+ weight[i]*( (Di1.p.j*(varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]))/(Di1.p.j+varv.p.j[i]+t(beta2.p.j[,i])%*%solve(Mi.p.j)%*%beta2.p.j[,i]) - g1.p[z] )
     
     muhat.p.j<-(1-fi*Di.p.j)*(ybar.ran[z]+t(beta1.p.j[,i])%*%(W.bar[,z]-w.i.bar.ran[,z]))+fi*Di.p.j*(beta0.p.j[i]+t(beta1.p.j[,i])%*%W.bar[,z]+t(beta2.p.j[,i])%*%XBARR.j[,i])+fi*Di.p.j*(t(beta2.p.j[,i])%*%varx.p.j[,,i]%*%solve(varx.p.j[,,i]+Ai.p.j)%*%(Xbar.ran[,z]-XBARR.j[,i]) )
     g.hat.p.j[z]<-g.hat.p.j[z]+ (muhat.p.j -muhat.p[z])^2     
     g.hat.wp.j[z]<-g.hat.wp.j[z]+weight[i]*((muhat.p.j -muhat.p[z] )^2)
     c1[i,z]<- muhat.p.j
   }
   G.hat.p.j[r,z]<-g.hat.p.j[z]
   G1.hat.p.j[r,z]<-g1.hat.p.j[z]
   #print("H2")
 }
 for (z in 1:m){
   #print("H3")
   mse.p.j[r,z]<-((m-1)/m)*g.hat.p.j[z]+g1.p[z]-((m-1)/m)*g1.hat.p.j[z]
   mse1.p.j[r,z]<-g1.p[z]-((m-1)/m)*g1.hat.p.j[z]
   mse2.p.j[r,z]<-((m-1)/m)*g.hat.p.j[z]
   
   mse.wp.j[r,z]<-g.hat.wp.j[z]+g1.p[z]-g1.hat.wp.j[z]
   mse1.wp.j[r,z]<-g1.p[z]-g1.hat.wp.j[z]
   mse2.wp.j[r,z]<-g.hat.wp.j[z]
   
   }
 #print("H4")  
 cat(r)
 }#loop for MC.
#-----------------------------------------------------------------------------------     
 RB.b0.p<-(1/MC)*( sum(B0.p)/beta0)-1
 RB.b1.1.p<-(1/MC)*(apply(B1.p,2,sum)[1]/beta1[1])-1
 RB.b1.2.p<-(1/MC)*(apply(B1.p,2,sum)[2]/beta1[2])-1
 RB.b2.1.p<-(1/MC)*(apply(B2.p,2,sum)[1]/beta2[1])-1
 RB.b2.2.p<-(1/MC)*(apply(B2.p,2,sum)[2]/beta2[2])-1
 RB.v.p<-(1/MC)*( sum(VARv.p)/varv0)-1
 RB.e.p<-(1/MC)*( sum(VARe.p)/vare0)-1
 RB.eta.1.p<-(1/MC)*(apply(VAReta.p,2,sum)[1]/vareta0[1,1])-1
 RB.eta.2.p<-(1/MC)*(apply(VAReta.p,2,sum)[2]/vareta0[2,2])-1
 RB.x.1.p<-(1/MC)*(apply(VARx.p,1:2,sum)[1]/varx0[1,1])-1
 RB.x.2.p<-(1/MC)*(apply(VARx.p,1:2,sum)[2]/varx0[2,2])-1
 RB.mu.1.p<-(1/MC)*(apply(MUx.p,2,sum)[1]/mux0[1])-1
 RB.mu.2.p<-(1/MC)*(apply(MUx.p,2,sum)[2]/mux0[2])-1
 
 RE.b0.p<-(1/MC)*( sum(abs( (B0.p/beta0)-1  ) ) )
 RE.b1.1.p<-(1/MC)*( sum(abs( (B1.p[,1]/beta1[1])-1 ) ) )
 RE.b1.2.p<-(1/MC)*( sum(abs( (B1.p[,2]/beta1[2])-1 ) ) )
 RE.b2.1.p<-(1/MC)*( sum(abs( (B2.p[,1]/beta2[1])-1 ) ) )
 RE.b2.2.p<-(1/MC)*( sum(abs( (B2.p[,2]/beta2[2])-1 ) ) )
 RE.v.p<-(1/MC)*(  sum(abs( (VARv.p/varv0)-1 ) ) )
 RE.e.p<-(1/MC)*(  sum(abs( (VARe.p/vare0)-1 ) ) )
 RE.eta.1.p<-(1/MC)*(  sum(abs( (VAReta.p[,1]/vareta0[1,1])-1 ) ) )
 RE.eta.2.p<-(1/MC)*(  sum(abs( (VAReta.p[,2]/vareta0[2,2])-1 ) ) )
 RE.x.1.p<-(1/MC)*(  sum(abs( (VARx.p[,,1]/varx0[1,1])-1 ) ) )
 RE.x.2.p<-(1/MC)*(  sum(abs( (VARx.p[,,2]/varx0[2,2])-1 ) ) )
 RE.mu.1.p<-(1/MC)*(  sum(abs( (MUx.p[,1]/mux0[1])-1 ) ) )
 RE.mu.2.p<-(1/MC)*(  sum(abs( (MUx.p[,2]/mux0[2])-1 ) ) )
 
 Abias.p<-rep(0,m)
 Amse.p.j<-rep(0,m)
 Amse.wp.j<-rep(0,m)
 Amse1.p.j<-rep(0,m)
 Amse1.wp.j<-rep(0,m)
 Amse2.p.j<-rep(0,m)
 Amse2.wp.j<-rep(0,m)
 
 
 EMSE.p<-rep(0,m)
 EMSE1.p<-rep(0,m)
 EMSE2.p<-rep(0,m)
 EMSEc.p<-rep(0,m)
 EG1.p<-rep(0,m)
 
 Bias.mse.p<-rep(0,m)
 Bias.mse.wp<-rep(0,m)
 
 for(i in 1:m){
   Amse.p.j[i]<-mean(mse.p.j[,i])
   Amse.wp.j[i]<-mean(mse.wp.j[,i])
   Amse1.p.j[i]<-mean(mse1.p.j[,i])
   Amse1.wp.j[i]<-mean(mse1.wp.j[,i])
   Amse2.p.j[i]<-mean(mse2.p.j[,i])
   Amse2.wp.j[i]<-mean(mse2.wp.j[,i])
    
   EMSE.p[i]<-mean(MSE.p[,i])
   EMSE1.p[i]<-mean(MSE1.p[,i])
   EMSE2.p[i]<-mean(MSE2.p[,i])
   EMSEc.p[i]<-mean(MSEc.p[,i])
   EG1.p[i]<-mean(G1.p[,i])
   Abias.p[i]<-mean(Bias.p[,i])
   
   Bias.mse.p[i]<-Amse.p.j[i]-EMSE.p[i]
   Bias.mse.wp[i]<-Amse.wp.j[i]-EMSE.p[i]            
 }
 
 
 Amse.p.j
 Amse.wp.j     
 Amse1.p.j
 Amse1.wp.j 
 Amse2.p.j
 Amse2.wp.j 
 
 EMSE.p
 EMSE1.p
 EMSE2.p
 EMSEc.p
 EG1.p
 
 mean(Amse.p.j)
 mean(Amse.wp.j)
 mean(Amse1.p.j)
 mean(Amse1.wp.j)
 mean(Amse2.p.j)
 mean(Amse2.wp.j)
 
 mean(EMSE.p)
 mean(EMSE1.p)
 mean(EMSE2.p)
 mean(EMSEc.p)
 mean(EG1.p)
 
 Abias.p
 Bias.mse.p
 Bias.mse.wp
 
 RB.v.p                             
 RB.e.p
 RB.eta.1.p
 RB.eta.2.p
 RB.b0.p
 RB.b1.1.p
 RB.b1.2.p
 RB.b2.1.p
 RB.b2.2.p
 RB.x.1.p
 RB.x.2.p
 RB.mu.1.p
 RB.mu.2.p
 
 RE.v.p                             
 RE.e.p
 RE.eta.1.p
 RE.eta.2.p
 RE.b0.p
 RE.b1.1.p
 RE.b1.2.p
 RE.b2.1.p
 RE.b2.2.p
 RE.x.1.p
 RE.x.2.p
 RE.mu.1.p
 RE.mu.2.p
 
 mean(abs(Abias.p))
 
 RB.mse.p.j<-rep(0,m)
 
 for(i in 1:m){
   RB.mse.p.j[i]<-(mean(mse.p.j[,i])/mean(MSE.p[,i]))-1}
 ARB.mse.p.j<-mean(RB.mse.p.j)

 ARB.mse.p.j
 RB.mse.p.j
 
 
 RB.mse.wp.j<-rep(0,m)
 
 for(i in 1:m){
   RB.mse.wp.j[i]<-(mean(mse.wp.j[,i])/mean(MSE.p[,i]))-1}
 
 ARB.mse.wp.j<-mean(RB.mse.wp.j)
 ARB.mse.wp.j
 RB.mse.wp.j
 
 
 CV.mse.p<-rep(0,m)
 for (i in 1:m){
   CV.mse.p[i]<-sqrt( mean(mse.p.j[,i]-MSE.p[,i])^2 )/ mean(MSE.p[,i])}
 ACV.mse.p<-mean(CV.mse.p)
 ACV.mse.p
 CV.mse.p
 
 
 CV.mse.wp<-rep(0,m)
 for (i in 1:m){
   CV.mse.wp[i]<-sqrt( mean(mse.wp.j[,i]-MSE.p[,i])^2 )/ mean(MSE.p[,i])}
 ACV.mse.wp<-mean(CV.mse.wp)
 ACV.mse.wp
 CV.mse.wp
 
 #Finish!
 #==========================================================
 #Table 1:
 g1<-rep(0,m)
 for(i in 1:m){
   Ai<-vareta0/to[i]
   Hi<-diag(p2)+Ai%*%solve(varx0)
   Di.b<-vare0/no[i]
   Mi<-solve(Ai)+solve(varx0)
   g1[i]<-(Di.b*(varv0+t(beta2)%*%solve(Mi)%*%beta2))/(Di.b+varv0+t(beta2)%*%solve(Mi)%*%beta2)
 }
 g1
 #-------------------------------------------------------------
 Area<-c(1:m)
 Result.RB<-cbind(Area,no,g1,100*RB.mse.p.j,100*RB.mse.wp.j)
 RESULT.RB<-round(Result.RB,digit=2)
 RESULT.RB
 
 Result.MSE<-cbind(Area,no,EMSE.p,EMSE1.p,EMSE2.p, EMSEc.p)
 RESULT.MSE<-round(Result.MSE,digit=2)
 RESULT.MSE
#-------------------------------------------------------------