c ##################################################################
c ##################################################################
c ###### ######
c ###### Copyright (c) 1993 ######
c ###### Northeast Parallel Architectures Center ######
c ###### Syracuse University. All rights reserved. ######
c ###### ######
c ##################################################################
c ##################################################################
c#######################################################################
c
c File: models.f
c
c AUTHORS: originally written by Thomas J. Finucane;
c rewritten and documented by
c Gang Cheng (gcheng@npac.syr.edu, (315)443-2083 )
c 06/07/91.
c
c MODIFICATION HISTORY:
c
c 06/06/92 (G. Cheng)
c Added documentation.
C
c#######################################################################
c *************** Option Price Models (Fortran 77 version) **************
c
c This Fortran 77 set contains 4 option pricing models with the
c subroutine names:
c
c bscall ---- Black-Scholes Model
c binomc ---- Constant Volatility American Call Binomial Model
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 Black-Scholes Model
c---------------------------------------------------------------------
c Routine name - bscall
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 (black-scholes model)
c---------------------------------------------------------------------
c Arguments:
c sp --- initial stock price sk --- exercise (striking) price
c r --- risk-free rate t --- time to maturity
c v --- initial volatility cval --- call price (output)
c---------------------------------------------------------------------
c Required routines - pnorml
c---------------------------------------------------------------------
subroutine bscall(sp,sk,r,t,v,cval)
real sp,sk,r,t,v,cval
implicit none
real hedge,h1,h2,p,tmp1,tmp2,pnorml
c --- calculate black-scholes call value
tmp1 = sk/exp(r*t)
tmp2 = v*sqrt(t)
c --- h = (log(sp/(sk/(exp(r*t)))))/(v*sqrt(t))+.5*(v*sqrt(t))
h1 = LOG(sp/tmp1)/tmp2 + 0.5*tmp2
h2 = h1 - tmp2
hedge = pnorml(h1)
p = pnorml(h2)
cval = sp*hedge - tmp1*p
return
end
c---------------------------------------------------------------------
c Routine name - pnorml (used by Black-Scholes Model)
c---------------------------------------------------------------------
c Purpose - This function supplies a polynomial
c approximation for the standard normal
c distribution.
c---------------------------------------------------------------------
c Required routines - none
c---------------------------------------------------------------------
function pnorml(z)
real z
real tmp1,tmp2,tmp3,tmp4
tmp1 = z/sqrt(2.)
tmp2 = abs(tmp1)
tmp3 = 1./(1. + 0.5*tmp2)
tmp4 = tmp3*exp(-tmp2*tmp2 - 1.26551223 + tmp3*(1.00002368 +
& tmp3*(.37409196 + tmp3*(.09678418 + tmp3*(-.18628806 +
& tmp3*(.27886807 + tmp3*(-1.13520398 + tmp3*(1.48851587 +
& tmp3*(-.82215223 + tmp3*.17087277)))))))))
if (tmp1.lt.0.) tmp4 = 2. - tmp4
pnorml = 1. - 0.5 * tmp4
return
end
c------------------------------------------------------------------------------
c Constant Volatility American Call Binomial Model
c------------------------------------------------------------------------------
c Routine name - binomc
c------------------------------------------------------------------------------
c Purpose - This subroutine calculates the value of an American
c call with 1 dividend using the binomial numerical
c procedure of Cox and Ross, and escrowed dividends.
c------------------------------------------------------------------------------
c Arguments:
c sp --- initial 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 reps --- number of discreted time steps
c cval --- call price (output)
c------------------------------------------------------------------------------
c Required routines - None
c------------------------------------------------------------------------------
subroutine binomc(sp,sk,r,t,v,tdiv,div,reps,cval)
real sp,sk,r,t,v,tdiv,div,cval
integer reps
parameter (MAXSIZE = 1000, P = 0.5)
implicit none
real r1,u,d,a1,a2,s0,escrow,xx,dt,tmp1,tmp2,pvdiv
real c(MAXSIZE)
integer i,j,n1,n2
c --- set formula parameters
dt = t/float(reps)
r1 = (r + 1.0)**dt
tmp1 = (r - 0.5*v*v)*dt
tmp2 = v*SQRT(dt)
u = tmp1 + tmp2
d = tmp1 - tmp2
n2 = INT(tdiv*365.0)
pvdiv = div/(r1**(tdiv*float(reps)/t))
s0 = LOG(sp - pvdiv)
c --- calculate binomial call value
do 950 i = 0, reps
c --- calculate reduction in stock price due to dividends
n1 = INT(FLOAT(reps-i)*t*365.0/FLOAT(reps))
if (n1 .ge. n2) then
escrow = 0.0
else
escrow = pvdiv*(r1**(reps-i))
endif
c --- evaluate binomial call formula recursively
do 900 j = 0, reps - i
xx = FLOAT(j)*u + FLOAT(reps-i-j)*d + s0
xx = AMIN1(xx,20.0)
a1 = EXP(xx) + escrow - sk
if (i .eq. 0) then
a2 = 0.0
else
a2 = (P*c(j+2) + (1.0-P) * c(j+1))/r1
endif
c(j+1) = AMAX1(a1,a2)
900 continue
950 continue
cval = c(1)
return
end
c---------------------------------------------------------------------------
c Stochastic Volatility American Call Binomial Model
c---------------------------------------------------------------------------
c Routine name - asvcorp2
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 - None
c---------------------------------------------------------------------------
subroutine asvcorp2(sp,sk,r,t,v,tdiv,div,psi,xmu,cor,reps,cval)
real sp,sk,r,t,v,tdiv,div,psi,xmu,cor,cval
integer reps
parameter (MAXSIZE = 262500)
implicit none
real sl(MAXSIZE),vl(MAXSIZE),eexerval(MAXSIZE)
real u,d,rx,xval,ux,dx,sig,dt,divadj,tmp1,tmp2
integer i,j,j1,iexer,reps1
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 --- create stock price lattice(sl) and volatility lattice(vl)
vl(1) = v*v
sl(1) = sp
reps1 = int(reps*tdiv/t) + 1
if (reps1.gt.reps) reps1 = reps
do 3000 i = 1, reps
if (i.eq.reps1) then
divadj = div
iexer = i
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 --- calculate call values at expiration
do 4000 i = 1, 2**reps
xval = sl(i) - sk
sl(i) = AMAX1(xval,0.0)
4000 continue
c --- set rx (periodic 1+rf) to risk free rate + 1
rx = (r+1.0)**(t/float(reps))
c --- calculate call price recursively
do 6000 i = reps, 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
c----------------------------------------------------------------------------
c Stochastic Volatility European Call Binomial Model
c----------------------------------------------------------------------------
c Routine name - svcorp2
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 - None
c----------------------------------------------------------------------------
subroutine svcorp2(sp,sk,r,t,v,psi,xmu,cor,reps,cval)
real sp,sk,r,t,v,psi,xmu,cor,cval
integer reps
parameter (MAXSIZE = 524288)
implicit none
real sl(MAXSIZE),vl(MAXSIZE)
real u,d,rx,xval,ux,dx,sig,dt,tmp1,tmp2
integer i,j,j1
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 --- create stock price lattice (sl)
vl(1) = v * v
sl(1) = sp
do 3000 i = 1, reps
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))
c --- calculate call price
tmp1 = rx**FLOAT(reps)
do 4000 i = 1, 2**reps
xval = sl(i) - sk
if (xval .gt. 0.0) cval = cval + xval/tmp1
4000 continue
cval = cval/float(2**reps)
return
end
Northeast Parallel Architectures Center, Syracuse University, npac@npac.syr.edu
This page is maintained by Gang Cheng, gcheng@npac.syr.edu