***------------------------------------------------------------***;
*Raw data are from 
Neumann, C.S. (1994). Structural equation modeling of symptoms of alcoholism and psychopathology.
Dissertation, University of Kansas, Lawrence. It is not publicly available
***------------------------------------------------------------***;
**input Neumann's data;
data neumann;
infile 'd:\data\alco.dat';
input v1-v10;
run;
*------------------------------------------------------------***;

options ls=100 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 by 1;
        do j=i to p by 1;
                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;
***--------------------------------------------------------------------***;
** Computing Mardia's multivariate skewness and kurtosis;
***--------------------------------------------------------------------***;

start Mardia(x,mskew0,mkurt0,mkurt_p,df);
N=nrow(x);
p=ncol(x);

x_t=x`;
x_bar=x_t*j(N,1,1)/N;
Scov=(x_t*x-N*x_bar*x_bar`)/N;
mkurt=0.0;
scov_in=inv(Scov);

do i=1 to N;
        x_i=x[i,]`;
        x_i0=x_i-x_bar;
		x_it0=x_i0`;
        d_i2=x_it0*scov_in*x_i0;
        mkurt=mkurt+d_i2*d_i2;
end;
mkurt=mkurt/N;
mkurt_p=mkurt/(p*(p+2));

mkurt0=(mkurt-p*(p+2))/sqrt(8*p*(p+2)/N);

mskew=0.0;
do i=1 to N;
        x_i=x[i,]`;
        x_i0=x_i-x_bar;
		x_it0=x_i0`;
        do j=1 to N;
                x_j=x[j,]`;
                x_j0=x_j-x_bar;

                d_ij2=x_it0*scov_in*x_j0;
                mskew=mskew+d_ij2*d_ij2*d_ij2;
        end;
end;
mskew=mskew/(N*N);

mskew0=N*mskew/6;
df=p*(p+1)*(p+2)/6;

finish;

***--------------------------------------------------------------------***;
*** Robustify data with Huber's type of weight;
***--------------------------------------------------------------------***;
start robmusig(x,H_c,H_beta,mu_1,sigma_1, xr);
ep=.00001;
n=nrow(x);
p=ncol(x);
xr=j(n,p,0);

**initial values;
mu_0=j(1,p,0); sigma_0=j(p,p,0);
do i=1 to n;
	x_i=x[i,];
	mu_0=mu_0+x_i;
	sigma_0=sigma_0+(x_i`)*x_i;
end;
mu_0=mu_0/n;
sigma_0=(sigma_0-n*(mu_0`)*mu_0)/(n-1);

dtheta=1;
n_div=1;
do jj=1 to 300 while (max(abs(dtheta))>ep);
	sig_in=inv(sigma_0);
	sumw1=0.0;
	sumx=j(1,p,0);
	sumxx=j(p,p,0);
	do i=1 to n;
		x_i=x[i,];
  		x_i0=x_i-mu_0;
  		sig_inx0=sig_in*(x_i0`);
  		d_i2=x_i0*sig_inx0;
  		d_i=sqrt(d_i2);
**Huber weight functions;
  		if d_i<=H_c then do;
     		w_i1=1.0;
     		w_i2=1.0/H_beta;
  		end;
  		else do;
     		w_i1=H_c/d_i;
     		w_i2=w_i1*w_i1/H_beta;
  		end;
  		sumw1=sumw1+w_i1;
  		sumx=sumx+w_i1*x_i;
  		sumxx=sumxx+w_i2*(x_i0`)*x_i0;
		**(ii);
	end;
	mu_1=sumx/sumw1;
	sigma_1=sumxx/n;
	dtheta=max(abs(mu_1-mu_0),abs(sigma_1-sigma_0));
	mu_0=mu_1;
	sigma_0=sigma_1;
end;

*print "robust mean=" mu_1;
*print "robust covariance=" sigma_1;
if jj<299 then do;
	do i=1 to n;
  		x_i=x[i,];
  		x_i0=x_i-mu_1;
		sig_inx0=sig_in*(x_i0`);
  		d_i2=x_i0*sig_inx0;
  		d_i=sqrt(d_i2);
**Huber weight functions;
  		if d_i<=H_c then do;
     		w_i1=1.0;
     		w_i2=1.0/H_beta;
  		end;
  		else do;
     		w_i1=H_c/d_i;
     		w_i2=w_i1*w_i1/H_beta;
  		end;
  		xr[i,]=sqrt(w_i2)*x_i0;
	end;
end;
finish;

***--------------------------------------------------------------------***;
**** Subroutine to evaluate sigma0 and (dsigma0/dtheta);
*Model as represented by Figure 7 of Deng & Yuan;
***--------------------------------------------------------------------***;
start sigdsig_sem(p,theta0,vsig,Sigma,vdsig);
ps=p*(p+1)/2;
p_x=2;
p_y=8;
lamb_x=theta0[1:2];
lamb_y=j(8,4,0);
lamb_y[1,1]=1.0;
lamb_y[2,1]=theta0[3];
lamb_y[3,2]=1.0;
lamb_y[4,2]=theta0[4];
lamb_y[5,3]=1.0;
lamb_y[6,3]=theta0[5];
lamb_y[7,4]=1.0;
lamb_y[8,4]=theta0[6];
lamb_xt=lamb_x`;
lamb_yt=lamb_y`;

gamma=j(4,1,0);
gamma[1]=theta0[7];
gamma_t=gamma`;

*AB=I-Beta in LISREL notation;
Beta=j(4,4,0);
Beta[2,1]=theta0[8];
Beta[3,1]=theta0[9];
Beta[3,2]=theta0[10];
Beta[4,2]=theta0[11];
Beta[4,3]=theta0[12];
A=i(4);
AB=A-Beta;

phi=1.0;

sig_zeta=diag(theta0[13:16]);

psi_x=diag(theta0[17:18]);
psi_y=diag(theta0[19:26]);

*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,2,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,4,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(2,2,0);

tt=j(8,4,0);
tt[2,1]=1;
tyy=tt*tly1+lyt1*(tt`);
txy=t2*(tt`);
ttt=(txx||txy)//((txy`)||tyy);
run vech(ttt,tttt);
dlamb_y[,1]=tttt;

tt=j(8,4,0);
tt[4,2]=1;
tyy=tt*tly1+lyt1*(tt`);
txy=t2*tt`;
ttt=(txx||txy)//((txy`)||tyy);
run vech(ttt,tttt);
dlamb_y[,2]=tttt;

tt=j(8,4,0);
tt[6,3]=1;
tyy=tt*tly1+lyt1*(tt`);
txy=t2*(tt`);
ttt=(txx||txy)//((txy`)||tyy);
run vech(ttt,tttt);
dlamb_y[,3]=tttt;

tt=j(8,4,0);
tt[8,4]=1;
tyy=tt*tly1+lyt1*(tt`);
txy=t2*(tt`);
ttt=(txx||txy)//((txy`)||tyy);
run vech(ttt,tttt);
dlamb_y[,4]=tttt;

*dgamma is the derivative with gamma;
*dgamma=j(ps,1,0);
t1=lamb_y*ABin;
t2=lamb_x*phi;
txx=j(2,2,0);

tt=j(4,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,5,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(2,2,0);
j=1;	
	tt=j(4,4,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;
do j=2 to 3;	
	tt=j(4,4,0);
	tt[3,j-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;
do j=4 to 5;	
	tt=j(4,4,0);
	tt[4,j-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;
end;


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

*dsig_zeta is the derivative with sig_zeta;
dsig_zeta=j(ps,4,0);

t1=lamb_y*ABin;
txx=j(p_x,p_x,0);
txy=j(p_x,p_y,0);
do j=1 to 4;
	tt=j(4,4,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;
*Model as represented by Figure 7 of Deng & Yuan;
***--------------------------------------------------------------------***;
start ML_SEM(n,p,dup,theta0,scov,T_ml,htheta_z,coef_xy,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;
	lamb=j(10,5,0);
	lamb_1=theta0[1:2];
	lamb_2=1//theta0[3];
	lamb_3=1//theta0[4];
	lamb_4=1//theta0[5];
	lamb_5=1//theta0[6];
	psi_1=diag(theta0[17:18]);
	psi_2=diag(theta0[19:20]);
	psi_3=diag(theta0[21:22]);
	psi_4=diag(theta0[23:24]);
	psi_5=diag(theta0[25:26]);
	psiL_1=inv(psi_1)*lamb_1;
	psiL_2=inv(psi_2)*lamb_2;
	psiL_3=inv(psi_3)*lamb_3;
	psiL_4=inv(psi_4)*lamb_4;
	psiL_5=inv(psi_5)*lamb_5;

	coef_1=inv(lamb_1`*psiL_1)*psiL_1;
	coef_2=inv(lamb_2`*psiL_2)*psiL_2;
	coef_3=inv(lamb_3`*psiL_3)*psiL_3;
	coef_4=inv(lamb_4`*psiL_4)*psiL_4;
	coef_5=inv(lamb_5`*psiL_5)*psiL_5;
	coef_xy=coef_1||coef_2||coef_3||coef_4||coef_5;

end;
finish;

*--------------------------------------------;
*LS regression using composite-scores;
*--------------------------------------------;
start regress(xy,reg_est, Rsq);
n=nrow(xy);
x=xy[,1];
y1=xy[,2];
y2=xy[,3];
y3=xy[,4];
y4=xy[,5];

x_t=x`;
s_xx=x_t*x; 
s_xy1=x_t*y1;
hgamma1_1=s_xy1/s_xx;
hsig2_11=ssq(y1-x*hgamma1_1)/(n-2);
Vargamma1_1=hsig2_11/s_xx; SEgamma1_1=sqrt(Vargamma1_1);
zgamma1_1=hgamma1_1/SEgamma1_1;
Rsq_1=1-ssq(y1-x*hgamma1_1)/ssq(y1);
*----------------;
y1_t=y1`;
y2_t=y2`;
s_y11=y1_t*y1; 
s_y21=y2_t*y1;
hbeta2_1=s_y21/s_y11;
hsig2_21=ssq(y2-y1*hbeta2_1)/(n-2);
Varbeta2_1=hsig2_21/s_y11; SEbeta2_1=sqrt(Varbeta2_1);
zbeta2_1=hbeta2_1/SEbeta2_1;
Rsq_2=1-ssq(y2-y1*hbeta2_1)/ssq(y2);
*----------------;
y3_t=y3`;
s_y22=y2_t*y2;
s_y31=y3_t*y1;
s_y32=y3_t*y2;
CCmat_3=(s_y11||s_y21)//(s_y21||s_y22);
cvec_3=s_y31//s_y32;
hbeta_3=inv(CCmat_3)*cvec_3;
yy_12=y1||y2;
hsig2_3=ssq(y3-yy_12*hbeta_3)/(n-3);
Covbeta_3=hsig2_3*inv(CCmat_3); 
SEbeta_3=sqrt(vecdiag(Covbeta_3));
zbeta_3=hbeta_3/SEbeta_3;
Rsq_3=1-ssq(y3-yy_12*hbeta_3)/ssq(y3);
*----------------;
y4_t=y4`;
s_y33=y3_t*y3;
s_y42=y4_t*y2;
s_y43=y4_t*y3;
CCmat_4=(s_y22||s_y32)//(s_y32||s_y33);
cvec_4=s_y42//s_y43;
hbeta_4=inv(CCmat_4)*cvec_4;
yy_23=y2||y3;
hsig2_4=ssq(y4-yy_23*hbeta_4)/(n-3);
Covbeta_4=hsig2_4*inv(CCmat_4); SEbeta_4=sqrt(vecdiag(Covbeta_4));
zbeta_4=hbeta_4/SEbeta_4;
Rsq_4=1-ssq(y4-yy_23*hbeta_4)/ssq(y4);

reg_est=j(6,3,0);
reg_est[1,]=hgamma1_1||SEgamma1_1||zgamma1_1;
reg_est[2,]=hbeta2_1||SEbeta2_1||zbeta2_1;
reg_est[3:4,]=(hbeta_3||SEbeta_3||zbeta_3);
reg_est[5:6,]=(hbeta_4||SEbeta_4||zbeta_4);
Rsq=Rsq_1||Rsq_2||Rsq_3||Rsq_4;
*print "Rsq=" Rsq[format=6.3];
*print "Alcol data";
*print reg_est;
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=" j psi_j;
		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 raw data (centralized);
*the model is as represented by Figure 7 of Deng & Yuan;
******************************************************************;
start PLS_A(x,y1,y2,y3,y4,reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
n=nrow(x); 
p_x=ncol(x);
p_y1=ncol(y1);
p_y2=ncol(y2);
p_y3=ncol(y3);
p_y4=ncol(y4);
x_t=x`;
y1_t=y1`;
y2_t=y2`;
y3_t=y3`;
y4_t=y4`;
s_xx=x_t*x/n;
s_y11=y1_t*y1/n;
s_y22=y2_t*y2/n;
s_y33=y3_t*y3/n;
s_y44=y4_t*y4/n;
*print "check s_y=" s_xx s_y11 s_y22 s_y33 s_y44;
*--------------------;
*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);

w_y4=j(p_y4,1,1);
c2_y4=w_y4`*s_y44*w_y4;
w_y4=w_y4/sqrt(c2_y4);

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


*environment variables;
r_xy1=(eta_x`*eta_y1)/n;
r_y12=(eta_y1`*eta_y2)/n;
r_y13=(eta_y1`*eta_y3)/n;
r_y23=(eta_y2`*eta_y3)/n;
r_y24=(eta_y2`*eta_y4)/n;
r_y34=(eta_y3`*eta_y4)/n;

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+sign(r_y24)*eta_y4;
ev_y3=sign(r_y13)*eta_y1+sign(r_y23)*eta_y2+sign(r_y34)*eta_y4;
ev_y4=sign(r_y24)*eta_y2+sign(r_y34)*eta_y3;

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);
ev_y4=ev_y4/sqrt(ssq(ev_y4)/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//w_y4;
*print "-------------mode A-------------";
dw=1; ep=.00001;
do j=1 to 50 while (max(abs(dw))>ep);
	w_x=x_t*eta_y1/n;
	*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;
	*w_2 is the vector of regression coefficients eta_y-->x2, var(eta_y)=1;
	c2_y1=w_y1`*s_y11*w_y1;
	w_y1=sign(w_y1[1])*w_y1/sqrt(c2_y1);
	eta_y1=y1*w_y1; *eta_y1 is the vector of predicted score;

	w_y2=y2_t*ev_y2/n;
	c2_y2=w_y2`*s_y22*w_y2;
	w_y2=sign(w_y2[1])*w_y2/sqrt(c2_y2);
	eta_y2=y2*w_y2; *eta_y2 is the vector of predicted score;
	
	w_y3=y3_t*ev_y3/n;
	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;

	w_y4=y4_t*ev_y4/n;
	c2_y4=w_y4`*s_y44*w_y4;
	w_y4=sign(w_y4[1])*w_y4/sqrt(c2_y4);
	eta_y4=y4*w_y4; *eta_y4 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=(eta_x`*eta_y1)/n;
	r_y12=(eta_y1`*eta_y2)/n;
	r_y13=(eta_y1`*eta_y3)/n;
	r_y23=(eta_y2`*eta_y3)/n;
	r_y24=(eta_y2`*eta_y4)/n;
	r_y34=(eta_y3`*eta_y4)/n;

	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+sign(r_y24)*eta_y4;
	ev_y3=sign(r_y13)*eta_y1+sign(r_y23)*eta_y2+sign(r_y34)*eta_y4;
	ev_y4=sign(r_y24)*eta_y2+sign(r_y34)*eta_y3;
	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);
	ev_y4=ev_y4/sqrt(ssq(ev_y4)/n);

*	check_y3=ev_y3`*ev_y3/n;
*print "check=1 (y3)=" check_y4;

*	check_x=ev_x`*ev_x/n;
*	print "check_j=1=" check_x;
	*print jj w_x1 w_x2;
	w_j=w_x//w_y1//w_y2//w_y3//w_y4;
	dw=w_j-w_0;
	w_0=w_j;
end;
print "mode A iteration time=" j;
w_a=w_j;
cc_xx=eta_x`*eta_x; 
cc_y11=eta_y1`*eta_y1;
cc_y21=eta_y2`*eta_y1;
cc_y22=eta_y2`*eta_y2;
cc_y23=eta_y2`*eta_y3;
cc_y33=eta_y3`*eta_y3;

c_xy1=eta_x`*eta_y1;
c_y21=eta_y2`*eta_y1;
c_y31=eta_y3`*eta_y1;
c_y32=eta_y3`*eta_y2;
c_y42=eta_y4`*eta_y2;
c_y43=eta_y4`*eta_y3;

cc_mat3=(cc_y11||cc_y21)//(cc_y21||cc_y22);
c_vec3=c_y31//c_y32;

cc_mat4=(cc_y22||cc_y23)//(cc_y23||cc_y33);
c_vec4=c_y42//c_y43;

hat_g11=c_xy1/cc_xx;
hat_b21=c_y21/cc_y11;
hat_b3=inv(cc_mat3)*c_vec3;
hat_b4=inv(cc_mat4)*c_vec4;

res_1=eta_y1-eta_x*hat_g11;
hsig2_1=ssq(res_1)/(n-2);
Var_g11=hsig2_1/cc_xx; 
SE_g11=sqrt(Var_g11);
z_g11=hat_g11/SE_g11;
Rsq_1=1-ssq(res_1)/ssq(eta_y1);

res_2=eta_y2-eta_y1*hat_b21;
hsig2_2=ssq(res_2)/(n-2);
Var_b21=hsig2_2/cc_y11;
SE_b21=sqrt(Var_b21);
z_b21=hat_b21/SE_b21;
Rsq_2=1-ssq(res_2)/ssq(eta_y2);

res_3=eta_y3-(eta_y1||eta_y2)*hat_b3;
hsig2_3=ssq(res_3)/(n-3);
Cov_b3=hsig2_3*inv(cc_mat3);
SE_b3=sqrt(vecdiag(Cov_b3));
z_b3=hat_b3/SE_b3;
Rsq_3=1-ssq(res_3)/ssq(eta_y3);

res_4=eta_y4-(eta_y2||eta_y3)*hat_b4;
hsig2_4=ssq(res_4)/(n-3);
Cov_b4=hsig2_4*inv(cc_mat4);
SE_b4=sqrt(vecdiag(Cov_b4));
z_b4=hat_b4/SE_b4;
Rsq_4=1-ssq(res_4)/ssq(eta_y4);

reg_plsA=j(6,3,0);
reg_plsA[1,]=hat_g11||SE_g11||z_g11;
reg_plsA[2,]=hat_b21||SE_b21||z_b21;
reg_plsA[3:4,]=hat_b3||SE_b3||z_b3;
reg_plsA[5:6,]=hat_b4||SE_b4||z_b4;
Rsq_A=Rsq_1||Rsq_2||Rsq_3||Rsq_4;


*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);
run ssw(w_y4,S_y44,w_sy4);

w_sb=w_sx//w_sy1//w_sy2//w_sy3//w_sy4;
*print "---------------------------------weights---------------------------------";
*print "mode A(w_x,w_y1,w_y2,w_y3,w_y4)--mode B_A(w_x,w_y1,w_y2,w_y3,w_y4)";
*print w_x w_y1 w_y2 w_y3 w_y4 w_sx w_sy1 w_sy2 w_sy3 w_sy4;

eta_sx=x*w_sx;
eta_sy1=y1*w_sy1;
eta_sy2=y2*w_sy2;
eta_sy3=y3*w_sy3;
eta_sy4=y4*w_sy4;

cc_sxx=eta_sx`*eta_sx; 
cc_sy11=eta_sy1`*eta_sy1;
cc_sy21=eta_sy2`*eta_sy1;
cc_sy22=eta_sy2`*eta_sy2;
cc_sy23=eta_sy2`*eta_sy3;
cc_sy33=eta_sy3`*eta_sy3;

c_sxy1=eta_sx`*eta_sy1;
c_sy21=eta_sy2`*eta_sy1;
c_sy31=eta_sy3`*eta_sy1;
c_sy32=eta_sy3`*eta_sy2;
c_sy42=eta_sy4`*eta_sy2;
c_sy43=eta_sy4`*eta_sy3;

cc_smat3=(cc_sy11||cc_sy21)//(cc_sy21||cc_sy22);
c_svec3=c_sy31//c_sy32;

cc_smat4=(cc_sy22||cc_sy23)//(cc_sy23||cc_sy33);
c_svec4=c_sy42//c_sy43;

hat_sg11=c_sxy1/cc_sxx;
hat_sb21=c_sy21/cc_sy11;
hat_sb3=inv(cc_smat3)*c_svec3;
hat_sb4=inv(cc_smat4)*c_svec4;

res_s1=eta_sy1-eta_sx*hat_sg11;
hsig2_s1=ssq(res_s1)/(n-2);
Var_sg11=hsig2_s1/cc_sxx; 
SE_sg11=sqrt(Var_sg11);
z_sg11=hat_sg11/SE_sg11;
Rsq_s1=1-ssq(res_s1)/ssq(eta_sy1);

res_s2=eta_sy2-eta_sy1*hat_sb21;
hsig2_s2=ssq(res_s2)/(n-2);
Var_sb21=hsig2_s2/cc_sy11;
SE_sb21=sqrt(Var_sb21);
z_sb21=hat_sb21/SE_sb21;
Rsq_s2=1-ssq(res_s2)/ssq(eta_sy2);

res_s3=eta_sy3-(eta_sy1||eta_sy2)*hat_sb3;
hsig2_s3=ssq(res_s3)/(n-3);
Cov_sb3=hsig2_s3*inv(cc_smat3);
SE_sb3=sqrt(vecdiag(Cov_sb3));
z_sb3=hat_sb3/SE_sb3;
Rsq_s3=1-ssq(res_s3)/ssq(eta_sy3);

res_s4=eta_sy4-(eta_sy2||eta_sy3)*hat_sb4;
hsig2_s4=ssq(res_s4)/(n-3);
Cov_sb4=hsig2_s4*inv(cc_smat4);
SE_sb4=sqrt(vecdiag(Cov_sb4));
z_sb4=hat_sb4/SE_sb4;
Rsq_s4=1-ssq(res_s4)/ssq(eta_sy4);

reg_plsBA=j(6,3,0);
reg_plsBA[1,]=hat_sg11||SE_sg11||z_sg11;
reg_plsBA[2,]=hat_sb21||SE_sb21||z_sb21;
reg_plsBA[3:4,]=hat_sb3||SE_sb3||z_sb3;
reg_plsBA[5:6,]=hat_sb4||SE_sb4||z_sb4;

Rsq_BA=Rsq_s1||Rsq_s2||Rsq_s3||Rsq_s4;
finish;
***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;
print "--------------------------------------------------------------------";
print "Neumann's data, 10-variable, 5-factor, factors f2-f5 further";
print "regress on factor 1, rescaled to have std less than 2";

print "---------------------------------------------------------";
use Neumann;
read all var _num_ into xmatr;

n=nrow(xmatr);
p=ncol(xmatr);
ps=p*(p+1)/2;

sc={  10.0
                                   55.0
                                   16.0
                                   14.0
                                   1.0
                                   2.0
                                   11.0
                                   10.0
                                   8.0
                                   5.0
};


rsc=diag(j(1,10,1)/sc);
xmatr=xmatr*rsc;

x1=xmatr[,5];
x2=xmatr[,6];
x3=xmatr[,8];
x4=xmatr[,7];
x5=xmatr[,4];
x6=xmatr[,3];
x7=xmatr[,2];
x8=xmatr[,1];
x9=xmatr[,10];
x10=xmatr[,9];

xf1=x1||x2;
Xf2=x3||x4;
xf3=x5||x6;
xf4=x7||x8;
xf5=x9||x10;
x=xf1||xf2||xf3||xf4||xf5;

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

run Mardia(x,mskew_0,mkurt_0,mkurt_p,df);
print "Mardia's measure of skew and kurt in raw data";
	print "df for mskew_0=" df;
	print "mskew_0=" mskew_0 [format=f10.3];
	print "mkurt_0=" mkurt_0 [format=f8.3];
	print "relative mkurt=" mkurt_p [format=f8.3];

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

prob=.90; *[rawdata, Mk=14.763, (.05, Mk=2.248), (.10, -0.110)];

qprob=1-prob;
*print "the q in Huber(q) =" qprob;
eta=p/2.0;
chip2=2*gaminv(prob,eta);
ck=sqrt(chip2);
cbeta=( p*probgam(chip2/2.0,eta+1)+ chip2*(1-prob) )/p;

run robmusig(x,ck,cbeta,mu1,sigma1, xr);
*run Mardia(xr,mskew_r,mkurt_r,mkurt_rp,df);
print "---------------------------------------";
*print xr;
*print "---------------------------------------";
*print "Mardia's measure of skew and kurt in robustified data";
*	print "df for mskew_0=" df;
*	print "mskew_r=" mskew_r [format=f10.3];
*	print "mkurt_r=" mkurt_r [format=f8.3];
*	print "relative mkurt=" mkurt_rp [format=f8.3];

*print "-----------------------------------------------------------";
*(SEM raw data)T_ml=  48.96027,
*(SEM robust data .90)T_ml= 40.984805;


*x=xr;
*xr=x;
Scov=xr`*(i(n)-j(n,n,1)/n)*xr/(n-1);
*print "scov=" scov;
run vech(Scov,vscov);

a0_sem={
 1.0260574 0.9310802 1.0994877 0.8090211  0.588969 0.476005 -0.0893442
  -0.6249662 0.225653 0.425656 -0.496624 -0.4852545 0.6798855  0.485086
  0.5025141 0.5137928  1.393926 0.3040502 0.2354705  0.435724 0.4249227
  0.2197304 0.2687321 0.5495969 0.4894514 0.7926654
};


a0=a0_sem`; *starting values for ML estimation;

q=nrow(a0);
theta0=a0;
df=ps-q;
run DP(p, dup);

run ML_sem(n,p,dup,theta0,scov,T_ml,htheta_z,coef_sem,nsem_div);
if nsem_div=1 then do;
	print "Fisher-scoring algorithm for SEM did not converge within 300 iterations";
end;
else do;
	RMSEA=sqrt((T_ml-df)/((n-1)*df));
	Sig_I=diag(Scov);
	SS_I=Scov*inv(Sig_I);
	Fml_I=trace(SS_I)-log(det(SS_I))-p;
	Tml_I=(n-1)*Fml_I;
	df_I=ps-p;
	CFI=1-max((T_ml-df),0)/max((Tml_I-df_I),(T_ml-df),0);
	pvalue=1-probchi(T_ml,df);
	print "(CB-SEM)T_ml, p_value, RMSEA, CFI =" T_ml pvalue RMSEA CFI;

	print "(CB-SEM)htheta--SE--Z";
	print htheta_z [format=f7.3];
	hgamma_z=htheta_z[7:12,];
*------------------R-square under CB-SEM------------------;
	gamma_11=theta0[7];
	beta_21=theta0[8];
	beta_31=theta0[9];
	beta_32=theta0[10];
	beta_42=theta0[11];
	beta_43=theta0[12];
	*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_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[14]+theta0[15];
	Rsq_cb3=1-theta0[15]/Var_eta3;
	*eta4=beta_42*eta_2+beta_43*eta_3+zeta_4
    =beta_42*(beta_21*eta_1+zeta_2)
    +beta_43*[(beta_31+beta_32*beta_21)*eta_1+beta_32*zeta_2+zeta_3]+zeta_4
    =beta_42*beta_21*eta_1+beta_42*zeta_2
    +beta_43*(beta_31+beta_32*beta_21)*eta_1+beta_43*beta_32*zeta_2+beta_43*zeta_3+zeta_4
    =[beta_42*beta_21+beta_43*(beta_31+beta_32*beta_21)]*eta_1
    +(beta_42+beta_43*beta_32)*zeta_2+beta_43*zeta_3+zeta_4;
	beta_4=beta_42*beta_21+beta_43*(beta_31+beta_32*beta_21);
	c_zeta2=beta_42+beta_43*beta_32;
	Var_eta4=beta_4*beta_4*Var_eta1
	+c_zeta2*c_zeta2*theta0[14]
	+beta_43*beta_43*theta0[15]+theta0[16];
	Rsq_cb4=1-theta0[16]/Var_eta4;
	Rsq_cb=Rsq_cb1||Rsq_cb2||Rsq_cb3||Rsq_cb4;
*----------------------FS reg----------------------;
	fx_1=xr[,1:2]*coef_sem[,1];
	fx_2=xr[,3:4]*coef_sem[,2];
	fx_3=xr[,5:6]*coef_sem[,3];
	fx_4=xr[,7:8]*coef_sem[,4];
	fx_5=xr[,9:10]*coef_sem[,5];
	xy_fs=fx_1||fx_2||fx_3||fx_4||fx_5;
	run regress(xy_fs,reg_fs,Rsq_fs);

	*----------EWC reg-------------;
	lamb_x=theta0[1:2];
	lamb_y1=1.0//theta0[3];
	lamb_y2=1.0//theta0[4];
	lamb_y3=1.0//theta0[5];
	lamb_y4=1.0//theta0[6];
	c_x=j(2,1,1)/(j(1,2,1)*lamb_x);
	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_y3=j(2,1,1)/(j(1,2,1)*lamb_y3);
	c_y4=j(2,1,1)/(j(1,2,1)*lamb_y4);
	x_c=xr[,1:2]*c_x;
	y_c1=xr[,3:4]*c_y1;
	y_c2=xr[,5:6]*c_y2;
	y_c3=xr[,7:8]*c_y3;
	y_c4=xr[,9:10]*c_y4;
	xy_c=x_c||y_c1||y_c2||y_c3||y_c4;
	run regress(xy_c,reg_ewc,Rsq_ewc);
	print "(gamma_11//b_21//b_31//b_32//b_42//b_43)(hat\gamma, SE, z)";
	print "---CB-SEM--- ---FSR--- ---EWCR---"; 
	print hgamma_z[format=f7.3] reg_fs[format=f7.3] reg_ewc[format=f7.3];

end;

x=xr[,1:2]; y1=xr[,3:4]; y2=xr[,5:6]; y3=xr[,7:8]; y4=xr[,9:10];
run PLS_A(x,y1,y2,y3,y4,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_fs//Rsq_ewc//Rsq_A//Rsq_BA)`;
*print "R-square(CB-SEM,FSR,EWCR,PLSA,PLSAB)";
*print Rsq_all[format=f6.3];


finish;
run main;