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     AUTHOR: Gang Cheng (gcheng@npac.syr.edu, (315)443-2083 )
c     06/07/92.   
c
c     MODIFICATION HISTORY:
c
c     08/06/92 (G. Cheng)
c       Added documentation.
C
c################################################################

c ***************  Option Price Models (Fortran 90 version) **************
c
c This Fortran 90 set contains 2 option pricing models with the 
c subroutine names:
c
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 total discreted time steps (reps = 17, reps >= 13)
c     number of discreted time steps at dividend time 
c      ( reps1 = MAX((INT(reps*tdiv1/t)+1),reps) )
c
c Note: cmf$ layout and CMPF MAP are used as MAPPING directives in 
c       CMfortran and MPfortran respectively to get best performance
c       on 8K CM2 and 8K DECmpp-12000
c ************************************************************************

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 --- total number of discreted time steps
c         reps1 --- number of discreted time steps at dividend time 
c                   ( reps1 = MAX((INT(reps*tdiv1/t)+1),reps) )
c         cval  --- call price (output)
c----------------------------------------------------------------------------
c       Required routines - None
c---------------------------------------------------------------------------
      SUBROUTINE asvcorp2(sp,sk,r,t,v,psi,xmu,cor,reps1,reps,div,cval)
      IMPLICIT NONE

      real sp,sk,r,t,v,psi,xmu,cor,div,cval
      integer reps,reps1

      real,array(1:(2**(reps-12)-1),1:8192) :: ss1,vv1
      real,array(1:2**(reps-13),1:8192) :: cearly_big
      real,array(1:8192) :: aatmp1,aatmp2,uux,ddx
      real,array(1:2**reps1) :: cearly,cback
      real u,d,dt,r1,tmp1,tmp2,rx,xreps
      integer i,j,n,nreps1,start,end,repsdiv,reps2,nreps2

cmf$  layout ss1(:serial,:news),vv1(:serial,:news)
cmf$  layout cearly_big(:serial,:news)
cmf$  layout uux(:news),ddx(:news),cearly(:news),cback(:news)
cmf$  layout aatmp1(:news),aatmp2(:news)

CMPF MAP ss1(MEMORY,ALLBITS),vv1(MEMORY,ALLBITS)
CMPF MAP cearly_big(MEMORY,ALLBITS)
CMPF MAP uux(ALLBITS),ddx(ALLBITS)
CMPF MAP aatmp1(ALLBITS),aatmp2(ALLBITS)
CMPF MAP cearly(ALLBITS),cback(ALLBITS)

      xreps = float(reps)
      dt = t/xreps
 
c 	--- SET RX (PERIODIC 1+RF) TO RISK FREE RATE + 1
      rx=(r+1.0)**dt

      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)
      vv1(1,:) = 0.0
      ss1(1,:) = 0.0
      vv1(1,1) = v*v
      ss1(1,1) = sp
      n = 1
      do i = 1,13
         aatmp1 = (r-0.5*vv1(1,:))*dt
         aatmp2 = sqrt(vv1(1,:)*dt)
	 uux = exp(aatmp1+aatmp2)
         ddx = exp(aatmp1-aatmp2)
         vv1(1,:) =  u*vv1(1,:) + 
     &      eoshift(d*vv1(1,:),dim = 1,shift = -n)
         ss1(1,:) = uux*ss1(1,:) + 
     &      eoshift(ddx*ss1(1,:),dim = 1,shift = -n)
         if (i.eq.reps1) then
             cearly = ss1(1,1:2**i) - sk
             ss1(1,1:2**i) = ss1(1,1:2**i) - div
         endif
         n = 2*n
      enddo
      if (reps.gt.13) then
         reps2 = reps - 13
         nreps1 = 2**(reps2+1) - 1
         repsdiv = reps1 - 14
         start = 2**repsdiv
         end = 2*start - 1
         
         do i = 1, nreps1/2
           aatmp1 = (r-0.5*vv1(i,:))*dt
           aatmp2 = sqrt(vv1(i,:)*dt)
	   uux = exp(aatmp1+aatmp2)
           ddx = exp(aatmp1-aatmp2)
           vv1(2*i,:) =  u*vv1(i,:)  
           vv1((2*i+1),:) = d*vv1(i,:)
           if (i.ge.start.and.i.le.end) then 
                 aatmp1 = uux*ss1(i,:)
                 aatmp2 = ddx*ss1(i,:)
                 cearly_big(2*(i-start)+1,:) = aatmp1 - sk
                 cearly_big(2*(i-start+1),:) = aatmp2 - sk 
                 ss1(2*i,:) =  aatmp1 - div
                 ss1((2*i+1),:) = aatmp2 - div
           else 
              ss1(2*i,:) = uux*ss1(i,:)
              ss1((2*i+1),:) = ddx*ss1(i,:)
           endif
           ss1(i,:) = 0.0
         enddo
         ss1 = ss1 - sk  
         where (ss1.lt.0.0) ss1 = 0.0
         if (reps1.le.13) then
             do i=nreps1/2+1,nreps1
               ss1(1,:) = ss1(1,:) + ss1(i,:) 
             enddo
             ss1(1,:) = ss1(1,:)/((2*rx)**(float(reps-13)))
             n = 1
             do i=1,13-reps1
               ss1(1,:) = ss1(1,:) + 
     &		eoshift(ss1(1,:),dim = 1,shift = n)  
               n = 2*n
             enddo
             nreps2 = 2**(13-reps1)
             cback = ss1(1,1:8192:nreps2)/
     &         ((2*rx)**float(13-reps1))
             where (cearly.gt.cback) cback = cearly
c 	--- calculate call price
          cval = SUM(cback,mask=cback.gt.0.0)/((2*rx)**(float(reps1)))
          else
             nreps2 = 2**(reps-reps1)
             cval = 0.0
             do i=1,2*start
                aatmp1 = 0.0
                n = nreps1/2 + (i-1)*nreps2 
                do j=1,nreps2
                  aatmp1 = aatmp1 + 
     &                     ss1(n+j,:)
                enddo
                aatmp1 = aatmp1/((2*rx)**(float(reps-reps1)))  
                where (cearly_big(i,:).gt.aatmp1) 
     &        		aatmp1 = cearly_big(i,:)
                cval = cval + SUM(aatmp1,mask=aatmp1.gt.0.0)
             enddo
c 	--- calculate call price
             cval = cval/((2*rx)**(float(reps1)))
          endif
      else
c 	--- calculate call price
          ss1 = ss1 - sk
          cval =  SUM(ss1(1,:),mask=ss1(1,:).gt.0.0)/((2*rx)**xreps)  
      endif

      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)
      implicit none

      real sp,sk,r,t,v,psi,xmu,cor,cval
      integer reps

      real,array(1:(2**(reps-12)-1),1:8192) :: ss1,vv1,uux,ddx,
     &                                     aatmp1,aatmp2
      integer,array(1:8192) :: index,ud
      real u,d,dt,r1,tmp1,tmp2,rx,xreps
      integer i,reps1,nreps1,n2

cmf$  layout ss1(:serial,:news),vv1(:serial,:news)
cmf$  layout uux(:serial,:news),ddx(:serial,:news)
cmf$  layout aatmp1(:serial,:news),aatmp2(:serial,:news)
cmf$  layout index(:news),ud(:news)

CMPF MAP ss1(ALLBITS,MEMORY),vv1(ALLBITS,MEMORY)
CMPF MAP uux(ALLBITS,MEMORY),ddx(ALLBITS,MEMORY)
CMPF MAP aatmp1(ALLBITS,MEMORY),aatmp2(ALLBITS,MEMORY)
CMPF MAP index(ALLBITS),ud(ALLBITS)

      xreps = float(reps)
      dt = t/xreps
      index = [0:8191]

c 	--- SET RX (PERIODIC 1+RF) TO RISK FREE RATE + 1
      rx=(r+1.0)**dt

      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)
        
      vv1(1,:) = v*v
      ss1(1,:) = sp
      n2 = 4096
      do i = 1,13
         aatmp1(1,:) = (r-0.5*vv1(1,:))*dt
         aatmp2(1,:) = sqrt(vv1(1,:)*dt)
	 uux(1,:) = exp(aatmp1(1,:)+aatmp2(1,:))
         ddx(1,:) = exp(aatmp1(1,:)-aatmp2(1,:))
         ud = index/n2
c         where (mod(ud,2) .ne. 0)
         where (BTEST(ud,POS=0))
           vv1(1,:) =  d*vv1(1,:)
           ss1(1,:) = ddx(1,:)*ss1(1,:)
         elsewhere
           vv1(1,:) =  u*vv1(1,:)
           ss1(1,:) = uux(1,:)*ss1(1,:)
         endwhere
         n2 = n2/2
      enddo

      if (reps.gt.13) then
         reps1 = reps - 13
         nreps1 = 2**(reps1+1) - 1
         
         do i = 1, nreps1/2
           aatmp1(i,:) = (r-0.5*vv1(i,:))*dt
           aatmp2(i,:) = sqrt(vv1(i,:)*dt)
	   uux(i,:) = exp(aatmp1(i,:)+aatmp2(i,:))
           ddx(i,:) = exp(aatmp1(i,:)-aatmp2(i,:))
           vv1(2*i,:) =  u*vv1(i,:)  
           vv1((2*i+1),:) = d*vv1(i,:)
           ss1(2*i,:) = uux(i,:)*ss1(i,:)
           ss1((2*i+1),:) = ddx(i,:)*ss1(i,:)
           ss1(i,:) = 0.0
         enddo
c 	--- calculate call price     
         ss1 = ss1 - sk  
         cval =  SUM(ss1,mask=ss1.gt.0.0)/((2*rx)**xreps)
      else
c 	--- calculate call price
         ss1 = ss1 - sk
         cval =  SUM(ss1(1,:),mask=ss1(1,:).gt.0.0)/((2*rx)**xreps)  
      endif

      return
      end
