***------------------------------------------------------------***;
*Kline 2016's data, the covariance matrix is obtined from the correlation matrix 
and standard deviations from Table 14.1 on page 342 of Kline SEM book 2016;
*p=12; *N=263;
***------------------------------------------------------------***;

options ls=125 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 5 of Deng & Yuan;
***--------------------------------------------------------------------***;
*The order of the factors in Kline 2016 (page342) is relabeled;
*new F1= old F4;
*new F2=old F3;
*new F3= old F2;
*new F4=old F1;
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]=theta0[10];
gamma_t=gamma`;

**AB=I-B in LISREL notation;
beta=j(3,3,0);
beta[2,1]=theta0[11];
beta[3,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;
t1=lamb_y*ABin;
t2=lamb_x*phi;
txx=j(3,3,0);
	tt=j(3,1,0);
	tt[1]=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=tttt;

*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 to evaluate sigma0 and (dsigma0/dtheta);
*Model2 as represented by the solid & dashed arrows of Figure 5 of Deng & Yuan;
***--------------------------------------------------------------------***;
*The order of the factors in Kline 2016 (page342) is relabeled;

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`;
*F2=gamma_21 F1+D2, F3=gamma_31 F1+ beta_32 F2+D3, 
F4=gamma_41 F1+beta_42 F2+beta_43 F3+D4;

**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 5 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 5 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 5 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];
*----------------;
hb_2=s_yy[2,1]/s_yy[1,1];
hsig2_2=n*(s_yy[2,2]-s_yy[2,1]*s_yy[1,2]/s_yy[1,1])/(n-2);
Varb_2=hsig2_2/(n*s_yy[1,1]); SEb_2=sqrt(Varb_2);
zb_2=hb_2/SEb_2;
Rsq_2=1-(s_yy[2,2]-s_yy[2,1]*s_yy[1,2]/s_yy[1,1])/s_yy[2,2];
*----------------;
hb_3=inv(s_yy[1:2,1:2])*s_yy[1:2,3];
hsig2_3=n*(s_yy[3,3]-s_yy[3,1:2]*inv(s_yy[1:2,1:2])*s_yy[1:2,3])/(n-3);
Covb_3=hsig2_3*inv(n*s_yy[1:2,1:2]); 
SEb_3=sqrt(vecdiag(Covb_3));
zb_3=hb_3/SEb_3;
Rsq_3=1-(s_yy[3,3]-s_yy[3,1:2]*inv(s_yy[1:2,1:2])*s_yy[1:2,3])/s_yy[3,3];

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

reg_est=j(4,3,0);

reg_est[1,]=hg_1||SEg_1||zg_1;
reg_est[2,]=(hb_2||SEb_2||zb_2);
reg_est[3: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 5 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 5 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_y12=w_y1`*s_y12*w_y2;
r_y13=w_y1`*s_y13*w_y3;
r_y23=w_y2`*s_y23*w_y3;

*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2+sign(r_y13)*eta_y3;
*ev_y2=sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
*ev_y3=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);

*check_y4=ev_y4`*ev_y4/n;
*print "check=1=" check_y4;

*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);
*eta_y1=y1*w_y1;
	w_x=s_xy1*w_y1;
	*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_y12)*eta_y1+sign(r_y23)*eta_y3;
	w_y2=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_y13)*eta_y1+sign(r_y23)*eta_y2;
	w_y3=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_y12=w_y1`*s_y12*w_y2;
	r_y13=w_y1`*s_y13*w_y3;
	r_y23=w_y2`*s_y23*w_y3;

	*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2+sign(r_y13)*eta_y3;
	*ev_y2=sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
	*ev_y3=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_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 5 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 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 regress2(n,w_sxt,w_syt,scov,reg_plsBA, Rsq_BA);

finish;

***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;
print "--------------------------------------------------------------------";
print "Kline's data, 12-variable, 4-factor, see Figures 14.1&14.2 on pages 346&348 of Kline 2016";
print "---------------------------------------------------------";

scov_0={
 0.881721 0.6379153 0.5587003 0.1387898 0.2069556 0.1018515 -0.113159 -0.160118 -0.177094 0.0534103 0.0643816 0.1062131,
0.6379153  1.034289 0.5708045 0.1491756 0.2434698 0.1305625  -0.10828 -0.149264 -0.181396 0.0694164  0.032007  0.177135,
0.5587003 0.5708045  0.877969 0.0863614 0.1758936 0.1134182 -0.106888 -0.135811 -0.126715 0.0626235 -0.036862 0.0553383,
0.1387898 0.1491756 0.0863614  0.315844 0.2075803 0.0739165  -0.10159 -0.117737  -0.10476 -0.006793 -0.036638 0.0354414,
0.2069556 0.2434698 0.1758936 0.2075803    0.5776 0.1788098 -0.118264 -0.141166 -0.141668 0.0815944 -0.043566 0.1049849,
0.1018515 0.1305625 0.1134182 0.0739165 0.1788098  0.274576 -0.043529 -0.073397 -0.082354 0.0525315 -0.001767 0.0230791,
-0.113159  -0.10828 -0.106888  -0.10159 -0.118264 -0.043529  0.342225 0.2682675 0.2369098 -0.030779 -0.026302 -0.069685,
-0.160118 -0.149264 -0.135811 -0.117737 -0.141166 -0.073397 0.2682675  0.370881 0.2613201 -0.048063 -0.027381 -0.044501,
-0.177094 -0.181396 -0.126715  -0.10476 -0.141668 -0.082354 0.2369098 0.2613201  0.534361 0.0083159  -0.01479 -0.061465,
0.0534103 0.0694164 0.0626235 -0.006793 0.0815944 0.0525315 -0.030779 -0.048063 0.0083159  0.505521 0.2269626 0.4006933,
0.0643816  0.032007 -0.036862 -0.036638 -0.043566 -0.001767 -0.026302 -0.027381  -0.01479 0.2269626  1.263376  0.426422,
0.1062131  0.177135 0.0553383 0.0354414 0.1049849 0.0230791 -0.069685 -0.044501 -0.061465 0.4006933  0.426422  1.002001 
};
set={10,11,12,7,8,9,4,5,6,1,2,3}; *re-order the variables;
scov_0=scov_0[set,set];
n=263;
p=ncol(scov_0);
ps=p*(p+1)/2;

a0_sem1={
  .464
 .492
 .864
1.126
 .991
1.766
 .812
1.031
 .892
-.065
-.332
-.259
 .907                             
 .233
 .091
 .473
 .290
1.021
 .256
 .105
 .070
 .301
 .199
 .213
 .198
 .260
 .374
 .383
};
a0_sem1=a0_sem1`; *starting values for ML estimation of model 1;
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;
	print "(CB-SEM model 1) T_ml=" T_ml;
	print "(CB-SEM model 1) htheta--SE--Z";
	print  htheta1_z[format=f8.3];
	hgamma1_z=htheta1_z[10:13,];

	*------------------R-square under CB-SEM------------------;
	gamma_11=theta0[10];
	beta_21=theta0[11];
	beta_31=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=beta_21*eta_1+zeta_2;
	Var_eta2=beta_21*beta_21*Var_eta1+theta0[15];
	Rsq_cb2=1-theta0[15]/Var_eta2;
	*eta3=beta_31*eta_1+beta_32*eta_2+zeta_3
	 =beta_31*eta_1+beta_32*(beta_21*eta_1+zeta_2)+zeta_3
	 =beta_31*eta_1+beta_32*beta_21*eta_1+beta_32*zeta_2+zeta_3
	 =(beta_31+beta_32*beta_21)*eta_1+beta_32*zeta_2+zeta_3;
	beta_3=beta_31+beta_32*beta_21;
	Var_eta3=beta_3*beta_3*Var_eta1
    +beta_32*beta_32*theta0[15]+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//b_21//b_31//b_32)(hat\gamma, SE, z)";
	print "---CB-SEM1--- ---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-------------------------";

a0_sem2={
  .464
 .492
 .864
1.126
 .991
1.766
 .812
1.031
 .892
 .1
 .1
-.065
-.332
-.259
 .907                             
 .233
 .091
 .473
 .290
1.021
 .256
 .105
 .070
 .301
 .199
 .213
 .198
 .260
 .374
 .383
};
a0_sem2=a0_sem2`; *starting values for ML estimation 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;
	print "(CB-SEM model 2) T_ml=" T_ml;
	print "(CB-SEM model 2)htheta--SE--Z";
	print  htheta2_z[format=f8.3];
	hgamma2_z=htheta2_z[10:15,];
	*------------------R-square under CB-SEM------------------;
	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_31//gamma_31//b_21//b_31//b_32)(hat\gamma, SE, z)";
	print "---CB-SEM2--- ---sem_FSR2--- ---EWC2---";
	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,PLSB_A)";
*	print Rsq_all2[format=f6.3];
end;

finish;
run main;