/* enter the following */
nburn=100; /* burn in iterations */
nsave=1; /* 1 if a machine readable output file is to be written, 0 otherwise */
apmout="momout1"; /* name of the machine readable output if nsave=1*/
hfile="sim"; /* input data file */
jf={1 2 3 4 5 11}; /* functions of interest */
outfile="momout"; /* name of the output file */

      /* do not write anything below */
let NG=100;
let ntaper={4 8 15};
cd=zeros(1,NG); cn=zeros(1,NG); cdd=zeros(1,NG); 
cdn=zeros(1,NG); cnn=zeros(1,NG); cvar=zeros(1,NG);
rnn=zeros(1,NG); rdd=zeros(1,NG); rnd=zeros(1,NG); rdn=zeros(1,NG);
if nburn < 0;
   nburn=0;
endif;


nrun=0;
output file=^outfile reset;

load in[]=^hfile;
niter=in[1];
nf=in[2];
if not (rows(in) == (niter*(4+nf) + 2));
   print "
          Error in simulation input file line 1 
          or unexpected end of file reached at line " (niter*nf +2) "
          or some of the data is not numeric \n";
   end;
endif;
g=zeros(1,nf);      /*initialize g*/

if nsave;
   egout=zeros(1,cols(jf));
   sdgout=zeros(1,cols(jf));
   sdnumout=zeros(4,cols(jf));
endif;
for i (1, cols(jf), 1);
   err=0;
   if jf[i] <= 0;
      print "function of interest negative; skipped \n";
      err=1;
   elseif jf[i]>nf;
      print "function of interest " jf[i] 
            "greater than limit of " nf "; skipped\n"; 
      err=1;
   elseif nburn > niter;
      print nburn "burn in requested but only " niter "available \n";
      end;
/* the following comment might be removed */
/* elseif nf > LG;
        print "file contains " nf "functions per record, "
              "exceeds storage capacity of " LG "\n";
*/
   endif; 
   if err /= 1;
      nrun=nrun+1;
/*
      for iburn (1, nburn,1)
          read data;
      endfor;
*/
      k=(4+nf)*nburn + 3;
      nuse=niter-nburn;
      if nuse < NG;
         print "mimimum" NG "analysis iterations are required;"
               nuse "avaliable";
         end;
      endif;
      ns=floor(nuse/NG);
      nuse=ns*NG;
      if nrun == 1;
         print "Posterior moment analysis";
         format /rdn 10, 0;
         print "Total iterations available:" niter;
         print "        Burn-in iterations:" nburn;
         print "       Analysis iterations:" nuse;
      endif;


/* form sufficiency statistics for all the subsequent work */
     td=0; tn=0; tdn=0; tdd=0; tnn=0; tvar=0;
     for ig (1,NG,1);
        gn=0; gd=0; gdn=0; gdd=0; gnn=0; gvar=0;
        for is (1,ns,1);
           iter=in[k]; k=k+1;
           wtlog=in[k]; k=k+3;
  /*           print "iter: " iter "wtlog" wtlog;
             print "g : "; */
           for ij (1,nf,1);
              g[ij]=in[k]; k=k+1;
              /* print  g[ij];*/
           endfor;
           ad=0;
           if wtlog > -709;
              ad=exp(wtlog);
           endif;
           an=ad*g[jf[i]];
           gd=gd+ad;
           gn=gn+an;
           gdn=gdn+ad*an;
           gdd=gdd+ad^2;
           gnn=gnn+an^2;
           gvar=gvar+an*g[jf[i]];
        endfor; /* is */
        td=td+gd;
        tn=tn+gn;
/* print "tn=" tn;*/
        tdn=tdn+gdn;
        tdd=tdd+gdd;
        tnn=tnn+gnn;
        tvar=tvar+gvar;
        cd[ig]=gd/ns;
        cn[ig]=gn/ns;
        cdn[ig]=gdn/ns;
        cdd[ig]=gdd/ns;
        cnn[ig]=gnn/ns;
        cvar[ig]=gvar/ns;
     endfor; /* ig */
/* form posterior moment and standard deviation, and print */
     eg=tn/td;
     varg=tvar/td-eg^2;
     sdg=-1;
     if varg > 0;
         sdg=sqrt(varg);
     endif;
     print;
     format /rdn 3,0; print  "Function" jf[i];;
     format /len 16,7;
     print " Posterior mean: " eg 
           " Posterior st dev:" sdg "\n"; 
     if nsave;
        egout[i]=eg;
        sdgout[i]=sdg;
     endif;
        
/* numerical standard error assuming no serial correlation */     
     varnum=(tnn-2*eg*tdn+tdd*eg^2)/(td^2);
     sdnum=-1;
     if varnum > 0;
        sdnum=sqrt(varnum);
     endif;
     print "Accuracy of approximation:                    Numerical   "
           "        Relative";
     print "                                              Standard    "
           "        Numerical";
     print "Method--                                      Error       "
           "        Efficiency";
     print "Assuming no serial correlation:             ";;
     print sdnum;;
     if nsave;
        sdnumout[1,i]=sdnum;
     endif;
     format /rd 15,6;
     print varg/(nuse*varnum);
/* get autocovariance of grouped means */
     barn=tn/nuse;
     bard=td/nuse;
     for ig (1,NG,1);
        cn[ig]=cn[ig]-barn;
        cd[ig]=cd[ig]-bard;
     endfor; /* ig*/
     for lag (0,NG-1,1);
        ann=0; 
        add=0; 
        adn=0;
        annd=0;
        for ig (lag+1,NG,1);
           ann=ann+cn[ig]*cn[ig-lag];
           add=add+cd[ig]*cd[ig-lag];
           annd=annd+cn[ig]*cd[ig-lag];
           adn=adn+cd[ig]*cn[ig-lag];
        endfor; /* ig */
        rnn[lag+1]=ann/NG;
        rdd[lag+1]=add/NG;
        rnd[lag+1]=annd/NG;
        rdn[lag+1]=adn/NG;
     endfor; /*lag*/

/* Numerical standard error with tapered autocovariance functions */
     for mm (1,3,1); 
        m=ntaper[mm];
        am=m;
        snn=rnn[1]; 
        sdd=rdd[1];
        snd=rnd[1];
        for lag (1,m-1,1);
           att=1-lag/am;
           snn=snn+2*att*rnn[lag+1];
           sdd=sdd+2*att*rdd[lag+1];
           snd=snd+att*(rdn[lag+1]+rnd[lag+1]);
        endfor; /*lag*/
        varnum=ns*nuse*(snn-2*eg*snd+sdd*eg^2)/(td^2);
        sdnum=-1;
        if varnum>0;
           sdnum=sqrt(varnum);
        endif;
        format /rdn 3,0; 
        print "Autocovariance function tapered to" m "%:" ;;
        format /re 19,7;
        print sdnum;;
        format /rd 16,6;
        print varg/(nuse*varnum);
        if nsave;
           sdnumout[1+mm,i]=sdnum;
        endif;
     endfor; /*mm*/
   endif; /* if err/=0 */ 
   if nrun== 0;
      print "error on input for function number, program stops";
      end;
   endif;
endfor;
if nsave;
   output file=^apmout reset;
   for i(1,cols(jf),1);
       format /rdn 8,0; print jf[i] nuse " ";;
       format /len 16,7; print egout[i] sdgout[i];
       for j(1,4,1);
          print sdnumout[j,i];;
       endfor;
       print;
   endfor;
endif;
end;