cccccccccccccccccccccccccccccccc   Date of last change: 4/2/92 cccccccccc
c     routine name              -bscall
c
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
c         required routines - pnorml
c----------------------------------------------------------------------

      subroutine bscall(s,e,r,t,v,cval)
      implicit none

      real s,e,r,t,v,h, cval,hedge,p1,p2,tmp1,tmp2
                                                                               
cccc  calculate black-scholes call value

      tmp1 = e/exp(r*t)
      tmp2 = v*sqrt(t)
c      h=(log(s/(e/(exp(r*t)))))/(v*sqrt(t))+.5*(v*sqrt(t))
      h= log(s/tmp1)/tmp2 + .5*tmp2
      p2 = h - tmp2
      call pnorml(h,hedge)
      call pnorml(p2,p1)
      cval=s*hedge - tmp1*p1
      return
      end
                                                                               
c     routine name          -pnorml                                           
c----------------------------------------------------------------------
c         purpose           -this function supplies a polynomial
c                            approximation for the standard normal
c                            distribution.
c         required routines - none
c----------------------------------------------------------------------
      subroutine pnorml(z,p)

      real z,p
      
      p = 1. - 0.5*ERFCC(z/sqrt(2.))
       
      RETURN 
      END
                                                                               

      FUNCTION ERFCC(X)

      real x,z,t

      Z=ABS(X)      
      T=1./(1.+0.5*Z)

      ERFCC = T*EXP(-Z*Z-1.26551223+T*(1.00002368+T*(.37409196+
     &    T*(.09678418+T*(-.18628806+T*(.27886807+T*(-1.13520398+
     &    T*(1.48851587+T*(-.82215223+T*.17087277)))))))))

      IF (X.LT.0.) ERFCC = 2. - ERFCC

      RETURN
      END

C
C       ROUTINE NAME        -BINOMC
C
C----------------------------------------------------------------------
C
C         PURPOSE           -THIS SUBROUTINE CALCULATES THE VALUE OF AN
C                            AMERICAN CALL WITH 1 DIV USING THE BINOMIAL
C                            NUMERICAL PROCEDURE OF COX AND ROSS, AND
C                            ESCROWED DIVIDENDS.
C
C----------------------------------------------------------------------
C
      SUBROUTINE BINOMC(S,E,R,T,V,TDIV,DIV,REPS,CVAL)

      REAL S,E,R,DIV,R1,U,D,T,P,A1,A2,C(1000)
      REAL V,CVAL,TDIV,S0,ESCROW,XX,DT,BIGV
      INTEGER REPS,I,J,N1,N2
C
CCCCC  CHECK FOR REASONABLE PARAMETERS
C
      CVAL=0.0

CCCC  SET FORMULA PARAMETERS
C
  500 CONTINUE
      DT=T/FLOAT(REPS)
      R1=(R+1.0)**(DT)
      P=0.50
      U=EXP((R-0.5*(V**2.0))*(DT)+V*SQRT(DT))
      D=EXP((R-0.5*(V**2.0))*(DT)-V*SQRT(DT))
      N2=INT(TDIV*365.0)
      PVDIV=DIV/(R1**(TDIV*FLOAT(REPS)/T))
      S0=S-PVDIV
C
CCCC  CALCULATE BINOMIAL CALL VALUE
C  
      L0=REPS
      I=-1
  900 I=I+1
      N1=INT(FLOAT(REPS-I)*T*365.0/FLOAT(REPS))
      IF (N1 .GE. N2) THEN
      ESCROW=0.0
      ELSE
      ESCROW=PVDIV*(R1**FLOAT(REPS-I))
      ENDIF
C
CCCC  CALCULATE REDUCTION IN STOCK PRICE DUE TO DIVIDENDS
C

C
CCCC  EVALUATE BINOMIAL CALL FORMULA RECURSIVELY
C
      L1=REPS-I
      J=-1
  950 J=J+1
      XX=(FLOAT(J)*LOG(U)+FLOAT(REPS-I-J)*LOG(D)+LOG(S0))

       IF (XX .GT. 20) XX=20.0
       A1=EXP(XX)+ESCROW-E

      IF (I. GT. 0) GOTO 1100
      A2=0.0
      GOTO 1200
 1100 A2=(P*C(J+2)+(1.0-P)*C(J+1))/R1
 1200 C(J+1)=AMAX1(A1,A2)

 1300 IF(J .LE. L1) GOTO 950
 1400 IF(I .LT. L0) GOTO 900
      CVAL=C(1)
   
      RETURN
      END

C
C
C     ROUTINE  NAME          - LPROB
C
C----------------------------------------------------------------------
C
C         PURPOSE            - THIS FUNCTION EVALUATES THE FUNCTION:
C
C               LN(P**REPS*N!/(J!(N-J)!)) , P=0.5
C
C         REQUIRED ROUTINES - NONE
C
C----------------------------------------------------------------------
C
C
      SUBROUTINE LPROB(N,J,LP)
      REAL N,J,XL,XH,NJ,X2,X,LF,LP,P
      P=0.5
      LF=0.0
      XH=N
      NJ=N-J
      IF ((J .LT. 0.5) .OR. (NJ .LT. 0.5)) GOTO 200
      IF (NJ .LT. J) GOTO 50
      XL=NJ+1
      GOTO 75
   50 XL=J+1
   75 CONTINUE
      IF (XL .GT. N) GOTO 200
      X=XL
      X2=(XH-XL+1.0)
  100 CONTINUE
      LF=LF+LOG(X)-LOG(X2)
      X=X+1.0
      X2=X2-1.0
      IF(X .LE. XH) GOTO 100
  200 CONTINUE
      LP=LF+N*LOG(P)
      RETURN
      END


C
C
C       ROUTINE NAME       - SVCORP2
C
C----------------------------------------------------------------------
C
C       PURPOSE           - THIS SUBROUTINE CALCULATES THE VALUE OF A
C                           EUROPEAN CALL WITH STOCHASTIC VOLATILITY
C                           AND A CORRELATION OF 1 OR -1 WITH THE STOCK
C                           PRICE BY USING A BINOMIAL APPROX. TO THE
C                           LOGNORMAL VOLATILITY DISTRIBUTION.
C                           (FULL LATTICE)
C       REQUIRED ROUTINES - FACT, BSCALL
C
C----------------------------------------------------------------------
C
      SUBROUTINE ASVCORP2(S,X,R,T,V,TDIV,DIV,PSI,XMU,COR,REPS,CVAL)
      REAL S,X,R,T,V0,TDIV,DIV,PSI,XMU,CVAL,V,U,D
      REAL RX,XVAL,COR,UX,DX,SIG,DT
      REAL SL(262500),VL(262500),EEXERVAL(262500)
      REAL TI,TITDIV,DIVADJ
      INTEGER REPS,J,I1,J1,I,IEXER,IR,JLAST,reps1
C
CCCC  NOTE: V AND VBAR ARE STD DEVIATIONS, V1 AND V2 ARE VARIANCES
C
      IEXER=0
C
CCCC  CHECK FOR REASONABLE PARAMETERS
C
      CVAL=0.0
      V0=V*V
      DT=T/FLOAT(REPS)
      U=EXP((XMU-0.5*(PSI**2.0))*(DT)+COR*PSI*SQRT(DT))
      D=EXP((XMU-0.5*(PSI**2.0))*(DT)-COR*PSI*SQRT(DT))
C
CCCC  SET RX (PERIODIC 1+RF) TO RISK FREE RATE + 1
C
      RX=(R+1.0)**(T/FLOAT(REPS))
C
CCCC  CREATE STOCK PRICE LATTICE (SL)
C

      VL(1)=V0
      SL(1)=S
c         reps1 = INT(reps*tdiv/t) + 1
c         if (reps1.gt.reps) reps1 = reps
      DO 3000 I=1,REPS
      TI=FLOAT(I)*DT
      TITDIV=TI-TDIV
      IF ((TI .GT. TDIV) .AND. (TDIV .LT. T) .AND. (TITDIV .LT. DT))
     1   THEN
c      if (i.eq.reps1) then
          DIVADJ=DIV
      ELSE
          DIVADJ=0.0
      ENDIF
      J=2**I
 1000 J1=J/2
      SIG=VL(J1)
      UX=EXP((R-0.5*(SIG))*(DT)+SQRT(SIG)*SQRT(DT))
      DX=EXP((R-0.5*(SIG))*(DT)-SQRT(SIG)*SQRT(DT))
    
      IF (DIVADJ .GT. 0.0) THEN
       XVAL=SL(J1)-X
       EEXERVAL(J1)=AMAX1(XVAL,0.0)
       IEXER=I
      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  
      J=J-2
      IF (J .GT. 0) GOTO 1000
 3000 CONTINUE
C
CCCC  CALCULATE CALL VALUES AT EXPIRATION

      I1=2**REPS
      DO 4000 I=1,I1
      XVAL=SL(I)-X
      SL(I)=AMAX1(XVAL,0.0)
 4000 CONTINUE

CCCC  CALCULATE CALL PRICE RECURSIVELYC
      IR=REPS
      DO WHILE (IR .GT. 0)
      J=2
      JLAST=2**IR
      DO WHILE (J .LE. JLAST)
        J1=J/2
        SL(J1)=(.5*(SL(J)+SL(J-1)))/RX
        IF (IR .EQ. IEXER) THEN
          SL(J1)=AMAX1(SL(J1),EEXERVAL(J1))
        ENDIF
        J=J+2
      END DO
      IR=IR-1
      END DO
      CVAL=SL(1)
      RETURN
      END


c       routine name       - svcorp2    
c ---------------------------------------------------------------------
c       purpose           - this subroutine calculates the value of a           
c                           european call with stochastic volatility            
c                           and a correlation of 1 or -1 with the stock         
c                           price by using a binomial approx. to the            
c                           lognormal volatility distribution.                  
c                           (full lattice)                                      
c       required routines - fact, bscall                                        
c                                                                               
c----------------------------------------------------------------------         
c                                                                               
      subroutine svcorp2(s,x,r,t,v,psi,xmu,cor,reps,cval)                       
      real s,x,r,t,v0,psi,xmu,cval,cx,xreps,xj,v,vx                 
      real u,d,p,r1,hedge,vbar,prob,ptot,fac                     
      real sx,rx,xx,dzero,px,probx,xval,cor                         
      real s1,s2,ux,dx,sig,dt,lp                                    
      real sl(524288),vl(524288)
      integer reps,j,n,i1,j1,i                                                  
c                                                                               
cccc  note: v and vbar are std deviations, v1 and v2 are variances              
c                                                                               
      dzero=0.0                                                                 
c                                                                               
cccc  check for reasonable parameters                                           
c                                                                               
      cval=0.0                                                                  
                                                                           
cccc  set r1 (periodic 1+rf) to xmu (volatility drift) + 1                      
c                                                                               
      r1=(xmu+1.0)**(t/float(reps))                                            
      v0=v*v                                                                    
      dt=t/float(reps)                                                         
      p=0.50                                                                    
      u=exp((xmu-0.5*(psi**2.0))*(dt)+cor*psi*sqrt(dt))                       
      d=exp((xmu-0.5*(psi**2.0))*(dt)-cor*psi*sqrt(dt))                       
c                                                                               
cccc  set rx (periodic 1+rf) to risk free rate + 1                              
c                                                                               
      rx=(r+1.0)**(t/float(reps))                                              
c                                                                               
cccc  create stock price lattice (sl)                                           
c                                                                               
      vl(1)=v0                                                                  
      sl(1)=s                                                                   
      do 3000 i=1,reps                                                          
c      j=2**i 
      do 3000 j = 2**i,0,-2                                                                   
 1000 j1=j/2                                                                    
      sig=vl(j1)                                                                
      ux=exp((r-0.5*(sig))*(dt)+sqrt(sig)*sqrt(dt))                        
      dx=exp((r-0.5*(sig))*(dt)-sqrt(sig)*sqrt(dt))
                        
      sl(j)=sl(j1)*dx                                                                                        
      sl(j-1)=sl(j1)*ux
                                                         
      vl(j)=vl(j1)*d                                                            
      vl(j-1)=vl(j1)*u                                                          
c     j=j-2                                                                     
c      if (j .gt. 0) goto 1000                                                   
 3000 continue                                                                  
c                                                                               
cccc  calculate call price                                                      
c                                                                               
      xreps=float(reps)                                                        
      i1=2**reps                                                                
      prob=1.0/float(i1)
        ii = 0
      do 4000 i=1,i1

      xval=sl(i)-x                                                              
      if(xval .le. dzero) goto 4000
                                             
      cval=cval+(xval)/(rx**xreps)
       ii = ii + 1                                  
 4000 continue
      cval = cval/float(i1)
                                                                  
      return                                                                    
      end    



