c*********************************************************************** c Program calculates prices of American and European calls / puts on c dividend paying stock. European value computed both by closed form c formula and by Crank-Nicholson method for accuracy comparison. c c Author: 29 Sept 2001 R. Jones c c*********************************************************************** implicit double precision (a-h,k-l,o-z) character*1 answer parameter (n = 600) dimension parm(15), ecall(0:n), eput(0:n), acall(0:n), aput(0:n), & xvcall(0:n), xvput(0:n), bscall(0:n), bsput(0:n), & arr(0:n,4) data smin, smax, kk / 0d0, 300d0, .02d0 / data nstart, nstop, nstep / 0, 240, 20 / h = (smax - smin) / n 100 write(*,9000) 9000 format('Enter riskless rate and strike price: '\) read (*,*) r, strike 102 write(*,9001) 9001 format('Enter maturity, volatility, div rate: '\) read (*,*) tmat, vol, div c Set exercise and maturity values do i = 1, n s = smin + i * h xvcall(i) = max( 0d0, s - strike ) xvput (i) = max( 0d0, strike - s ) acall (i) = xvcall(i) ecall (i) = xvcall(i) aput (i) = xvput (i) eput (i) = xvput (i) enddo c Calculate Black-Scholes formula prices ... do i = nstart, nstop, nstep s = smin + i * h call BS(0,s,strike,tmat,r,vol,div,bscall(i)) call BS(1,s,strike,tmat,r,vol,div,bsput (i)) enddo c Calculate Crank-Nicholson prices ... if ( tmat .le. 0d0 ) goto 200 k = tmat / max( 1, nint( tmat / kk ) ) ! bend timestep to fit imin = 1 ! forces value at s=0 imax = 0 parm(1) = vol parm(2) = r parm(3) = div call CNSET(n,smin,smax,k,ifn,0,imin,imax,parm,arr) do t = 0d0, tmat - k/2, k vmin = 0d0 call CNSTEP(vmin,acall,arr) call CNSTEP(vmin,ecall,arr) vmin = strike call CNSTEP(vmin,aput ,arr) vmin = strike * exp( -r * (t + k) ) call CNSTEP(vmin,eput ,arr) do i = 0, n acall(i) = max( acall(i), xvcall(i) ) aput (i) = max( aput (i), xvput (i) ) enddo enddo c Print to screen and file ... 200 write(*,9500) 9500 format(/4x,'S',6x,'BScall',4x,'Ecall',4x,'Acall',4x, & 'BSput ',4x,'Eput ',4x,'Aput ') 9502 format(f6.1 ,2x, 6f9.3) 9004 format(/'Do another? (y/n): '\) 9006 format(a1) do i = nstart, nstop, nstep s = smin + i * h write(*,9502) s, bscall(i), ecall(i), acall(i), & bsput (i), eput (i), aput (i) enddo write(*,9004) read (*,9006) answer if ( answer .eq. 'y' .or. answer .eq. 'Y' ) goto 100 stop end c*********************************************************************** C C The coefficients for the Black-Scholes option pricing model, with C S being the stock price, and R and D the assumed constant interest C rate and proportional dividend rate respectively, are C FNA(S) = SIGS * SIGS * S * S / 2.0 C FNB(S) = (R - div ) * S C FNC(S) = -R C C Coefficients for pde satisfied by pdf of S given below C*********************************************************************** double precision function COEFF() implicit double precision(a-h,k-l,o-z) dimension parm(15) entry FNA(s,ifn,parm) sigma = parm(1) r = parm(2) div = parm(3) FNA = (sigma * s) ** 2 / 2d0 return entry FNB(s,ifn,parm) FNB = (r - div) * s return entry FNC(s,ifn,parm) FNC = -r return entry FMIN(t,ifn,parm) ! value at smin passed as t FMIN = t return entry FMAX(t,ifn,parm) FMAX = t return end C**** Routine to calculate default-free BS call value ****************** subroutine BS (IFLAG, S, K, T, R, SIG, D, P) C*********************************************************************** C Subroutine calculates BLACK-SCHOLES value of european option C C INPUTS: S current stock price C K exercise price C T time to expirty C R interest rate C D proportional dividend rate on stock C SIG stock volatility C IFLAG 0 for CALL 1 for PUT C C OUTPUTS: P option value C C CALLS: subroutine NDTR() for cumulative normal distribution fcn. C C 9 Sept 2001 ... several fixes R. Jones C*********************************************************************** implicit double precision (A-H,K-L,O-Z) C Calculate option values KDISC = K * exp( -R * T ) SDISC = S * exp( -D * T ) if ( S .le. 0d0 ) then if ( IFLAG .EQ. 1 ) then P = KDISC else P = 0D0 endif return endif if ( T .le. 0D0 ) then if ( IFLAG .eq. 1 ) then P = max( 0D0, K - S ) else P = max( 0D0, S - K ) endif return endif SKDIS = max( 1D-20, SDISC/max( 1D-20, KDISC ) ) TROOT = sqrt( T ) SIG2 = SIG * SIG XPLUS = ( log( SKDIS ) + SIG2 * T / 2D0 ) / (SIG * TROOT) XMINUS= XPLUS - SIG * TROOT call NDTR( XPLUS , P , D1 ) call NDTR( XMINUS, P2, D2 ) C = SDISC * P - KDISC * P2 if ( IFLAG .eq. 1) then P = C - SDISC + K * exp( -R * T) else P = C endif return end subroutine NDTR(X,P,D) C Calculates cumulative normal distribution and density functions C X = value for which CDF is calculated C P = return value of normal CDF C D = return value of normal density function implicit double precision (A-H,O-Z) AX = dmin1(1.0D20, dabs(X)) T = 1.0D0/(1.0D0 + .2316419D0 * AX) D = 0.3989423D0 * dexp(-AX*AX/2.0D0) P = 1.0D0 - D*T*((((1.330274D0 * T - 1.821256D0) * T 1 + 1.781478D0) * T - 0.3565638D0) * T + 0.3193815D0) if (X) 1,2,2 1 P = 1.0D0 - P 2 return end