***------------------------------------------------------------***;
*Data are the sample covariance matrix in Table 2 of Weston & Gore 2006;
*p=12; *N=403;
***------------------------------------------------------------***;

options ls=95 nodate nonumber;
title;
proc iml;
reset noname;
***--------------------------------------------------------------------***;
**** The following subroutine is to perform vech(.) function;
***--------------------------------------------------------------------***;
start vech(A,Va);
l=0;
p=nrow(A);
pstar=p*(p+1)/2;
Va=j(pstar,1,0);
do i=1 to p;
        do j=i to p;
                l=l+1;
                Va[l]=A[j,i];
        end;
end;
finish;
***--------------------------------------------------------------------***;
* creates duplication matrix as defined in Magnus&Neudecker*;
***--------------------------------------------------------------------***;
start DP(p, dup);
  Dup=j(p*p,p*(p+1)/2,0);
  count=0;
  do j=1 to p;
    do i=j to p;
      count=count+1;
        if i=j then do;
                Dup[(j-1)*p+j, count]=1;
        end;
        else do;
                Dup[(j-1)*p+i, count]=1;
                Dup[(i-1)*p+j, count]=1;
        end;
    end;
  end;
finish;
***--------------------------------------------------------------------***;
**** Subroutine to evaluate sigma0 and (dsigma0/dtheta);
*Model1 as represented by the solid arrows of Figure 6 of Deng & Yuan;
***--------------------------------------------------------------------***;
start sigdsig_sem1(p,theta0,vsig,Sigma,vdsig);
p_x=3;
p_y=9;
ps=p*(p+1)/2;
lamb_x=theta0[1:3];
lamb_xt=lamb_x`;

lamb_y=j(p_y,3,0);
lamb_y[1,1]=1.0;
lamb_y[2:3,1]=theta0[4:5];
lamb_y[4,2]=1.0;
lamb_y[5:6,2]=theta0[6:7];
lamb_y[7,3]=1.0;
lamb_y[8:9,3]=theta0[8:9];
lamb_yt=lamb_y`;

gamma=j(3,1,0);
gamma[1:2]=theta0[10:11];
gamma_t=gamma`;

**AB=I-B in LISREL notation;
beta=j(3,3,0);
beta[2,1]=theta0[12];
beta[3,2]=theta0[13];
A=i(3);
AB=A-beta;
*phi=theta0[14];
sig_zeta=diag(theta0[14:16]);
psi_x=diag(theta0[17:19]);
psi_y=diag(theta0[20:28]);

phi=1;
*computing Sigma based on structured parameters;
ABin=inv(AB);
ABin_t=ABin`;
sigyy=lamb_y*ABin*(gamma*phi*gamma_t+sig_zeta)*ABin_t*lamb_yt+psi_y;
sigxy=lamb_x*phi*gamma_t*ABin_t*lamb_yt;
sigxx=lamb_x*phi*lamb_xt+psi_x;
sigma=j(p,p,0);
sigma[1:p_x,1:p_x]=sigxx;
sigma[1:p_x,(p_x+1):p]=sigxy;
sigma[(p_x+1):p,1:p_x]=sigxy`;
sigma[(p_x+1):p,(p_x+1):p]=sigyy;

run vech(sigma,vsig);

***** The following is to evaluate (dsigma0/dtheta);

*dlamb_x is the derivative with Lamb_x;
dlamb_x=j(ps,p_x,0);
t1=phi*gamma_t*ABin_t*lamb_yt;
phl=phi*lamb_xt;
lph=lamb_x*phi;
tyy=j(p_y,p_y,0);
do j=1 to p_x;
	tt=j(3,1,0);
	tt[j]=1;
	txx=tt*phl+lph*tt`;
	txy=tt*t1;
	ttt=(txx||txy)//(txy`||tyy);
	run vech(ttt,tttt);
	dlamb_x[,j]=tttt;
end;


*dlamb_y is the derivative with Lamb_y;
dlamb_y=j(ps,6,0);
t1=ABin*(gamma*phi*gamma_t+sig_zeta)*ABin_t;
t2=lamb_x*phi*gamma_t*ABin_t;
tly1=t1*lamb_yt;
lyt1=lamb_y*t1;
txx=j(3,3,0);
do j=1 to 2;
	tt=j(p_y,3,0);
	tt[j+1,1]=1;
	tyy=tt*tly1+lyt1*(tt`);
	txy=t2*tt`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dlamb_y[,j]=tttt;
end;
do j=3 to 4;
	tt=j(p_y,3,0);
	tt[j+2,2]=1;
	tyy=tt*tly1+lyt1*(tt`);
	txy=t2*tt`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dlamb_y[,j]=tttt;
end;
do j=5 to 6;
	tt=j(p_y,3,0);
	tt[j+3,3]=1;
	tyy=tt*tly1+lyt1*(tt`);
	txy=t2*tt`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dlamb_y[,j]=tttt;
end;

*dgamma is the derivative with gamma;
dgamma=j(ps,2,0);
t1=lamb_y*ABin;
t2=lamb_x*phi;
txx=j(3,3,0);
do j=1 to 2;
	tt=j(3,1,0);
	tt[j]=1;
	tyy=t1*(tt*phi*gamma_t+gamma*phi*(tt`))*(t1`);
	txy=t2*tt`*t1`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
    dgamma[,j]=tttt;
end;

*dbeta is the derivative with beta;
dbeta=j(ps,2,0);
t1=ABin*(gamma*phi*gamma_t+sig_zeta)*ABin_t;
t2=lamb_x*phi*gamma_t*ABin_t;
t3=lamb_y*ABin;
txx=j(3,3,0);
do j=1 to 2;
	tt=j(3,3,0);
	tt[j+1,j]=1;
	tyy=t3*tt*t1*lamb_yt+lamb_y*t1*(tt`)*(t3`);
	txy=t2*(tt`)*(t3`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dbeta[,j]=tttt;
end;
*dphi is the derivative with phi;
*txx=lamb_x*lamb_t;
*t1=lamb_y*ABin*gamma; *tyy=t1*(t1`);
*txy=lamb_x*(t1`);
*ttt=(txx||txy)//((txy`)||tyy);
*run vech(ttt,tttt);
*dphi=tttt;

*dsiz_zeta is the derivative with siz_zeta;
dsig_zeta=j(ps,3,0);
t1=lamb_y*ABin;
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
do j=1 to 3;
	tt=j(3,3,0);
	tt[j,j]=1;
	tyy=t1*tt*(t1`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dsig_zeta[,j]=tttt;
end;

*dpsi_x is the derivative with psi_x;
dpsi_x=j(ps,p_x,0);
tyy=j(p_y,p_y,0);
txy=j(p_x,p_y,0);
do j=1 to p_x;
	tt=j(p_x,p_x,0);
	tt[j,j]=1;
	txx=tt;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dpsi_x[,j]=tttt;
end;

*dpsi_y is the derivative with psi_y;
dpsi_y=j(ps,p_y,0);
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
do j=1 to p_y;
	tt=j(p_y,p_y,0);
	tt[j,j]=1;
	tyy=tt;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dpsi_y[,j]=tttt;
end;

vdsig=dlamb_x||dlamb_y||dgamma||dbeta||dsig_zeta||dpsi_x||dpsi_y;

finish;
***--------------------------------------------------------------------***;
**** Subroutine to evaluate sigma0 and (dsigma0/dtheta);
*Model2 as represented by the solid & dashed arrows of Figure 6 of Deng & Yuan;
*\eta_1=gamma_11 \xi_1+zeta_1, 
\eta_2=gamma_21 \xi_1+ beta_21 \eta_1+zeta_2,
\eta_3=gamma_31 \xi_1+ beta_31 \eta_1+ beta_32 \eta_2+zeta_3;
***--------------------------------------------------------------------***;
start sigdsig_sem2(p,theta0,vsig,Sigma,vdsig);
p_x=3;
p_y=9;
ps=p*(p+1)/2;
lamb_x=theta0[1:3];
lamb_xt=lamb_x`;

lamb_y=j(p_y,3,0);
lamb_y[1,1]=1.0;
lamb_y[2:3,1]=theta0[4:5];
lamb_y[4,2]=1.0;
lamb_y[5:6,2]=theta0[6:7];
lamb_y[7,3]=1.0;
lamb_y[8:9,3]=theta0[8:9];
lamb_yt=lamb_y`;

gamma=theta0[10:12];
gamma_t=gamma`;
*
\eta_1=gamma_11 \xi_1+zeta_1, 
\eta_2=gamma_21 \xi_1+ beta_21 \eta_1+zeta_2,
\eta_3=gamma_31 \xi_1+ beta_31 \eta_1+ beta_32 \eta_2+zeta_3;

**AB=I-B in LISREL notation;
beta=j(3,3,0);
beta[2,1]=theta0[13];
beta[3,1]=theta0[14];
beta[3,2]=theta0[15];
A=i(3);
AB=A-beta;
*phi=theta0[14];
sig_zeta=diag(theta0[16:18]);
psi_x=diag(theta0[19:21]);
psi_y=diag(theta0[22:30]);

phi=1;
*computing Sigma based on structured parameters;
ABin=inv(AB);
ABin_t=ABin`;
sigyy=lamb_y*ABin*(gamma*phi*gamma_t+sig_zeta)*ABin_t*lamb_yt+psi_y;
sigxy=lamb_x*phi*gamma_t*ABin_t*lamb_yt;
sigxx=lamb_x*phi*lamb_xt+psi_x;
sigma=j(p,p,0);
sigma[1:p_x,1:p_x]=sigxx;
sigma[1:p_x,(p_x+1):p]=sigxy;
sigma[(p_x+1):p,1:p_x]=sigxy`;
sigma[(p_x+1):p,(p_x+1):p]=sigyy;

run vech(sigma,vsig);

***** The following is to evaluate (dsigma0/dtheta);

*dlamb_x is the derivative with Lamb_x;
dlamb_x=j(ps,p_x,0);
t1=phi*gamma_t*ABin_t*lamb_yt;
phl=phi*lamb_xt;
lph=lamb_x*phi;
tyy=j(p_y,p_y,0);
do j=1 to p_x;
	tt=j(3,1,0);
	tt[j]=1;
	txx=tt*phl+lph*tt`;
	txy=tt*t1;
	ttt=(txx||txy)//(txy`||tyy);
	run vech(ttt,tttt);
	dlamb_x[,j]=tttt;
end;


*dlamb_y is the derivative with Lamb_y;
dlamb_y=j(ps,6,0);
t1=ABin*(gamma*phi*gamma_t+sig_zeta)*ABin_t;
t2=lamb_x*phi*gamma_t*ABin_t;
tly1=t1*lamb_yt;
lyt1=lamb_y*t1;
txx=j(3,3,0);
do j=1 to 2;
	tt=j(p_y,3,0);
	tt[j+1,1]=1;
	tyy=tt*tly1+lyt1*(tt`);
	txy=t2*tt`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dlamb_y[,j]=tttt;
end;
do j=3 to 4;
	tt=j(p_y,3,0);
	tt[j+2,2]=1;
	tyy=tt*tly1+lyt1*(tt`);
	txy=t2*tt`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dlamb_y[,j]=tttt;
end;
do j=5 to 6;
	tt=j(p_y,3,0);
	tt[j+3,3]=1;
	tyy=tt*tly1+lyt1*(tt`);
	txy=t2*tt`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dlamb_y[,j]=tttt;
end;

*dgamma is the derivative with gamma;
dgamma=j(ps,3,0);
t1=lamb_y*ABin;
t2=lamb_x*phi;
txx=j(3,3,0);
do j=1 to 3;
	tt=j(3,1,0);
	tt[j]=1;
	tyy=t1*(tt*phi*gamma_t+gamma*phi*(tt`))*(t1`);
	txy=t2*tt`*t1`;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
    dgamma[,j]=tttt;
end;

*dbeta is the derivative with beta;
dbeta=j(ps,3,0);
t1=ABin*(gamma*phi*gamma_t+sig_zeta)*ABin_t;
t2=lamb_x*phi*gamma_t*ABin_t;
t3=lamb_y*ABin;
txx=j(3,3,0);
do j=1 to 2;
	tt=j(3,3,0);
	tt[j+1,1]=1;
	tyy=t3*tt*t1*lamb_yt+lamb_y*t1*(tt`)*(t3`);
	txy=t2*(tt`)*(t3`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dbeta[,j]=tttt;
end;
j=3;
	tt=j(3,3,0);
	tt[3,2]=1;
	tyy=t3*tt*t1*lamb_yt+lamb_y*t1*(tt`)*(t3`);
	txy=t2*(tt`)*(t3`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dbeta[,j]=tttt;
*dphi is the derivative with phi;
*txx=lamb_x*lamb_t;
*t1=lamb_y*ABin*gamma; 
*tyy=t1*(t1`);
*txy=lamb_x*(t1`);
*ttt=(txx||txy)//((txy`)||tyy);
*run vech(ttt,tttt);
*dphi=tttt;

*dsiz_zeta is the derivative with siz_zeta;
dsig_zeta=j(ps,3,0);
t1=lamb_y*ABin;
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
do j=1 to 3;
	tt=j(3,3,0);
	tt[j,j]=1;
	tyy=t1*tt*(t1`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dsig_zeta[,j]=tttt;
end;

*dpsi_x is the derivative with psi_x;
dpsi_x=j(ps,p_x,0);
tyy=j(p_y,p_y,0);
txy=j(p_x,p_y,0);
do j=1 to p_x;
	tt=j(p_x,p_x,0);
	tt[j,j]=1;
	txx=tt;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dpsi_x[,j]=tttt;
end;

*dpsi_y is the derivative with psi_y;
dpsi_y=j(ps,p_y,0);
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
do j=1 to p_y;
	tt=j(p_y,p_y,0);
	tt[j,j]=1;
	tyy=tt;
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dpsi_y[,j]=tttt;
end;

vdsig=dlamb_x||dlamb_y||dgamma||dbeta||dsig_zeta||dpsi_x||dpsi_y;

finish;

***--------------------------------------------------------------------***;
**Subroutine for parameter estimation by minimizing the F_ml function using Fisher-scoring method**;
**coefficients of factor-scores are subsequently obtained;
*Model1 as represented by the solid arrows of Figure 6 of Deng & Yuan;
***--------------------------------------------------------------------***;
start ML_SEM1(n,p,dup,theta0,scov,T_ml,htheta_z,coef_xt,coef_yt,n_div);
ep=.00001;
dtheta=1;
n_div=1;
do jj=1 to 300 while (max(abs(dtheta))>ep);
	run sigdsig_sem1(p,theta0,vsig0,Sigma0, vdsig);
  	sig_in=inv(Sigma0);
	** weight given by normal theory ;
	weight=0.5*dup`*(sig_in@sig_in)*dup;
	vdsig_t=vdsig`;
	dsw=vdsig_t*weight;
	dswds=dsw*vdsig;
	dwd_in=inv(dswds);
	cov_res=scov-sigma0;
	run vech(cov_res,vcov_res);
	dtheta=dwd_in*dsw*vcov_res;
	theta0=theta0+dtheta;
end;
***** Test start here;
if jj<300 then do;
	n_div=0;
	Ssig_in=Scov*sig_in;
	f_ml=trace(ssig_in)-log(det(ssig_in))-p;
	T_ml=(n-1)*f_ml;
	SE_hat=sqrt(vecdiag(dwd_in)/(n-1));
	*SE_hat=sqrt(vecdiag(dwd_in)/n);
	z_score=theta0/SE_hat;
	htheta_z=theta0||SE_hat||z_score;
*R-square for the structural model;
	*gamma=theta0[6:11];
	*print "gamma_s=" gamma_s;
	*print "T_ml=" T_ml;
*computation of factor-score;
	p_x=3;p_y=9;
	lamb_x=theta0[1:3];
	lamb_xt=lamb_x`;
	lamb_y=j(p_y,3,0);
	lamb_y[1,1]=1.0;
	lamb_y[2:3,1]=theta0[4:5];
	lamb_y[4,2]=1.0;
	lamb_y[5:6,2]=theta0[6:7];
	lamb_y[7,3]=1.0;
	lamb_y[8:9,3]=theta0[8:9];
	lamb_yt=lamb_y`;
	psi_x=diag(theta0[17:19]);
	psi_y=diag(theta0[20:28]);
	psi_xin=inv(psi_x);
	psi_yin=inv(psi_y);
	coef_xt=inv(lamb_xt*psi_xin*lamb_x)*lamb_xt*psi_xin;
	coef_yt=inv(lamb_yt*psi_yin*lamb_y)*lamb_yt*psi_yin;
end;
finish;
***--------------------------------------------------------------------***;
**Subroutine for parameter estimation by minimizing the F_ml function using Fisher-scoring method**;
**coefficients of factor-scores are subsequently obtained;
*Model2 as represented by the solid & dashed arrows of Figure 6 of Deng & Yuan;
***--------------------------------------------------------------------***;
start ML_SEM2(n,p,dup,theta0,scov,T_ml,htheta_z,coef_xt,coef_yt,n_div);
ep=.00001;
dtheta=1;
n_div=1;
do jj=1 to 300 while (max(abs(dtheta))>ep);
	run sigdsig_sem2(p,theta0,vsig0,Sigma0, vdsig);
  	sig_in=inv(Sigma0);
	** weight given by normal theory ;
	weight=0.5*dup`*(sig_in@sig_in)*dup;
	vdsig_t=vdsig`;
	dsw=vdsig_t*weight;
	dswds=dsw*vdsig;
	dwd_in=inv(dswds);
	cov_res=scov-sigma0;
	run vech(cov_res,vcov_res);
	dtheta=dwd_in*dsw*vcov_res;
	theta0=theta0+dtheta;
end;
***** Test start here;
if jj<300 then do;
	n_div=0;
	Ssig_in=Scov*sig_in;
	f_ml=trace(ssig_in)-log(det(ssig_in))-p;
	T_ml=(n-1)*f_ml;
	SE_hat=sqrt(vecdiag(dwd_in)/(n-1));
	*SE_hat=sqrt(vecdiag(dwd_in)/n);
	z_score=theta0/SE_hat;
	htheta_z=theta0||SE_hat||z_score;
*R-square for the structural model;
	*gamma=theta0[6:11];
	*print "gamma_s=" gamma_s;
	*print "T_ml=" T_ml;
*computation of factor-score;
	p_x=3;p_y=9;
*computation of factor-score;
	p_x=3;p_y=9;
	lamb_x=theta0[1:3];
	lamb_xt=lamb_x`;
	lamb_y=j(p_y,3,0);
	lamb_y[1,1]=1.0;
	lamb_y[2:3,1]=theta0[4:5];
	lamb_y[4,2]=1.0;
	lamb_y[5:6,2]=theta0[6:7];
	lamb_y[7,3]=1.0;
	lamb_y[8:9,3]=theta0[8:9];
	lamb_yt=lamb_y`;
	psi_x=diag(theta0[19:21]);
	psi_y=diag(theta0[22:30]);
	psi_xin=inv(psi_x);
	psi_yin=inv(psi_y);
	coef_xt=inv(lamb_xt*psi_xin*lamb_x)*lamb_xt*psi_xin;
	coef_yt=inv(lamb_yt*psi_yin*lamb_y)*lamb_yt*psi_yin;
end;
finish;
*--------------------------------------------;
*LS regression using the sample covariance matrix;
*Model1 as represented by the solid arrows of Figure 6 of Deng & Yuan;
*--------------------------------------------;
start regress1(n,c_xt,c_yt,scov,reg_est, Rsq);
p_x=3;p_y=9;
p=p_x+p_y;
scov=(n-1)*scov/n;
scov_xx=scov[1:p_x,1:p_x];
scov_xy=scov[1:p_x,(p_x+1):p];
scov_yx=scov_xy`;
scov_yy=scov[(p_x+1):p,(p_x+1):p];
c_x=c_xt`; c_y=c_yt`;
s_xx=c_xt*scov_xx*c_x;
s_xy=c_xt*scov_xy*c_y;
s_yy=c_yt*scov_yy*c_y;

hg_1=s_xy[1]/s_xx;
hsig2_1=n*(s_yy[1,1]-s_xy[1]*s_xy[1]/s_xx)/(n-2);
Varg_1=hsig2_1/(n*s_xx); SEg_1=sqrt(Varg_1);
zg_1=hg_1/SEg_1;
Rsq_1=1-(s_yy[1,1]-s_xy[1]*s_xy[1]/s_xx)/s_yy[1,1];
*----------------;
Cmat_2=(s_xx||s_xy[1])//(s_xy[1]||s_yy[1,1]);
cvec_2=s_xy[2]//s_yy[2,1];
cvec_2t=cvec_2`;
hgb_2=inv(Cmat_2)*cvec_2;
hsig2_2=n*(s_yy[2,2]-cvec_2t*inv(Cmat_2)*cvec_2)/(n-3);
Covgb_2=hsig2_2*inv(n*Cmat_2); 
SEgb_2=sqrt(vecdiag(Covgb_2));
zgb_2=hgb_2/SEgb_2;
Rsq_2=1-(s_yy[2,2]-cvec_2t*inv(Cmat_2)*cvec_2)/s_yy[2,2];
*----------------;
hb_3=s_yy[2,3]/s_yy[2,2];
hsig2_3=n*(s_yy[3,3]-s_yy[3,2]*s_yy[2,3]/s_yy[2,2])/(n-2);
Varb_3=hsig2_3/(n*s_yy[2,2]); 
SEb_3=sqrt(Varb_3);
zb_3=hb_3/SEb_3;
Rsq_3=1-(s_yy[3,3]-s_yy[3,2]*s_yy[2,3]/s_yy[2,2])/s_yy[3,3];

*----------------;

reg_est=j(4,3,0);

reg_est[1,]=hg_1||SEg_1||zg_1;
reg_est[2:3,]=(hgb_2||SEgb_2||zgb_2);
reg_est[4,]=hb_3||SEb_3||zb_3;
*print reg_est;
Rsq=Rsq_1||Rsq_2||Rsq_3;
finish;
*--------------------------------------------;
*LS regression using the sample covariance matrix;
*Model2 as represented by the solid & dashed arrows of Figure 6 of Deng & Yuan;
***--------------------------------------------------------------------***;
start regress2(n,c_xt,c_yt,scov,reg_est, Rsq);
p_x=3;p_y=9;
p=p_x+p_y;
scov=(n-1)*scov/n;
scov_xx=scov[1:p_x,1:p_x];
scov_xy=scov[1:p_x,(p_x+1):p];
scov_yx=scov_xy`;
scov_yy=scov[(p_x+1):p,(p_x+1):p];
c_x=c_xt`; c_y=c_yt`;
s_xx=c_xt*scov_xx*c_x;
s_xy=c_xt*scov_xy*c_y;
s_yx=s_xy`;
s_yy=c_yt*scov_yy*c_y;
*print "sxy syx=" s_xy s_yx;

hg_1=s_xy[1]/s_xx;
hsig2_1=n*(s_yy[1,1]-s_xy[1]*s_xy[1]/s_xx)/(n-2);
Varg_1=hsig2_1/(n*s_xx); SEg_1=sqrt(Varg_1);
zg_1=hg_1/SEg_1;
Rsq_1=1-(s_yy[1,1]-s_xy[1]*s_xy[1]/s_xx)/s_yy[1,1];
*----------------;
Cmat_2=(s_xx||s_xy[1])//(s_xy[1]||s_yy[1,1]);
cvec_2=s_xy[2]//s_yy[2,1];
cvec_2t=cvec_2`;
hgb_2=inv(Cmat_2)*cvec_2;
hsig2_2=n*(s_yy[2,2]-cvec_2t*inv(Cmat_2)*cvec_2)/(n-3);
Covgb_2=hsig2_2*inv(n*Cmat_2); 
SEgb_2=sqrt(vecdiag(Covgb_2));
zgb_2=hgb_2/SEgb_2;
Rsq_2=1-(s_yy[2,2]-cvec_2t*inv(Cmat_2)*cvec_2)/s_yy[2,2];
*----------------;
Cmat_3=(s_xx||s_xy[1,1:2])//(s_yx[1:2,1]||s_yy[1:2,1:2]);
cvec_3=s_xy[3]//s_yy[1:2,3];
cvec_3t=cvec_3`;
hgb_3=inv(Cmat_3)*cvec_3;
hsig2_3=n*(s_yy[3,3]-cvec_3t*inv(Cmat_3)*cvec_3)/(n-4);
Covgb_3=hsig2_3*inv(n*Cmat_3); 
SEgb_3=sqrt(vecdiag(Covgb_3));
zgb_3=hgb_3/SEgb_3;
Rsq_3=1-(s_yy[3,3]-cvec_3t*inv(Cmat_3)*cvec_3)/s_yy[3,3];

reg_est=j(6,3,0);
reg_est[1,]=hg_1||SEg_1||zg_1;
reg_est[2,]=(hgb_2[1]||SEgb_2[1]||zgb_2[1]);
reg_est[3,]=(hgb_3[1]||SEgb_3[1]||zgb_3[1]);
reg_est[4,]=(hgb_2[2]||SEgb_2[2]||zgb_2[2]);
reg_est[5:6,]=(hgb_3[2:3]||SEgb_3[2:3]||zgb_3[2:3]);

Rsq=Rsq_1||Rsq_2||Rsq_3;
finish;
*--------------------------------------------;
*Transforming weight of model A to B using structured 
Sigma=h*w*w'+Psi;
*--------------------------------------------;
start ssw(w,S,w_ss);
p=nrow(S);
w_t=w`;
c2=w_t*S*w;
d_s=vecdiag(S);
t1=c2-sum(w#w#d_s);
t2=ssq(w)*ssq(w)-sum(w#w#w#w);
h=t1/t2; 
psi=d_s-h*w#w;
*print "p=" p;
do j=1 to p;
	psi_j=psi[j];
	if psi_j<0 then do;
		print "Heywood case, j, psi_j, w=" j psi_j w;
		psi[j]=.05;
	end;
end;
w_s=w/psi;
w_ss=w_s/sqrt(w_s`*S*w_s);
*hlamb=sqrt(h)*w;
*hlamb2=hlamb#hlamb;
*maximal reliability coefficient;
*a_1=sum(hlamb2);
*a_12=a_1*a_1;
*a_2=sum(hlamb2#psi);
*rho_A=a_12/(a_12+a_2);
*b_1=sum(hlamb2/psi);
*rho_B=b_1/(1+b_1);

*print "w psi w_s, D_s, lamb, lamb2=" w psi w_s d_s hlamb hlamb2;
*print "PLS: p_j, rho_A, rho_B=" p rho_A rho_B;
finish;
******************************************************************;
*Subroutine conducting PLS-SEM mode A, and B_A;
*PLS-SEM mode A based on the sample covariance matrix;
*Model1 as represented by the solid arrows of Figure 6 of Deng & Yuan;
*--------------------------------------------;
start PLS_A1(n,scov0, reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
scov=(n-1)*scov0/n;
p_x=3;
p_y1=3;
p_y2=3;
p_y3=3;
s_xx=scov[1:3,1:3];s_xy1=scov[1:3,4:6];s_xy2=scov[1:3,7:9];s_xy3=scov[1:3,10:12];
s_y1x=s_xy1`;      s_y11=scov[4:6,4:6];s_y12=scov[4:6,7:9];s_y13=scov[4:6,10:12];
s_y2x=s_xy2`;      s_y21=s_y12`;       s_y22=scov[7:9,7:9];s_y23=scov[7:9,10:12];
s_y3x=s_xy3`;      s_y31=s_y13`;       s_y32=s_y23`;       s_y33=scov[10:12,10:12];

*--------------------;
*starting value of weights;
w_x=j(p_x,1,1);
c2_x=w_x`*s_xx*w_x;
w_x=w_x/sqrt(c2_x);

w_y1=j(p_y1,1,1);
c2_y1=w_y1`*s_y11*w_y1;
w_y1=w_y1/sqrt(c2_y1);

w_y2=j(p_y2,1,1);
c2_y2=w_y2`*s_y22*w_y2;
w_y2=w_y2/sqrt(c2_y2);

w_y3=j(p_y3,1,1);
c2_y3=w_y3`*s_y33*w_y3;
w_y3=w_y3/sqrt(c2_y3);

*initial composites;
*eta_x=x*w_x;
*eta_y1=y1*w_y1;
*eta_y2=y2*w_y2;
*eta_y3=y3*w_y3;


*environment variables;
r_xy1=w_x`*s_xy1*w_y1;
r_xy2=w_x`*s_xy2*w_y2;
r_y12=w_y1`*s_y12*w_y2;
r_y23=w_y2`*s_y23*w_y3;

*ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2;
*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2;
*ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
*ev_y3=sign(r_y23)*eta_y2;

*ev_x=ev_x/sqrt(ssq(ev_x)/n);
*ev_y1=ev_y1/sqrt(ssq(ev_y1)/n);
*ev_y2=ev_y2/sqrt(ssq(ev_y2)/n);
*ev_y3=ev_y3/sqrt(ssq(ev_y3)/n);

*w_x0=w_x; *w_y10=w_y1; *w_y20=w_y2;
w_0=w_x//w_y1//w_y2//w_y3;
*print "-------------mode A-------------";
dw=1; ep=.00001;
do j=1 to 50 while (max(abs(dw))>ep);
	*ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2;
	w_x=sign(r_xy1)*s_xy1*w_y1+sign(r_xy2)*s_xy2*w_y2;
	*w_x is the vector of regression coefficients eta_x-->x, var(eta_x)=1;
	c2_x=w_x`*s_xx*w_x;
	w_x=sign(w_x[1])*w_x/sqrt(c2_x);
	*eta_x=x*w_x; *eta_x is the vector of predicted score;

	*w_y1=y1_t*ev_y1/n, *ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2;
	w_y1=sign(r_xy1)*(s_y1x)*w_x+sign(r_y12)*s_y12*w_y2;
	c2_y1=w_y1`*s_y11*w_y1;
	w_y1=sign(w_y1[1])*w_y1/sqrt(c2_y1);
	
	*w_y2=y2_t*ev_y2/n, ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
	w_y2=sign(r_xy2)*s_y2x*w_x+sign(r_y12)*s_y21*w_y1+sign(r_y23)*s_y23*w_y3;
	c2_y2=w_y2`*s_y22*w_y2;
	w_y2=sign(w_y2[1])*w_y2/sqrt(c2_y2);

	*w_y3=y3_t*ev_y3/n; 
	w_y3=sign(r_y23)*s_y32*w_y2;
	c2_y3=w_y3`*s_y33*w_y3;
	w_y3=sign(w_y3[1])*w_y3/sqrt(c2_y3);
	*eta_y3=y3*w_y3; *eta_y3 is the vector of predicted score;

	*Note that for the corrected model, 
	the estimates do not depend the sign infront of eta_y or ev_y,
	how about misspecified models?;

	*environment variables;
	r_xy1=w_x`*s_xy1*w_y1;
	r_xy2=w_x`*s_xy2*w_y2;
	r_y12=w_y1`*s_y12*w_y2;
	r_y23=w_y2`*s_y23*w_y3;

	*ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2;
	*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2;
	*ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
	*ev_y3=sign(r_y23)*eta_y2;
	w_j=w_x//w_y1//w_y2//w_y3;
	dw=w_j-w_0;
	w_0=w_j;
end;
print "mode A model 1 iteration time=" j;
w_a=w_j;
w_xt=w_x`;
w_y=(w_y1||j(3,2,0))//(j(3,1,0)||w_y2||j(3,1,0))//(j(3,2,0)||w_y3);
w_yt=w_y`;	
scov=scov0;
run regress1(n,w_xt,w_yt,scov,reg_plsA, Rsq_A);

*print "--------mode A structured to mode B--------";
*print "--------block-x---------";
run ssw(w_x,S_xx,w_sx); 
*print "--------block-y1---------";
run ssw(w_y1,S_y11,w_sy1);
*print "--------block-y2---------";
run ssw(w_y2,S_y22,w_sy2);
run ssw(w_y3,S_y33,w_sy3);

w_sb=w_sx//w_sy1//w_sy2//w_sy3;
print "---------------------------------weights model 1---------------------------------";
print "model 1[mode A(w_x,w_y1,w_y2,w_y3)--mode AB(w_x,w_y1,w_y2,w_y3)]";
print w_x w_y1 w_y2 w_y3 w_sx w_sy1 w_sy2 w_sy3;
w_sxt=w_sx`;
w_sy=(w_sy1||j(3,2,0))//(j(3,1,0)||w_sy2||j(3,1,0))//(j(3,2,0)||w_sy3);
w_syt=w_sy`;	
scov=scov0;
run regress1(n,w_sxt,w_syt,scov,reg_plsBA, Rsq_BA);

finish;

******************************************************************;
*Subroutine conducting PLS-SEM mode A, and B_A;
*PLS-SEM mode A based on the sample covariance matrix;
*Model2 as represented by the solid & dashed arrows of Figure 6 of Deng & Yuan;
*--------------------------------------------;
start PLS_A2(n,scov0, reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
scov=(n-1)*scov0/n;
p_x=3;
p_y1=3;
p_y2=3;
p_y3=3;
s_xx=scov[1:3,1:3];s_xy1=scov[1:3,4:6];s_xy2=scov[1:3,7:9];s_xy3=scov[1:3,10:12];
s_y1x=s_xy1`;      s_y11=scov[4:6,4:6];s_y12=scov[4:6,7:9];s_y13=scov[4:6,10:12];
s_y2x=s_xy2`;      s_y21=s_y12`;       s_y22=scov[7:9,7:9];s_y23=scov[7:9,10:12];
s_y3x=s_xy3`;      s_y31=s_y13`;       s_y32=s_y23`;       s_y33=scov[10:12,10:12];

*--------------------;
*starting value of weights;
w_x=j(p_x,1,1);
c2_x=w_x`*s_xx*w_x;
w_x=w_x/sqrt(c2_x);

w_y1=j(p_y1,1,1);
c2_y1=w_y1`*s_y11*w_y1;
w_y1=w_y1/sqrt(c2_y1);

w_y2=j(p_y2,1,1);
c2_y2=w_y2`*s_y22*w_y2;
w_y2=w_y2/sqrt(c2_y2);

w_y3=j(p_y3,1,1);
c2_y3=w_y3`*s_y33*w_y3;
w_y3=w_y3/sqrt(c2_y3);

*initial composites;
*eta_x=x*w_x;
*eta_y1=y1*w_y1;
*eta_y2=y2*w_y2;
*eta_y3=y3*w_y3;


*environment variables;
r_xy1=w_x`*s_xy1*w_y1;
r_xy2=w_x`*s_xy2*w_y2;
r_xy3=w_x`*s_xy3*w_y3;
r_y12=w_y1`*s_y12*w_y2;
r_y13=w_y1`*s_y13*w_y3;
r_y23=w_y2`*s_y23*w_y3;

*ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2+sign(r_xy3)*eta_y3;
*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2+sign(r_y13)*eta_y3;
*ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
*ev_y3=sign(r_xy3)*eta_x+sign(r_y13)*eta_y1+sign(r_y23)*eta_y2;

*ev_y1=ev_y1/sqrt(ssq(ev_y1)/n);
*ev_y2=ev_y2/sqrt(ssq(ev_y2)/n);
*ev_y3=ev_y3/sqrt(ssq(ev_y3)/n);


*w_x0=w_x; *w_y10=w_y1; *w_y20=w_y2;
w_0=w_x//w_y1//w_y2//w_y3;
*print "-------------mode A-------------";
dw=1; ep=.00001;
do j=1 to 50 while (max(abs(dw))>ep);
	*ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2+sign(r_xy3)*eta_y3;
	w_x=sign(r_xy1)*s_xy1*w_y1+sign(r_xy2)*s_xy2*w_y2+sign(r_xy3)*s_xy3*w_y3;
	*w_x is the vector of regression coefficients eta_y1-->x, var(eta_y1)=1;
	c2_x=w_x`*s_xx*w_x;
	w_x=sign(w_x[1])*w_x/sqrt(c2_x);
	*eta_x=x*w_x; *eta_x is the vector of predicted score;

	*w_y1=y1_t*ev_y1/n, ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2+sign(r_y13)*eta_y3;
	w_y1=sign(r_xy1)*(s_y1x)*w_x+sign(r_y12)*s_y12*w_y2+sign(r_y13)*s_y13*w_y3;
	c2_y1=w_y1`*s_y11*w_y1;
	w_y1=sign(w_y1[1])*w_y1/sqrt(c2_y1);
	
	*w_y2=y2_t*ev_y2/n, ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
	w_y2=sign(r_xy2)*s_y2x*w_x+sign(r_y12)*s_y21*w_y1+sign(r_y23)*s_y23*w_y3;
	c2_y2=w_y2`*s_y22*w_y2;
	w_y2=sign(w_y2[1])*w_y2/sqrt(c2_y2);

	*w_y3=y3_t*ev_y3/n, ev_y3=sign(r_xy3)*eta_x+sign(r_y13)*eta_y1+sign(r_y23)*eta_y2;
	w_y3=sign(r_xy3)*s_y3x*w_x+sign(r_y13)*s_y31*w_y1+sign(r_y23)*s_y32*w_y2;
	c2_y3=w_y3`*s_y33*w_y3;
	w_y3=sign(w_y3[1])*w_y3/sqrt(c2_y3);
	*eta_y3=y3*w_y3; *eta_y3 is the vector of predicted score;

	*Note that for the corrected model, 
	the estimates do not depend the sign infront of eta_y or ev_y,
	how about misspecified models?;

	*environment variables;
	r_xy1=w_x`*s_xy1*w_y1;
	r_xy2=w_x`*s_xy2*w_y2;
	r_xy3=w_x`*s_xy3*w_y3;
	r_y12=w_y1`*s_y12*w_y2;
	r_y13=w_y1`*s_y13*w_y3;
	r_y23=w_y2`*s_y23*w_y3;

	*ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2+sign(r_xy3)*eta_y3;
	*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2+sign(r_y13)*eta_y3;
	*ev_y2=sign(r_xy2)*eta_x+sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
	*ev_y3=sign(r_xy3)*eta_x+sign(r_y13)*eta_y1+sign(r_y23)*eta_y2;

	w_j=w_x//w_y1//w_y2//w_y3;
	dw=w_j-w_0;
	w_0=w_j;
end;
print "mode A model 2 iteration time=" j;
w_a=w_j;
w_xt=w_x`;
w_y=(w_y1||j(3,2,0))//(j(3,1,0)||w_y2||j(3,1,0))//(j(3,2,0)||w_y3);
w_yt=w_y`;	

scov=scov0;
run regress2(n,w_xt,w_yt,scov,reg_plsA, Rsq_A);

*print "--------mode A structured to mode B--------";
*print "--------block-x---------";
run ssw(w_x,S_xx,w_sx); 
*print "--------block-y1---------";
run ssw(w_y1,S_y11,w_sy1);
*print "--------block-y2---------";
run ssw(w_y2,S_y22,w_sy2);
run ssw(w_y3,S_y33,w_sy3);

w_sb=w_sx//w_sy1//w_sy2//w_sy3;
*print "---------------------------------weights model 2---------------------------------";
*print "model 2 [mode A(w_x,w_y1,w_y2,w_y3)--mode B_A(w_x,w_y1,w_y2,w_y3)]";
*print w_x w_y1 w_y2 w_y3 w_sx w_sy1 w_sy2 w_sy3;
w_sxt=w_sx`;
w_sy=(w_sy1||j(3,2,0))//(j(3,1,0)||w_sy2||j(3,1,0))//(j(3,2,0)||w_sy3);
w_syt=w_sy`;	
scov=scov0;
run regress2(n,w_sxt,w_syt,scov,reg_plsBA, Rsq_BA);

finish;

***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;
print "--------------------------------------------------------------------";
print "Weston & Gore Jr 2006, 12-variable, 4-factor, see Figure 2 and Figure 4 of the article";
print "---------------------------------------------------------";
p=12; N=403;
ps=p*(p+1)/2;
scov_L={
4.811 
3.909 5.074 
3.179 3.041 3.503 
2.976 2.931 2.402 4.590 
2.743 2.867 2.418 3.626 3.860 
2.764 2.821 2.410 3.668 3.597 4.237 
3.260 3.592 2.279 3.413 2.926 2.852 7.511 
3.070 3.665 2.394 2.518 2.840 2.910 5.505 7.220 
2.596 2.817 2.553 2.630 2.671 2.618 5.773 5.716 6.958
0.257 0.224 0.236 0.260 0.225 0.260 0.243 0.281 0.256 0.072 
0.336 0.297 0.269 0.288 0.278 0.286 0.299 0.299 0.270 0.037 0.064 
0.280 0.239 0.198 0.230 0.239 0.280 0.177 0.296 0.209 0.039 0.038 0.078
};
Scov0=j(p,p,0);
count=0;
do i=1 to p;
	do j=1 to i;
 		count=count+1;
 		Scov0[i,j]=scov_L[count];
		Scov0[j,i]=scov_L[count];
	end;
end;
print scov0[format=f6.3];
*print Scov0;
print "----------------------------";

*re-order the variables;
set={7,8,9,4,5,6,10,11,12,1,2,3};
scov_0=scov0[set,set];
*print "-----------------------------";
*print scov_0[format=f6.3];
*print "-----------------------------";
**starting values for ML estimation;
lambda={
2.381,
2.365,
2.403, 
 .976,
 .993,
1.144,
1.011,
 .963,
 .795};
gamma={1.186, .046};
beta={.057, 10.361};
sig_zeta={
2.302,
 .008,
 .967
};
 
psi={ 
1.840,
1.628,
1.186,
 .882,
 .325,
 .580,
 .044,
 .027,
 .049,
 .794,
1.350,
 .962};

a0_sem1=lambda//gamma//beta//sig_zeta//psi; **starting values for ML estimation model 1;
a0_sem2=lambda//gamma//{.209}//{.378}//beta//sig_zeta//psi; **starting values for ML estimation model 2;

run Dp(p,dup);

theta0=a0_sem1;
q1=nrow(a0_sem1);
df1=ps-q1;
print "N,p=" n p;
print "df_1=" df1;
scov=scov_0;
run ML_sem1(n,p,dup,theta0,scov,T_ml,htheta1_z,coef1_semx,coef1_semy,n1_div);
if n1_div=1 then do;
	print "Fisher-scoring did not converge within 300 iterations under ML-SEM 1";
end;
else do;
	p_value=1-probchi(T_ml,df1);
	print "(SEM model 1) T_ml, df, p_value=" T_ml df1 p_value;
	print "htheta--SE--Z";
	print  htheta1_z[format=f8.3];
	hgamma1_z=htheta1_z[10:13,];

	*------------------R-square under CB-SEM------------------;
	gamma_11=theta0[10];
	gamma_21=theta0[11];
	beta_21=theta0[12];
	beta_32=theta0[13];
	*sig_zeta=diag(theta0[14:16]);

	*eta1=gamma_11*xi_1+zeta_1;
	Var_eta1=gamma_11*gamma_11+theta0[14];
	Rsq_cb1=1-theta0[14]/Var_eta1;
	*eta2=gamma_21*xi_1+beta_21*eta_1+zeta_2
	=gamma_21*xi_1+beta_21*(gamma_11*xi_1+zeta_1)+zeta_2
	=gamma_21*xi_1+beta_21*gamma_11*xi_1+beta_21*zeta_1+zeta_2
	=(gamma_21+beta_21*gamma_11)*xi_1+beta_21*zeta_1+zeta_2;
	beta_2=gamma_21+beta_21*gamma_11;
	Var_eta2=beta_2*beta_2+beta_21*beta_21*theta0[14]+theta0[15];
	Rsq_cb2=1-theta0[15]/Var_eta2;
	*eta3=beta_32*eta_2+zeta_3;
	Var_eta3=beta_32*beta_32*Var_eta2+theta0[16];
	Rsq_cb3=1-theta0[16]/Var_eta3;
	Rsq_cbm1=Rsq_cb1||Rsq_cb2||Rsq_cb3;
*---------------------------------------------------------;
	scov=scov_0;
	run regress1(n,coef1_semx,coef1_semy,scov,reg_fsr1, Rsq_fsr1);
	*----------EWC reg-------------;
	lamb_x=theta0[1:3];
	lamb_y1=1.0//theta0[4:5];
	lamb_y2=1.0//theta0[6:7];
	lamb_y3=1.0//theta0[8:9];
	c_x=j(3,1,1)/(j(1,3,1)*lamb_x);
	c_y1=j(3,1,1)/(j(1,3,1)*lamb_y1);
	c_y2=j(3,1,1)/(j(1,3,1)*lamb_y2);
	c_y3=j(3,1,1)/(j(1,3,1)*lamb_y3);
	c_xt=c_x`;
	c_y=(c_y1||j(3,2,0))//(j(3,1,0)||c_y2||j(3,1,0))//(j(3,2,0)||c_y3);
	c_yt=c_y`;
	scov=scov_0;
	run regress1(n,c_xt,c_yt,scov,reg_ewc1,Rsq_ewc1);
	print "(gamma_11//gamma_21//b_21//b_32)(hat\gamma, SE, z)";
	print "---CB-SEM1--- ---sem_FSR1--- ---EWCR1---";
	print hgamma1_z[format=f7.3] reg_fsr1[format=f7.3] reg_ewc1[format=f7.3];
*-------------PLS-SEM---------------;
	scov=scov_0;
	run PLS_A1(n,scov, reg_plsA1,reg_plsBA1, Rsq_A1, Rsq_BA1);
	print "---plsA1------plsB_A1---";
	print reg_plsA1[format=f7.3] reg_plsBA1[format=f7.3];

	Rsq_all1=(Rsq_cbm1//Rsq_fsr1//Rsq_ewc1//Rsq_A1//Rsq_BA1)`;
	print "R-square(CB-SEM,FSR,EWCR,PLSA,PLSBA)";
	print Rsq_all1[format=f6.3];

end;
print "------------------------------------------------";
print "------------------------------------------------";

print "-------------------------Model 2-------------------------";


theta0=a0_sem2;
q2=nrow(a0_sem2);
df2=ps-q2;
print "df_2=" df2;
scov=scov_0;
run ML_SEM2(n,p,dup,theta0,scov,T_ml,htheta2_z,coef2_semx,coef2_semy,n2_div);
if n2_div=1 then do;
	print "Fisher-scoring did not converge within 300 iterations under ML-SEM 2";
end;
else do;
	p_value=1-probchi(T_ml,df2);
	print "(SEM model 2) T_ml, df, p_value=" T_ml df2 p_value;
	
	print "htheta--SE--Z";
	print  htheta2_z[format=f8.3];
	hgamma2_z=htheta2_z[10:15,];
	*------------------R-square under CB-SEM------------------;
*
\eta_1=gamma_11 \xi_1+zeta_1, 
\eta_2=gamma_21 \xi_1+ beta_21 \eta_1+zeta_2,
\eta_3=gamma_31 \xi_1+ beta_31 \eta_1+ beta_32 \eta_2+zeta_3;

	gamma_11=theta0[10];
	gamma_21=theta0[11];
	gamma_31=theta0[12];
	beta_21=theta0[13];
	beta_31=theta0[14];
	beta_32=theta0[15];
	*sig_zeta=diag(theta0[16:18]);

	*eta1=gamma_11*xi_1+zeta_1;
	Var_eta1=gamma_11*gamma_11+theta0[16];
	Rsq_cb1=1-theta0[16]/Var_eta1;
	*eta2=gamma_21*xi_1+beta_21*eta_1+zeta_2
	=gamma_21*xi_1+beta_21*(gamma_11*xi_1+zeta_1)+zeta_2
	=gamma_21*xi_1+beta_21*gamma_11*xi_1+beta_21*zeta_1+zeta_2
	=(gamma_21+beta_21*gamma_11)*xi_1+beta_21*zeta_1+zeta_2;
	beta_2=gamma_21+beta_21*gamma_11;
	Var_eta2=beta_2*beta_2+beta_21*beta_21*theta0[16]+theta0[17];
	Rsq_cb2=1-theta0[17]/Var_eta2;

	*eta3=gamma_31*xi_1+beta_31*eta_1+beta_32*eta_2+zeta_3
	=gamma_31*xi_1+beta_31*(gamma_11*xi_1+zeta_1)
	+beta_32[(gamma_21+beta_21*gamma_11)*xi_1+beta_21*zeta_1+zeta_2]+zeta_3
	=gamma_31*xi_1+beta_31*gamma_11*xi_1+beta_31*zeta_1
	+beta_32*(gamma_21+beta_21*gamma_11)*xi_1
	+beta_32*beta_21*zeta_1+beta_32*zeta_2+zeta_3
	=[gamma_31+beta_31*gamma_11+beta_32*(gamma_21+beta_21*gamma_11)]*xi_1
	+(beta_31+beta_32*beta_21)*zeta_1+beta_32*zeta_2+zeta_3;

	beta_3=gamma_31+beta_31*gamma_11+beta_32*(gamma_21+beta_21*gamma_11);
	gamma_3=(beta_31+beta_32*beta_21);
	Var_eta3=beta_3*beta_3+gamma_3*gamma_3*theta0[16]
	+beta_32*beta_32*theta0[17]+theta0[18];
	Rsq_cb3=1-theta0[18]/Var_eta3;
	Rsq_cbm2=Rsq_cb1||Rsq_cb2||Rsq_cb3;
*---------------------------------------------------------;
	scov=scov_0;
	run regress2(n,coef2_semx,coef2_semy,scov,reg_fsr2, Rsq_fsr2);
	*----------EWC reg-------------;
	lamb_x=theta0[1:3];
	lamb_y1=1.0//theta0[4:5];
	lamb_y2=1.0//theta0[6:7];
	lamb_y3=1.0//theta0[8:9];
	c_x=j(3,1,1)/(j(1,3,1)*lamb_x);
	c_y1=j(3,1,1)/(j(1,3,1)*lamb_y1);
	c_y2=j(3,1,1)/(j(1,3,1)*lamb_y2);
	c_y3=j(3,1,1)/(j(1,3,1)*lamb_y3);
	c_xt=c_x`;
	c_y=(c_y1||j(3,2,0))//(j(3,1,0)||c_y2||j(3,1,0))//(j(3,2,0)||c_y3);
	c_yt=c_y`;
	scov=scov_0;
	run regress2(n,c_xt,c_yt,scov,reg_ewc2, Rsq_ewc2);
	print "(gamma_11//gamma_21//gamma_31//b_21//b_31//b_32)(hat\gamma, SE, z)";
	print "---CB-SEM2--- ---sem_FSR2--- ---EWCR2---";
	print hgamma2_z[format=f7.3] reg_fsr2[format=f7.3] reg_ewc2[format=f7.3];
	scov=scov_0;
	run PLS_A2(n,scov, reg_plsA2,reg_plsBA2,Rsq_A2, Rsq_BA2);
	print "---plsA2------plsB_A2---";
	print reg_plsA2[format=f7.3] reg_plsBA2[format=f7.3];

	Rsq_all2=(Rsq_cbm2//Rsq_fsr2//Rsq_ewc2//Rsq_A2//Rsq_BA2)`;
	print "R-square(CB-SEM,FSR,EWCR,PLSA,PLSBA)";
	print Rsq_all2[format=f6.3];
	
end;



finish;
run main;