***------------------------------------------------------------***;
*Data are from Table 10.1 on page 202 of Shumacker and Lomax 2010;
*p=9, four factors;
*The model is as in Figure 10.1 on page 196 of the book;
***------------------------------------------------------------***;
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;
***--------------------------------------------------------------------***;
**** The following subroutine is to evaluate sigma0 and (dsigma0/dtheta);
***--------------------------------------------------------------------***;
*The order of the factors in Kline 2011 is relabeled (page274);
*new F1= old F4;
*new F2=old F3;
*new F3= old F2;
*new F4=old F1;
*Kline model 1, Figure 10.3 without the dashed line;
*New model: F1-->F2, F2-->F3, F3-->F4, F2-->F4;
*F2=gamma_21 F1+D2, F3= beta_32 F2+D3, F4=beta_42 F2+beta_43 F3+D4;
*fy1        fx1     fy2        fy1     fy3        fy1        fy2;   

start sigdsig_sem(p,theta0,vsig,Sigma,vdsig);
p_x=5;
p_y=4;
ps=p*(p+1)/2;
lamb_x=j(p_x,2,0);
lamb_x[1:3,1]=theta0[1:3];
lamb_x[4:5,2]=theta0[4:5];
lamb_xt=lamb_x`;

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

gamma=j(2,2,0);
gamma=theta0[8:9]||theta0[10:11];
gamma_t=gamma`;

**AB=I-B in LISREL notation;
beta=j(2,2,0);
beta[2,1]=theta0[12];
A=i(2);
AB=A-beta;
phi=(1||theta0[13])//(theta0[13]||1);
sig_zeta=diag(theta0[14:15]);
psi_x=diag(theta0[16:20]);
psi_y=diag(theta0[21:24]);
*T_ml=57.167, df=21;
*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 3;
	tt=j(p_x,2,0);
	tt[j,1]=1;
	txx=tt*phl+lph*tt`;
	txy=tt*t1;
	ttt=(txx||txy)//(txy`||tyy);
	run vech(ttt,tttt);
	dlamb_x[,j]=tttt;
end;
do j=4 to 5;
	tt=j(p_x,2,0);
	tt[j,2]=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,2,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(p_x,p_x,0);
j=1;
	tt=j(p_y,2,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;
j=2;
	tt=j(p_y,2,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;

*dgamma is the derivative with gamma;
dgamma=j(ps,4,0);
t1=lamb_y*ABin;
t2=lamb_x*phi;
txx=j(p_x,p_x,0);
c=0;
do j=1 to 2;
do k=1 to 2;
	c=c+1;
	tt=j(2,2,0);
	tt[k,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[,c]=tttt;
end;
end;
*dbeta is the derivative with beta;
dbeta=j(ps,1,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(p_x,p_x,0);
j=1;
	tt=j(2,2,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;

*dphi is the derivative with phi;
t1=lamb_y*ABin*gamma;
	tt=j(2,2,0);
	tt[1,2]=1;
	tt[2,1]=1;
	txx=lamb_x*tt*lamb_xt;
	txy=lamb_x*tt*(t1`);
	tyy=t1*tt*(t1`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dphi=tttt;

*dsiz_zeta is the derivative with siz_zeta;
dsig_zeta=j(ps,2,0);
t1=lamb_y*ABin;
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
do j=1 to 2;
	tt=j(2,2,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||dphi||dsig_zeta||dpsi_x||dpsi_y;
vdsig=dlamb_x||dlamb_y||dgamma||dbeta||dphi||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;
***--------------------------------------------------------------------***;
start ML_SEM(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_sem(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=5;p_y=4;
	lamb_x=j(p_x,2,0);
	lamb_x[1:3,1]=theta0[1:3];
	lamb_x[4:5,2]=theta0[4:5];
	lamb_xt=lamb_x`;
	lamb_y=j(p_y,2,0);
	lamb_y[1,1]=1.0;
	lamb_y[2,1]=theta0[6];
	lamb_y[3,2]=1.0;
	lamb_y[4,2]=theta0[7];
	lamb_yt=lamb_y`;
	psi_x=diag(theta0[16:20]);
	psi_y=diag(theta0[21:24]);
	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;
*--------------------------------------------;
start regress(n,c_xt,c_yt,scov0,reg_est, Rsq);
*xy;
p_x=5;p_y=4;
p=p_x+p_y;
scov=(n-1)*scov0/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=inv(s_xx)*s_xy[,1]; *gamma_{11},gamma_{12};
hsig2_1=n*(s_yy[1,1]-(s_xy[,1]`)*inv(s_xx)*s_xy[,1])/(n-3);
Varg_1=hsig2_1*inv(n*s_xx); SEg_1=sqrt(vecdiag(Varg_1));
zg_1=hg_1/SEg_1;

Rsq_1=1-(s_yy[1,1]-(s_xy[,1]`)*inv(s_xx)*s_xy[,1])/s_yy[1,1];
*----------------;
Cmat_2=(s_xx||s_xy[,1])//(s_yx[1,]||s_yy[1,1]);
cvec_2=s_xy[,2]//s_yy[1,2];
cvec_2t=cvec_2`;
hgb_2=inv(Cmat_2)*cvec_2; *gamma_{21},gamma_{22},beta_{21};
hsig2_2=n*(s_yy[2,2]-cvec_2t*inv(Cmat_2)*cvec_2)/(n-4);
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];

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

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

Rsq=Rsq_1||Rsq_2;

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;
*--------------------------------------------;
start PLS_A(n,scov0, reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
scov=(n-1)*scov0/n;
p_x1=3;
p_x2=2;
p_y1=2;
p_y2=2;
s_x11=scov[1:3,1:3];s_x12=scov[1:3,4:5];s_x1y1=scov[1:3,6:7];s_x1y2=scov[1:3,8:9];
s_x21=s_x12`;       s_x22=scov[4:5,4:5];s_x2y1=scov[4:5,6:7];s_x2y2=scov[4:5,8:9];
s_y1x1=s_x1y1`;     s_y1x2=s_x2y1`;     s_y11=scov[6:7,6:7]; s_y12=scov[6:7,8:9];
s_y2x1=s_x1y2`;     s_y2x2=s_x2y2`;     s_y21=s_y12`;        s_y22=scov[8:9,8:9];

*--------------------;
*starting value of weights;
w_x1=j(p_x1,1,1);
c2_x1=w_x1`*s_x11*w_x1;
w_x1=w_x1/sqrt(c2_x1);

w_x2=j(p_x2,1,1);
c2_x2=w_x2`*s_x22*w_x2;
w_x2=w_x2/sqrt(c2_x2);

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);

*initial composites;
*eta_x1=x1*w_x1;
*eta_x2=x2*w_x2;
*eta_y1=y1*w_y1;
*eta_y2=y2*w_y2;


*environment variables;
r_x1y1=w_x1`*s_x1y1*w_y1;
r_x1y2=w_x1`*s_x1y2*w_y2;
r_x2y1=w_x2`*s_x2y1*w_y1;
r_x2y2=w_x2`*s_x2y2*w_y2;
r_y12=w_y1`*s_y12*w_y2;

*ev_x1=sign(r_x1y1)*eta_y1+sign(r_x1y2)*eta_y2;
*ev_x2=sign(r_x2y1)*eta_y1+sign(r_x2y2)*eta_y2;
*ev_y1=sign(r_x1y1)*eta_x1+sign(r_x2y1)*eta_x2+sign(r_y12)*eta_y2;
*ev_y2=sign(r_x1y2)*eta_x1+sign(r_x2y2)*eta_x2+sign(r_y12)*eta_y1;

*w_x0=w_x; *w_y10=w_y1; *w_y20=w_y2;
w_0=w_x1//w_x2//w_y1//w_y2;
*print "-------------mode A-------------";
dw=1; ep=.00001;
do j=1 to 50 while (max(abs(dw))>ep);
	*ev_x1=sign(r_x1y1)*eta_y1+sign(r_x1y2)*eta_y2;
	w_x1=sign(r_x1y1)*s_x1y1*w_y1+sign(r_x1y2)*s_x1y2*w_y2;
	*w_x1 is the vector of regression coefficients eta_x1-->x1, var(eta_x1)=1;
	c2_x1=w_x1`*s_x11*w_x1;
	w_x1=sign(w_x1[1])*w_x1/sqrt(c2_x1);
	*eta_x1=x1*w_x1; *eta_x1 is the vector of predicted score;

	*ev_x2=sign(r_x2y1)*eta_y1+sign(r_x2y2)*eta_y2;
	w_x2=sign(r_x2y1)*s_x2y1*w_y1+sign(r_x2y2)*s_x2y2*w_y2;
	*w_x2 is the vector of regression coefficients eta_x2-->x2, var(eta_x2)=1;
	c2_x2=w_x2`*s_x22*w_x2;
	w_x2=sign(w_x2[1])*w_x2/sqrt(c2_x2);
	*eta_x2=x2*w_x2; *eta_x2 is the vector of predicted score;

	*w_y1=y1_t*ev_y1/n,*ev_y1=sign(r_x1y1)*eta_x1+sign(r_x2y1)*eta_x2+sign(r_y12)*eta_y2;
	w_y1=sign(r_x1y1)*(s_y1x1)*w_x1+sign(r_x2y1)*s_y1x2*w_x2+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_x1y2)*eta_x1+sign(r_x2y2)*eta_x2+sign(r_y12)*eta_y1;
	w_y2=sign(r_x1y2)*s_y2x1*w_x1+sign(r_x2y2)*s_y2x2*w_x2+sign(r_y12)*s_y21*w_y1;
	c2_y2=w_y2`*s_y22*w_y2;
	w_y2=sign(w_y2[1])*w_y2/sqrt(c2_y2);

	*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_x1y1=w_x1`*s_x1y1*w_y1;
	r_x1y2=w_x1`*s_x1y2*w_y2;
	r_x2y1=w_x2`*s_x2y1*w_y1;
	r_x2y2=w_x2`*s_x2y2*w_y2;
	r_y12=w_y1`*s_y12*w_y2;

	*ev_x1=sign(r_x1y1)*eta_y1+sign(r_x1y2)*eta_y2;
	*ev_x2=sign(r_x2y1)*eta_y1+sign(r_x2y2)*eta_y2;
	*ev_y1=sign(r_x1y1)*eta_x1+sign(r_x2y1)*eta_x2+sign(r_y12)*eta_y2;
	*ev_y2=sign(r_x1y2)*eta_x1+sign(r_x2y2)*eta_x2+sign(r_y12)*eta_y1;

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

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

*print "--------mode A structured to mode B--------";
*print "--------block-x1---------";
run ssw(w_x1,S_x11,w_sx1);
*print "--------block-x2---------";
run ssw(w_x2,S_x22,w_sx2); 
*print "--------block-y1---------";
run ssw(w_y1,S_y11,w_sy1);
*print "--------block-y2---------";
run ssw(w_y2,S_y22,w_sy2);

w_sb=w_sx1//w_sx2//w_sy1//w_sy2;
*print "---------------------------------weights---------------------------------";
*print "mode A(w_x1,w_x2,w_y1,w_y2)--mode AB(w_x1,w_x2,w_y1,w_y2)";
*print w_x1 w_x2 w_y1 w_y2 w_sx1 w_sx2 w_sy1 w_sy2;
w_sx=(w_sx1||j(3,1,0))//(j(2,1,0)||w_sx2);
w_sxt=w_sx`;
w_sy=(w_sy1||j(2,1,0))//(j(2,1,0)||w_sy2);
w_syt=w_sy`;	
scov=scov0;
run regress(n,w_sxt,w_syt,scov,reg_plsBA, Rsq_BA);

finish;

***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;
print "--------------------------------------------------------------------";
print "The model is as in Figure 10.1 on page 196 of the book";
*Table 10.1 on page 202 of Shumacker and Lomax 2010;
*9-variable, 4-factor, see Figure 10.1
*Educational achievement data from Table 10.1 on page 202 of Shumacker and Lomax 2010 3rd ED.;
print "---------------------------------------------------------";
n=200;
p=9;
ps=p*(p+1)/2;
print "N, p=" N p;
p_x=5;
p_y=4;
scov=j(p,p,0);
scov_vec={
1.024
 .792 1.077
1.027  .919 1.844
 .756  .697 1.244 1.286
 .567  .537  .876  .632 .852
 .445  .424  .677  .526 .518 .670
 .434  .389  .635  .498 .475 .545 .716
 .580  .564  .893  .716 .546 .422 .373 .851
 .491  .499  .888  .646 .508 .389 .339 .629 .871};
 c=0;
do i=1 to p;
do j=1 to i;
	c=c+1;
	scov[i,j]=scov_vec[c];
	scov[j,i]=scov_vec[c];
end;
end;
*print "scov=" scov;

set={5,6,7,8,9,1,2,3,4}; *re-order the variables;
scov_a=scov[set,set];
*print "sample cov";
*print scov_a;
**get starting value using estimate from MLE;
a0={.730 .735 .703 .814 .772
.917 .759
.299 .177 .480 .611
.548
.726
.335 .225
.320 .130 .222 .188 .274
.160 .351 .205 .342};
a0=a0`; *starting values for ML estimation;
q=nrow(a0);
df=ps-q;
*print "df=" df;

run Dp(p,dup);

scov=scov_a;
theta0=a0;
run ML_sem(n,p,dup,theta0,scov,T_ml,htheta_z,coef_semx,coef_semy,n_div);

if n_div=1 then do;
	print "Fisher-scoring did not converge within 300 iterations under ML-SEM";
end;
else do;
	p_v=1-probchi(T_ml,df);
	print "T_ml, df, p_value=" T_ml df p_v;
	print "htheta--SE--Z";
	print htheta_z[format=f8.3];
	hgamma_z=htheta_z[8:12,];
	*------------------R-square under CB-SEM------------------;
*sig_zeta=diag(theta0[14:15]);

	gamma_11=theta0[8];
	gamma_21=theta0[9];
	gamma_12=theta0[10];
	gamma_22=theta0[11];
	beta_21=theta0[12];
	phi_21=theta0[13];

	*eta1=gamma_11*xi_1+gamma_12*xi_2+zeta_1;
	Var_eta1=gamma_11*gamma_11+gamma_12*gamma_12
		+2*gamma_11*gamma_12*phi_21+theta0[14];
	Rsq_cb1=1-theta0[14]/Var_eta1;
	*eta2=gamma_21*xi_1+gamma_22*xi_2+beta_21*eta_1+zeta_2
	=gamma_21*xi_1+gamma_22*xi_2
	+beta_21*(gamma_11*xi_1+gamma_12*xi_2+zeta_1)+zeta_2
	=gamma_21*xi_1+gamma_22*xi_2
	+beta_21*gamma_11*xi_1+beta_21*gamma_12*xi_2
	+beta_21*zeta_1+zeta_2
	=(gamma_21+beta_21*gamma_11)*xi_1
	+(gamma_22+beta_21*gamma_12)*xi_2
	+beta_21*zeta_1+zeta_2;
	gb_1=gamma_21+beta_21*gamma_11;
	gb_2=gamma_22+beta_21*gamma_12;
	Var_eta2=gb_1*gb_1+gb_2*gb_2
		+2*gb_1*gb_2*phi_21+beta_21*beta_21*theta0[14]+theta0[15];
	Rsq_cb2=1-theta0[15]/Var_eta2;
	Rsq_cb=Rsq_cb1||Rsq_cb2;
*---------------------------------------------------------;
	scov=scov_a;
	run regress(n,coef_semx,coef_semy,scov,reg_fsr, Rsq_fsr);
	*----------EWC reg-------------;
	lamb_x1=theta0[1:3];
	lamb_x2=theta0[4:5];
	lamb_y1=1.0//theta0[6];
	lamb_y2=1.0//theta0[7];
	c_x1=j(3,1,1)/(j(1,3,1)*lamb_x1);
	c_x2=j(2,1,1)/(j(1,2,1)*lamb_x2);
	c_y1=j(2,1,1)/(j(1,2,1)*lamb_y1);
	c_y2=j(2,1,1)/(j(1,2,1)*lamb_y2);
	c_x=(c_x1||j(3,1,0))//(j(2,1,0)||c_x2);
	c_xt=c_x`;
	c_y=(c_y1||j(2,1,0))//(j(2,1,0)||c_y2);
	c_yt=c_y`;
	scov=scov_a;
	run regress(n,c_xt,c_yt,scov,reg_ewc, Rsq_ewc);
	print "(gamma_11//gamma_21//gamma_12//gamma_22//beta_21)(hat\gamma, SE, z)";
	print "---CB-SEM--- ---FSR--- ---EWCR---"; 
	print hgamma_z[format=f7.3] reg_fsr[format=f7.3] reg_ewc[format=f7.3];
*-------------PLS-SEM---------------;
	scov=scov_a;
	run PLS_A(n,scov, reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
	print "---plsA------plsB_A---";
	print reg_plsA[format=f7.3] reg_plsBA[format=f7.3];
	Rsq_all=(Rsq_cb//Rsq_fsr//Rsq_ewc//Rsq_A//Rsq_BA)`;
	print "R-square(CB-SEM,FSR,EWCR,PLSA,PLSB_A)";
	print Rsq_all[format=f6.3];
end;
*------------------------------------------------;

finish;
run main;