/* enter the following */
nf=3; /* number of files to be read */
ng=1; /*number of functions of interest in each file */
ff={"try1" "try2" "try3"}; /* names of files to be read */
outfile="out"; /* output file name */

      /* do not write anything below */
let LF=250;
let LG=20;
let ntaper={4 8 15};
g=zeros(nf,ng); 
sdnum1=zeros(nf,ng); sdnum2=zeros(nf,ng); 
sdnum3=zeros(nf,ng); sdnum4=zeros(nf,ng); 

output file=^outfile reset;

for j (1, nf, 1);
   hfile=ff[j];
   load in[]=^hfile;
   if not (rows(in) == (ng*8));
      print "
          Error in input file " hfile ": " ng "functions needed.\n"; 
      end;
   endif;
   in=reshape(in,ng,8);
   g[j,.]=(in[.,3])';
   sdnum1[j,.]=(in[.,5])'; sdnum2[j,.]=(in[.,6])';
   sdnum3[j,.]=(in[.,7])'; sdnum4[j,.]=(in[.,8])';
endfor;

for i (1, ng, 1);
   format /rdn 3, 0;
   print "\nMoment " i;
   print "N.S.E. estimate            Mean            ";;
   print "N.S.E.           Equality chi sq";
   for k (1,4,1);
      if k==1; print "i.i.d.       ";; sdnum=sdnum1;
      else;
         format /rdn 2,0; print ntaper[k-1] "% taper     ";;
         if k==2; sdnum=sdnum2; elseif k==3; sdnum=sdnum3;
         else; sdnum=sdnum4;
         endif;
      endif;
      eg=0; wtsum=0;
      for j (1,nf,1);
         eg=eg+g[j,i]/(sdnum[j,i])^2;
         wtsum=wtsum+1/(sdnum[j,i])^2;
      endfor;
      eg=eg/wtsum; nse=1/wtsum^(1/2);
      format /len 16,7;
      print "   " eg " " nse "     ";;
      gvar=zeros(nf-1,nf-1);
      for j (1,nf-2,1);
         gvar[j,j]=(sdnum[j,i])^2+(sdnum[j+1,i])^2;
         gvar[j,j+1]=-(sdnum[j+1,i])^2;
         gvar[j+1,j]=gvar[j,j+1];
      endfor; 
      j=nf-1;
      gvar[j,j]=(sdnum[j,i])^2+(sdnum[j+1,i])^2;
      ginv=invpd(gvar);
      g1=g[1:nf-1,i]; g2=g[2:nf,i];
      cstat=(g2-g1)'*ginv*(g2-g1);
      p=cdfchic(cstat,nf-1);
      print p;
   endfor;
endfor;
end;