c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                Copyright (c) 1993                    ######
c     ######      Northeast Parallel Architectures Center         ######
c     ######     Syracuse University.  All rights reserved.       ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################

c#######################################################################
c
c     File: models.f
c
c     AUTHORS: originally written by Thomas J. Finucane;
c              rewritten and documented by
               Gang Cheng (gcheng@npac.syr.edu, (315)443-2083 )
c     06/07/91.   
c
c     MODIFICATION HISTORY:
c
c     06/06/92 (G. Cheng)
c       Added documentation.
C
c#######################################################################


c ***************  Option Price Models (Fortran 77 version) **************
c
c This Fortran 77 set contains 4 option pricing models with the 
c subroutine names:
c
c     bscall   ----  Black-Scholes Model
c     binomc   ----  Constant Volatility American Call Binomial Model
c     asvcorp2 ----  Stochastic Volatility American Call Binomial Model
c     svcorp2  ----  Stochastic Volatility European Call Binomial Model
c     
c ------------------------------------------------------------------------
c Model parameters:  meaning (name = typical value)
c
c     initial stock price (sp = 90)    
c     exercise (striking) price (sk = 100)   
c     risk-free rate (r = 0.15)   
c     time to maturity (t = 0.25)   
c     time to dividend (early exercise) (tdiv = 0.125)
c     initial volatility (v = estimated)    
c     variance of stock volatility (psi = estimated)
c     correlation between stock price and its volatility (cor = estimated)
c     drift of the variance (xmu = 0)
c     dividend amount (div = 1.1)
c     number of discreted time steps (reps = 17)
c
c ************************************************************************


c---------------------------------------------------------------------
c                         Black-Scholes Model
c---------------------------------------------------------------------
c     Routine name           - bscall 
c---------------------------------------------------------------------
c         Purpose            - This subroutine calculates the value of
c                              european call option using the formula
c                              derived by black and scholes for a
c                              payout protected european call. no
c                              adjustment is made for dividends.
c                              (black-scholes model)
c---------------------------------------------------------------------
c     Arguments: 
c        sp --- initial stock price  sk --- exercise (striking) price
c        r  --- risk-free rate       t  --- time to maturity
c        v  --- initial volatility   cval --- call price (output)
c---------------------------------------------------------------------
c     Required routines - pnorml
c---------------------------------------------------------------------
      subroutine bscall(sp,sk,r,t,v,cval)
      real sp,sk,r,t,v,cval

      implicit none

      real hedge,h1,h2,p,tmp1,tmp2,pnorml
                                                                           
c     --- calculate black-scholes call value
      tmp1 = sk/exp(r*t)
      tmp2 = v*sqrt(t)
c     --- h = (log(sp/(sk/(exp(r*t)))))/(v*sqrt(t))+.5*(v*sqrt(t))
      h1 = LOG(sp/tmp1)/tmp2 + 0.5*tmp2
      h2 = h1 - tmp2
      hedge = pnorml(h1)
      p =  pnorml(h2)

      cval = sp*hedge - tmp1*p

      return
      end

c---------------------------------------------------------------------    
c     Routine name         - pnorml (used by Black-Scholes Model)         
c---------------------------------------------------------------------
c         Purpose          - This function supplies a polynomial
c                            approximation for the standard normal
c                            distribution.
c---------------------------------------------------------------------
c     Required routines - none
c---------------------------------------------------------------------
      function pnorml(z)
      real z

      real tmp1,tmp2,tmp3,tmp4
      
      tmp1 = z/sqrt(2.)
      tmp2 = abs(tmp1)
      tmp3 = 1./(1. + 0.5*tmp2)
      tmp4 = tmp3*exp(-tmp2*tmp2 - 1.26551223 + tmp3*(1.00002368 + 
     &     tmp3*(.37409196 + tmp3*(.09678418 + tmp3*(-.18628806 +
     &     tmp3*(.27886807 + tmp3*(-1.13520398 + tmp3*(1.48851587 + 
     &     tmp3*(-.82215223 + tmp3*.17087277)))))))))

      if (tmp1.lt.0.) tmp4 = 2. - tmp4

      pnorml = 1. - 0.5 * tmp4
       
      return 
      end
                                                                               
c------------------------------------------------------------------------------
c                Constant Volatility American Call Binomial Model
c------------------------------------------------------------------------------
c     Routine name    - binomc 
c------------------------------------------------------------------------------
c     Purpose         - This subroutine calculates the value of an American 
c                       call with 1 dividend using the binomial numerical
c                       procedure of Cox and Ross, and escrowed dividends. 
c------------------------------------------------------------------------------
c     Arguments: 
c         sp --- initial stock price        sk --- exercise(striking) price
c         r  --- risk-free rate             t  --- time to maturity
c         v  --- initial volatility       tdiv --- time to dividend 
c         div --- dividend amount   reps --- number of discreted time steps
c         cval  --- call price (output)
c------------------------------------------------------------------------------
c     Required routines - None
c------------------------------------------------------------------------------
      subroutine binomc(sp,sk,r,t,v,tdiv,div,reps,cval)
      real sp,sk,r,t,v,tdiv,div,cval
      integer reps

      parameter (MAXSIZE = 1000, P = 0.5)
      implicit none

      real r1,u,d,a1,a2,s0,escrow,xx,dt,tmp1,tmp2,pvdiv
      real c(MAXSIZE)
      integer i,j,n1,n2

c     ---  set formula parameters
      dt    = t/float(reps)
      r1    = (r + 1.0)**dt
      tmp1  = (r - 0.5*v*v)*dt
      tmp2  = v*SQRT(dt)
      u     = tmp1 + tmp2
      d     = tmp1 - tmp2
      n2    = INT(tdiv*365.0)
      pvdiv = div/(r1**(tdiv*float(reps)/t))
      s0    = LOG(sp - pvdiv)

c     ---  calculate binomial call value
      do 950 i = 0, reps
         
c     ---  calculate reduction in stock price due to dividends
         n1 = INT(FLOAT(reps-i)*t*365.0/FLOAT(reps))
         if (n1 .ge. n2) then
            escrow = 0.0
         else
            escrow = pvdiv*(r1**(reps-i))
         endif
         
c     ---  evaluate binomial call formula recursively
         do 900 j = 0, reps - i
            xx = FLOAT(j)*u + FLOAT(reps-i-j)*d + s0
            xx = AMIN1(xx,20.0)
            a1 = EXP(xx) + escrow - sk
            
            if (i .eq. 0) then 
               a2 = 0.0
            else
               a2 = (P*c(j+2) + (1.0-P) * c(j+1))/r1
            endif
            c(j+1) = AMAX1(a1,a2)
 900     continue
         
 950  continue

      cval = c(1)
   
      return
      end
      
c---------------------------------------------------------------------------
c               Stochastic Volatility American Call Binomial Model
c---------------------------------------------------------------------------
c       Routine name      - asvcorp2 
c---------------------------------------------------------------------------   
c       Purpose           - This subroutine calculates the value of a
c                           one-dividend American call with stochastic 
c                           volatility and a correlation with the stock
c                           price by using a binomial approx. to the
c                           lognormal volatility distribution. (full lattice)
c---------------------------------------------------------------------------- 
c       Arguments: 
c         sp --- initia stock price    sk --- exercise (striking) price
c         r  --- risk-free rate        t  --- time to maturity
c         v  --- initial volatility    tdiv --- time to dividend
c         div --- dividend amount      psi --- variance of stock volatility
c         cor --- correlation between stock price(sp) and its volatility(v)
c         xmu --- drift of the variance  
c         reps --- number of discreted time steps
c         cval  --- call price (output)
c----------------------------------------------------------------------------
c       Required routines - None
c---------------------------------------------------------------------------
      subroutine asvcorp2(sp,sk,r,t,v,tdiv,div,psi,xmu,cor,reps,cval)
      real sp,sk,r,t,v,tdiv,div,psi,xmu,cor,cval
      integer reps
 
      parameter (MAXSIZE = 262500)
      implicit none

      real sl(MAXSIZE),vl(MAXSIZE),eexerval(MAXSIZE)
      real u,d,rx,xval,ux,dx,sig,dt,divadj,tmp1,tmp2
      integer i,j,j1,iexer,reps1

c     ---  initilize parameters
      dt   = t/float(reps)
      tmp1 = (xmu-0.5*(psi*psi))*dt
      tmp2 = cor*psi*SQRT(dt)
      u    = exp(tmp1 + tmp2)
      d    = exp(tmp1 - tmp2)

c     ---  create stock price lattice(sl) and volatility lattice(vl)
      vl(1) = v*v
      sl(1) = sp
      reps1 = int(reps*tdiv/t) + 1
      if (reps1.gt.reps) reps1 = reps
      do 3000 i = 1, reps
         if (i.eq.reps1) then
            divadj = div
            iexer = i
         else
            divadj = 0.0
         endif
         do 2000 j = 2**i, 1, -2 
            j1 = j/2
            sig = vl(j1)
            tmp1 = (r-0.5*sig)*dt
            tmp2 = SQRT(sig*dt)
            ux = EXP(tmp1 + tmp2)
            dx = EXP(tmp1 - tmp2)
            if (divadj .gt. 0.0) then
               xval = sl(j1) - sk
               eexerval(j1) = AMAX1(xval,0.0)
            endif            
            sl(j)   = sl(j1)*dx - divadj
            sl(j-1) = sl(j1)*ux - divadj
            vl(j)   = vl(j1)*d
            vl(j-1) = vl(j1)*u  
 2000    continue
 3000 continue
      
c     ---  calculate call values at expiration
      
      do 4000 i = 1, 2**reps
         xval = sl(i) - sk
         sl(i) = AMAX1(xval,0.0)
 4000 continue

c     ---  set rx (periodic 1+rf) to risk free rate + 1
      rx = (r+1.0)**(t/float(reps))

c     ---  calculate call price recursively
      do 6000 i = reps, 1, -1
         do 5000 j = 2, 2**i, 2
            j1 = j/2
            sl(j1) = (0.5*(sl(j) + sl(j-1)))/rx
            if (i .eq. iexer) sl(j1) = AMAX1(sl(j1),eexerval(j1))
 5000    continue
 6000 continue

      cval = sl(1)

      return
      end

c----------------------------------------------------------------------------
c             Stochastic Volatility European Call Binomial Model
c----------------------------------------------------------------------------
c       Routine name       - svcorp2 
c----------------------------------------------------------------------------
c       Purpose           - This subroutine calculates the value of a          
c                           european call with stochastic volatility     
c                           and a correlation with the stock price by  
c                           using a binomial approx. to the lognormal   
c                           volatility distribution. (full lattice)       
c---------------------------------------------------------------------------- 
c       Arguments: 
c         sp --- initia stock price    sk --- exercise (striking) price
c         r  --- risk-free rate        t  --- time to maturity
c         v  --- initial volatility    psi --- variance of stock volatility
c         cor --- correlation between stock price(sp) and its volatility(v)
c         xmu --- drift of the variance  
c         reps --- number of discreted time steps
c         cval  --- call price (output)
c----------------------------------------------------------------------------
c       Required routines - None
c----------------------------------------------------------------------------
      subroutine svcorp2(sp,sk,r,t,v,psi,xmu,cor,reps,cval)                    
      real sp,sk,r,t,v,psi,xmu,cor,cval
      integer reps

      parameter (MAXSIZE = 524288)
      implicit none
   
      real sl(MAXSIZE),vl(MAXSIZE)
      real u,d,rx,xval,ux,dx,sig,dt,tmp1,tmp2
      integer i,j,j1

c     ---  initilize parameters
      dt   = t/float(reps)
      tmp1 = (xmu-0.5*(psi*psi))*dt
      tmp2 = cor*psi*SQRT(dt)
      u    = exp(tmp1 + tmp2)
      d    = exp(tmp1 - tmp2)
                                                                           
c     ---  create stock price lattice (sl)     
      vl(1) = v * v                           
      sl(1) = sp             
      do 3000 i = 1, reps                 
         do 2000 j = 2**i,1,-2                  
            j1 = j/2                             
            sig = vl(j1)  
            tmp1 = (r-0.5*sig)*dt
            tmp2 = SQRT(sig*dt)
            ux = EXP(tmp1 + tmp2)
            dx = EXP(tmp1 - tmp2)

            sl(j)   = sl(j1)*dx             
            sl(j-1) = sl(j1)*ux            
            vl(j)   = vl(j1)*d                       
            vl(j-1) = vl(j1)*u                      
 2000    continue
 3000 continue                                                                  

c     ---  set rx (periodic 1+rf) to risk free rate + 1   
      rx = (r + 1.0)**(t/float(reps))   

c     ---  calculate call price                            
      tmp1 = rx**FLOAT(reps)
      do 4000 i = 1, 2**reps
         xval = sl(i) - sk
         if (xval .gt. 0.0) cval = cval + xval/tmp1
 4000 continue
      
      cval = cval/float(2**reps)
      
      return                                                                    
      end    
