*************************************************
This program performs moderation analysis by NML,
as described in Yuan, Cheng and Maxwell (2012).
The first few lines need to be modified to read data
correctly. Then, for any model in the form of (A5) 
and (A6) of the paper only the last part of the 
program needs to be 
modified.
*************************************************;
options linesize=110 nodate nonumber ps=2000;
title;
data yxu;
infile 'd:\moderation\simudata.dat';
input x1-x6;
run;

proc IML;
reset noname;

*==========================================*;
*Sandwich-type covariance matrix for LS analysis;
start sandwich_LS(hatee,CC_in,C_mat,
 SE_SW0,SE_SW1,SE_SW2,SE_SW3,SE_SW4);

n=nrow(C_mat);
p_c=ncol(C_mat);
B0_mat=j(p_c,p_c,0);
B2_mat=j(p_c,p_c,0);
B3_mat=j(p_c,p_c,0);
B4_mat=j(p_c,p_c,0);
do i=1 to n;
	c_it=C_mat[i,];
	c_i=c_it`;
	cc_i=c_i*c_it;	
	t_ci=c_it*CC_in*c_i; *for B2_mat;
	vc_i=1/(1-t_ci);
	vc_i2=vc_i*vc_i; *for B3_mat;

	r_ci=min(4,(n*t_ci/p_c)); *for B4_mat;
	vc_ir=vc_i**r_ci;

	delta_i2=hatee[i];
***********;
	B0_i=delta_i2*cc_i;
	B2_i=vc_i*B0_i;
	B3_i=vc_i2*B0_i; 
	B4_i=vc_ir*B0_i; 
	
	B0_mat=B0_mat+B0_i;	
	B2_mat=B2_mat+B2_i;	
	B3_mat=B3_mat+B3_i;	
	B4_mat=B4_mat+B4_i;
end;

Cov_SW0=CC_in*B0_mat*CC_in;
Var_SW0=vecdiag(Cov_SW0);
SE_SW0=sqrt(Var_SW0);

SE_SW1=sqrt(n/(n-p_c))*SE_SW0;

Cov_SW2=CC_in*B2_mat*CC_in;
Var_SW2=vecdiag(Cov_SW2);
SE_SW2=sqrt(Var_SW2);

Cov_SW3=CC_in*B3_mat*CC_in;
Var_SW3=vecdiag(Cov_SW3);
SE_SW3=sqrt(Var_SW3);

Cov_SW4=CC_in*B4_mat*CC_in;
Var_SW4=vecdiag(Cov_SW4);
SE_SW4=sqrt(Var_SW4);

finish;

*==========================================*;
*Sandwich-type covariance matrix for ML analysis
of MMR model with heterogeneous variances;
start sandwich_ML(hate,tau_vec,CC_in,C_mat,H_mat, 
 SE_SW0,SE_SW1,SE_SW2,SE_SW3,SE_SW4);
n=nrow(C_mat);
p_c=ncol(C_mat);
p_h=ncol(H_mat);
p=p_c+p_h;
HH_mat=(H_mat`)*H_mat;
HH_in=inv(HH_mat);
A_mat=j(p,p,0);
B_mat0=j(p,p,0);
*B_mat1=j(p,p,0);
B_mat2=j(p,p,0);
B_mat3=j(p,p,0);
B_mat4=j(p,p,0);

do i=1 to n;
	c_it=C_mat[i,];
	c_i=c_it`;
	cc_i=c_i*c_it;	
	h_it=H_mat[i,];
	h_i=h_it`;
	ch_i=c_i*h_it;
	hh_i=h_i*h_it;

	t_ci=c_it*CC_in*c_i; *for B_mat2;
	t_hi=h_it*HH_in*h_i;
	vc_i=1/(1-t_ci);
	vh_i=1/(1-t_hi);
	vch_i=vc_i*vh_i;
	svch_i=sqrt(vch_i);

	vc_i2=vc_i*vc_i; *for B_mat3;
	vh_i2=vh_i*vh_i;

	r_ci=min(4,(n*t_ci/p_c)); *for B_mat4;
	r_hi=min(4,(n*t_hi/p_h));
	vc_ir=vc_i**r_ci;
	vh_ir=vh_i**r_hi;
	svch_ir=sqrt(vc_ir*vh_ir);

	delta_i=hate[i];
	delta_i2=delta_i*delta_i;
	tau_i2=tau_vec[i];
	tau_i4=tau_i2*tau_i2;
	tau_i6=tau_i4*tau_i2;
	tau_i8=tau_i6*tau_i2;
	dt_i2=delta_i2-tau_i2;
	A_i11=cc_i/tau_i2;
	A_i12=delta_i*ch_i/tau_i4;
	A_i22=(2*delta_i2-tau_i2)*hh_i/(2*tau_i6);
	*print "--================";
	*print A_i11 A_i12 A_i22;
	A_mat=A_mat+( (A_i11||A_i12)//((A_i12`)||A_i22) );

	B_i11=delta_i2*cc_i/tau_i4;
	B_i12=delta_i*dt_i2*ch_i/(2*tau_i6);
	B_i22=dt_i2*dt_i2*hh_i/(4*tau_i8);
	B_i21=B_i12`;
	B_mat0=B_mat0+( (B_i11||B_i12)//((B_i21)||B_i22) );
	B_i11_2=B_i11*vc_i; B_i12_2=B_i12*svch_i;
	B_i22_2=B_i22*vh_i; B_i21_2=B_i12_2`;
	B_mat2=B_mat2+( (B_i11_2||B_i12_2)//((B_i21_2)||B_i22_2) );

	B_i11_3=B_i11*vc_i2; B_i12_3=B_i12*vch_i;
	B_i22_3=B_i22*vh_i2; B_i21_3=B_i12_3`;
	B_mat3=B_mat3+( (B_i11_3||B_i12_3)//((B_i21_3)||B_i22_3) );

	B_i11_4=B_i11*vc_ir; B_i12_4=B_i12*svch_ir;
	B_i22_4=B_i22*vh_ir; B_i21_4=B_i12_4`;
	B_mat4=B_mat4+( (B_i11_4||B_i12_4)//((B_i21_4)||B_i22_4) );

end;
A_in=inv(A_mat);

Cov_SW0=A_in*B_mat0*A_in;
Var_SW0=vecdiag(Cov_SW0);
SE_SW0=sqrt(Var_SW0);

SE_SW1=sqrt(n/(n-p))*SE_SW0;

Cov_SW2=A_in*B_mat2*A_in;
Var_SW2=vecdiag(Cov_SW2);
SE_SW2=sqrt(Var_SW2);

Cov_SW3=A_in*B_mat3*A_in;
Var_SW3=vecdiag(Cov_SW3);
SE_SW3=sqrt(Var_SW3);

Cov_SW4=A_in*B_mat4*A_in;
Var_SW4=vecdiag(Cov_SW4);
SE_SW4=sqrt(Var_SW4);

finish;
*========================================================*;
*          evaluateing the R-square of beta predicted by u;
start Rsquare(u,gamma,sig, R2);
	n=nrow(u); 
	scov_u=u`*(i(n)-j(n,n,1)/n)*u/(n-1);
	
	sig_model=gamma`*scov_u*gamma;
	R2=sig_model/(sig_model+sig);
print "R2=" R2;
finish;
*========================================================*;
*========================================================*;
*========================================================*;
start NML(y, C_mat,H_mat,gamma_1,sig_1);
n=nrow(C_mat); 
p_c=ncol(C_mat);
p_h=ncol(H_mat);

C_matt=C_mat`;
CC_mat=C_matt*C_mat;
CC_in=inv(CC_mat);
H_matt=H_mat`;
epsilon=.0001;

print "==========================================================";
print "                   Least Square regression with SEs";
Cy=C_mat`*y;
gamma_LS=CC_in*Cy;
hy=C_mat*gamma_LS;
hate=y-hy;
hatee=hate#hate;
RSS=sum(hatee);
hsig2=Rss/(n-p_c);
hsig2_LS=RSS/n;
BIC_LS=n+n*log(hsig2_LS)+(p_c+1)*log(n);
print "BIC_LS=" BIC_LS;
Cov_LS=hsig2*CC_in;
Var_LS=vecdiag(Cov_LS);
SE_LS=sqrt(Var_LS);
run sandwich_LS(hatee,CC_in,C_mat,SE_SW0,SE_SW1,SE_SW2,SE_SW3,SE_SW4);
			z_ls=gamma_LS/SE_LS;
			z_sw0=gamma_LS/SE_SW0;
			z_sw1=gamma_LS/SE_SW1;
			z_sw2=gamma_LS/SE_SW2;
			z_sw3=gamma_LS/SE_SW3;
			z_sw4=gamma_LS/SE_SW4;
print "---------------------------------LS----------------------------";
print "hgamma, SE_LS, z_LS, SE_w0, z_sw0, SE_w1, z_sw1,SE_w2, z_sw2,SE_w3, z_sw3,SE_w4, z_sw4";
tt_LS=gamma_LS||SE_LS||z_LS||SE_sw0||z_sw0||SE_sw1||z_sw1||SE_sw2||z_sw2||SE_sw3||z_sw3||SE_sw4||z_sw4;
print tt_LS [format=F7.3];
print "hat.sigma^2=" hsig2;


print "==========================================================";
print "                   NML with SEs";

gamma_0=gamma_LS;
sumH=j(1,n,1)*H_mat;
sig_0t=hsig2*sumH/sum(sumH); *starting value;
sig_0=sig_0t`;
*ML by IRLS;
  n_it=0;	
  delta=.01;
  do i=1 to 300 while (delta>epsilon);
  		n_it=n_it+1;
			tau_vec=H_mat*sig_0;
			W_mat1=diag(j(n,1,1)/tau_vec);
			CW_mat=C_matt*W_mat1;
			CWC=CW_mat*C_mat;
			CWC_in=inv(CWC);
			CWy=CW_mat*y;
			gamma_1=CWC_in*CWy;
			tau_vec2=tau_vec#tau_vec;
			W_mat2=diag(j(n,1,1)/(2*tau_vec2));
			hatmu=C_mat*gamma_1;
			hate=y-hatmu;
			hatee=hate#hate;
			HW=H_matt*W_mat2;
			HWH=HW*H_mat;
			HWH_in=inv(HWH);
			HWee=HW*hatee;
			sig_1=HWH_in*HWee;
			delta=max(abs(gamma_1-gamma_0),abs(sig_1-sig_0));
			gamma_0=gamma_1;
			sig_0=sig_1;
	end; *(iterations for IRLS);
  if n_it<300 then do;
			print "number of iterations in IRLS=" n_it;

		sum1=sum(log(tau_vec));
		sum2=sum(hatee/tau_vec);
		BIC_ML=sum1+sum2+(p_c+p_h)*log(n);
		print "BIC_ML=" BIC_ML;

	  	htheta=gamma_1//sig_1;
  		SE_gamma=sqrt(vecdiag(CWC_in));
  		SE_sig=sqrt(vecdiag(HWH_in));
  		SE_ML=SE_gamma//SE_sig;
  		run sandwich_ML(hate,tau_vec,CC_in,C_mat,H_mat, 
                 SE_SW0,SE_SW1,SE_SW2,SE_SW3,SE_SW4);

			z_ml=htheta/SE_ML;
			z_sw0=htheta/SE_SW0;
			z_sw1=htheta/SE_SW1;
			z_sw2=htheta/SE_SW2;
			z_sw3=htheta/SE_SW3;
			z_sw4=htheta/SE_SW4;
   end;

print "================================================================";
print "---------------------------------ML----------------------------";

print "htheta, SE_ml, z_ml, SE_w0, z_sw0, SE_w1, z_sw1,SE_w2, z_sw2,SE_w3, z_sw3,SE_w4, z_sw4";
tt_ML=htheta||SE_ml||z_ml||SE_sw0||z_sw0||SE_sw1||z_sw1||SE_sw2||z_sw2||SE_sw3||z_sw3||SE_sw4||z_sw4;
print tt_ML [format=F7.3];
print "================================================================";

finish;
*========================================================*;

****only this section needs respecification****;
use yxu;
read all var _num_ into v;
n=nrow(v); 

y=v[,1];x1=v[,2];x2=v[,3]; 
x0=j(n,1,1);u0=j(n,1,1);
u1=v[,4];u2=v[,5];u3=v[,6];

*potential predictors of y;
xu_00=x0#u0; xu_01=x0#u1; xu_02=x0#u2; xu_03=x0#u3;
xu_10=x1#u0; xu_11=x1#u1; xu_12=x1#u2; xu_13=x1#u3;
xu_20=x2#u0; xu_21=x2#u1; xu_22=x2#u2; xu_23=x2#u3;

*selected predictors for y;
C_mat=
  xu_00||xu_01||xu_02 
||xu_10||xu_12||xu_13
||xu_20||xu_21||xu_22||xu_23;

*selected predictors for variance;
H_mat=(x1#x1)||(x2#x2)||(2*x1#x2)||(2*x1#x0)||(2*x0#x2)||j(n,1,1);
run NML(y, C_mat,H_mat,hgamma,hsig);

*R-square for beta_i1;
u=u2||u3; gamma=hgamma[5:6]; sig=hsig[1];
run Rsquare(u,gamma,sig, R2);

*R-square for beta_i2;
u=u1||u2||u3; gamma=hgamma[8:10]; sig=hsig[2];
run Rsquare(u,gamma,sig, R2);