## Dec 14, 2008
## Estimating the parameters
est.g34.inv<-function(par, data, cor=TRUE){
x1<-data$x1
x2<-data$x2
## factor correlation matrix
phi<-matrix(c(1,      par[1],  par[2],
              par[1], 1,       par[3],
              par[2], par[3],  1    ), nrow=3, byrow=T) 
 
## factor residual variance
psi1<-diag(1,9,9)
for (i in 1:9){
  psi1[i,i]<-par[3+i]
}
 
## factor loading matrix
lam11<-diag(1,3,3)
for (i in 1:3){
  lam11[i,i]<-par[12+i]
}
 
lam12<-matrix(par[16:33],nrow=6,byrow=T)
 
lambda1<-rbind(lam11,lam12)
 
## the population covariance matrix
sigma1<-lambda1%*%phi%*%t(lambda1)+psi1
 
## the observed covariance/correlation matrix
if (cor){
  S1<-cor(x1)
}else{
  S1<-cov(x1)
}
 
## the likelihood function
n1<-dim(x1)[1]
p1<-dim(x2)[2]
 
l1<-n1*(log(det(sigma1))+sum(diag(S1%*%solve(sigma1)))-log(det(S1))-p1)/2
 
## for the second group
## factor residual variance
psi2<-diag(1,9,9)
for (i in 1:9){
  psi2[i,i]<-par[33+i]
}
 
## factor loading matrix
lam21<-diag(1,3,3)
for (i in 1:3){
  lam21[i,i]<-par[43+i]
}
 
lam22<-matrix(par[46:63],nrow=6,byrow=T)
 
lambda2<-rbind(lam21,lam22)
 
## the population covariance matrix
sigma2<-lambda2%*%phi%*%t(lambda2)+psi2
 
## the observed covariance/correlation matrix
if (cor){
  S2<-cor(x2)
}else{
  S2<-cov(x2)
}
 
## the likelihood function
n2<-dim(x2)[1]
p2<-dim(x2)[2]
 
l2<-n2*(log(det(sigma2))+sum(diag(S2%*%solve(sigma2)))-log(det(S2))-p2)/2
 
l1+l2
}
 
## estimate the model with different factor correlation matrices
 
est.g34.v<-function(par, data, cor=TRUE){
x1<-data$x1
x2<-data$x2
## factor correlation matrix
phi1<-matrix(c(1,      par[1],  par[2],
              par[1], 1,       par[3],
              par[2], par[3],  1    ), nrow=3, byrow=T) 
 
## factor residual variance
psi1<-diag(1,9,9)
for (i in 1:9){
  psi1[i,i]<-par[3+i]
}
 
## factor loading matrix
lam11<-diag(1,3,3)
for (i in 1:3){
  lam11[i,i]<-par[12+i]
}
 
lam12<-matrix(par[16:33],nrow=6,byrow=T)
 
lambda1<-rbind(lam11,lam12)
 
## the population covariance matrix
sigma1<-lambda1%*%phi1%*%t(lambda1)+psi1
 
## the observed covariance/correlation matrix
if (cor){
  S1<-cor(x1)
}else{
  S1<-cov(x1)
}
 
## the likelihood function
n1<-dim(x1)[1]
p1<-dim(x2)[2]
 
l1<-n1*(log(det(sigma1))+sum(diag(S1%*%solve(sigma1)))-log(det(S1))-p1)/2
 
## for the second group
## factor correlation matrix
phi2<-matrix(c(1,      par[64],  par[65],
              par[64], 1,       par[66],
              par[65], par[66],  1    ), nrow=3, byrow=T) 
 
## factor residual variance
psi2<-diag(1,9,9)
for (i in 1:9){
  psi2[i,i]<-par[33+i]
}
 
## factor loading matrix
lam21<-diag(1,3,3)
for (i in 1:3){
  lam21[i,i]<-par[43+i]
}
 
lam22<-matrix(par[46:63],nrow=6,byrow=T)
 
lambda2<-rbind(lam21,lam22)
 
## the population covariance matrix
sigma2<-lambda2%*%phi2%*%t(lambda2)+psi2
 
## the observed covariance/correlation matrix
if (cor){
  S2<-cor(x2)
}else{
  S2<-cov(x2)
}
 
## the likelihood function
n2<-dim(x2)[1]
p2<-dim(x2)[2]
 
l2<-n2*(log(det(sigma2))+sum(diag(S2%*%solve(sigma2)))-log(det(S2))-p2)/2
 
l1+l2
}
 
 
## data
#x<-read.table("test.txt")
#res<-nlm(hof, par, x=x, cor=T, hessian=T)
 
## data generation
library(mvtnorm)
sim.data.g3<-function(par,n1,n2){
 
## factor correlation matrix
phi<-matrix(c(1,      par[1],  par[2],
              par[1], 1,       par[3],
              par[2], par[3],  1    ), nrow=3, byrow=T) 
 
## factor residual variance
psi<-diag(1,9,9)
for (i in 1:9){
  psi[i,i]<-par[3+i]
}
 
## factor loading matrix
lam1<-diag(1,3,3)
for (i in 1:3){
  lam1[i,i]<-par[12+i]
}
 
lam2<-matrix(par[16:33],nrow=6,byrow=T)
 
lambda<-rbind(lam1,lam2)
 
x<-array(NA, dim=c(n1,9))
 
for (i in 1:n1){
  f<-rmvnorm(1,c(0,0,0),phi)
  x[i,]<-rmvnorm(1,lambda%*%t(f),psi)
}
 
## for the second group
## factor residual variance
psi2<-diag(1,9,9)
for (i in 1:9){
  psi2[i,i]<-par[33+i]
}
 
## factor loading matrix
lam21<-diag(1,3,3)
for (i in 1:3){
  lam21[i,i]<-par[43+i]
}
 
lam22<-matrix(par[46:63],nrow=6,byrow=T)
 
lambda2<-rbind(lam21,lam22)
 
y<-array(NA, dim=c(n2,9))
 
for (i in 1:n2){
  f<-rmvnorm(1,c(0,0,0),phi)
  y[i,]<-rmvnorm(1,lambda2%*%t(f),psi2)
}
 
list(x1=x,x2=y)
}
 
## starting values
par<-c(.5,.4,.3,rep(.5,9),rep(.707,3), rep(c(.41,.41,0,0,.439,.439,.422,0,.422),2),
       rep(.5,9),rep(.707,3), rep(c(.707,0,0,0,.707,0,0,0,.707),2))
 
par1<-c(.5,.4,.3,rep(.45,9),rep(.707,3), rep(c(.41,.41,0,0,.439,.439,.422,0,.422),2),
       rep(.45,9),rep(.707,3), rep(c(.707,0,0,0,.707,0,0,0,.707),2))
par2<-c(.5,.4,.3,rep(.45,9),rep(.707,3), rep(c(.41,.41,0,0,.439,.439,.422,0,.422),2),
       rep(.45,9),rep(.707,3), rep(c(.707,0,0,0,.707,0,0,0,.707),2),.5,.4,.3)
 
n1<-500
n2<-500
 
R<-1000
 sim.res<-NULL
 chi.res<-NULL
 i<-1
 
  while (i<=R){
     data.i<-sim.data.g3(par,n1,n2)
     res1.i<-nlm(est.g34.inv,par1,data=data.i,cor=T)
     res2.i<-nlm(est.g34.v,par2,data=data.i,cor=T)
     if (res1.i$code<=2 && res2.i$code<=2 ){
        sim.res<-rbind(sim.res,c(res1.i$estimate,res2.i$estimate))
        chi.res<-rbind(chi.res,c(res1.i$minimum,res2.i$minimum))
        i<-i+1
     }
   }
 
 
 
dim(sim.res)