PROGRAM MOMENT
C  This is a program for analysis of posterior simulator output.
      IMPLICIT REAL*8 (A-H,O-Z)
C     REAL*4 CPSEC
      PARAMETER(LG=1000,NG=100)
      LOGICAL LSAVE
      DIMENSION CD(NG),CN(NG),CDN(NG),CDD(NG),CNN(NG),CVAR(NG)
      DIMENSION RNN(0:NG),RDD(0:NG),RND(0:NG),RDN(0:NG)
      DIMENSION G(LG),NTAPER(3)
      DIMENSION SDOUT(4)
      DATA NTAPER/4,8,15/
      CALL NOWPR(6)
C     TIME=CPSEC()
      READ(5,*)NBURN,NSAVE
      LSAVE=.TRUE.
      IF (NSAVE.EQ.0)LSAVE=.FALSE.
      IF(LSAVE)CALL FOPEN('Moment output file',2,.TRUE.)
      CALL FOPEN('Data',1,.FALSE.)
      NRUN=0
   10 READ(5,*,END=1010,ERR=1020)JF
      REWIND(1)
      IF(JF.LE.0)GO TO 1030
      READ(1,*,END=1110,ERR=1120)NITER,NF
      IF(JF.GT.NF)GO TO 1040
      IF(NBURN.GT.NITER)GO TO 1050
      IF(NBURN.LT.0)NBURN=0
      IF(NF.GT.LG)GO TO 1060
      NRUN=NRUN+1
      DO 20 IBURN=1,NBURN
      READ(1,*,END=1130,ERR=1140)ITER,WTLOG
   20 READ(1,*,END=1150,ERR=1160)(G(IF),IF=1,NF)
      NUSE=NITER-NBURN
      IF(NUSE.LT.NG)GO TO 1170
      NS=NUSE/NG
      NUSE=NS*NG
      IF(NRUN.EQ.1)WRITE(6,30)NITER,NBURN,NUSE
   30 FORMAT(/,' Posterior moment analysis',/,
     1  '  Total iterations available:',I10,/,
     2  '          Burn-in iterations:',I10,/,
     3  '         Analysis iterations:',I10)
C  Form sufficient summary statistics for all subsequent work.
  100 TD=0.0D0
      TN=0.0D0
      TDN=0.0D0
      TDD=0.0D0
      TNN=0.0D0
      TVAR=0.0D0
      DO 120 IG=1,NG
      GN=0.0D0
      GD=0.0D0
      GDN=0.0D0
      GDD=0.0D0
      GNN=0.0D0
      GVAR=0.0D0
      DO 110 IS=1,NS
      READ(1,*,END=1130,ERR=1140)ITER,WTLOG
      READ(1,*,END=1150,ERR=1160)(G(IF),IF=1,NF)
      AD=0.0D0
      IF(WTLOG.GT.-7.0D2)AD=DEXP(WTLOG)
      AN=AD*G(JF)
      GD=GD+AD
      GN=GN+AN
      GDN=GDN+AD*AN
      GDD=GDD+AD**2
      GNN=GNN+AN**2
  110 GVAR=GVAR+AN*G(JF)
      TD=TD+GD
      TN=TN+GN
      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
  120 CVAR(IG)=GVAR/NS
C  Form posterior moment and standard deviation, and print.
  200 EG=TN/TD
      VARG=TVAR/TD-EG**2
      SDG=-1.0D0
      IF(VARG.GT.0.0D0)SDG=DSQRT(VARG)
      WRITE(6,210)JF,EG,SDG
  210 FORMAT(/,' Function',I4,' Posterior mean:',1PE16.7,
     1  ' Posterior st dev:',1PE16.7) 
      IF(LSAVE)WRITE(2,220)JF,NUSE,EG,SDG
  220 FORMAT(I4,I6,1PE16.7,1PE16.7)
C  Numerical standard error assuming no serial correlation
  300 VARNUM=(TNN-2.0D0*EG*TDN+TDD*EG**2)/TD**2
      SDNUM=-1.0D0
      IF(VARNUM.GT.0.0D0)SDNUM=DSQRT(VARNUM)
      WRITE(6,310)SDNUM,VARG/(NUSE*VARNUM)
  310 FORMAT(/,14X,'Accuracy of approximation:',9X,'Numerical',8X,
     1  'Relative',/,49X,'Standard',9X,'Numerical',/,' Method--',
     2  40X,'Error',12X,'Efficiency',/,9X,
     3  'Assuming no serial correlation:',1PE20.7,0PF15.6)
      IF(LSAVE)SDOUT(1)=SDNUM
C  Get autocovariance function of grouped means.
      BARN=TN/NUSE
      BARD=TD/NUSE
  400 DO 410 IG=1,NG
      CN(IG)=CN(IG)-BARN
  410 CD(IG)=CD(IG)-BARD
      DO 430 LAG=0,NG-1
      ANN=0.0D0
      ADD=0.0D0
      AND=0.0D0
      ADN=0.0D0
      DO 420 IG=LAG+1,NG
      ANN=ANN+CN(IG)*CN(IG-LAG)
      ADD=ADD+CD(IG)*CD(IG-LAG)
      AND=AND+CN(IG)*CD(IG-LAG)
  420 ADN=ADN+CD(IG)*CN(IG-LAG)
      RNN(LAG)=ANN/NG
      RDD(LAG)=ADD/NG
      RND(LAG)=AND/NG
  430 RDN(LAG)=ADN/NG
C  Numerical standard error with tapered autocovariance functions
  500 DO 540 MM=1,3
      M=NTAPER(MM)
      AM=M
      SNN=RNN(0)
      SDD=RDD(0)
      SND=RND(0)
      DO 510 LAG=1,M-1
      ATT=1.0D0-LAG/AM
      SNN=SNN+2.0D0*ATT*RNN(LAG)
      SDD=SDD+2.0D0*ATT*RDD(LAG)
  510 SND=SND+ATT*(RDN(LAG)+RND(LAG))
      VARNUM=NS*NUSE*(SNN-2.0D0*EG*SND+SDD*EG**2)/TD**2
      SDNUM=-1.0D0
      IF(VARNUM.GT.0.0D0)SDNUM=DSQRT(VARNUM)
  520 WRITE(6,530)M,SDNUM,VARG/(NUSE*VARNUM)
  530 FORMAT(' Autocovariance function tapered to',
     1  I3'%:',1PE20.7,0PF15.6)
  540 IF(LSAVE)SDOUT(MM+1)=SDNUM
      IF(LSAVE)WRITE(2,550)(SDOUT(I),I=1,4)
  550 FORMAT((4(1PE16.7)))
      GO TO 10
 1010 IF(NRUN.EQ.0)GO TO 1020
C     TIME=CPSEC()
C     WRITE(6,1015)NRUN,TIME,TIME/NRUN
C1015 FORMAT(/,'   Moments evaluated:',I12,/,
C    1         '            CPU time:',F12.4,/,
C    2         ' CPU time per moment:',F12.4)
      CALL NOWPR(6)
      STOP
 1020 WRITE(6,1025)
 1025 FORMAT(/,' Error on input for function number, program stops.')
      CALL TERM
 1030 WRITE(6,1035)JF
 1035 FORMAT(/,' Function of interest',I4,' negative; skipped')
      GO TO 10
 1040 WRITE(6,1045)JF,NF
 1045 FORMAT(/,' Function of interest',I4,' greater than limit of',I4,
     1  '; skipped.')
      GO TO 10
 1050 WRITE(6,1055)NBURN,NITER
 1055 FORMAT(/,I10 ' burn in requested but only',I10, ' available.')
      CALL TERM
 1060 WRITE(6,1065)NF,LG
 1065 FORMAT(/,' File contains',I5,' functions per record,'
     1  ' exceeds storeage capacity of',I5)
      CALL TERM
 1110 WRITE(6,1115)
 1115 FORMAT(/,' *** Unexpected end of simulation input file, line 1')
      CALL TERM
 1120 WRITE(6,1125)
 1125 FORMAT(/,' *** Read error in simulation input file, line 1')
      CALL TERM
 1130 WRITE(6,1135)ITER
 1135 FORMAT(/,' *** Unexpected end of simulation input file',
     1  ' at or after iteration record',I16)
      CALL TERM
 1140 WRITE(6,1145)ITER
 1145 FORMAT(/,' *** Read error at or after iteration record',I16)
      CALL TERM
 1150 WRITE(6,1155)ITER
 1155 FORMAT(/,' *** Unexpected end of simulation input file in'
     1  ,' iteration record',I16)
      CALL TERM
 1160 WRITE(6,1165)ITER
 1165 FORMAT(/,' *** Read error in simulation input file iteration'
     1  ,' record',I16)
      CALL TERM
 1170 WRITE(6,1175)NG,NUSE
 1175 FORMAT(/,' *** Minimum',I4,' analysis iterations required;'
     1  ,/,I10,' available') 
      CALL TERM
      END