***------------------------------------------------------------***;
** Data FROM Table 7_4 of the LISREL8 manual on page 193, 
and the model is as in Figure 7_2 on page 194
of J\"oreskog, K.G., \& S\"orbom, D. (1993)
{\em LISREL 8: Structural equation modeling with the simplis command language}.
Lincolnwood, IL: Scientific Software International.
*** All tests are based on normal theory estimates;

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

options ls=75 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;

***--------------------------------------------------------------------***;
**** subroutine to evaluate sigma0 and (dsigma0/dtheta);
***--------------------------------------------------------------------***;

start sigdsig_sem(p,theta0,vsig,Sigma,vdsig);
p_x=4; p_y=8;
ps=p*(p+1)/2;

lamb_x=j(p_x,2,0);
lamb_x[1:2,1]=theta0[1:2];
lamb_x[3:4,2]=theta0[3:4];
lamb_xt=lamb_x`;

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

gamma=j(4,2,0);
gamma[1,1]=theta0[9];
gamma[2,2]=theta0[10];
gamma_t=gamma`;

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

phi=j(2,2,1);
phi[2,1]=theta0[13];
phi[1,2]=theta0[13];
sig_zeta=diag(theta0[14:17]);
psi_x=diag(theta0[18:21]);
psi_y=diag(theta0[22:29]);



*computing Sigma based on structured parameters;
*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 2;
	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=3 to 4;
	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,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(p_x,p_x,0);
j=1;
	tt=j(p_y,4,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,4,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;
j=3;
	tt=j(p_y,4,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;
j=4;
	tt=j(p_y,4,0);
	tt[j+4,4]=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,2,0);
t1=lamb_y*ABin;
t2=lamb_x*phi;
txx=j(p_x,p_x,0);
do j=1 to 2;
	tt=j(4,2,0);
	tt[j,j]=1;
	tyy=t1*(tt*phi*gamma_t+gamma*phi*(tt`))*(t1`);
	txy=t2*(tt`)*(t1`);
	ttt=(txx||txy)//((txy`)||tyy);
	run vech(ttt,tttt);
	dgamma[,j]=tttt;
end;	


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

*dphi is the derivative with phi;
tt=j(2,2,0);
tt[2,1]=1;
tt[1,2]=1;
txx=lamb_x*tt*lamb_xt;
t1=lamb_y*ABin*gamma; tyy=t1*tt*(t1`);
txy=lamb_x*tt*(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||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);

**if idiv=1 then model does not converge;
p_x=4; p_y=8;
ep=.00001;
**ep=E-5;
ep1=.00000000000000000001;
**ep1=E-20;
ps=p*(p+1)/2;
q=nrow(theta0);
df=ps-q;
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;
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;

*computation of factor-score;
	lamb_x=j(p_x,2,0);
	lamb_x[1:2,1]=theta0[1:2];
	lamb_x[3:4,2]=theta0[3:4];
	lamb_xt=lamb_x`;

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

	psi_x=diag(theta0[18:21]);
	psi_y=diag(theta0[22:29]);
	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=4;p_y=8;
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;

*----------------;
*y1=g1 x1+e_1;
hat_g1=s_xy[1,1]/s_xx[1,1];
hsig2_g1=n*(s_yy[1,1]-s_xy[1,1]*s_xy[1,1]/s_xx[1,1])/(n-2);
Var_g1=hsig2_g1/(n*s_xx[1,1]); SE_g1=sqrt(Var_g1);
z_g1=hat_g1/SE_g1;
Rsq_1=1-(s_yy[1,1]-s_xy[1,1]*s_xy[1,1]/s_xx[1,1])/s_yy[1,1];

*----------------;
*y2=g2 x2+e_2;
hat_g2=s_xy[2,2]/s_xx[2,2];
hsig2_g2=n*(s_yy[2,2]-s_xy[2,2]*s_xy[2,2]/s_xx[2,2])/(n-2);
Var_g2=hsig2_g2/(n*s_xx[2,2]); SE_g2=sqrt(Var_g2);
z_g2=hat_g2/SE_g2;
Rsq_2=1-(s_yy[2,2]-s_xy[2,2]*s_xy[2,2]/s_xx[2,2])/s_yy[2,2];

*----------------;
*y3=b1 y1+e_3;
hat_b1=s_yy[3,1]/s_yy[1,1];
hsig2_b1=n*(s_yy[3,3]-s_yy[3,1]*s_yy[1,3]/s_yy[1,1])/(n-2);
Var_b1=hsig2_b1/(n*s_yy[1,1]); SE_b1=sqrt(Var_b1);
z_b1=hat_b1/SE_b1;
Rsq_3=1-(s_yy[3,3]-s_yy[3,1]*s_yy[1,3]/s_yy[1,1])/s_yy[3,3];

*----------------;
*y4=b2 y2+e_4;
hat_b2=s_yy[4,2]/s_yy[2,2];
hsig2_b2=n*(s_yy[4,4]-s_yy[4,2]*s_yy[2,4]/s_yy[2,2])/(n-2);
Var_b2=hsig2_b2/(n*s_yy[2,2]); SE_b2=sqrt(Var_b2);
z_b2=hat_b2/SE_b2;
Rsq_4=1-(s_yy[4,4]-s_yy[4,2]*s_yy[2,4]/s_yy[2,2])/s_yy[4,4];

*----------------;
reg_est=j(4,3,0);
reg_est[1,]=hat_g1||SE_g1||z_g1;
reg_est[2,]=hat_g2||SE_g2||z_g2;
reg_est[3,]=hat_b1||SE_b1||z_b1;
reg_est[4,]=hat_b2||SE_b2||z_b2;

Rsq=Rsq_1||Rsq_2||Rsq_3||Rsq_4;
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;
*--------------------------------------------;
start PLS_A(n,scov0, reg_plsA,reg_plsBA, Rsq_A, Rsq_BA);
scov=(n-1)*scov0/n;
p_x1=2;
p_x2=2;
p_y1=2;
p_y2=2;
p_y3=2;
p_y4=2;

s_x11=scov[1:2,1:2];s_x12=scov[1:2,3:4];s_x1y1=scov[1:2,5:6];s_x1y2=scov[1:2,7:8];
s_x1y3=scov[1:2,9:10];s_x1y4=scov[1:2,11:12];

s_x21=s_x12`;s_x22=scov[3:4,3:4];s_x2y1=scov[3:4,5:6];s_x2y2=scov[3:4,7:8];
s_x2y3=scov[3:4,9:10];s_x2y4=scov[3:4,11:12];

s_y1x1=s_x1y1`;s_y1x2=s_x2y1`;s_y11=scov[5:6,5:6];s_y12=scov[5:6,7:8];
s_y13=scov[5:6,9:10];s_y14=scov[5:6,11:12];

s_y2x1=s_x1y2`;s_y2x2=s_x2y2`;s_y21=s_y12`;s_y22=scov[7:8,7:8];
s_y23=scov[7:8,9:10];s_y24=scov[7:8,11:12];

s_y3x1=s_x1y3`;s_y3x2=s_x2y3`;s_y31=s_y13`;s_y32=s_y23`;
s_y33=scov[9:10,9:10];s_y34=scov[9:10,11:12];

s_y4x1=s_x1y4`;s_y4x2=s_x2y4`;s_y41=s_y14`;s_y42=s_y24`;
s_y43=s_y34`;s_y44=scov[11:12,11:12];

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

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;


*environment variables;
r_x1y1=w_x1`*s_x1y1*w_y1;
r_x2y2=w_x2`*s_x2y2*w_y2;
r_y13=w_y1`*s_y13*w_y3;
r_y24=w_y2`*s_y24*w_y4;

*ev_y1=sign(r_x1y1)*eta_x1+sign(r_y13)*eta_y3;
*ev_y2=sign(r_x2y2)*eta_x2+sign(r_y24)*eta_y4;

w_0=w_x1//w_x2//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*ev_y1/n, ev_x=sign(r_xy1)*eta_y1+sign(r_xy2)*eta_y2;
	w_x1=sign(r_x1y1)*s_x1y1*w_y1;
	*w_x1 is the vector of regression coefficients eta_y1-->x1, var(ev_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;

	w_x2=sign(r_x2y2)*s_x2y2*w_y2;
	*w_x2 is the vector of regression coefficients eta_y2-->x2, var(ev_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_y13)*eta_y3;
	w_y1=sign(r_x1y1)*(s_y1x1)*w_x1+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_x2y2)*eta_x2+sign(r_y24)*eta_y4;
	w_y2=sign(r_x2y2)*s_y2x2*w_x2+sign(r_y24)*s_y24*w_y4;
	c2_y2=w_y2`*s_y22*w_y2;
	w_y2=sign(w_y2[1])*w_y2/sqrt(c2_y2);

	w_y3=sign(r_y13)*s_y31*w_y1;
	c2_y3=w_y3`*s_y33*w_y3;
	w_y3=sign(w_y3[1])*w_y3/sqrt(c2_y3);

	w_y4=sign(r_y24)*s_y42*w_y2;
	c2_y4=w_y4`*s_y44*w_y4;
	w_y4=sign(w_y4[1])*w_y4/sqrt(c2_y4);

	*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_x2y2=w_x2`*s_x2y2*w_y2;
	r_y13=w_y1`*s_y13*w_y3;
	r_y24=w_y2`*s_y24*w_y4;
	
	w_j=w_x1//w_x2//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;
w_x=(w_x1||j(2,1,0))//(j(2,1,0)||w_x2);
w_xt=w_x`;
w_y=j(8,4,0);
w_y[1:2,1]=w_y1;
w_y[3:4,2]=w_y2;
w_y[5:6,3]=w_y3;
w_y[7:8,4]=w_y4;
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_x1,S_x11,w_sx1); 
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);
*print "--------block-y1---------";
run ssw(w_y3,S_y33,w_sy3);
*print "--------block-y2---------";
run ssw(w_y4,S_y44,w_sy4);

w_sb=w_sx1//w_sx2//w_sy1//w_sy2//w_sy3//w_sy4;
*print "---------------------------------weights model 1---------------------------------";
*print "mode A(w_x1,w_x2,w_y1,w_y2,w_y3,w_y4)";
*print w_x1 w_x2 w_y1 w_y2 w_y3 w_y4;
*print "mode AB(w_x1,w_x2,w_y1,w_y2,w_y3,w_y4)";
*print w_sx1 w_sx2 w_sy1 w_sy2 w_sy3 w_sy4;
w_sx=(w_sx1||j(2,1,0))//(j(2,1,0)||w_sx2);
w_sxt=w_sx`;
w_sy=j(8,4,0);
w_sy[1:2,1]=w_sy1;
w_sy[3:4,2]=w_sy2;
w_sy[5:6,3]=w_sy3;
w_sy[7:8,4]=w_sy4;
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 "--------------------------------------------------------------------";
print "Data FROM Table 7_4 of the LISREL8 manual on page 193, 
and the model is as in Figure 7_2 on page 194 of Joreskog and Sorbom, D. (1993)";
print "---------------------------------------------------------";

n=189;
p=12;
ps=p*(p+1)/2;
print "N, p=" N p;
p_x=4;
p_y=8;
scov=j(p,p,0);
scov_vec={
.247
.227 .248
.479 .443 2.313
.312 .288 1.225 2.264
.198 .196  .432  .283 .250
.201 .205  .437  .299 .232 .250
.402 .375 1.459 1.053 .468 .467 2.453
.212 .185  .906 1.212 .267 .271 1.028 1.749
.185 .176  .375  .250 .203 .202  .397  .221 .250
.176 .172  .351  .247 .196 .199  .381  .210 .233 .250
.376 .330 1.308 1.020 .418 .419 1.517  .906 .502 .491 2.476
.231 .204  .907 1.186 .274 .273 1.004 1.107 .319 .312 1.326 2.281
};
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;
scov_a=scov;
a0={
.479  .473 1.261 1.009
1.014  .701  .976  .796
.421 1.123 .881 .937
.723
.052 .272 .061 .320
.018 .024 .724 1.245
.021  .015  .919  .995  .011  .023  .811 1.225};
a0=a0`; *starting values for ML estimation;

q=nrow(a0);
df0=ps-q;
*print "df0=" df0;
run Dp(p,dup);
theta0=a0;
scov=scov_a;
run ML_SEM(n, p, dup, theta0,scov,  T_ml,htheta_z,coef_xt,coef_yt,n_div);
if n_div=1 then do;
	print "Fisher-scoring algorithm did not converge within 300 iterations";
end;
else do;
	p_v=1-probchi(T_ml,df0);
	print "T_ml, df, p_value=" T_ml df0 p_v;

	print "(SEM model) htheta--SE--Z";
	print  htheta_z[format=f8.3];
	hgamma1_z=htheta_z[9:12,];
	*------------------R-square under CB-SEM------------------;
	gamma_11=theta0[9];
	gamma_22=theta0[10];
	beta_31=theta0[11];
	beta_42=theta0[12];
	phi_21=theta0[13];

	*eta1=gamma_11*xi_1+zeta_1;
	Var_eta1=gamma_11*gamma_11+theta0[14];
	Rsq_cb1=1-theta0[14]/Var_eta1;
	*eta2=gamma_22*xi_2+zeta_2;
	Var_eta2=gamma_22*gamma_22+theta0[15];
	Rsq_cb2=1-theta0[15]/Var_eta2;
	*eta3=beta_31*eta_1+zeta_3;
	Var_eta3=beta_31*beta_31*Var_eta1+theta0[16];
	Rsq_cb3=1-theta0[16]/Var_eta3;
	*eta4=beta_42*eta_2+zeta_4;
	Var_eta4=beta_42*beta_42*Var_eta2+theta0[17];
	Rsq_cb4=1-theta0[17]/Var_eta4;
	Rsq_cb=Rsq_cb1||Rsq_cb2||Rsq_cb3||Rsq_cb4;
*---------------------------------------------------------;
	scov=scov_a;
	run regress(n,coef_xt,coef_yt,scov,reg_fsr, Rsq_fsr);
	*----------EWC reg-------------;
	lamb_x1=theta0[1:2];
	lamb_x2=theta0[3:4];
	lamb_y1=1.0//theta0[5];
	lamb_y2=1.0//theta0[6];
	lamb_y3=1.0//theta0[7];
	lamb_y4=1.0//theta0[8];
	c_x1=j(2,1,1)/(j(1,2,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_y3=j(2,1,1)/(j(1,2,1)*lamb_y3);
	c_y4=j(2,1,1)/(j(1,2,1)*lamb_y4);
	c_x=(c_x1||j(2,1,0))//(j(2,1,0)||c_x2);
	c_xt=c_x`;
	c_y=j(8,4,0);
	c_y[1:2,1]=c_y1;
	c_y[3:4,2]=c_y2;
	c_y[5:6,3]=c_y3;
	c_y[7:8,4]=c_y4;
	c_yt=c_y`;
	scov=scov_a;
	run regress(n,c_xt,c_yt,scov,reg_ewc, Rsq_ewc);
	print "[gamma_11//gamma_22//b_31//b_42](hat\gamma, SE, z)";
	print "---CB-SEM--- ---FSR--- ---EWCR---"; 
	print hgamma1_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;