c **********************************************************************
c
c MCarlo valuation of World Bank NASDAQ 100 capped/floored return bonds:
c    These bonds have a single payment at maturity that is linked to the
c    performance of the NASDAQ 100 stock index over the bond's life.
c    Final payment is computed as principal plus the sum of the returns
c    on the NASDAQ index over reset intervals of length DT, subject to 
c    a CAP c    on the return over each interval.  Overall maturity 
c    value is also subject to a specified FLOOR value.
c
c    For example: maturity 10 yrs, reset interval 0.25 yr,
c                 cap rate .20 (20%/yr continuously compounded),
c                 floor 1.05 (minimum payout per $1 of par)
c
c Author:  Apr 2001 R. Jones
c
c          Uses Numerical Recipes routine RAN1 to generate uniform rv's
c               and GASDEV (Box-Mueller method) to get Normal rv's
c          Model for stock index process in Black-Scholes (lognormal)
c               with parameters annualized volatility, dividend rate,
c               risk-free rate  (e.g., .25, .01, .05 cc/yr)
c
c **********************************************************************

      implicit double precision (a-h,k-l,o-z)
      parameter ( ntrials = 20000 )
      real*4 GASDEV     

 100  write(*,9002)
9002  format('Enter bond maturity (yrs) [negative to quit]    : ',\)
      read(*,*) TMAT
      if ( TMAT .le. 0 ) stop

      write(*,9000)
9000  format('Enter index volatility, div. rate, riskfree rate: ',\)
      read(*,*) sigma, div, rf
      write(*,9001)
9001  format('Enter reset interval, cap rate, lifetime floor  : ',\)
      read(*,*) dt, cap, floor

c     Transform to per reset interval values:
      s2    = sigma * sigma
      disc  = exp( - rf * TMAT )

      stub  = mod( TMAT, DT )
      nmat  = int( TMAT / DT )

c     Do Monte Carlo simulation
      idum   = -1
      ATOT   = 0d0
      ATOT2  = 0d0

      do II = 1, ntrials
         TOT = 0d0   

c        Handle full reset intervals
         do IT = 1, nmat
            z = GASDEV(idum)
            y = sqrt( s2 * DT )
            x = z * y
            x = x + ( rf - div - s2 * .5d0 ) * DT
            w = min( x, log( 1 + cap * DT ) )
            R = exp( w ) - 1d0
            TOT = TOT + R
         enddo

c        Handle stub period         
         z = GASDEV(idum)
         x = z * sqrt( s2 * stub ) 
         x = x + ( rf - div - s2 * .5d0 ) * stub 
         w = min( x, log( 1 + cap * stub ) )
         R = exp( w ) - 1d0
         TOT = TOT + R

c        Value for this trial
         TOT = TOT + 1d0
         VAL = disc * max( floor, TOT )
         ATOT = ATOT + VAL
         ATOT2 = ATOT2 + VAL * VAL
      enddo

c     Report result
      value = ATOT / ntrials
      stdev = sqrt( max( 0d0, ATOT2 / ntrials - value * value ) )
      serr  = stdev / sqrt( dble(ntrials) )
      write(*,9500) value*1d2, stdev*1d2, serr*1d2
9500  format('Value:',f8.3,3x,'Stdev:',f8.3,3x,'St.error:',f7.3)
      write(*,*)
      goto 100      
      end


c***********************************************************************
c     Function returns standard normal rv. cf. Numerical Recipes p.203
 
      FUNCTION GASDEV (IDUM)
      DATA ISET /0/
      SAVE GSET,ISET
      IF (ISET.EQ.0) THEN
   1     V1=2.*RAN1(IDUM)-1.
         V2=2.*RAN1(IDUM)-1.
         R=V1*V1+V2*V2
         IF (R.GE.1.) GO TO 1
         FAC=SQRT(-2.*LOG(R)/R)
         GSET=V1*FAC
         GASDEV=V2*FAC
         ISET=1
      ELSE
         GASDEV=GSET
         ISET=0
      ENDIF
      RETURN
      END
 
 
c**** PORTABLE RANDOM NUMBER GENERATOR FROM 'NUMERICAL RECIPES' P.196 **
c     Returns a uniform random number between 0.0 and 1.0 .  Set IDUM
c     to any negative integer to initialize the sequence.
 
      FUNCTION RAN1 (IDUM)
      DIMENSION R(97)
      PARAMETER (M1=259200, IA1=7141, IC1=54773, RM1=1./M1)
      PARAMETER (M2=134456, IA2=8121, IC2=28411, RM2=1./M2)
      PARAMETER (M3=243000, IA3=4561, IC3=51349)
      DATA IFF /0/
 
      IF (IDUM .LT. 0 .OR. IFF .EQ. 0) THEN
         IFF=1
         IX1=MOD(IC1-IDUM,M1)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX2=MOD(IX1,M2)
         IX1=MOD(IA1*IX1+IC1,M1)
         IX3=MOD(IX1,M3)
         DO 100 J=1,97
            IX1=MOD(IA1*IX1+IC1,M1)
            IX2=MOD(IA2*IX2+IC2,M2)
            R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
100      CONTINUE
         IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
      IF (J.GT.97 .OR. J.LT.1) PAUSE
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
      END