nobs <- 500
nsim <- 1000

for (i in 1:nsim) {
    x1 <- (rchisq(nobs, df=1) - 1)/sqrt(2)
    x2 <- (rchisq(nobs, df=1) - 1)/sqrt(2)
    
    covadj <- matrix(c(0.085,0.093,0.093,0.907), ncol=2)
    
    x.matrix <- cbind(x1,x2)%*%chol(covadj)
    f1 <- x.matrix[,1]
    f2 <- x.matrix[,2]
    
    y <- 1.074*f1 + .109*f2 + (rchisq(nobs, df=1)-1)*sqrt(0.871)/sqrt(2) + 1.554
    
    ex1 <- 1*f1 + (rchisq(nobs, df=1)-1)*sqrt(0.207)/sqrt(2) + 1.146
    ex2 <- 1.086*f1 + (rchisq(nobs, df=1)-1)*sqrt(0.203)/sqrt(2) + 1.122
    ex3 <- 1.624*f1 + (rchisq(nobs, df=1)-1)*sqrt(0.553)/sqrt(2) + 1.393
    ex4 <- 1.262*f1 + (rchisq(nobs, df=1)-1)*sqrt(0.298)/sqrt(2) + 1.203
    in1 <- 1*f2 + (rchisq(nobs, df=1)-1)*sqrt(0.945)/sqrt(2) + 2.286
    in2 <- 0.747*f2 + (rchisq(nobs, df=1)-1)*sqrt(0.860)/sqrt(2) + 1.956
    in3 <- 0.888*f2 + (rchisq(nobs, df=1)-1)*sqrt(0.865)/sqrt(2) + 1.983
    in4 <- 0.675*f2 + (rchisq(nobs, df=1)-1)*sqrt(0.631)/sqrt(2) + 1.636
    
    obs.matrix <- cbind(ex1,ex2,ex3,ex4,in1,in2,in3,in4,y)
    
    print(i)
    write.table(obs.matrix, 
                file=paste("z:/mplus/mvchi/mvchi",i,".dat",sep=""), 
                row.names=FALSE, col.names=FALSE)
}

write.table(paste("mvchi",1:nsim,".dat", sep=""), 
            file="z:/MPlus/MVChi/mvchilist.dat",
            quote=FALSE, col.names=FALSE, row.names=FALSE)