(* 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="mout";  (* name of the output file *) 


(* DO NOT MODIFY ANYTHING BELOW *)

                 (* DO NOT MODIFY ANYTHING BELOW *)
Off[General::spell];
Off[General::spell1];
NG=100;
LG=1000;
<<for.m;
ntaper={4, 8, 15};
Array[cn, NG]; Array[cd, NG]; Array[cdn, NG];
Array[cdd, NG]; Array[cnn, NG]; Array[cvar, NG];
Array[rnn, NG]; Array[rdd, NG];
Array[rnd, NG]; Array[rdn, NG];
If[nsave==1, outa=OpenWrite[apmout]];
If[nburn<0, nburn=0];
nrun=0;
OpenRead::noopen="Cannot open input file";
Check[h=OpenRead[hfile], Quit[1]];
out=OpenWrite[outfile];                               

For[n=1, n <= Length[jf], n++, 
    If[NumberQ[jf[[n]]]==False,
       Write[out, "invalid function of interest\n"];
       Quit[1];
      ];   
    err=0;
    If[jf[[n]]<=0,
       WriteString[out,"function of interest negative, skipped \n"];
       err=1;
      ];             
    SetStreamPosition[hfile,0];
    Check[niter=Read[h, Number],
          WriteString[out, "read error in simulation input file, line 1 \n"];
          Quit[1]];
    If[niter==EndOfFile,
       WriteString[out,"unexpected end of simulation input file, line 1 \n"];
       Quit[1];
      ];
    Check[nf=Read[h,Number],
          WriteString[out,"read error in simulation input file, line 1 \n"];
          Quit[1]];
    If[nf==EndOfFile,
       WriteString[out,"unexpected end of simulation input file, line 1 \n"];
       Quit[1];
      ];    

    If[jf[[n]]>nf,
       WriteString[out,"function of interest greater than limit of", nf,
                   " skipped \n"];
       err=1
      ];
    If[nburn > niter,
       WriteString[out,nburn," burn in requested but only ",
                   niter, " available\n"];
       Quit[1];
      ];
    If[nf > LG,
       WriteString[out,"file contains ", nf, " functions per record, exceeds
                   storage capacity of ", LG, "\n"];
       Quit[1];
      ];

    If[err==0, 
       nrun++;
       For[iburn=1, iburn <= nburn, iburn++,
           Check[iter=Read[h, Number],
                 WriteString[out,"read error before, at or after iteration
                             record ", iter, "\n"];
                 Quit[1]];
           If[iter==EndOfFile,
              WriteString[out,"unexpected end of simulation input file before,
                          at or after iteration record ", iter, "\n"];
              Quit[1];
             ];
           Check[wtlog=Read[h, Number],
                 WriteString[out,"read error before, at or after iteration
                             record ", iter, "\n"];
                 Quit[1]];
           If[wtlog==EndOfFile,
              WriteString[out,"unexpected end of simulation input file before,
                          at or after iteration record ", iter, "\n"];
              Quit[1];
             ];
           Read[h, {Number, Number}];
           Check[g=ReadList[ h, Number, nf ], 
                 WriteString[out,"read error in simulation input file iteration
                             record ", iter, "\n"];
                 Quit[1]];
           If[Length[g] < nf,
              WriteString[out,"unexpected end of simulation input file in",
                          iter, "\n"];
              Quit[1];
             ];
          ];       

       nuse=niter-nburn;
       If[nuse < NG,
          WriteString[out, "minimum ", NG,
                " analysis iterations are required; ",
                nuse, " avaliable \n" ];
          Quit[1];
         ];  (* If nuse < NG *)
       ns=Floor[nuse / NG];

(* Format the output header *)

       If[nrun == 1,
          WriteString [out, "Posterior Moment Analysis\n"];
          WriteString [out, "Total iterations available: ", niter, "\n"];
          WriteString [out, "        Burn-in iterations: ", nburn, "\n"];
          WriteString [out, "       Analysis iterations: ", nuse, "\n\n"];
         ];

(* form sufficiency statistics for all subsequent work *)

       td=0; tn=0; tdn=0; tdd=0; tnn=0; tvar=0; 
       For[ig=1, ig <= NG, ig++,
           gd=0; gn=0; gdn=0; gdd=0; gnn=0; gvar=0;
           For[is=1, is <= ns, is++,
           Check[iter=Read[h, Number],
                 WriteString[out,"read error before, at or after iteration
                             record ", iter, "\n"];
                 Quit[1]];
           If[iter==EndOfFile,
              WriteString[out,"unexpected end of simulation input file before,
                          at or after iteration record ", iter, "\n"];
              Quit[1];
             ];
           Check[wtlog=Read[h, Number],
                 WriteString[out,"read error before, at or after iteration
                             record ", iter, "\n"];
                 Quit[1]];
           If[wtlog==EndOfFile,
              WriteString[out,"unexpected end of simulation input file before,
                          at or after iteration record ", iter, "\n"];
              Quit[1];
             ];
           Read[h, {Number, Number}];
           Check[g=ReadList[ h, Number, nf ], 
                 WriteString[out,"read error in simulation input file iteration
                             record ", iter, "\n"];
                 Quit[1]];
           If[Length[g] < nf,
              WriteString[out,"unexpected end of simulation input file in",
                          iter, "\n"];
              Quit[1];
             ];    
               ad=0;
               If[wtlog > -700, ad=Exp[wtlog]];
                  an=ad*g[[ jf[[n]] ]];
                  gd=gd+ad;
                  gn=gn+an;
                  gdn=gdn + ad*an;
                  gdd=gdd + ad*ad;
                  gnn=gnn + an*an;
                  gvar=gvar + an*g[[ jf[[n]] ]];
                 ];  (* For is *)
               td=td+gd;
               tn=tn+gn;
               tdn=tdn+gdn;
               tdd=tdd+gdd;
               tnn=tnn+gnn;
               tvar=tvar+gvar;

               cn[ig]=gn/ns;
               cd[ig]=gd/ns;
               cdn[ig]=gdn/ns;
               cdd[ig]=gdd/ns;
               cnn[ig]=gnn/ns;
               cvar[ig]=gvar/ns;

             ]; (* For ig *)
          eg=tn/td;
          varg=tvar/td - eg^2;
          sdg=-1;
          If[varg>0, sdg=Sqrt[varg]];
          WriteString[out, "\nFunction    ", jf[[n]],
                 "  Posterior mean: ", e[eg,16,7],
                 "  Posterior standard deviation ", e[sdg,16,7], "\n\n"];
          If[nsave==1, WriteString[outa,jf[[n]]," ",nuse," ",e[eg,16,7],"   ",
                                        e[sdg,16,7],"\n"]];

 (* 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]];
          WriteString[out, "               Accuracy of approximation ",
                      "     Numerical ",
                      "     Relative", "\n"];
          WriteString[out, "                                         ",
                      "     Standard  ",
                      "     Numerical  ", "\n"];
          WriteString[out, "Method--                                 ",
                      "     Error     ",
                      "     Efficiency  ", "\n"];
          WriteString[out, "        Assuming no serial correlation: ",
                      "      ",e[sdnum,16,7],
                      "      ", (varg/(nuse*varnum)), "\n"];
          If[nsave==1, WriteString[outa,e[sdnum,16,7]," "]];

(* get autocovariance of grouped means *)

          barn=tn/nuse;
          bard=td/nuse;
          For[ig=1, ig <=NG, ig++,
              cn[ig]=cn[ig]-barn;
              cd[ig]=cd[ig]-bard;
          ];
          For[lag=0, lag<=NG-1, lag++,
              ann=0; add=0; and=0; adn=0;
              For[ig=lag+1, ig <= NG, ig++,
                  ann=ann+cn[ig]*cn[ig-lag];
                  add=add+cd[ig]*cd[ig-lag];
                  and=and+cn[ig]*cd[ig-lag];
                  adn=adn+cd[ig]*cn[ig-lag];
                 ]; (* For ig *)
          (* lag+1 is for lag below *)
              rnn[lag+1]=ann/NG;    
              rdd[lag+1]=add/NG;
              rnd[lag+1]=and/NG;
              rdn[lag+1]=adn/NG;
             ]; (* For lag *)

(* numerical standard error with tapered autocovariance functions *)

          For[mm=1, mm<=3, mm++,
              m=ntaper[[mm]];
              am=m;
              snn=rnn[1]; sdd=rdd[1]; snd=rnd[1];
              For[lag=1, lag<=m-1, lag++,
                  att=1-lag/am;
                  snn=snn+2*att*rnn[lag+1];
                  sdd=sdd+2*att*rdd[lag+1];
                  snd=snd+att*(rnd[lag+1]+rnd[lag+1]);
                 ]; (* For lag *)
              varnum=ns*nuse*(snn-2*eg*snd+sdd*eg^2)/(td^2);
              sdnum=-1;
              If[varnum > 0, sdnum=Sqrt[varnum]];
              WriteString[out, "Autocovariance function tapered",
                          " to  ", PaddedForm[m,4], "%    ", e[sdnum,16,7], "      ",
                          (varg/(nuse*varnum)), "\n"];
              If[nsave==1, WriteString[outa,e[sdnum,16,7]," "]];
                                                      
             ]; (* For mm *)
          If[nsave==1, WriteString[outa,"\n"]];
      ];  (* If err==0 *)
    If[nrun==0, WriteString[out, "Error on input for function number,
                      program stops"];
      ];
   ];   (* For jf *)
Close[hfile];
Quit[1];