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 = .80;

         gam0= -1 ; gam1= 1 ; gam2= -1 ; gam3 =1;
         b10 = 2 ; b20 = 3 ;
         b21 = .8 ; b31 = .6 ;
         ps1= .5; ps2=.5 ; ps3=.5; ps4=.5 ;
         mf1 = 0;
         ** 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
                     + gam3*gam3*15*ph1*ph1*ph1;
              *phf2=7;
              ps4 = sqrt(phf2*(1-reliab)/reliab);
         run;


     run;

 proc iml;
   itn=1000;
   n=500;
              **********This itn =  1000 and n = 500 is creating data for
                        the simulation to run 1000 times and use data with
                        sample size =  500;
   
   f1=j(itn*n,1,0);
   hold1=j(n,1,0);

   do j= 1 to n;
    hold1[j,]= sqrt(12)*(ranuni(125)-.5);
   end;
   do i= 1 to itn;
    f1[(i-1)*n+1:i*n,]=hold1;
   end;
   f=f1;
     create fdata from f;
     append from f;
     close fdata;
   quit;
   run;



proc datasets nolist;
     delete saveparm save;
run;
     %macro docalis(sam,obs);
     data start;
      do i=1 to %eval(&sam*&obs);
      a: %do j=1 %to 4;
         x&j=rannor(125);                     /*x&j=0.70710678*(x1*x1-1);*/
         %end;
         /*y1=rannor(125);*/

      output;
      end;
    data a; merge start fdata; y1=col1; drop col1;
     data a;  set a;
     if _n_=1 then set param;

     f1= mf1 + ph1*y1;
     y4= gam0 + gam1*f1 + gam2*f1*f1 + gam3*f1*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 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 telling
                        the simulation to run 1000 times and use 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{PH1} into phic;
     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 telling
                        the simulation to run 1000 times and use data with
                        sample size =  500;
                        


 reliab = .80;

         b0= -1 ; b1= 1 ; b2= -1 ; b3=1;
         a0 = 2 ; c0 = 3 ;
         a1 = .8 ; c1 = .6;
         ps1= .5; ps2=.5 ; ps3=.5; ps4=.5 ;
         muf = 0;
         ** 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
      + b3*b3*15*phisd*phisd*phisd;

      *phf2 = 7;

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




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

do it=1 to itn;

use fdata;
  read point(1:n) var{col1} into fs;


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=fs[i,1];

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+b3#f#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);
phiout=phic[it,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;

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

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

 fhatsq=fhat#fhat;

 fhat3 = fhat#fhat#fhat;
 fhat4 = fhat#fhat#fhat#fhat;
 fhat5 = fhat#fhat#fhat#fhat#fhat;


 hh=(int||fhat||(fhatsq-int#sr2hat)||(fhat3-3#fhat#sr2hat));

 mff = sum(fhatsq)/n;
 mff2= sum(fhatsq-sr2hat)/n;
 mff3 = sum(fhat3 -3#fhat#sr2hat)/n;
 mff4 = sum(fhat4 - 6#(fhatsq-sr2hat)#sr2hat - 3#sr2hat#sr2hat)/n;

 mhat=hh`*hh/n;
 mhyhat=hh`*newy/n;

meanf=sum(fhat)/n;

 sigman=j(4,4,0);
 sigman[2,2]=sr2hat;
 sigman[2,3]=2#meanf#sr2hat;  sigman[3,2]=2#meanf#sr2hat;
 sigman[3,3]=4#mff2#sr2hat + 2#sr2hat#sr2hat;
 sigman[4,2]=3#mff2#sr2hat; sigman[2,4]=sigman[4,2];
 sigman[4,3]= 6#mff3#sr2hat + 6#meanf#sr2hat#sr2hat; sigman[3,4]=sigman[4,3];
 sigman[4,4] = 6#sr2hat#sr2hat#sr2hat+ 18#mff2#sr2hat#sr2hat + 9#mff4#sr2hat ;

 bhatn=inv(mhat-sigman)*mhyhat;

/*calculating the standard errors */

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

   subn = j(4,n,0);
   do t = 1 to n;

   sigmat3=j(4,4,0);
  sigmat3[2,2]=sr2hat;
  sigmat3[2,3]=2#fhat[t,1]#sr2hat;
  sigmat3[3,2]=2#fhat[t,1]#sr2hat;
  sigmat3[3,3]=4#sr2hat#(fhatsq[t,1]-sr2hat) + 2#sr2hat#sr2hat;
  sigmat3[4,2]=3#(fhatsq[t,1]-sr2hat)#sr2hat; sigmat3[2,4]=sigmat3[4,2];
   sigmat3[4,3]= 6#(fhat3[t,1]-3#fhat[t,1]#sr2hat)#sr2hat + 6#fhat[t,1]#sr2hat#sr2hat; sigmat3[3,4]=sigmat3[4,3];
   sigmat3[4,4] = 6#sr2hat#sr2hat#sr2hat+ 18#(fhatsq[t,1]-sr2hat)#sr2hat#sr2hat +
                  9#(fhat4[t,1] - 6#(fhatsq[t,1]-sr2hat)#sr2hat -3#sr2hat#sr2hat)#sr2hat ;

  bighn=(hh[t,])`*hh[t,] - sigmat3;
    subn[,t]=bighn*bhatn;
   end;
***********************;
dtn=littlem - subn;

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

part1=j(4,1,0); part2=j(4,4,0);
part1[3,1]=sum(-newy)/n; part1[4,1]=sum(-3#fhat#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;
part2[4,1]=part2[2,3]; part2[1,4]=part2[2,3];
part2[4,2]=part2[3,3]; part2[2,4]=part2[3,3];
part2[4,3]=-10*sum(fhat3)/n+30*meanf*sr2hat; part2[3,4]=part2[4,3];
part2[4,4]=-15*sum(fhat4)/n+90*sum(fhatsq)/n*sr2hat-45*sr2hat*sr2hat;

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

do t=1 to n;

   part1=j(4,1,0); part2=j(4,4,0);
   part1[2,1]=newy[t,];  part1[3,1]=2*fhat[t,]#newy[t,]; part1[4,1]=(3*fhatsq[t,]-3*sr2hat)#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;

   part2[4,1]=part2[2,3]; part2[1,4]=part2[2,3]; part2[4,2]=part2[3,3];   part2[2,4]=part2[3,3];

   part2[4,3]=5*fhat4[t,]-30*fhatsq[t,]#sr2hat+15*sr2hat*sr2hat; part2[3,4]=part2[4,3];
   part2[4,4]=6*fhat5[t,]-60*fhat3[t,]#sr2hat+90*fhat[t,]*sr2hat*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 = inv(mhat-sigman)*omegahn*inv(mhat-sigman);

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


myy=newy`*newy/n;
cm=(myy||mhyhat`)//(mhyhat||mhat);
rcm=(inv(root(cm)))`;
sigma0=j(1,5,0)//(j(4,1,0)||sigman);
msm=rcm*sigma0*rcm`;
ll=eigval(msm);
l=1/ll[1,];
lk=1><(l-(1/n));
a3=lk-4/n;
bhatn3=inv(mhat-sigman#a3)*mhyhat;
a8=lk-9/n;
bhatn8=inv(mhat-sigman#a8)*mhyhat;
a11=lk-1/n;
bhatn1=inv(mhat-sigman#a11)*mhyhat;
a00=lk;
bhat00=inv(mhat-sigman#a00)*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; bhat00=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; bhat00=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; bhat00=bhat;

   end;


result=result//(a1hat||c1hat||suuhat||saahat||seehat||bhatn`||segamn0||segamn1||segamn2||
                segamn3||bhatn8`||mean1||mean2||mean3||phiout||bhatn3`||bhatn1`||bhat00`||lk);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 col6=gamn0
                col7=gamn1 col8=gamn2 col9=gamn3
                col10=segamn0 col11=segamn1 col12=segamn2 col13= segamn3
                col14=bhatn80 col15=bhatn81 col16=bhatn82 col17=bhatn83
                col18=mean1 col19=mean2 col20=mean3 col21=phi col22=bhatn30 col23=bhatn31 col24=bhatn32 col25=bhatn33
                col26=bhatn10 col27=bhatn11 col28=bhatn12 col29=bhatn13 col30=bhatn00 col31=bhatn01 col32=bhatn02 col33=bhatn03
                col34=modif;
           run;