c ************  Option Price Models (Fortran 77 + MP version) ************
c
c This Fortran 77 + MP(API) set contains 2 option pricing models with 
c the (host) 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 discreted time steps (reps = 17)
c
c ************************************************************************

c----------------------------------------------------------------------------
c             Stochastic Volatility European Call Binomial Model
c----------------------------------------------------------------------------
c       Routine name       - svcorp2 (HOST program)
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 - Node Program
c----------------------------------------------------------------------------

      SUBROUTINE svcorp2(sp,sk,r,t,v,psi,xmu,cor,reps,reps1,cval)
      implicit none
      
      integer info, nodes, n
      integer mytid, tids(0:32)
      COMMON /APIENV/ nodes, mytid, tids
      
      integer msgtype1,msgtype2,msgtype3
      data msgtype1/1/ msgtype2/2/ msgtype3/3/
      
      real*4 sp,sk,r,t,v,psi,xmu,cor,cval,cval_tmp, xval
      integer reps, reps1, reps2, MAXSIZE
      parameter (MAXSIZE = 524288)
      
      real*4 sl(MAXSIZE),vl(MAXSIZE)
      real*4 u,d,rx,ux,dx,sig,dt,tmp1,tmp2
      integer i,j,j1 
      integer dest, msglen
      external FADD 
      integer nbuf(4), allgrp

c     ---  initilize parameters
      call mp_task_query(nbuf,4,3)
      allgrp = nbuf(4)
      
      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, reps1
         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  

      if (reps1 .ne. 0) then
         sl(1) = sl(mytid+1)
         vl(1) = vl(mytid+1)
      endif

      reps2 = reps - reps1
      do 3001 i = 1, reps2
         do 2001 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                      
 2001    continue
 3001 continue  

c     --- sum final values in the binomial lattice
      cval_tmp = 0.0
      do 4000 i = 1, 2**reps2
         xval = sl(i) - sk
         if (xval .gt. 0.0) cval_tmp = cval_tmp + xval
 4000 continue
    
c     ---  calculate call price 
      dest = 0
      msglen = 4
      call MP_REDUCE(cval_tmp, cval, msglen, dest, FADD, allgrp)

      if (mytid .eq. 0) then
c     ---  set rx (periodic 1+rf) to risk free rate + 1   
         rx = (r + 1.0)**(t/float(reps))   
         tmp1 = (2*rx)**reps
         cval = cval/tmp1
      endif

      return
      END

      real FUNCTION FADD(rin1, rin2, rout, msglen)
      real*4 rin1, rin2, rout
      integer msglen
      
      rout = rin1 + rin2
      
      return
      END

c---------------------------------------------------------------------------
c               Stochastic Volatility American Call Binomial Model
c---------------------------------------------------------------------------
c       Routine name      - asvcorp2 (HOST program)
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 - Node Program
c---------------------------------------------------------------------------
      subroutine asvcorp2(sp,sk,r,t,v,tdiv,div,psi,
     &     xmu,cor,reps,reps1,cval)
      
      implicit none
      
      integer info, nodes
      integer mytid, tids(0:32)
      COMMON /APIENV/ nodes, mytid, tids
       
      integer msgtype1,msgtype2,msgtype3
      data msgtype1/1/ msgtype2/2/ msgtype3/3/
      
      real*4 sp,sk,r,t,v,tdiv,div,psi,xmu,cor,cval,xval
      integer reps, reps1, reps2, MAXSIZE
      parameter (MAXSIZE = 524288)
      
      real*4 sl(MAXSIZE),vl(MAXSIZE)
      real*4 eexerval(MAXSIZE),eexerval1(MAXSIZE)
      real*4 u,d,rx,ux,dx,sig,dt,divadj,tmp1,tmp2
      integer i,j,j1,iexer,iexer1
      integer dest, msglen
      integer nbuf(4), allgrp

c     ---  initilize parameters
      call mp_task_query(nbuf,4,3)
      allgrp = nbuf(4)

      dt   = t/float(reps)
      tmp1 = (xmu-0.5*(psi*psi))*dt
      tmp2 = cor*psi*SQRT(dt)
      u    = exp(tmp1 + tmp2)
      d    = exp(tmp1 - tmp2)
      iexer = int(reps*tdiv/t) + 1
      if (iexer.gt.reps) iexer = reps
      
cc     ---  create stock price lattice (sl)     
      vl(1) = v * v                           
      sl(1) = sp             
      do 3000 i = 1, reps1
         if (i.eq.iexer) then
            divadj = div
         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  
  
      if (reps1 .ne. 0) then
         sl(1) = sl(mytid+1)
         vl(1) = vl(mytid+1)
      endif

      reps2 = reps-reps1
      iexer1 = iexer-reps1
      
      do 3001 i = 1, reps2
         if (i.eq.iexer1) then
            divadj = div
         else
            divadj = 0.0
         endif
         do 2001 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
               eexerval1(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                      
 2001    continue
 3001 continue  

c     ---  calculate call values at expiration
      do 4001 i = 1, 2**reps2
         xval = sl(i) - sk
         sl(i) = AMAX1(xval,0.0)
 4001 continue
      
c     ---  set rx (periodic 1+rf) to risk free rate + 1   
      rx = (r + 1.0)**(t/float(reps))   
      
c     --- sum final values in the binomial lattice
      do 6001 i = reps2, 1, -1
         do 5001 j = 2, 2**i, 2
            j1 = j/2
            sl(j1) = (0.5*(sl(j) + sl(j-1)))/rx
            if (i .eq. iexer1) sl(j1) = AMAX1(sl(j1),eexerval1(j1))
 5001    continue
 6001 continue
      
      xval = sl(1)
      dest = 0
      msglen = 4
      call MP_GATHER(xval,sl,msglen,dest,allgrp)
      
      if (mytid .eq. 0) then
c     --- calculate call price 
         do 6000 i = reps1, 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)  
      endif
      
      return
      END
      
      
