***------------------------------------------------------------***;
*Kline (1998, p.254) data, only a correlation matrix is available;
*p=12; *N=158;
***------------------------------------------------------------***;

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 model is the same as in Kline 1998 Figure 8.4 on
page 255 (structured SEM model);
*Model: F1-->F2, F2-->F3, F3-->F4;

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

lamb_y=j(p_y,3,0);
lamb_y[1,1]=1.0;
lamb_y[2:3,1]=theta0[3:4];
lamb_y[4,2]=1.0;
lamb_y[5:6,2]=theta0[5:6];
lamb_y[7,3]=1.0;
lamb_y[8:10,3]=theta0[7: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,2]=theta0[12];
A=i(3);
AB=A-beta;
sig_zeta=diag(theta0[13:15]);
psi_x=diag(theta0[16:17]);
psi_y=diag(theta0[18:27]);

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(p_x,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,7,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);
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 7;
	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(p_x,p_x,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,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(p_x,p_x,0);
j=1;
	tt=j(3,3,0);
	tt[2,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;
j=2;
	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||dphi||dsig_zeta||dpsi_x||dpsi_y;
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;
***--------------------------------------------------------------------***;
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=2;p_y=10;
	lamb_x=theta0[1:2];
	lamb_xt=lamb_x`;
	lamb_y=j(p_y,3,0);
	lamb_y[1,1]=1.0;
	lamb_y[2:3,1]=theta0[3:4];
	lamb_y[4,2]=1.0;
	lamb_y[5:6,2]=theta0[5:6];
	lamb_y[7,3]=1.0;
	lamb_y[8:10,3]=theta0[7:9];
	lamb_yt=lamb_y`;
	psi_x=diag(theta0[16:17]);
	psi_y=diag(theta0[18:27]);
	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,scov,reg_est, Rsq);
*xy;
p_x=2;p_y=10;
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=s_yy[3,2]/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(3,3,0);

reg_est[1,]=hg_1||SEg_1||zg_1;
reg_est[2,]=(hb_2||SEb_2||zb_2);
reg_est[3,]=hb_3||SEb_3||zb_3;
*print reg_est;
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;
*the model is as represented by Figure 4 of Deng & Yuan;
*originally from page 255 of Kline 1998;
*--------------------------------------------;
start PLS_A(n,scov0, reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
scov=(n-1)*scov0/n;
p_x=2;
p_y1=3;
p_y2=3;
p_y3=4;
s_xx=scov[1:2,1:2];s_xy1=scov[1:2,3:5];s_xy2=scov[1:2,6:8];s_xy3=scov[1:2,9:12];
s_y1x=s_xy1`;      s_y11=scov[3:5,3:5];s_y12=scov[3:5,6:8];s_y13=scov[3:5,9:12];
s_y2x=s_xy2`;      s_y21=s_y12`;       s_y22=scov[6:8,6:8];s_y23=scov[6:8,9:12];
s_y3x=s_xy3`;      s_y31=s_y13`;       s_y32=s_y23`;       s_y33=scov[9:12,9: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_y23=w_y2`*s_y23*w_y3;

*ev_y1=sign(r_xy1)*eta_x+sign(r_y12)*eta_y2;
*ev_y2=sign(r_y12)*eta_y1+sign(r_y23)*eta_y3;
*ev_y3=sign(r_y23)*eta_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;
	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_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_y23)*eta_y2;
	w_y3=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_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(4,2,0)||w_y3);
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-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 "[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(4,2,0)||w_sy3);
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 "Kline's data, 12-variable, 4-factor, see Table 8.1 on page 254 of Kline 1998";
print "---------------------------------------------------------";


scov_0={
 1.00   0.42  -0.43  -0.40  -0.35  -0.39  -0.24  -0.31  -0.25  -0.14  -0.25  -0.16,
 0.42   1.00  -0.50  -0.40  -0.38  -0.43  -0.37  -0.33  -0.25  -0.17  -0.26  -0.18,
-0.43  -0.50   1.00   0.66   0.67   0.78   0.69   0.63   0.49   0.18   0.42   0.33,
-0.40  -0.40   0.66   1.00   0.60   0.56   0.49   0.49   0.32   0.09   0.25   0.27,
-0.35  -0.38   0.67   0.60   1.00   0.73   0.70   0.72   0.58   0.17   0.46   0.35,
-0.39  -0.43   0.78   0.56   0.73   1.00   0.73   0.87   0.53   0.14   0.42   0.36,
-0.24  -0.37   0.69   0.49   0.70   0.73   1.00   0.72   0.60   0.15   0.44   0.38,
-0.31  -0.33   0.63   0.49   0.72   0.87   0.72   1.00   0.59   0.15   0.45   0.38,
-0.25  -0.25   0.49   0.32   0.58   0.53   0.60   0.59   1.00   0.25   0.77   0.59,
-0.14  -0.17   0.18   0.09   0.17   0.14   0.15   0.15   0.25   1.00   0.19  -0.29,
-0.25  -0.26   0.42   0.25   0.46   0.42   0.44   0.45   0.77   0.19   1.00   0.58,
-0.16  -0.18   0.33   0.27   0.35   0.36   0.38   0.38   0.59  -0.29   0.58   1.00
};
n=158;
p=ncol(scov_0);
ps=p*(p+1)/2;

a0_sem={
 .614
 .684
 .816
 .950
 .845
 .954
 .233
 .854
 .659
-.650
 .991
 .642                             
 .332
 .155
 .536
 .256
 .505
 .50
 .301
 .199
 .213
 .198
 .260
 .374
 .383
 .50
 .50
};
a0_sem=a0_sem`; *starting values for ML estimation;
run Dp(p,dup);

theta0=a0_sem;
q=nrow(a0_sem);
df=ps-q;
print "N,p=" n p;
print "df=" df;
scov=scov_0;
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 CB-SEM";
end;
else do;
	print "(CB-SEM) T_ml=" T_ml;
	print "(CB-SEM) htheta--SE--Z";
	print  htheta_z[format=f8.3];
	hgamma_z=htheta_z[10:12,];
	*------------------R-square under CB-SEM------------------;
	gamma_11=theta0[10];
	beta_21=theta0[11];
	beta_32=theta0[12];
	*sig_zeta=diag(theta0[13:15]);

	*eta1=gamma_11*xi_1+zeta_1;
	Var_eta1=gamma_11*gamma_11+theta0[13];
	Rsq_cb1=1-theta0[13]/Var_eta1;
	*eta2=beta_21*eta_1+zeta_2;
	Var_eta2=beta_21*beta_21*Var_eta1+theta0[14];
	Rsq_cb2=1-theta0[14]/Var_eta2;
	*eta3=beta_32*eta_2+zeta_3;
	Var_eta3=beta_32*beta_32*Var_eta2+theta0[15];
	Rsq_cb3=1-theta0[15]/Var_eta3;

	Rsq_cb=Rsq_cb1||Rsq_cb2||Rsq_cb3;
*---------------------------------------------------------;

	scov=scov_0;
	run regress(n,coef_semx,coef_semy,scov,reg_fsr, Rsq_fsr);
	*----------EWC reg-------------;
	lamb_x=theta0[1:2];
	lamb_y1=1.0//theta0[3:4];
	lamb_y2=1.0//theta0[5:6];
	lamb_y3=1.0//theta0[7:9];
	c_x=j(2,1,1)/(j(1,2,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(4,1,1)/(j(1,4,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(4,2,0)||c_y3);
	c_yt=c_y`;
	scov=scov_0;
	run regress(n,c_xt,c_yt,scov,reg_ewc, Rsq_ewc);
	print "(gamma_11//b_21//b_32)(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_0;
	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;