options nonotes;
 
options validvarname = v6;  ***This option is needed if you are
                               running the program in SAS8.  Since I wrote
                               the program using SAS6.12, I used the CALIS 
                               variable intercep (instead of intercept).
                               Without this options statement, SAS8 will not
                               recognize intercep.  Note, if you are running
                               the program using SAS 6.12, this option 
                               can be deleted;
                               


     data param;
    reliab = .75;

         gam0= 3 ; gam1= 10 ; gam2= -2 ;
         b10 = 2 ; b20 = 3 ;
         b21 = .8 ; b31 = .6 ;
         ps1= .5; ps2=.5 ; ps3=.5; ps4=.5 ;
         mf1 = 3;
         ** ph1 should be the standard deviation of phi;  ph1= 1;
          ps1 = sqrt(ph1*(1-reliab)/reliab);
              ps2 = sqrt(b21*b21*ph1*(1-reliab)/reliab);
              ps3 = sqrt(b31*b31*ph1*(1-reliab)/reliab);
              phf2 = gam1*gam1*ph1 + gam2*gam2*(2*ph1*ph1+4*ph1*mf1*mf1) + 4*gam1*gam2*mf1*ph1;
              phf2 = 11;

              ps4 = sqrt(phf2*(1-reliab)/reliab);
         run;


     run;

proc datasets nolist;
     delete saveparm save;
run;
     %macro docalis(sam,obs);
     data start;
      do i=1 to %eval(&sam*&obs);
       %do j=1 %to 4;
                 x&j=rannor(125);
                 %end;
                yf1=ranuni(125); y1=sqrt(12)*(yf1-.5);

      output;
      end;

     data a;  set start;
     if _n_=1 then set param;

     f1= mf1 + ph1*y1;
     y4= gam0 + gam1*f1 + gam2*f1*f1 + ps4*x1;


     Y1   =          f1   + ps1*x2;
     Y2   = b10+ b21*f1   + ps2*x3;
     Y3   = b20+ b31*f1   + ps3*x4;


     keep y1-y4;
     run;

     data a; set a;
            %do p=1 %to &sam;
               if %eval(1+(&p-1)*&obs)<=_n_<=%eval(&p*&obs) then iter=&p;

            %end;

     proc means data=a noprint; by iter;
     var y1-y4;
     output out=outmean mean=mean1-mean4 std = std1-std4;
     run;

     data a; set a;
      merge a outmean; by iter;
     run;

     data a; set a;
     y1c=y1-mean1;
     y2c=y2-mean2;
     y3c=y3-mean3;
     y4c=y4-mean4;

     output;
     keep iter Y1c Y2c Y3c Y4c y1 y2 y3 y4;
     run;


     PROC CALIS METHOD=ML aug uCOV  DATA=a MAXITER=5000 maxfunc=2000 outram=modelone outest=covpar noprint;
     by iter;
     var y1 y2 y3 y4;
     LINEQS
     Y1   = m1 intercep  +  F1   + E1,
     Y2   = b10 intercep + B21 F1   + E2,
     Y3   = b20 intercep + B31 F1   + E3,
     y4   = b30 intercep + b41 f1 + e4;

     STD
     E1-E4 = TH1-TH4,
     F1=PH1;

     BOUNDS
     0<=PH1,
     0<=TH1-TH4;

     Parameters m2, m3;
     b10=b21*m1 + m2;
     b20=b31*m1+ m3;


 data saveparm(keep = iter m2 m3 B21 B31 th1 th2 th3 PH1
             sem2 sem3 seB21 seB31 seth1 seth2 seth3 seph1);
     set modelone;

     if _NAME_="M2" then do; m2 = _estim_; sem2=_stderr_; end;
     if _NAME_="M3" then do; m3 = _estim_; sem3=_stderr_; end;
     if _NAME_="B21" then do; B21 = _estim_; seb21=_stderr_; end;
     if _NAME_="B31" then do; B31 = _estim_; seb31=_stderr_; end;
     if _NAME_="TH1" then do; th1 = _estim_; seth1=_stderr_; end;
     if _NAME_="TH2" then do; th2 = _estim_; seth2=_stderr_; end;
     if _NAME_="TH3" then do; th3 = _estim_; seth3=_stderr_; end;
     if _NAME_="PH1" then do; PH1 = _estim_; seph1=_stderr_; end;


 proc means noprint; by iter;
         var  m2 m3 B21 B31 th1 th2 th3 PH1 sem2 sem3 seB21 seB31 seth1 seth2 seth3 seph1;
         output out= new(keep = iter  m2 m3 B21 B31 th1 th2 th3 PH1  sem2 sem3 seB21 seB31 seth1 seth2 seth3 sePH1)
         mean= m2 m3 B21 B31 th1 th2 th3 PH1 sem2 sem3 seB21 seB31 seth1 seth2 seth3 seph1;
         run;

  data savecov(keep = iter m2 m3 b21 b31 TH1 TH2 TH3) ; set covpar;

       if (_TYPE_ = "COV") then do;
          if (_NAME_= "M2") then output;
          if (_NAME_= "M3") then output;
          if (_NAME_= "B21") then output;
          if (_NAME_= "B31") then output;
          if (_NAME_= "TH1") then output;
          if (_NAME_= "TH2") then output;
          if (_NAME_= "TH3") then output;

          end;

  data save; set new;
       merge new outmean; by iter;
       run;


   %mend docalis;

   %docalis(1000,500);  **********This (1000,500) is for running
                        the simulation 1000 times and using data with
                        sample size =  500;

      run;

 data cov; set savecov;

run;

proc iml;
 use cov;
 read all var{m2 m3 b21 b31 th1 th2 th3} into covdata;


  use save;
     read all var{TH1 TH2 TH3} into thetac;
     read all var{B21 B31} into betac;
     read all var{mean2 mean4 mean1 mean3} into means;

itn=1000;
n=500;        **********This itn =  1000 and n = 500 is for running
                        the simulation 1000 times and using data with
                        sample size =  500;

 reliab = .75;

         b0= 3 ; b1= 10 ; b2= -2 ;
         a0 = 2 ; c0 = 3 ;
         a1 = .8 ; c1 = .6;
         ps1= .5; ps2=.5 ; ps3=.5; ps4=.5 ;
         muf = 3;
         ** phisd should be the standard deviation of phi;  phisd= 1;

ps1 = sqrt(phisd*(1-reliab)/reliab);
     ps2 = sqrt(a1*a1*phisd*(1-reliab)/reliab);
     ps3 = sqrt(c1*c1*phisd*(1-reliab)/reliab);
     phf2 = b1*b1*phisd + b2*b2*(2*phisd*phisd+4*phisd*muf*muf) + 4*b1*b2*muf*phisd;
     phf2=11;

     ps4 = sqrt(phf2*(1-reliab)/reliab);

se=ps3;
su=ps1;
sa=ps2;
zeta=ps4;

do it=1 to itn;

y=j(n,1);
x=y;
z=y;

newy=y;
int=y;
do i=1 to n;

kx1 = rannor(125); kx2 = rannor(125);  kx3 = rannor(125);  kx4 = rannor(125);

y1=sqrt(12)*(ranuni(125)-.5);




f=muf+phisd#y1;

y[i,]=c0+c1#f+se#kx4;
x[i,]=f+su#kx2;
z[i,]=a0+a1#f+sa#kx3;
newy[i,]=b0+b1#f+b2#f#f+zeta#kx1;
end;
junk=y||x||z||newy;

bar=(y||x||z)`*int/n; mean1=bar[2,1]; mean2=bar[3,1]; mean3=bar[1,1];
fbar=bar[2,];


covarian=j(7,7,0);


count=it-1;
fcount=7*count;
row1=fcount+6;
row2=fcount+7;
row3=fcount+1;
row4=fcount+2;
row5=fcount+3;
row6=fcount+4;
row7=fcount+5;

covarian[1,]=covdata[row1,];
covarian[2,]=covdata[row2,];
covarian[3,]=covdata[row3,];
covarian[4,]=covdata[row4,];
covarian[5,]=covdata[row5,];
covarian[6,]=covdata[row6,];
covarian[7,]=covdata[row7,];

mm=(y||x||z)`*(y||x||z)/(n-1)-bar*bar`#n/(n-1);
a1hat=betac[it,1];
c1hat=betac[it,2];
a0hat=bar[3,]-a1hat#bar[2,];
c0hat=bar[1,]-c1hat#bar[2,];

suuhat=thetac[it,1];
saahat=thetac[it,2];
seehat=thetac[it,3];

check=(suuhat<0.0005)+(saahat<0.0005)+(seehat<0.0005);

if check=0 then do;
/* The following is used to estimate 3 and 4th moments of rho */
v1=z-a0hat-a1hat#x;
v2=y-c0hat-c1hat#x;

dep3=j(4,1,0);
dep3[1,]=sum(v1##3)/n;
dep3[2,]=sum((v1##2)#v2)/n;
dep3[3,]=sum(v1#(v2##2))/n;
dep3[4,]=sum(v2##3)/n;

projx3=j(4,3,0);
projx3[1,1]=1; projx3[1,3]=0; projx3[2,3]=0; projx3[3,3]=0;
projx3[1,2]=-a1hat**3;
projx3[2,2]=(-a1hat**2)*(c1hat);
projx3[3,2]=(-c1hat**2)*(a1hat);
projx3[4,1]=0;
projx3[4,2]=-c1hat**3; projx3[4,3]=1;

expa3u3=ginv(projx3`*projx3)*projx3`*dep3;
sa3hat=expa3u3[1,1];
su3hat=expa3u3[2,1];
se3hat=expa3u3[3,1];

dep4=j(5,1,0);
dep4[1,]=sum(v1##4)/n;
dep4[2,]=sum((v1##3)#v2)/n;
dep4[3,]=sum((v1##2)#(v2##2))/n;
dep4[4,]=sum((v2##3)#v1)/n;
dep4[5,]=sum(v2##4)/n;

inter4=j(5,1,0);
inter4[1,1]=6*a1hat*a1hat*saahat*suuhat;
inter4[2,1]=3*a1hat*c1hat*saahat*suuhat;
inter4[3,1]=saahat*seehat+a1hat*a1hat*suuhat*seehat+c1hat*c1hat*saahat*suuhat;
inter4[4,1]=3*a1hat*c1hat*seehat*suuhat;
inter4[5,1]=6*c1hat*c1hat*seehat*suuhat;

projx4=j(5,3,0);
projx4[1,1]=1; projx4[1,3]=0; projx4[2,3]=0; projx4[3,3]=0; projx4[4,3]=0;
projx4[1,2]=a1hat**4;
projx4[2,2]=(a1hat**3)*c1hat;
projx4[3,2]=(a1hat**2)*(c1hat**2);
projx4[4,2]=(c1hat**3)*a1hat;
projx4[5,2]=c1hat**4; projx4[5,1]=0; projx4[5,3]=1;

expa4u4=ginv(projx4`*projx4)*projx4`*(dep4-inter4);

sa4hat=max(sa3hat**2/saahat,saahat**2,expa4u4[1,1]);
su4hat=max(su3hat**2/suuhat,suuhat**2,expa4u4[2,1]);
se4hat=max(se3hat**2/seehat,seehat**2,expa4u4[3,1]);



c=1/(suuhat#seehat#a1hat#a1hat+suuhat#saahat#c1hat#c1hat+saahat#seehat);

sr2hat=(c##2)#(((a1hat#seehat#suuhat)##2)#saahat + ((saahat#suuhat#c1hat)##2)#seehat + ((saahat#seehat)##2)#suuhat);

sr3hat=(c##3)#(((a1hat#seehat#suuhat)##3)#sa3hat + ((saahat#suuhat#c1hat)##3)#se3hat + ((saahat#seehat)##3)#su3hat);

sr4hat=(c##4)#(((a1hat#seehat#suuhat)##4)#sa4hat + ((saahat#suuhat#c1hat)##4)#se4hat + ((saahat#seehat)##4)#su4hat
       + 6#seehat##2#saahat##2#suuhat##4#a1hat##2#c1hat##2#saahat#seehat+
         6#saahat##2#suuhat##2#seehat##4#a1hat##2#saahat#suuhat+
         6#suuhat##2#seehat##2#saahat##4#c1hat##2#seehat#suuhat);




fhat=c#(seehat#suuhat#a1hat#(z-a0hat) + saahat#suuhat#c1hat#(y-c0hat) + saahat#seehat#x);

fhatsq=fhat#fhat;

hh=(int||fhat||(fhatsq-int#sr2hat));
mff= sum(fhatsq)/n;
mff2=sum(fhatsq-int#sr2hat)/n;
mhat=hh`*hh/n;
mhyhat=hh`*newy/n;

meanf=sum(fhat)/n;


sigma=j(3,3,0);
sigma[2,2]=sr2hat;
sigma[2,3]=2#meanf#sr2hat+sr3hat; sigma[3,2]=sigma[2,3];
sigma[3,3]=4#mff2#sr2hat+sr4hat + 4#meanf#sr3hat - sr2hat#sr2hat;

sigman=j(3,3,0);
sigman[2,2]=sr2hat;
sigman[2,3]=2#meanf#sr2hat; sigman[3,2]=sigman[2,3];
sigman[3,3]=4#sr2hat#mff - 2#sr2hat#sr2hat;

bhat=ginv(mhat-sigma)*mhyhat;
bhatn=ginv(mhat-sigman)*mhyhat;

/*calculating the standard errors */

littlem = hh`# (newy`//newy`//newy`);

 subn = j(3,n,0);
 do t=1 to n;
  sigmat =j(3,3,0);
  sigmat[2,2]=sr2hat;
  sigmat[2,3]=2#fhat[t,1]#sr2hat;
  sigmat[3,2]=2#fhat[t,1]#sr2hat;
  sigmat[3,3]=4#sr2hat#fhatsq[t,1]- 2#sr2hat#sr2hat;
  *sigmat[3,3]=4#sr2hat#(phihat+(bar[2,]##2))+2#sr2hat#sr2hat;
  bighn=(hh[t,])`*hh[t,] - sigmat;
  subn[,t]=bighn*bhatn;
  end;

dtn=littlem - subn;

dldf=j(3,1,0);
dldsig=j(3,1,0);

part1=j(3,1,0); part2=j(3,3,0);
part1[3,1]=sum(-newy)/n;

part2[1,3]=-1; part2[2,2]=-1; part2[2,3]=-3*sum(fhat)/n; part2[3,3]=-6*sum(fhatsq)/n+6*sr2hat;
part2[3,1]=-1; part2[3,2]=-3*sum(fhat)/n;

dldsig=part1-part2*bhatn;

sigvv=j(2,2,0);
sigvv[1,1]=saahat+a1hat*a1hat*suuhat;
sigvv[1,2]=a1hat*c1hat*suuhat;   sigvv[2,1]=a1hat*c1hat*suuhat;
sigvv[2,2]=seehat+ c1hat*c1hat*suuhat;

invsig=inv(sigvv);

modbet=(a1hat||c1hat);

dsigth1= 1 - modbet*invsig*suuhat*modbet`
           - suuhat*modbet*invsig*modbet`
           + suuhat*modbet*invsig*(modbet`*modbet)*invsig*modbet`*suuhat;


dsigth2 = suuhat*modbet*invsig*diag((1||0))*invsig`*modbet`*suuhat;
dsigth3 = suuhat*modbet*invsig*diag((0||1))*invsig`*modbet`*suuhat;

dsigb11 = -suuhat*(1||0)*invsig*suuhat*modbet` - suuhat*modbet*invsig*(1//0)*suuhat
          +suuhat*modbet*invsig*((1//0)*modbet + modbet`*(1||0))*suuhat*invsig*modbet`*suuhat;

dsigb12 = -suuhat*(0||1)*invsig*suuhat*modbet` - suuhat*modbet*invsig*(0//1)*suuhat
          +suuhat*modbet*invsig*((0//1)*modbet + modbet`*(0||1))*suuhat*invsig*modbet`*suuhat;


dsigdthe = (0||0||dsigb11||dsigb12||dsigth1||dsigth2||dsigth3);

coeff = (i(4)||j(4,3,0))//dsigdthe;
covnew = coeff*covarian*coeff`;

biggam=-suuhat*(a1hat||c1hat)*inv(diag((saahat||seehat))+suuhat*(a1hat||c1hat)`*(a1hat||c1hat));

betat=j(1,4,0);
average=j(3,4,0);

do t=1 to n;

   part1=j(3,1,0); part2=j(3,3,0);
   part1[2,1]=newy[t,];  part1[3,1]=2*fhat[t,]#newy[t,];

   part2[1,2]=1;  part2[1,3]=2*fhat[t,];  part2[2,2]= 2*fhat[t,]; part2[2,3]=3*fhatsq[t,]-3*sr2hat;
   part2[3,3]=4*fhatsq[t,]#fhat[t,]-12*fhat[t,]#sr2hat;

   part2[2,1]=1;  part2[3,1]=2*fhat[t,]; part2[3,2]=3*fhatsq[t,]-3*sr2hat;

   dldf = part1-part2*bhatn;

   betat=(biggam||biggam#x[t,]);

   dldfbeta=dldf*betat;


   average=average + dldfbeta#(1/n);

end;



climit=average||dldsig;

addition= climit*covnew*climit`;

omegahn = dtn* dtn` # (1/n)#(1/n) + addition;
vhatn = ginv(mhat-sigman)*omegahn*ginv(mhat-sigman);

segamn0=sqrt(vhatn[1,1]);
segamn1=sqrt(vhatn[2,2]);
segamn2=sqrt(vhatn[3,3]);

myy=newy`*newy/n;
cm=(myy||mhyhat`)//(mhyhat||mhat);
rcm=(ginv(root(cm)))`;
sigma0=j(1,4,0)//(j(3,1,0)||sigman);
msm=rcm*sigma0*rcm`;
ll=eigval(msm);
l=1/ll[1,];
lk=1><(l-(1/n));
a3=lk-3/n;
bhatn3=ginv(mhat-sigman#a3)*mhyhat;
a8=lk-8/n;
bhatn8=ginv(mhat-sigman#a8)*mhyhat;
a11=lk-1/n;
bhatn1=ginv(mhat-sigman#a11)*mhyhat;

sigma0=j(1,4,0)//(j(3,1,0)||sigma);
msm=rcm*sigma0*rcm`;
ll=eigval(msm);
l=1/ll[1,];
lk=1><(l-(1/n));
a8=lk-8/n;
bhat8=ginv(mhat-sigma#a8)*mhyhat;

end;

if suuhat=0 then do;
   fhat=x;
   hh=(int||fhat||(fhat##2));
   bhat=ginv(hh`*hh)*hh`*newy; biv=0//0//0;
   projhh=hh*ginv(hh`*hh)*hh`;
   sig=(newy`*(I(n)-projhh)*newy)#(1/n);
   vhat=inv(hh`*hh)#sig;
   segam0=sqrt(vhat[1,1]);
   segam1=sqrt(vhat[2,2]);
   segam2=sqrt(vhat[3,3]);
   bhatn=bhat; bhat3=bhat; bhat1=bhat; bhat8=bhat; sigma33=0; sigma23=0;
   segamn0=segam0; segamn1=segam1; segamn2=segam2; sa3hat=0; sa4hat=0; su3hat=0; su4hat=0; se3hat=0; se4hat=0;
   sr2hat=0; sr3hat=0; sr4hat=0; bhatn3=bhat; bhatn1=bhat; bhatn8=bhat;
   end;

if saahat=0 then do;
   fhat=(z-a0hat)#(1/a1hat);
   hh=(int||fhat||(fhat##2));
   bhat=ginv(hh`*hh)*hh`*newy;biv=0//0//0;
   projhh=hh*ginv(hh`*hh)*hh`;
   sig=(newy`*(I(n)-projhh)*newy)#(1/n);
   vhat=inv(hh`*hh)#sig;
   segam0=sqrt(vhat[1,1]);
   segam1=sqrt(vhat[2,2]);
   segam2=sqrt(vhat[3,3]);
   bhatn=bhat; bhat3=bhat; bhat1=bhat; bhat8=bhat; sigma33=0; sigma23=0;
   segamn0=segam0; segamn1=segam1; segamn2=segam2; sa3hat=0; sa4hat=0; su3hat=0; su4hat=0; se3hat=0; se4hat=0;
   sr2hat=0; sr3hat=0; sr4hat=0; bhatn3=bhat; bhatn1=bhat; bhatn8=bhat;

   end;

 if seehat=0 then do;
   fhat=(y-c0hat)#(1/c1hat);
   hh=(int||fhat||(fhat##2));
   bhat=ginv(hh`*hh)*hh`*newy;biv=0//0//0;
   projhh=hh*ginv(hh`*hh)*hh`;
   sig=(newy`*(I(n)-projhh)*newy)#(1/n);
   vhat=inv(hh`*hh)#sig;
   segam0=sqrt(vhat[1,1]);
   segam1=sqrt(vhat[2,2]);
   segam2=sqrt(vhat[3,3]);
   bhatn=bhat; bhat3=bhat; bhat1=bhat; bhat8=bhat; sigma33=0; sigma23=0;
   segamn0=segam0; segamn1=segam1; segamn2=segam2; sa3hat=0; sa4hat=0; su3hat=0; su4hat=0; se3hat=0; se4hat=0;
   sr2hat=0; sr3hat=0; sr4hat=0; bhatn3=bhat; bhatn1=bhat; bhatn8=bhat;

   end;


result=result//(a1hat||c1hat||suuhat||saahat||seehat||seehat||bhatn`||segamn0||segamn1||segamn2||
                sr2hat||bhatn8`||bhat`||bhat8`||mean1||mean2||mean3);end;

create res from result;
append from result;
close res;



quit;
run;


data finale;
set res; rename
                col1=a1hat col2=c1hat col3=psiuu col4=psiaa col5=psiee col16=psiee col7=gamn0 col8=gamn1 col9=gamn2
                col10=segamn0 col11=segamn1 col12=segamn2 col13= sr2hat
                col14=bhatn80 col15=bhatn81 col16=bhatn82
                col17=bhat0 col18=bhat1 col19=bhat2
                col20=bhat80 col21=bhat81 col22=bhat82 col23=mean1 col24=mean2 col25=mean3;
           run;