program cashacc
      implicit real*8(a-h,o-z)
      double precision cumnorm(625,2),firm(250,3)
      common /normdta/ cumnorm,ncum
      open(unit=15,file='FirmOut1.out')
      open(unit=19,file='CumNorm.dta')
C----------------------------------------------------------------------C
C cumnorm(625,2) - Table of the standard normal distribution [input]   C
C firm(250,3) - Firm outputs: (1) mean return, (2) correlation with    C
C               market variation, (3) risk aversion coefficient        C
C               [output]                                               C
C nfirm - Number of Firms to be simulate [parameter]                   C
C ncum - Number of elements in the cumulative table [parameter]        C
C znorm - Normal random variable [internal]                            C
      nfirm=250
      ncum=607
      rk1 = 0.25*0.01153
      rk2 = 0.50*0.01153
      rk3 = 2.062
      do 10 i=1,ncum
        read(19,9995) idum,(cumnorm(i,j),j=1,2)
   10 continue
C Draw firm level variables
      do 50 i=1,250
   20   call rndnorm(znorm1)
        if(znorm1.lt.-3.49.or.znorm1.gt.3.49) goto 20
        firm(i,1) = 0.03261 + rk1*znorm1
   30   call rndnorm(znorm2)
        if(znorm2.lt.-3.49.or.znorm2.gt.3.49) goto 30
        firm(i,2) = 1.00 + rk2*znorm2
   40   call rndnorm(znorm3)
        if(znorm3.lt.-3.49.or.znorm3.gt.3.49) goto 40
        firm(i,3) = -89.865 + rk3*znorm3
   50 continue
C Write Output to a FirmOut1.out
      do 60 i=1,250
        write(15,9999) (firm(i,j),j=1,3)
   60 continue
 9994 format(x,10f10.4)
 9995 format(x,i5,f8.4,x,f8.4)
 9996 format(x,f10.6,x,f10.6)
 9997 format(x,'Summary Results')
 9998 format(x,'Firm Level Results')
 9999 format(x,25(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