cc compare.for          2/18/92  "Compare Models"

      program draw        
      implicit none
      PARAMETER (INFILE0=1,INFILE1=2)
      PARAMETER (OUTFILE1=9)
      PARAMETER (IOfile = 'file.io')      
                                       
c      include 'models.h'

      real psi,rho
      real k,mcall1,mcall2,sz,rate1,rate2,rate,oldrate,div1,div2,div3
      integer mtime,flag,ifdiv1,ifdiv2,ifdiv3,reps1
      character*22 string1
      character*12 string2
      real t,mcall,tdiv1,szdiv,tmp
      real biasbscall,biaseur,biasamc,biasams
      real*8 mbiasbscall,mbiaseur,mbiasamc,mbiasams
      real*8 vbiasbscall,vbiaseur,vbiasamc,vbiasams
      real ceur,cams,camc,cbscall,xmu
      character*40 fname,fname1
      integer*4 total,records,recordf,nn(4),sknn(5),totalsk
      real v(4),units,ends,tmp,m,xunits,ratio
      integer reps(4),n,itmp,j,i
      real*8 me(4,4),rms(4,4),tme(4,4),trms(4,4),fileme(4),filerms(4),
     &       tfileme(4),tfilerms(4),skme(4,5),skrms(4,5)

      open(UNIT=INFILE0,FILE=IOfile,STATUS='OLD',
     &     ACCESS='SEQUENTIAL',ERR=998)
      read(INFILE0,55,ERR=999,END=234) fname1

      do 123 i = 1,4
          nn(i) = 0
          sknn(i) = 0
        do 124 j = 1,4
          tme(i,j) = 0.0
          trms(i,j) = 0.0
          skme(i,j) = 0.0
          skrms(i,j) = 0.0
 124    continue
          skme(i,5) = 0.0
          skrms(i,5) = 0.0
          sknn(5) = 0
 123    continue

      total = 0
      records = 0
      totalsk = 0

 1111  continue

      open(UNIT=INFILE1,FILE=fname1,STATUS='OLD',
     &     ACCESS='SEQUENTIAL',ERR=998)

      recordf = 0

      do 125 i = 1,4
        do 125 j = 1,4
          me(i,j) = 0.0
          rms(i,j) = 0.0
 125   continue

      do 222 i=1, 30000
777    read(INFILE1,'(A22,I3,F7.2,2(F7.2),F7.2,A12,I2,F6.2,
     & F8.4,3(I4,F8.4),4(1X,F9.4))',ERR=999,END=333) string1,mtime,k,
     & mcall1,mcall2,sz,string2,flag,rate1,rate2,ifdiv1,div1,
     & ifdiv2,div2,ifdiv3,div3,biasbscall,biasamc,biaseur,biasams
 
       total = total + 1
c       print *,'records=',records,'bias=',biasbscall
       if (mtime.eq.0.or.biasbscall.eq.0.0) goto 222
       
       if (mtime.le.15) then
         me(1,1) = me(1,1) + biasbscall
         me(2,1) = me(2,1) + biasamc
         me(3,1) = me(3,1) + biaseur
         me(4,1) = me(4,1) + biasams
         rms(1,1) = rms(1,1) + biasbscall*biasbscall
         rms(2,1) = rms(2,1) + biasamc*biasamc
         rms(3,1) = rms(3,1) + biaseur*biaseur
         rms(4,1) = rms(4,1) + biasams*biasams
         nn(1) = nn(1) + 1
       else if (mtime.le.30) then 
         me(1,2) = me(1,2) + biasbscall
         me(2,2) = me(2,2) + biasamc
         me(3,2) = me(3,2) + biaseur
         me(4,2) = me(4,2) + biasams
         rms(1,2) = rms(1,2) + biasbscall*biasbscall
         rms(2,2) = rms(2,2) + biasamc*biasamc
         rms(3,2) = rms(3,2) + biaseur*biaseur
         rms(4,2) = rms(4,2) + biasams*biasams
         nn(2) = nn(2) + 1
       else if (mtime.le.60) then
         me(1,3) = me(1,2) + biasbscall
         me(2,3) = me(2,2) + biasamc
         me(3,3) = me(3,2) + biaseur
         me(4,3) = me(4,2) + biasams
         rms(1,3) = rms(1,3) + biasbscall*biasbscall
         rms(2,3) = rms(2,3) + biasamc*biasamc
         rms(3,3) = rms(3,3) + biaseur*biaseur
         rms(4,3) = rms(4,3) + biasams*biasams
         nn(3) = nn(3) + 1
       else if (mtime.ge.61) then
         me(1,4) = me(1,4) + biasbscall
         me(2,4) = me(2,4) + biasamc
         me(3,4) = me(3,4) + biaseur
         me(4,4) = me(4,4) + biasams
         rms(1,4) = rms(1,4) + biasbscall*biasbscall
         rms(2,4) = rms(2,4) + biasamc*biasamc
         rms(3,4) = rms(3,4) + biaseur*biaseur
         rms(4,4) = rms(4,4) + biasams*biasams
         nn(4) = nn(4) + 1
         totalsk = totalsk + 1
         ratio = sz/k
         if (ratio.lt.0.85) then
           skme(1,1) = skme(1,1) + biasbscall
           skme(2,1) = skme(2,1) + biasamc
           skme(3,1) = skme(3,1) + biaseur
           skme(4,1) = skme(4,1) + biasams
           skrms(1,1) = skrms(1,1) + biasbscall*biasbscall
           skrms(2,1) = skrms(2,1) + biasamc*biasamc
           skrms(3,1) = skrms(3,1) + biaseur*biaseur
           skrms(4,1) = skrms(4,1) + biasams*biasams
           sknn(1) = sknn(1) + 1
         else if (ratio.lt.0.95) then 
           skme(1,2) = skme(1,2) + biasbscall
           skme(2,2) = skme(2,2) + biasamc
           skme(3,2) = skme(3,2) + biaseur
           skme(4,2) = skme(4,2) + biasams
           skrms(1,2) = skrms(1,2) + biasbscall*biasbscall
           skrms(2,2) = skrms(2,2) + biasamc*biasamc
           skrms(3,2) = skrms(3,2) + biaseur*biaseur
           skrms(4,2) = skrms(4,2) + biasams*biasams
           sknn(2) = sknn(2) + 1
         else if (ratio.lt.1.05) then
           skme(1,3) = skme(1,2) + biasbscall
           skme(2,3) = skme(2,2) + biasamc
           skme(3,3) = skme(3,2) + biaseur
           skme(4,3) = skme(4,2) + biasams
           skrms(1,3) = skrms(1,3) + biasbscall*biasbscall
           skrms(2,3) = skrms(2,3) + biasamc*biasamc
           skrms(3,3) = skrms(3,3) + biaseur*biaseur
           skrms(4,3) = skrms(4,3) + biasams*biasams
           sknn(3) = sknn(3) + 1
         else if (ratio.lt.1.15) then
           skme(1,4) = skme(1,4) + biasbscall
           skme(2,4) = skme(2,4) + biasamc
           skme(3,4) = skme(3,4) + biaseur
           skme(4,4) = skme(4,4) + biasams
           skrms(1,4) = skrms(1,4) + biasbscall*biasbscall
           skrms(2,4) = skrms(2,4) + biasamc*biasamc
           skrms(3,4) = skrms(3,4) + biaseur*biaseur
           skrms(4,4) = skrms(4,4) + biasams*biasams
           sknn(4) = sknn(4) + 1
         else 
           skme(1,5) = skme(1,5) + biasbscall
           skme(2,5) = skme(2,5) + biasamc
           skme(3,5) = skme(3,5) + biaseur
           skme(4,5) = skme(4,5) + biasams
           skrms(1,5) = skrms(1,5) + biasbscall*biasbscall
           skrms(2,5) = skrms(2,5) + Biasamc*biasamc
           skrms(3,5) = skrms(3,5) + biaseur*biaseur
           skrms(4,5) = skrms(4,5) + biasams*biasams
           sknn(5) = sknn(5) + 1
         endif
       endif
       records = records + 1
       recordf = recordf + 1
222   continue
 
333   continue
      close(INFILE1)


      do 225 i = 1,4
         fileme(i) = 0.0
         filerms(i) = 0.0
       do 226 j = 1,4
         fileme(i) = fileme(i) + me(i,j)
         filerms(i) = filerms(i) + rms(i,j)
         tme(i,j) = tme(i,j) + me(i,j)
         trms(i,j) = trms(i,j) + rms(i,j)
 226  continue
         fileme(i) = fileme(i)/float(recordf)
         filerms(i) = sqrt(filerms(i))/float(recordf)
 225  continue
       
       print *,'file:', fname1
       print *,'number of records processed:',recordf
       print *,'             model 1   model 2   model 3   model 4'
       write(*,10) fileme(1),fileme(2),fileme(3),fileme(4)
       write(*,11) filerms(1),filerms(2),filerms(3),filerms(4)
 10     format('Mean error ',4(F9.6,1X))
 11     format('RMS  error ',4(F9.6,1X))
       print *
    
       read(INFILE0,55,ERR=999,END=234) fname1
       goto 1111


 234   close(INFILE0)

      do 127 i = 1,4
         tfileme(i) = 0.0
         tfilerms(i) = 0.0
       do 128 j = 1,4
         tfileme(i) = tfileme(i) + tme(i,j)
         tfilerms(i) = tfilerms(i) + trms(i,j)
        if (nn(j).ne.0) then
         tme(i,j) = tme(i,j)/float(nn(j))
         trms(i,j) = sqrt(trms(i,j)/float(nn(j)))
        endif
 128  continue
         tfileme(i) = tfileme(i)/float(records)
         tfilerms(i) = sqrt(tfilerms(i))/float(records)
 127  continue

      print *,'------------------------------------------'
      print *,'Total number of records:',total
      print *,'total number of records processed:',records
      print *
      
      print *,'Overall errors:'
      print *
      print *,'Model    Mean err     RMS err'
      write(*,21) tfileme(1),tfilerms(1)
 21   format('  1     ',5(F9.6,3X))
      write(*,22) tfileme(2),tfilerms(2)
 22   format('  2     ',5(F9.6,3X))
      write(*,23) tfileme(3),tfilerms(3)
 23   format('  3     ',5(F9.6,3X))
      write(*,24) tfileme(4),tfilerms(4)
 24   format('  4     ',5(F9.6,3X))

      print *
      print *,'Breakdown by expiration time:'
      print *
      print *,'Mean errors:'
      print *,'Model   0-15 days   16-30 days  ',
     & '31-60 days   61+ days'
      print *,' (N)     ',nn(1),'        ',nn(2),
     & '      ',nn(3),'       ',nn(4)
      write(*,21) tme(1,1),tme(1,2),tme(1,3),tme(1,4)
      write(*,22) tme(2,1),tme(2,2),tme(2,3),tme(2,4)
      write(*,23) tme(3,1),tme(3,2),tme(3,3),tme(3,4)
      write(*,24) tme(4,1),tme(4,2),tme(4,3),tme(4,4)

      print *
      print *,'RMS errors:'
      print *,'Model   0-15 days   16-30 days  ',
     & '31-60 days   61+ days'
      print *,' (N)     ',nn(1),'        ',nn(2),
     & '      ',nn(3),'       ',nn(4)
      write(*,21) trms(1,1),trms(1,2),trms(1,3),trms(1,4)
      write(*,22) trms(2,1),trms(2,2),trms(2,3),trms(2,4)
      write(*,23) trms(3,1),trms(3,2),trms(3,3),trms(3,4)
      write(*,24) trms(4,1),trms(4,2),trms(4,3),trms(4,4)

      do 227 i = 1,4
       do 227 j = 1,5
        if (sknn(j).ne.0) then
         skme(i,j) = skme(i,j)/float(sknn(j))
         skrms(i,j) = sqrt(skrms(i,j)/float(sknn(j)))
        endif
 227  continue

      print *
      print *,'Breakdown by S/K ratio for long ',
     &         'maturaty time (61+ days):'
      print *,'Total number of records:',totalsk
      print *
      print *,'Mean errors:'
      print *,'Model    < 0.85    0.85-0.95   ',
     & '0.95-1.05    1.05-1.15    > 1.15'
      print *,' (N)     ',sknn(1),'       ',sknn(2),
     & '       ',sknn(3),'       ',sknn(4),'       ',sknn(5)
      write(*,21) skme(1,1),skme(1,2),skme(1,3),
     &              skme(1,4),skme(1,5)
      write(*,22) skme(2,1),skme(2,2),skme(2,3),
     &              skme(2,4),skme(2,5)
      write(*,23) skme(3,1),skme(3,2),skme(3,3),
     &              skme(3,4),skme(3,5)
      write(*,24) skme(4,1),skme(4,2),skme(4,3),
     &              skme(4,4),skme(4,5)

      print *
      print *,'RMS errors:'
      print *,'Model    < 0.85    0.85-0.95   ',
     & '0.95-1.05    1.05-1.15    > 1.15'
      print *,' (N)     ',sknn(1),'       ',sknn(2),
     & '       ',sknn(3),'       ',sknn(4),'       ',sknn(5)
      write(*,21) skrms(1,1),skrms(1,2),skrms(1,3),
     &              skrms(1,4),skrms(1,5)
      write(*,22) skrms(2,1),skrms(2,2),skrms(2,3),
     &              skrms(2,4),skrms(2,5)
      write(*,23) skrms(3,1),skrms(3,2),skrms(3,3),
     &              skrms(3,4),skrms(3,5)
      write(*,24) skrms(4,1),skrms(4,2),skrms(4,3),
     &              skrms(4,4),skrms(4,5)

c      print *, '---- Program finished ----'
c     print *, 'Total records =',total
      stop 
       
999   print *, 'read or write file error!'
      total = total + 1
c      goto 777

998   stop 'open or read file error!'

55    format (A40)
66    format (F9.6,F9.6)

      end





