c **********************************************************
c *  Copyright (c) Northeast Parallel Architecture Center  *
c *  Syracuse University, Syracuse, NY 13244-4110          *
c *  All rights reserved                                   *
c **********************************************************

c Please be informed that this Fortran 90 version of the
c option pricing models is not well-documented
c

c  15 Feb 92, subroutine to compute dividend adjusted, implied
c  volatilities. This version returns an array of four volatilities,
c  as follows:
c
c     vol(1) -- volatility from dividend adjusted, BS model (bscalld)
c     vol(2) -- volatility from dividend adjusted, constant volatility
c                binomial method (binomc)
c     vol(3) -- volatility from dividend adjusted, stochastic volatility
c                European full binomial (svcorp2e)
c     vol(4) -- volatility from dividend adjusted, stochastic volatility
c                American full binomial (svcorp2)


      subroutine price(s,x,r,t,v,tdiv,div,psi,xmu,rho,rep,index,pr)
      implicit none
      real s,x,r,t,v,tdiv,div,psi,xmu,rho,pr,szdiv
      integer rep,index,reps1


      if (index .eq. 1) then
         szdiv = div*exp(-r*tdiv)    
         call  bscall(s-szdiv,x,r,t,v,pr)
      else if (index .eq. 2) then
         call BINOMC(s,x,r,t,v,tdiv,div,rep,pr)
      else if (index .eq. 3) then
         szdiv = div*exp(-r*tdiv)         
         call svcorp2(s-szdiv,x,r,t,v,psi,xmu,rho,rep,pr)
      else 
c         reps1 =INT(rep*tdiv/t) + 1
c         if (reps1.gt.rep) reps1 = rep
         call asvcorp2(s,x,r,t,v,tdiv,div,psi,xmu,rho,rep,pr)
      endif
      
      return
      end
      

      subroutine implied_vol(s,x,r,t,c,tdiv,div,psi,xmu,rho,reps,vol)
      implicit none

c      Use bisection method to find volatility
c      reps contains the number of intervals for methods 2,3, and 4

      real s,x,r,t,c,div,tdiv,vol(4),psi,xmu,rho
      integer i,nmax,reps(4),j,kk
      real fl,fh,f,v1,v2,v0,dv,fmid,eps,vmid
      parameter(nmax = 20, eps = 1.0e-6)

CMPF  ONFE vol
CMPF  ONFE reps


       if (div.eq.0.0) then 
          kk = 3
       else
          kk = 4
       endif

c      do j = 1,4
       do j = 1,kk
         v1 = 0.0000001
         v2 = 1.4

         call price(s,x,r,t,v2,tdiv,div,psi,xmu,rho,reps(j),j,fmid)
         fmid = fmid - c
         call price(s,x,r,t,v1,tdiv,div,psi,xmu,rho,reps(j),j,f)
         f = f - c

         if (f*fmid .gt. 0.0) then
c            print *,'implied_vol: root not bracketed'
            vol(1) = -1.0
            return

         else
            if (f .lt. 0.0) then
               v0 = v1
               dv = v2 - v1
            else
               v0 = v2
               dv = v1 - v2
            endif
         
            do i = 1, nmax
               dv = dv/2.0
               vmid = v0 + dv
         call price(s,x,r,t,vmid,tdiv,div,psi,xmu,rho,reps(j),j,fmid)
               fmid = fmid - c
               if (fmid .lt. 0.0) v0 = vmid
            enddo
            
            vol(j) = v0
         endif
      enddo
         if (div.eq.0.0) vol(4) = vol(3) 

      return
      end
