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 AUTHOR: Gang Cheng (gcheng@npac.syr.edu, (315)443-2083 )
c 06/07/92.
c
c MODIFICATION HISTORY:
c
c 08/06/92 (G. Cheng)
c Added documentation.
C
c################################################################
c *************** Option Price Models (Fortran 90 version) **************
c
c This Fortran 90 set contains 2 option pricing models with the
c subroutine names:
c
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 total discreted time steps (reps = 17, reps >= 13)
c number of discreted time steps at dividend time
c ( reps1 = MAX((INT(reps*tdiv1/t)+1),reps) )
c
c Note: cmf$ layout and CMPF MAP are used as MAPPING directives in
c CMfortran and MPfortran respectively to get best performance
c on 8K CM2 and 8K DECmpp-12000
c ************************************************************************
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 --- total number of discreted time steps
c reps1 --- number of discreted time steps at dividend time
c ( reps1 = MAX((INT(reps*tdiv1/t)+1),reps) )
c cval --- call price (output)
c----------------------------------------------------------------------------
c Required routines - None
c---------------------------------------------------------------------------
SUBROUTINE asvcorp2(sp,sk,r,t,v,psi,xmu,cor,reps1,reps,div,cval)
IMPLICIT NONE
real sp,sk,r,t,v,psi,xmu,cor,div,cval
integer reps,reps1
real,array(1:(2**(reps-12)-1),1:8192) :: ss1,vv1
real,array(1:2**(reps-13),1:8192) :: cearly_big
real,array(1:8192) :: aatmp1,aatmp2,uux,ddx
real,array(1:2**reps1) :: cearly,cback
real u,d,dt,r1,tmp1,tmp2,rx,xreps
integer i,j,n,nreps1,start,end,repsdiv,reps2,nreps2
cmf$ layout ss1(:serial,:news),vv1(:serial,:news)
cmf$ layout cearly_big(:serial,:news)
cmf$ layout uux(:news),ddx(:news),cearly(:news),cback(:news)
cmf$ layout aatmp1(:news),aatmp2(:news)
CMPF MAP ss1(MEMORY,ALLBITS),vv1(MEMORY,ALLBITS)
CMPF MAP cearly_big(MEMORY,ALLBITS)
CMPF MAP uux(ALLBITS),ddx(ALLBITS)
CMPF MAP aatmp1(ALLBITS),aatmp2(ALLBITS)
CMPF MAP cearly(ALLBITS),cback(ALLBITS)
xreps = float(reps)
dt = t/xreps
c --- SET RX (PERIODIC 1+RF) TO RISK FREE RATE + 1
rx=(r+1.0)**dt
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)
vv1(1,:) = 0.0
ss1(1,:) = 0.0
vv1(1,1) = v*v
ss1(1,1) = sp
n = 1
do i = 1,13
aatmp1 = (r-0.5*vv1(1,:))*dt
aatmp2 = sqrt(vv1(1,:)*dt)
uux = exp(aatmp1+aatmp2)
ddx = exp(aatmp1-aatmp2)
vv1(1,:) = u*vv1(1,:) +
& eoshift(d*vv1(1,:),dim = 1,shift = -n)
ss1(1,:) = uux*ss1(1,:) +
& eoshift(ddx*ss1(1,:),dim = 1,shift = -n)
if (i.eq.reps1) then
cearly = ss1(1,1:2**i) - sk
ss1(1,1:2**i) = ss1(1,1:2**i) - div
endif
n = 2*n
enddo
if (reps.gt.13) then
reps2 = reps - 13
nreps1 = 2**(reps2+1) - 1
repsdiv = reps1 - 14
start = 2**repsdiv
end = 2*start - 1
do i = 1, nreps1/2
aatmp1 = (r-0.5*vv1(i,:))*dt
aatmp2 = sqrt(vv1(i,:)*dt)
uux = exp(aatmp1+aatmp2)
ddx = exp(aatmp1-aatmp2)
vv1(2*i,:) = u*vv1(i,:)
vv1((2*i+1),:) = d*vv1(i,:)
if (i.ge.start.and.i.le.end) then
aatmp1 = uux*ss1(i,:)
aatmp2 = ddx*ss1(i,:)
cearly_big(2*(i-start)+1,:) = aatmp1 - sk
cearly_big(2*(i-start+1),:) = aatmp2 - sk
ss1(2*i,:) = aatmp1 - div
ss1((2*i+1),:) = aatmp2 - div
else
ss1(2*i,:) = uux*ss1(i,:)
ss1((2*i+1),:) = ddx*ss1(i,:)
endif
ss1(i,:) = 0.0
enddo
ss1 = ss1 - sk
where (ss1.lt.0.0) ss1 = 0.0
if (reps1.le.13) then
do i=nreps1/2+1,nreps1
ss1(1,:) = ss1(1,:) + ss1(i,:)
enddo
ss1(1,:) = ss1(1,:)/((2*rx)**(float(reps-13)))
n = 1
do i=1,13-reps1
ss1(1,:) = ss1(1,:) +
& eoshift(ss1(1,:),dim = 1,shift = n)
n = 2*n
enddo
nreps2 = 2**(13-reps1)
cback = ss1(1,1:8192:nreps2)/
& ((2*rx)**float(13-reps1))
where (cearly.gt.cback) cback = cearly
c --- calculate call price
cval = SUM(cback,mask=cback.gt.0.0)/((2*rx)**(float(reps1)))
else
nreps2 = 2**(reps-reps1)
cval = 0.0
do i=1,2*start
aatmp1 = 0.0
n = nreps1/2 + (i-1)*nreps2
do j=1,nreps2
aatmp1 = aatmp1 +
& ss1(n+j,:)
enddo
aatmp1 = aatmp1/((2*rx)**(float(reps-reps1)))
where (cearly_big(i,:).gt.aatmp1)
& aatmp1 = cearly_big(i,:)
cval = cval + SUM(aatmp1,mask=aatmp1.gt.0.0)
enddo
c --- calculate call price
cval = cval/((2*rx)**(float(reps1)))
endif
else
c --- calculate call price
ss1 = ss1 - sk
cval = SUM(ss1(1,:),mask=ss1(1,:).gt.0.0)/((2*rx)**xreps)
endif
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)
implicit none
real sp,sk,r,t,v,psi,xmu,cor,cval
integer reps
real,array(1:(2**(reps-12)-1),1:8192) :: ss1,vv1,uux,ddx,
& aatmp1,aatmp2
integer,array(1:8192) :: index,ud
real u,d,dt,r1,tmp1,tmp2,rx,xreps
integer i,reps1,nreps1,n2
cmf$ layout ss1(:serial,:news),vv1(:serial,:news)
cmf$ layout uux(:serial,:news),ddx(:serial,:news)
cmf$ layout aatmp1(:serial,:news),aatmp2(:serial,:news)
cmf$ layout index(:news),ud(:news)
CMPF MAP ss1(ALLBITS,MEMORY),vv1(ALLBITS,MEMORY)
CMPF MAP uux(ALLBITS,MEMORY),ddx(ALLBITS,MEMORY)
CMPF MAP aatmp1(ALLBITS,MEMORY),aatmp2(ALLBITS,MEMORY)
CMPF MAP index(ALLBITS),ud(ALLBITS)
xreps = float(reps)
dt = t/xreps
index = [0:8191]
c --- SET RX (PERIODIC 1+RF) TO RISK FREE RATE + 1
rx=(r+1.0)**dt
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)
vv1(1,:) = v*v
ss1(1,:) = sp
n2 = 4096
do i = 1,13
aatmp1(1,:) = (r-0.5*vv1(1,:))*dt
aatmp2(1,:) = sqrt(vv1(1,:)*dt)
uux(1,:) = exp(aatmp1(1,:)+aatmp2(1,:))
ddx(1,:) = exp(aatmp1(1,:)-aatmp2(1,:))
ud = index/n2
c where (mod(ud,2) .ne. 0)
where (BTEST(ud,POS=0))
vv1(1,:) = d*vv1(1,:)
ss1(1,:) = ddx(1,:)*ss1(1,:)
elsewhere
vv1(1,:) = u*vv1(1,:)
ss1(1,:) = uux(1,:)*ss1(1,:)
endwhere
n2 = n2/2
enddo
if (reps.gt.13) then
reps1 = reps - 13
nreps1 = 2**(reps1+1) - 1
do i = 1, nreps1/2
aatmp1(i,:) = (r-0.5*vv1(i,:))*dt
aatmp2(i,:) = sqrt(vv1(i,:)*dt)
uux(i,:) = exp(aatmp1(i,:)+aatmp2(i,:))
ddx(i,:) = exp(aatmp1(i,:)-aatmp2(i,:))
vv1(2*i,:) = u*vv1(i,:)
vv1((2*i+1),:) = d*vv1(i,:)
ss1(2*i,:) = uux(i,:)*ss1(i,:)
ss1((2*i+1),:) = ddx(i,:)*ss1(i,:)
ss1(i,:) = 0.0
enddo
c --- calculate call price
ss1 = ss1 - sk
cval = SUM(ss1,mask=ss1.gt.0.0)/((2*rx)**xreps)
else
c --- calculate call price
ss1 = ss1 - sk
cval = SUM(ss1(1,:),mask=ss1(1,:).gt.0.0)/((2*rx)**xreps)
endif
return
end
Northeast Parallel Architectures Center, Syracuse University, npac@npac.syr.edu
This page is maintained by Gang Cheng, gcheng@npac.syr.edu