      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)
      real timediff,timearray(2)

      integer reps,errst,nt,srep,rep0,repn,repp,reps1       
      integer repp2,repn2,bmethd,smethd,j1,j2,j3,j4                                                                
                                                       
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=100
      repn=200
      repp=200                                                  
      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
       sz = 90.


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


  200 continue                                         

ccccc generate tables                  


      do 88888 sz=90., 110., 5.
      print *
      print *,'sz= ',sz    
     
c      print *,'odd # (5,7,9,11,13,15)'

      do 66666 repn2 = 17, 19
         repp2 = repn2
c         print *,'# = ',repn2
                    
      timediff = dtime(timearray)

      CALL SVCORP2(SZ,K,R,T,V,PSI,XMU,1.0,repn2,CLn)
      timediff = dtime(timearray)
      print *,'#=',repn2,' running time = ',timediff,' seconds'
c      print *,'#= ',repn2,'  cval(1)=',cln

66666  continue

       stop
c      CALL SVCORP2(SZ,K,R,T,V,PSI,XMU,1.0,REPP2,CLP)      
c      print *,'call value(+1)=',clp
     
        
c99999 continue
88888 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


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((xmu-0.5*(sig))*(dt)+sqrt(sig)*sqrt(dt))                        
      dx=exp((xmu-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)
c      print *,'prob = ',prob,'i1 = ',i1                                                       
        ii = 0
      do 4000 i=1,i1
c      print *,'i= ',i,'  sl(i)= ',sl(i)                                                            
      xval=sl(i)-x                                                              
      if(xval .le. dzero) goto 4000
                                             
      cval=cval+(xval)/(rx**xreps)
       ii = ii + 1                                  
 4000 continue
c       print *,'ii = ',ii
      cval = cval/float(i1)
                                                                  
      return                                                                    
      end                                                                       
