c***********************************************************************
c                                                            
c  Monte Carlo valuation of European put options on non-dividend
c  paying security.  Trinomial tree simulation on fixed grid. Transition
c  probabilities from explicit finite difference formulas.  Simplest
c  linear congruential random number generator (use MULT = 65539 for 
c  old IBM rand, or 630360016 from Marse & Roberts (1983) if desired).
c
c  Remarks: Jun 2005 R. Jones
c           Note exercise price fixed at XPRICE = 50.00 here
c
c***********************************************************************

      implicit  double precision (A-H,K-L,O-Z)
      integer   SEED, MULT, MODLUS
      parameter ( N = 101 )
      dimension U(N), DIS(N), PARM(15), INO(N), IUP(N)

C**** SET MODEL PARAMETERS *********************************************

C     Constants needed for Random Number Generator:

      SEED   = 1973272912              ! not too important
      MULT   = 452807053               ! from Dudewicz (1980): very impt
      MODLUS = 2147483647              ! largest non-neg 32 bit integer
      RNMAX  = 2d0 ** 31 - 1d0         ! real version of same

C     Grid parameters:

      data SMIN, SMAX, K / 8.0, 108.0, .025 /
      H = ( SMAX - SMIN ) / ( N - 1 )

C     Financial model parameters: Black-Scholes model, annualized parms

      data R, SIGMA, XPRICE / .10, .20, 50.0 /
      PARM(1) = SIGMA
      PARM(2) = R

C**** Setup transition probability arrays ******************************

 50   do 100 I = 1, N
         S    = SMIN + ( SMAX - SMIN ) * dble(I-1) / dble(N-1)
         A    = FNA(S,IFN,PARM)
         B    = FNB(S,IFN,PARM)
         C    = FNC(S,IFN,PARM)

C        Probabilities of moves UP, DOWN, or NO change.
         PRUP = K * (2 * A + H * B) / ( 2 * H**2 * (1 + K * C ))
         PRDN = K * (2 * A - H * B) / ( 2 * H**2 * (1 + K * C ))

C           If any negative probabilities send message
            if ( PRUP .lt. 0.0 .or. PRDN .lt. 0.0 ) then
               print*, ' *** Need finer state grid ***'
               stop
            endif
            if ( 1.0 - PRUP - PRDN .lt. 0.0 ) then
               K = K / 2
               print*, ' Changing time step to ', K
               go to 50
            endif

C           Adjust probabilities at boundaries to stay within them
            if ( I .eq. 1 ) PRDN = 0.0
            if ( I .eq. N ) PRUP = 0.0

         PRNO = 1.0 - PRUP - PRDN

C        Let's store the critical random numbers as integers for speed
         INO(I) = INT ( PRNO * RNMAX )
         IUP(I) = INT ( (PRNO + PRUP) * RNMAX )

C        This line stores discount factor for state I in DIS(I)
         DIS(I) = 1 + K * C

C        While here, let us also store terminal option payoff in U()
         U(I) = max ( 0.0D0 , XPRICE - S )
100   continue

C**** Generate random walks and accumulate outcomes ********************

150   print*, ' Enter maturity (yrs), sample size, current stock price',
     &        ' (0''s to quit):'
      read *, TMAT, NWALKS, S0
      if ( TMAT .le. 0.0 ) stop

      I0    = 1 + NINT ( (S0 - SMIN) / H )
      S0    = SMIN + H * ( I0 - 1 )
      TOTAL = 0.0
      TOT2  = 0.0

      do 250 J = 1, NWALKS
         I      = I0
         DISFAC = 1.0

C        Take one random walk out to TMAT
         do 200 T = 0.0, TMAT - K/2, K

C           Accumulate discount factor
            DISFAC = DISFAC * DIS(I)

C           Generate next random integer from the last
            SEED = SEED * MULT
            if ( SEED .lt. 0 ) SEED = SEED + MODLUS + 1

C           See which way the state changed
            if ( SEED .gt. INO(I) ) then
               if ( SEED .le. IUP(I) ) then
                  I = I + 1
                  else
                  I = I - 1
               endif
            endif
200      continue

C        Determine payoff for this walk and accumulate result
         PAY   = DISFAC * U(I)
         TOTAL = TOTAL + PAY
         TOT2  = TOT2  + PAY * PAY
250   continue

C     Calculate average payoff and standard deviation
      VALUE = TOTAL / NWALKS
      STDEV = sqrt ( TOT2 - NWALKS * VALUE ** 2 ) / NWALKS

C**** Print results to screen ******************************************

      print*
      print '('' Initial stock P: '',F8.3)', S0
      print '('' Option maturity: '',F8.3)', TMAT
      print '('' Estimated value: '',F8.3)', VALUE
      print '('' Standard deviat: '',F8.3)', STDEV
      print*
      go to 150

      end

c**** FUNCTION DEFINITIONS REQUIRED ************************************
c     The coefficients for the Black-Scholes option pricing model, with
c     S being the stock price, and R and D the assumed constant interest
c     rate and proportional dividend rate respectively, are            
c     FNA()  = SIGS * SIGS * S * S / 2.0                               
c     FNB()  = (R - D) * S                                             
c     FNC()  = -R                                                      
c***********************************************************************

      double precision function COEFF()

      implicit double precision (A-H,K-L,O-Z)
      dimension PARM(15)

      entry FNA(S,IFN,PARM)
         SIGMA  = PARM(1)
         FNA    = (SIGMA * S) ** 2 / 2.0
      return
      
      entry FNB(S,IFN,PARM)
         R      = PARM(2)
         FNB    = R * S
      return

      entry FNC(S,IFN,PARM)
         R      = PARM(2)
         FNC    = -R
      return
      end