      program main
   
      real a1,a2,xmu,rho              
      real s,k,r,t,v,vic(3),hedge                          
      real cl(4),clp,cl0,cln,biasp,bias0,biasn  
      real sz,kz,tz,cz                           
      real er(4),clb,bias,biasx                  
      real emax(4),emin(4),ser(4),varm(75,4),mean(4)  
      real sd(4),en,sz2,k2                               
      real ksh,ksl,vh,vl,cdvh(3),cdvl(3),cv0(3),kl,kh
      real clba,va,ksa,ka,ks               
      real pdvh(3),pdvl(3)                    
      real clbl,clbh,psi,none,pone,cor
      real early,earlyn,earlyp,sig,sign,sigp

      real vv(200),hh(200),ccl(200)
      
      integer reps,errst,nt,srep,rep0,repn,repp,reps1,repb       
      integer repp2,repn2,bmethd,smethd,j1,j2,j3,j4               
            real timediff,timearray(2)
                                                 
                                                       
c     bmethd=0 - forced recombination binomial (repp,repn)    
c     bmethd=1 - full binomial                 (repp2,pepn2)
c     smethd=0 - control variate simulation                     
c     smethd=1 - no control variates in simulation          
                                                         
c   set initial values                          


      bmethd=1
      smethd=0                                                   
      s=100.00                                                     
      sz=s-10.00                                                          
      sz2=s+10.00                                   
      sz = 110
      rep0=200
      repn=200
      repp=200
      repb=100                                                  
      repp2=17
      repn2=17                                                      
      srep=1024
      nt=90                                                   
      psi=1.0
      xmu=0.0                                                    

      rho=-1.0                                                     
      if (sz .lt. 106.0) rho = 0.6
      rho = -1.0
      none=-1.0                                                
      pone=1.0                                                        
      iflg=0                                                   
      k=100.00                                                   
      r=0.0                                                             
      t=.25                                                      
      v=0.15            
      div = 1.0
      tdiv = 0.125

cccc repeat for a range of values of stock price (sz)


  200 continue                                         

ccccc generate tables                  

      do 99999 sz=90., 110., 5.

ccc    print headers                      

      write(*,9050)                                
      write(*,9100)                                                         
      write(*,9120)                                                         
      write(*,9150) sz                                                      
      write(*,9160) k                                                       
      write(*,9180) r                                                       
      write(*,9200) t                                                       
      write(*,9220) v                                                       
      write(*,9225) psi                                                     
      write(*,9228) xmu                                                     
      write(*,9233) tdiv
      write(*,9277) div
      write(*,9229) rep0,repp,repn                 
      write(*,9235) repp2,repn2  
      write(*,9230) srep,nt    
      write(*,9320)             
      write(*,9340)    
      write(*,9360)     

          timediff = dtime(timearray)

      CALL SVCOR0(SZ,K,R,T,V,TDIV,DIV,PSI,XMU,REP0,REPB,CL0)

      print *,'call value(0)=',cl0  

          timediff = dtime(timearray)
          write(6,40) timediff

40        format(//,'time took:',f6.3,'seconds'/)

99999 continue

 9050 format(1h /)
 9100 format(1h ,'American Call Option Pricing Errors by Correlation')
 9120 format(1h ,61(1h-))
 9140 format(1h ,'base values  : ')
 9150 format(1h ,19x,'Stock price       :  ',f7.3)
 9160 format(1h ,19x,'Striking price    :  ',f7.3)
 9180 format(1h ,19x,'Interest rate     :  ',f7.3)
 9200 format(1h ,19x,'Time to maturity  :  ',f7.3)
 9220 format(1h ,19x,'Volatility        :  ',f7.3)
 9225 format(1h ,19x,'Psi               :  ',f7.3)
 9227 format(1h ,19x,'Rho               :  ',f7.3)
 9228 format(1h ,19x,'Mu                :  ',f7.3)
 9229 format(1h ,19x,'rep0,repp,repn    :  ',3(i5,','))
 9230 format(1h ,19x,'Srep, nt          :  ',i5,',',i5)
 9235 format(1h ,19x,'repp2,repn2       :  ',2(i5,','))
 9233 format(1h ,19x,'time to dividend  :  ',f7.3)   
 9277 format(1h ,19x,'dividend amount   :  ',f7.3)   
 9234 format(1h ,19x,'simulation method :  ',i3)   
 9280 format(1h ,19x,'option value      :  ',f9.5/)
 9300 format(1h ///1h ,'         ',3x,'stoch vol',6x,' b-s  ',
     1       8x,'stoch vol')
 9320 format(1h ///1h ,'  rho',4x,'simulated',5x,'F-B',5x,
     1       'standard',4x,'% early')
 9340 format(1h ,11x,'price',6x,'price',6x,'error',5x,
     1 'exercise')
 9360 format(1h ,7(1h-),2x,4(9(1h-),2x)/)
 9380 format(1h ,1x,f5.2,2x,f9.5,2x,f9.5,2x,f8.6,2x,f9.2)
 9400 FORMAT(1H ,'Mean Error :',7X,F9.5,2X)   
 9420 FORMAT(1H ,'Std. Deviation :',3X,F9.5,2X)
 9440 FORMAT(1H ,'Minimum Error :',4X,F9.5,2X)
 9460 FORMAT(1H ,'Maximum Error :',4X,F9.5,2X)
 9480 FORMAT(1H ,19X,2(9(1H-),2X)/)


      stop
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
C
C
C
C       ROUTINE NAME       - SVCOR0
C
C----------------------------------------------------------------------
C
C       PURPOSE           - THIS SUBROUTINE CALCULATES THE VALUEOF A
C                           EUROPEAN CALL WITH STOCHASTIC VOLATILITY
C                           AND A CORRELATION OF 0 WITH THE STOCK
C                           PRICE BY USING A BINOMIAL APPROX. TO THE
C                           LOGNORMAL VOLATILITY DISTRIBUTION.
C
C       REQUIRED ROUTINES - FACT, BSCALL
C
C----------------------------------------------------------------------
C
      SUBROUTINE SVCOR0(S,X,R,T,V,TDIV,DIV,PSI,XMU,REPS,REPB,CVAL)

      REAL S,X,R,T,V0,TDIV,DIV,PSI,XMU,CVAL,CX,XREPS,XJ,V
      REAL U,D,VBAR,V1,PROB,PTOT,DT,LP
      INTEGER REPS,J,REPB

      CVAL=0.0
CCCC  SET U,D,P
C
      V0=V*V
      DT=T/FLOAT(REPS)
      U=EXP((XMU-0.5*(PSI**2.0))*(DT)+PSI*SQRT(DT))
      D=EXP((XMU-0.5*(PSI**2.0))*(DT)-PSI*SQRT(DT))
C
CCCC  CALCULATE TERMINAL VOLATILITIES AND PROBABILITIES
C
      J=0
      PTOT=0.0
 1000 CONTINUE
      XREPS=FLOAT(REPS)
      XJ=FLOAT(J)
      CALL LPROB(XREPS,XJ,LP)
      PROB=0.0
      IF (LP .LT. -174.0) GOTO 2000
      PROB=EXP(LP)
 2000 CONTINUE
      PTOT=PTOT+PROB
      V1=V0*(U**FLOAT(J))*(D**FLOAT(REPS-J))
C
CCCC  CALCULATE ARITHMETIC AVERAGE VOLATILITY
C
      VBAR=SQRT((V1+V0)/2.0)     

C
CCCC  CALCULATE CALL VALUE USING AVERAGE VOLATILITY VBAR
C
C     IF (VBAR .GT. BIGV) THEN VBAR=BIGV
      CALL BINOMC(S,X,R,T,VBAR,TDIV,DIV,REPB,CX)
C
CCCC  WEIGHT BINOMIAL CALL VALUE BY PROB. AND ADD TO TOTAL
C
      CVAL=CVAL+(CX*PROB)     

      J=J+1
      IF (J .LE. REPS) GOTO 1000
 
      RETURN
      END
C
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



 American Call Option Pricing Errors by Correlation
 -------------------------------------------------------------
                    Stock price       :   90.000
                    Striking price    :  100.000
                    Interest rate     :    0.000
                    Time to maturity  :    0.250
                    Volatility        :    0.150
                    Psi               :    1.000
                    Mu                :    0.000
                    time to dividend  :    0.125
                    dividend amount   :    1.000
                    rep0,repp,repn    :    200,  200,  200,
                    repp2,repn2       :     17,   17,
                    Srep, nt          :   1024,   90
 


   rho    simulated     F-B     standard    % early
            price      price      error     exercise
 -------  ---------  ---------  ---------  ---------

call value(0)=   0.191497


time took:50.240seconds

 

 American Call Option Pricing Errors by Correlation
 -------------------------------------------------------------
                    Stock price       :   95.000
                    Striking price    :  100.000
                    Interest rate     :    0.000
                    Time to maturity  :    0.250
                    Volatility        :    0.150
                    Psi               :    1.000
                    Mu                :    0.000
                    time to dividend  :    0.125
                    dividend amount   :    1.000
                    rep0,repp,repn    :    200,  200,  200,
                    repp2,repn2       :     17,   17,
                    Srep, nt          :   1024,   90
 


   rho    simulated     F-B     standard    % early
            price      price      error     exercise
 -------  ---------  ---------  ---------  ---------

call value(0)=   0.846107


time took:49.920seconds

 

 American Call Option Pricing Errors by Correlation
 -------------------------------------------------------------
                    Stock price       :  100.000
                    Striking price    :  100.000
                    Interest rate     :    0.000
                    Time to maturity  :    0.250
                    Volatility        :    0.150
                    Psi               :    1.000
                    Mu                :    0.000
                    time to dividend  :    0.125
                    dividend amount   :    1.000
                    rep0,repp,repn    :    200,  200,  200,
                    repp2,repn2       :     17,   17,
                    Srep, nt          :   1024,   90
 


   rho    simulated     F-B     standard    % early
            price      price      error     exercise
 -------  ---------  ---------  ---------  ---------

call value(0)=    2.59615


time took:50.200seconds

 

 American Call Option Pricing Errors by Correlation
 -------------------------------------------------------------
                    Stock price       :  105.000
                    Striking price    :  100.000
                    Interest rate     :    0.000
                    Time to maturity  :    0.250
                    Volatility        :    0.150
                    Psi               :    1.000
                    Mu                :    0.000
                    time to dividend  :    0.125
                    dividend amount   :    1.000
                    rep0,repp,repn    :    200,  200,  200,
                    repp2,repn2       :     17,   17,
                    Srep, nt          :   1024,   90
 


   rho    simulated     F-B     standard    % early
            price      price      error     exercise
 -------  ---------  ---------  ---------  ---------

call value(0)=    5.82778


time took:50.090seconds

 

 American Call Option Pricing Errors by Correlation
 -------------------------------------------------------------
                    Stock price       :  110.000
                    Striking price    :  100.000
                    Interest rate     :    0.000
                    Time to maturity  :    0.250
                    Volatility        :    0.150
                    Psi               :    1.000
                    Mu                :    0.000
                    time to dividend  :    0.125
                    dividend amount   :    1.000
                    rep0,repp,repn    :    200,  200,  200,
                    repp2,repn2       :     17,   17,
                    Srep, nt          :   1024,   90
 


   rho    simulated     F-B     standard    % early
            price      price      error     exercise
 -------  ---------  ---------  ---------  ---------

call value(0)=   10.18582


time took:51.280seconds



