program cashacc
      implicit real*8(a-h,o-z)
      double precision cumnorm(625,2),firm(250,3),firmret(51,250),
     1        aggret(51),asset(4,51,250),dasset(4,51),a0(4),
     2        a1(4),targ(4),alevel(4,51),sasset(1000,4,51),
     3        vasset(1000,4,51)
C----------------------------------------------------------------------C
C cumnorm( ) - Cumulative Normal Distribution                          C
C firm( )    - Firm Level Rate of Return                               C
C firmret( ) - Firm Level Returns (Rate of Return Times Asset Level)   C
C aggret( )  - Aggregate Rate of Return                                C
C assetd( )  - Aggregate Rate Change in Asset Values                   C
C assetl( )  - Aggregate Asset Value                                   C
C fasstd( )  - Firm Level Change in Asset Level (Failure)              C
C----------------------------------------------------------------------C
      common /normdta/ cumnorm,ncum
      data targ /1.00,0.99,0.85,0.50/
      data a0 /0.01,0.00,0.00,0.00/
      data a1 /0.50,0.15,0.10,0.25/
      open(unit=15,file='FirmOut2.out')
      open(unit=16,file='Summary.out')
      open(unit=18,file='FirmOut1.out')
      open(unit=19,file='CumNorm.dta')
      nfirm=250
      nyear=51
      ndraw=1000
      ncum=607
C----------------------------------------------------------------------C
C Load initial values - Firm Parameters and Normal Cumulative Dist.    C
C----------------------------------------------------------------------C
      do 10 i=1,nfirm
        read(18,9999) (firm(i,j),j=1,3)
   10 continue
      do 20 i=1,ncum
        read(19,9995) idum,(cumnorm(i,j),j=1,2)
   20 continue
      do 300 idraw=1,ndraw
C----------------------------------------------------------------------C
C Draw aggregate return variability                                    C
C----------------------------------------------------------------------C
      do 50 i=1,nyear
   25   call rndnorm(znorm)
        if(znorm.lt.-3.50.or.znorm.gt.3.50) goto 25
        aggret(i) = (0.25*0.01153)*znorm
   50 continue
C----------------------------------------------------------------------C
C Draw returns for each firm                                           C
C----------------------------------------------------------------------C
      do 100 i=1,nyear
        do 75 j=1,nfirm
   60      call rndnorm(znorm)
           if(znorm.lt.-3.49.or.znorm.gt.3.49) goto 60
           firmret(i,j) = firm(i,1)+aggret(i)*firm(i,2)+
     1              (0.75*0.01153)*znorm
   75   continue
  100 continue
      do 125 i=1,nyear
        do 115 j=1,250
          asset(1,i,j) = 0.0
          do 110 k=2,4
            asset(k,i,j) = 0.0
            if(i.eq.1) then
              draw = rand()
              if(draw.gt.targ(k)) asset(k,1,j) = 1.0
            else
              if(asset(k,i-1,j).eq.1.0) then
                 asset(k,i,j) = 1.0
              else
                draw = rand()
                if(draw.gt.targ(k)) asset(k,i,j) = 1.0
              endif
            endif
  110     continue
  115   continue
  125 continue
      do 150 i=1,nyear
        aggret(i) = 0.0
        do 140 j=1,nfirm
          aggret(i) = aggret(i) + firmret(i,j)
  140   continue
        aggret(i) = aggret(i) / real(nfirm)
        do 145 k=1,4
  142     call rndnorm(znorm)
          if(znorm.lt.-3.49.or.znorm.gt.3.49) goto 142
          dasset(k,i) = a0(k)+a1(k)*(aggret(i)-0.02) + 0.03*znorm
  145   continue
  150 continue
      do 158 k=1,4
        alevel(k,1) = 100.0
        do 156 i=1,nyear-1
           alevel(k,i+1) = alevel(k,i)*(1.0+dasset(k,i))
  156   continue
  158 continue
C---------------------------------------------------------------------C
C Survival Statistics                                                 C
C---------------------------------------------------------------------C
      do 180 k=1,4
        do 175 i=1,nyear
          sasset(idraw,k,i) = 0.0
          do 160 j=1,nfirm
            sasset(idraw,k,i) = sasset(idraw,k,i) + (1.0-asset(k,i,j))
  160     continue
          sasset(idraw,k,i) = sasset(idraw,k,i)/real(nfirm)
  175   continue
  180 continue
      do 250 k=1,4
        do 230 j=1,nfirm
          do 210 i=1,nyear
            asset(k,i,j) = (1.0-asset(k,i,j))*alevel(k,i)
  210     continue
  230   continue
  250 continue
      do 290 k=1,4
        do 275 i=1,nyear
          vasset(idraw,k,i) = 0.0
          do 260 j=1,nfirm
            vasset(idraw,k,i) = vasset(idraw,k,i) + asset(k,i,j)
  260     continue
          vasset(idraw,k,i) = vasset(idraw,k,i)/real(nfirm)
  275   continue
  290 continue
  300 continue
      do 350 k=1,4
        do 325 idraw=1,ndraw
          write(16,9993) k,idraw,(sasset(idraw,k,i),i=1,nyear)
  325   continue
        write(16,9992)
  350 continue
      do 400 k=1,4
        do 375 idraw=1,ndraw
          write(16,9993) k,idraw,(vasset(idraw,k,i),i=1,nyear)
  375   continue
        write(16,9992)
  400 continue
 9992 format(x)
 9993 format(x,i5,x,i5,x,51(e15.8))
 9994 format(x,10f8.4)
 9995 format(x,i5,f8.4,x,f8.4)
 9996 format(x,26(e15.8))
 9997 format(x,'Asset',i5)
 9998 format(x,51(f3.1,x))
 9999 format(x,51(e15.8,x))
      end
      subroutine rndnorm(znorm)
C Subroutine to draw normal variate
      implicit real*8(a-h,o-z)
      double precision cumnorm(625,2)
      common /normdta/ cumnorm,ncum
C Draw a uniform variate
    5 draw = rand()
C If the variate is below our table redraw
      if (draw.lt.0.0002.or.draw.gt.0.9998) goto 5
C Lookup and interpolation with the cumlative normal
      do 10 i=2,ncum-1
         if (draw.lt.cumnorm(i,1)) then
           znorm = cumnorm(i-1,2) + (draw-cumnorm(i-1,1))/(cumnorm(i,1)-
     1             cumnorm(i-1,1))*(cumnorm(i,2)-cumnorm(i-1,2))
           goto 20
         endif 
   10 continue
      goto 5
   20 continue
      end