c ************  Option Price Models (Fortran 77 + MP version) ************
c
c This Fortran 77 + MP(PVM) 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
      include '/usr/local/include/fpvm3.h'
      
      integer info, nodes, n
      integer mytid, tids(0:32)
      COMMON /PVMENV/ nodes, mytid, tids
      real*4 data(6)
      
      integer msgtype1,msgtype2,msgtype3
      data msgtype1/1/ msgtype2/2/ msgtype3/3/
      
      real sp,sk,r,t,v,psi,xmu,cor,cval,cval_tmp
      integer reps, reps1, MAXSIZE
      parameter (MAXSIZE = 32)
      
      real sl(MAXSIZE),vl(MAXSIZE)
      real u,d,rx,ux,dx,sig,dt,tmp1,tmp2
      integer i,j,j1 

c     --- Interface part to PVM
      
c     --- Computation part 
     
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     --- prepare data to be sent
      data(1) = u
      data(2) = d
      data(3) = r
      data(4) = dt
      data(5) = float(reps-reps1)
      data(6) = sk
      n = 6	

c     broadcast data to all node programs 
      call  pvmfinitsend( 0, info )
      call pvmfpack( REAL4,    data, n, 1, info ) 
      call pvmfmcast( nodes, tids, msgtype1, info )
     
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  

c     ---  set rx (periodic 1+rf) to risk free rate + 1   
      rx = (r + 1.0)**(t/float(reps))   
      tmp1 = (2*rx)**reps
      
c     --- send initial stock and volatility values to nodes
      DO 11 i = 0,nodes-1
         call pvmfinitsend( 0, info )
         call pvmfpack( REAL4, sl(i+1), 1, 1, info )
         call pvmfpack( REAL4, vl(i+1), 1, 1, info )
         call pvmfsend( tids(i), msgtype2, info ) 
 11   CONTINUE
      
c     ---  calculate call price 
      cval = 0.0
      do 30 i=1,nodes 
         call pvmfrecv( -1, msgtype3, info ) 
         call pvmfunpack( REAL4, cval_tmp, 1, 1, info )
         cval = cval + cval_tmp
 30   continue 
      cval = cval/tmp1

      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
      include '/usr/local/include/fpvm3.h'
      
      integer info, nodes, n
      integer mytid, tids(0:32)
      COMMON /PVMENV/ nodes, mytid, tids
      real*4 data(9)
      
      integer msgtype1,msgtype2,msgtype3
      data msgtype1/1/ msgtype2/2/ msgtype3/3/
      
      real sp,sk,r,t,v,tdiv,div,psi,xmu,cor,cval,xval
      integer reps, reps1, MAXSIZE
      parameter (MAXSIZE = 32)
      
      real sl(MAXSIZE),vl(MAXSIZE),eexerval(MAXSIZE)
      real u,d,rx,ux,dx,sig,dt,divadj,tmp1,tmp2
      integer i,j,j1,iexer

c     --- Interface part to PVM
      
c     --- Computation part 
      
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)
      iexer = int(reps*tdiv/t) + 1
      if (iexer.gt.reps) iexer = reps
      
c     ---  set rx (periodic 1+rf) to risk free rate + 1   
      rx = (r + 1.0)**(t/float(reps))   
      
c     --- prepare data to be sent
      data(1) = u
      data(2) = d
      data(3) = r
      data(4) = dt
      data(5) = float(reps-reps1)
      data(6) = sk
      data(7) = float(iexer-reps1)
      data(8) = rx
      data(9) = div
      n = 9	

c     broadcast data to all node programs 
      call pvmfinitsend( 0, info )
      call pvmfpack( REAL4,    data, n, 1, info ) 
      call pvmfmcast( nodes, tids, msgtype1, info )

c     ---  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  

c     --- send initial stock and volatility values to nodes
      DO 11 i = 0,nodes-1
         call pvmfinitsend( 0, info )
         call pvmfpack( REAL4, sl(i+1), 1, 1, info )
         call pvmfpack( REAL4, vl(i+1), 1, 1, info )
         call pvmfsend( tids(i), msgtype2, info ) 
 11   CONTINUE

      do 30 i=0, nodes-1 
         call pvmfrecv( tids(i), msgtype3, info ) 
         call pvmfunpack( REAL4, sl(i+1), 1, 1, info )
 30   continue 
      
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)  

      return
      END


      subroutine pvminit(who,pnodes)
      implicit none
      include '/usr/local/include/fpvm3.h'
      integer who, pnodes
      
      integer info, nodes, msgtype
      integer mytid, tids(0:32)
      character*8 arch
      COMMON /PVMENV/ nodes, mytid, tids
      
      CHARACTER*80 PRGNAM_eu,PRGNAM_am
      PARAMETER ( PRGNAM_eu='node_eu', PRGNAM_am='node_am')
      
c     ------------ Starting up all the tasks ---------------------------

c     Enroll this program in PVM 
      call pvmfmytid( mytid )
      
c     Initiate nodes instances of node program 
c      print *,'How many slave programs (1-32)?'
c      read *, nodes

      nodes = pnodes

c     If arch is set to '*' then ANY configured machine is acceptable
c     otherwise arch should be set to architecture type you wish to use.

      arch = '*'
c	for modules (node program) run on alpha nodes, this tag must be set
      arch = 'ALPHA'
      if (who.eq.1) then
         call pvmfspawn( PRGNAM_am, PVMARCH, arch, nodes, tids, info )
      else
         call pvmfspawn( PRGNAM_eu, PVMARCH, arch, nodes, tids, info )
      endif

      if (info .lt. 0) then
c         WRITE(6,*) 'Failed to load program'
         stop
      endif

      return
      end


      subroutine pvmdestroy()  

      implicit none
      include '/usr/local/include/fpvm3.h'
      
      integer info, nodes, i
      integer mytid, tids(0:32)
      COMMON /PVMENV/ nodes, mytid, tids
  
c     ---- kill the node processes
      do i = 0, nodes-1
         call pvmfkill(tids(i),info)
      enddo

c     program finished leave PVM before exiting 
      call pvmfexit(info) 
      
      return
      end

