MathFinance - Financial Functions
Library in Fortran 90
Here you find some sample source code written in FORTRAN 90.
It covers the Cox-Ross-Rubinstein binomial model, the analytic
Black-Scholes formula as well as Monte Carlo Simulation. These
functions can easily we converted into a dynamic link library, such
that you can use them in other applications such as Microsoft Excel.
!-----------------------------------------------------------------------------------------
!wy_US produces values and greeks of European and American puts and calls
!following the Cox, Ross, Rubinstein binomial model.
!-----------------------------------------------------------------------------------------
recursive double precision function wy_US(spot, strike, vol, rd, rf, tau, phi, ea, &
convention, quotation, greek, N) result (f)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_US
parameter percent = 0.01
parameter nmax = 1002
double precision :: spot, strike, vol , rd , rf, tau
integer ::ea, convention, quotation, greek , N, phi
double precision, dimension(-nmax: nmax):: s
double precision, dimension(-1: nmax):: V
double precision :: dt, DF, u, p, q, relstrike, rforward, disc, valuefactor
integer :: i, j, k
select case (greek)
case (0) !value
f = 0. !initialize, then do two binomial trees, one with N steps and one with N+1 steps
do k = 0, 1
! set parameters:
dt = tau / (N + k)
u = exp(vol * Sqrt(dt))
select case (convention)
case (0) !linear for for mat less than one year, annually compounded for higher maturities
if (tau .le. 1.) then
rforward = (1. + rd * tau * 365. / 360.) / (1. + rf * tau * 365. / 360.)
disc = 1. / (1. + rd * tau * 365. / 360.)
else
rforward = exp((log(1. + rd) - log(1. + rf)) * tau)
disc = exp(-log(1 + rd) * tau)
end if
case (1) !annually compounded
rforward = exp((log(1. + rd) - log(1. + rf)) * tau)
disc = exp(-log(1. + rd) * tau)
case (2) !continuously compounded
rforward = exp((rd - rf) * tau)
disc = exp(-rd * tau)
end select
DF = exp(log(disc) / (N+k))
p = (exp(log(rforward) / (N+k)) - 1. / u) / (u - 1. / u) * DF
q = DF - p
relstrike = strike / spot
! build the stock price tree, assuming spot = 1:
do j = -(N + k), N + k
s(j) = u**j
end do
! set final time value:
do j = 0, N + k
V(j) = wy_payoff((s(2 * j - N - k)), relstrike, phi)
end do
! compute present value by stepping backwards in the value tree:
do i = N + k, 1,-1
do j = 0, i - 1
V(j) = (p * V(j + 1) + q * V(j))
if (ea .eq. 2) then
V(j) = max(V(j), wy_payoff((s(2 * j - i + 1)), relstrike, phi))
end if
end do
end do
f = f + V(0)
end do
f = f * spot / 2.
case (1) !delta
f = (wy_US(spot*(1.+percent) , strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- wy_US(spot*(1.-percent) , strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (2. * percent * spot)
valuefactor = 1 / spot
case (2) !gamma
f = (wy_US(spot*(1.+percent) , strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
+ wy_US(spot*(1.-percent) , strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- 2. * wy_US(spot, strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (percent*percent*spot*spot)
case (3) !theta
f = (wy_US(spot, strike, vol, rd, rf, tau*(1.-percent) , phi, ea, convention ,0 , 0, N) &
- wy_US(spot, strike, vol, rd, rf, tau*(1.+percent) , phi, ea, convention ,0 , 0, N)) / (2. * percent * tau)
case (4) !vega
f = (wy_US(spot, strike, vol*(1.+percent) , rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- wy_US(spot, strike, vol*(1.-percent) , rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (200. * percent * vol)
case (5) !rhod
f = (wy_US(spot, strike, vol, rd*(1.+percent) , rf, tau, phi, ea, convention ,0 , 0, N) &
- wy_US(spot, strike, vol, rd*(1.-percent) , rf, tau, phi, ea, convention ,0 , 0, N)) / (2. * percent * rd)
case (6) !rhof
f = (wy_US(spot, strike, vol, rd, rf*(1.+percent) , tau, phi, ea, convention ,0 , 0, N) &
- wy_US(spot, strike, vol, rd, rf*(1.-percent) , tau, phi, ea, convention ,0 , 0, N)) / (2. * percent * rf)
case (7) !dstrike
f = (wy_US(spot, strike*(1.+percent) , vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- wy_US(spot, strike*(1.-percent) , vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (2. * percent * strike)
case (8) !d2strike
f = (wy_US(spot, strike*(1.+percent) , vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
+ wy_US(spot, strike*(1.-percent) , vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- 2. * wy_US(spot, strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (percent*percent*strike*strike)
case (9) !dduration
f = (wy_US(spot, strike, vol, rd, rf, tau*(1.+percent) , phi, ea, convention ,0 , 0, N) &
- wy_US(spot, strike, vol, rd, rf, tau*(1.-percent) , phi, ea, convention ,0 , 0, N)) / (2. * percent * tau)
case (10)!leverage
f = spot * (wy_US(spot*(1.+percent) , strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- wy_US(spot*(1.-percent) , strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (2. * percent * spot) &
/ wy_US(spot, strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)
case (11)!vomma
f = (wy_US(spot, strike, vol*(1.+3.*percent) , rd, rf, tau, phi, ea, convention ,0 , 0, N) &
+ wy_US(spot, strike, vol*(1.-3.*percent) , rd, rf, tau, phi, ea, convention ,0 , 0, N) &
- 2. * wy_US(spot, strike, vol, rd, rf, tau, phi, ea, convention ,0 , 0, N)) / (9.*percent*percent*vol*vol)
end select
!quotation adjustment:
select case (greek)
case (0, 3, 4, 5, 6, 9)
select case (quotation)
case (0) !domestic
case (1) !%domestic
f = f / strike * 100.
case (2) !foreign
f = f / strike / spot
case (3) !%foreign
f = f / spot * 100.
end select
case (1, 7)
select case (quotation)
case (0) !plain derivative
case (1) !plain derivative in %
f = f * 100.
case (2) !spot * plain derivative of %foreign quoted value
f = f - wy_US(spot, strike, vol, rd, rf, tau, phi, ea, convention, 0, 0, N) * valuefactor
case (3) !spot * plain derivative of %foreign quoted value in%
f = 100. * (f - wy_US(spot, strike, vol, rd, rf, tau, phi, ea, convention, 0, 0, N) * valuefactor)
end select
case (2)
select case (quotation)
case (0) !plain derivative
case (1) !change of delta when spot changes by 1%
f = f * spot
end select
case (8)
select case (quotation)
case (0) !plain derivative
case (1) !change of delta when strike changes by 1%
f = f * strike
end select
end select
end function
!---------------------------------------------------------------------------------------------------------------
double precision function wy_payoff(spot, strike, phi)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_PAYOFF
double precision spot, strike
integer :: phi
wy_payoff = max(0.D00, phi * (spot - strike))
end function
!---------------------------------------------------------------------------------------------------------------
!******************************************************************************************************************
!*** vanilla function follows
!*** analytic valuation and greeks for European style puts and calls
!*** allows different daycount and interest rate conventions
!*** author: Uwe Wystup
!*** Date : March 1999
!*** reference: Aspects of Symmetry, Homogeneity and Duality of the Black-Scholes pricing formula for European
!*** put and call options, by Uwe Wystup
!***
!*** List of Greeks:
!*** ---------------
!*** 0 :value
!*** 1 :spotdelta
!*** 2 :gamma
!*** 3 :annual theta
!*** 4 :vega
!*** 5 :rho domestic
!*** 6 :rho foreign
!*** 7 :dstrike !derivative wrt strike
!*** 8 :d2strike !second derivative wrt strike
!*** 9 :dduration !derivative wrt to T(time to expiration) = - annual theta
!*** 10:leverage !spot * spotdelta / value
!*** 11:vomma !second derivative wrt vol
!*** 12:retrieved_vola !volatility implied by the option value, enter value for vol
!*** 13:forwarddelta !derivative wrt forward
!*** 14:probdelta !common probability factor in spot delta and forward delta, N(d1)
!*** 15:daily theta
!*** 16:forward
!*** 17:swaprate
!*** 18:d1
!*** 19:d2
!*** 20:disc !domestic discount factor
!*** 21:strike_of_delta
!******************************************************************************************************************
recursive double precision Function wy_vanilla(spot , strike , vol , rd , rf , tau , phi , convention, &
quotation , greek) result(f)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_VANILLA
double precision :: spot, strike, vol, rd, rf, tau
integer :: phi, convention, quotation, greek
double precision :: d1, d2, forward, disc, convrd, convrf, dconvrd, dconvrf, valuefactor
Select Case (convention)
Case (0) !linear for for mat less than one year, annually compounded for higher maturities
If (tau .lt. 1.) Then
forward = spot * (1. + rd * tau * 365. / 360.) / (1. + rf * tau * 365. / 360.)
disc = 1. / (1. + rd * tau * 365. / 360.)
convrd = (rd * 365. / 360.) / (1. + rd * tau * 365. / 360.)
convrf = (rf * 365. / 360.) / (1. + rf * tau * 365. / 360.)
dconvrd = (365. / 360.) / (1. + rd * tau * 365. / 360.)
dconvrf = (365. / 360.) / (1. + rf * tau * 365. / 360.)
Else
forward = spot * Exp((Log(1. + rd) - Log(1. + rf)) * tau)
disc = Exp(-Log(1. + rd) * tau)
convrd = Log(1. + rd)
convrf = Log(1. + rf)
dconvrd = 1. / (1. + rd)
dconvrf = 1. / (1. + rf)
End If
Case (1) !annually compounded
forward = spot * Exp((Log(1. + rd) - Log(1. + rf)) * tau)
disc = Exp(-Log(1. + rd) * tau)
convrd = Log(1. + rd)
convrf = Log(1. + rf)
dconvrd = 1. / (1. + rd)
dconvrf = 1. / (1. + rf)
Case (2) !continuously compounded
forward = spot * Exp((rd - rf) * tau)
disc = Exp(-rd * tau)
convrd = rd
convrf = rf
dconvrd = 1.
dconvrf = 1.
End Select
If (greek .ne. 21) Then
d1 = (Log(forward / strike) + vol * vol * 0.5 * tau) / (vol * Sqrt(tau))
d2 = d1 - vol * Sqrt(tau)
End If
Select Case (greek)
Case (0) !value
f = phi * disc * (forward * wy_nc(phi * d1) - strike * wy_nc(phi * d2))
Case (1) !spotdelta
f = phi * disc * forward / spot * wy_nc(phi * d1)
valuefactor = 1. / spot
Case (2) !gamma
f = disc * forward / spot * wy_ndf(d1) / (vol * Sqrt(tau) * spot)
Select Case (quotation)
Case (0) !plain derivative
Case (1) !change in delta when spot changes by 1% in % (trader's Gamma)
f = f * spot
End Select
Case (3) !annual theta
f = -disc * (wy_ndf(d1) * vol * forward / (2. * Sqrt(tau)) &
+ phi * (convrd * strike * wy_nc(phi * d2) - convrf * forward * wy_nc(phi * d1)))
Case (4) !vega
f = disc * forward * wy_ndf(d1) * Sqrt(tau) / 100.
Case (5) !domestic rho
f = phi * strike * tau * disc * wy_nc(phi * d2) * dconvrd
Case (6) !foreign rho
f = -phi * forward * tau * disc * wy_nc(phi * d1) * dconvrf
Case (7) !dstrike
f = -phi * disc * wy_nc(phi * d2)
valuefactor = 0.
Case (8) !d2strike
f = disc * wy_ndf(d2) / (strike * vol * Sqrt(tau))
Case (9) !dduration
f = disc * (wy_ndf(d1) * vol * forward / (2. * Sqrt(tau)) &
+ phi * (convrd * strike * wy_nc(phi * d2) - convrf * forward * wy_nc(phi * d1)))
Case (10) !leverage
f = forward * wy_nc(phi * d1) / (forward * wy_nc(phi * d1) - strike * wy_nc(phi * d2))
Case (11) !vomma
f = disc * forward * wy_ndf(d1) * Sqrt(tau) * d1 * d2 / vol
Case (12) !retrieved_vola
f = wy_find_vol(spot, strike, vol, rd, rf, tau, phi, convention, quotation, 0.1D0)
Select Case (quotation)
Case (0) !plain volatility
Case (1) !volatility in %
f = f * 100.
End Select
Case (13) !forwarddelta
f = phi * wy_nc(phi * d1) * disc
valuefactor = 1. / forward
Case (14) !probdelta
f = phi * wy_nc(phi * d1)
valuefactor = 1. / (forward * disc)
Case (15) !daily theta
f = (-disc * (wy_ndf(d1) * vol * forward / (2. * Sqrt(tau)) &
+ phi * (convrd * strike * wy_nc(phi * d2) - convrf * forward * wy_nc(phi * d1)))) / 365.
Case (16) !forward
f = forward
Case (17) !swaprate
f = (forward - spot) * 10000.
Case (18) !d1
f = d1
Case (19) !d2
f = d2
Case (20) !domestic discount factor
f = disc
Case (21) !strike_of_delta (spot delta is given, strike is wanted)
Select Case (quotation)
Case (0) !plain derivative
f = forward * Exp(-phi * wy_ncinv(phi * strike / disc / forward * spot) &
* vol * Sqrt(tau) + 0.5D0 * vol * vol * tau)
Case (1) !plain derivative in %
f = forward * Exp(-phi * wy_ncinv(phi * strike / 100.D0/ disc / forward * spot) &
* vol * Sqrt(tau) + 0.5D0 * vol * vol * tau)
Case (2) !spot * plain derivative of %foreign quoted value
f = wy_strike_of_pfdelta(spot, strike, vol, disc, forward, tau, phi)
Case (3) !spot * plain derivative of %foreign quoted value in%
f = wy_strike_of_pfdelta(spot, strike / 100.D0, vol, disc, forward, tau, phi)
End Select
End Select
Select Case (greek)
Case (0, 3, 4, 5, 6, 9, 15)
Select Case (quotation)
Case (0) !domestic
Case (1) !%domestic
f = f / strike * 100.
Case (2) !foreign
f = f / strike / spot
Case (3) !%foreign
f = f / spot * 100.
End Select
Case (1, 7, 13, 14)
Select Case (quotation)
Case (0) !plain derivative
Case (1) !plain derivative in %
f = f * 100.
Case (2) !spot * plain derivative of %foreign quoted value
f = f - phi * disc * (forward * wy_nc(phi * d1) - strike * wy_nc(phi * d2)) * valuefactor
Case (3) !spot * plain derivative of %foreign quoted value in%
f = 100. * (f - phi * disc * (forward * wy_nc(phi * d1) - strike * wy_nc(phi * d2)) * valuefactor)
End Select
End Select
End Function
!------------------------------------------------------------------------------------------------------------------
! find_vol finds the volatility given a value of a vanilla option based on Newton's method
!------------------------------------------------------------------------------------------------------------------
recursive double precision Function wy_find_vol(spot, strike, vol, rd, rf, tau, phi, convention, quotation, &
startvol) result(f)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_FIND_VOL
double precision :: spot, strike, vol, rd, rf, tau, startvol
integer :: phi, convention, quotation
double precision :: func, dfunc
integer :: counter
Select Case (quotation) !invert the quotation:
Case (0) !domestic
Case (1) !%domestic
vol = vol * strike / 100.D0
Case (2) !foreign
vol = vol * strike * spot
Case (3) !%foreign
vol = vol * spot / 100.D0
End Select
!start Newton's method:
f = startvol
func = wy_vanilla(spot, strike, startvol, rd, rf, tau, phi, convention, 0, 0) - vol
dfunc = wy_vanilla(spot, strike, startvol, rd, rf, tau, phi, convention, 0, 4) * 100.D0
counter = 0
Do while (Abs(func / dfunc) .ge. 1.0D-7 .and. counter .le. 20)
f = f - func / dfunc
func = wy_vanilla(spot, strike, f, rd, rf, tau, phi, convention, 0, 0) - vol
dfunc = wy_vanilla(spot, strike, f, rd, rf, tau, phi, convention, 0, 4) * 100.D0
counter = counter + 1
end do
End Function
!---------------------------------------------------------------------------------------------------------------
! nc returns the cumulative distribution function of a standard normal random variable, see Hull
!---------------------------------------------------------------------------------------------------------------
recursive double precision function wy_nc(x) result (f)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_NC
double precision :: x
double precision, dimension(1: 5):: a
If (x .lt. -7.D0) Then
f = wy_ndf(x) / sqrt(1.D0 + x * x)
ElseIf (x .gt. 7.D0) Then
f = 1.D0 - wy_nc(-x)
Else
f = 0.2316419D0
a(1) = 0.31938153D0
a(2) = -0.356563782D0
a(3) = 1.781477937D0
a(4) = -1.821255978D0
a(5) = 1.330274429D0
f = 1.D0 / (1.D0 + f * Abs(x))
f = 1.D0 - wy_ndf(x) * (a(1) * f + a(2) * f ** 2.D0 + a(3) * f ** 3.D0 + a(4) * f ** 4.D0 + a(5) * f ** 5.D0)
If (x .le. 0.D0) Then
f = 1.D0 - f
end if
End If
End Function
!---------------------------------------------------------------------------------------------------------------
! ndf returns the density function of a standard normal random variable
double precision function wy_ndf(x)
double precision :: x
!DEC$ ATTRIBUTES DLLEXPORT :: WY_NDF
wy_ndf = 0.398942280401433D0 * exp(-x * x * 0.5D0) !0.398942280401433 = 1/sqareroot(2*pi)
end function
!---------------------------------------------------------------------------------------------------------------
!---------------------------------------------------------------------------------------------------------------
! ncinv returns the inverse of a standard normal distribution function (nc in this library)
!---------------------------------------------------------------------------------------------------------------
double precision function wy_ncinv(x)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_NCINV
double precision :: x
double precision :: func, dfunc
if (x .le. 0.D0) then
wy_ncinv = -1.5D1
elseif (x .ge. 1.D0) then
wy_ncinv = 1.5D1
else
!start Newton's method:
wy_ncinv = 0.D0
func = 0.5D0 - x
dfunc = 0.398942280401433D0
Do while (Abs(func / dfunc) .ge. 1.D-7)
wy_ncinv = wy_ncinv - func / dfunc
func = wy_nc(wy_ncinv) - x
dfunc = wy_ndf(wy_ncinv)
end do
endif
end function
!---------------------------------------------------------------------------------------------------------------
!CARDANO produces the roots of the equation x³+ax²+bx+c=0.
!output: 0 : discriminant D which classifies the solution:
! D>0: one real y(1) and two complex conjugate solutions
! y(2)+y(3)i and y(2)-y(3)i
! D=0: three real solutions including one double solution (y(3)=0)
! D<0: three distinct real solutions y(1)<y(2)<y(3)
!output: 1,2,3 : y(output)
!reference: Harris/Stocker, Handbook of Mathematics and Computational Science,
! Springer 1998, Chapter 2.5 Cubic Equations
!-----------------------------------------------------------------------------
double precision function wy_cardano(a,b,c,output)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_CARDANO
parameter epsilon = 0.000001
double precision a,b,c
double precision p, q, D, phi, u, v
double precision x(3)
double precision y(3)
double precision radicand
integer output
p=(3.*b-a*a)/3.
q=c+2.*a*a*a/27.-a*b/3.
D=(P/3.)**3.+q*q/4.
if (D .LT. 0.) then !the irreducible case has a trigonometric form:
phi=acos(-q/(2*sqrt((abs(p)/3)**3)))
x(1)=-a/3.+2*sqrt(abs(p)/3)*cos(phi/3.)
x(2)=-a/3.-2*sqrt(abs(p)/3)*cos((phi-3.141592654)/3.)
x(3)=-a/3.-2*sqrt(abs(p)/3)*cos((phi+3.141592654)/3.)
! sort solutions according to order:
y(1) = min(x(1),x(2),x(3))
y(3) = max(x(1),x(2),x(3))
y(2) = x(1)+x(2)+x(3)-y(1)-y(3)
else
radicand =abs(-q/2.+sqrt(D))
if (radicand .ge. epsilon) then
u=sign(1.,-q/2.+sqrt(D))*exp(log(radicand)/3.)
else
u=0
end if
radicand =abs(-q/2.-sqrt(D))
if (radicand .ge. epsilon) then
v=sign(1.,-q/2.-sqrt(D))*exp(log(radicand)/3.)
else
v=0
end if
y(1)=-a/3.+u+v
y(2)=-a/3.-(u+v)/2.
y(3)=sqrt(3.)*(u-v)/2.
end if
select case (output)
case (0)
wy_cardano = D
case (1:3)
wy_cardano = y(output)
end select
end function
!------------------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------------------
double precision function wy_normrv
!DEC$ ATTRIBUTES DLLEXPORT :: WY_NORMRV
double precision :: x
!call random_seed()
call random_number(x)
wy_normrv = wy_ncinv(x)
end function
!------------------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------
! Monte Carlo Simulator for the geometric Brownian motion model
! illustrates the techniques 1) antithetic variables and 2) control variate
! Prices Vanilla and fixed strike Asian options
!-------------------------------------------------------------------------------------------
recursive double precision function wy_mc(spot, strike, vol, rd, rf, tau, phi, convention, &
quotation, payofftype, para1, para2, Nsim, Nsteps, greek) result(f)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_MC
integer, parameter :: Nstepsmax = 1000
double precision, parameter :: days = 365.0D0
double precision :: spot, strike, vol, rd, rf, tau, para1, para2
integer :: phi, convention, quotation, payofftype, Nsim, Nsteps, greek, i, j
double precision :: forward, disc, convrd, convrf, dconvrd, dconvrf, fc1, fc2, beta
double precision :: my, volrt, delta_t, sfinal1, sfinal2, zeta, f1, f2, geoavg1, geoavg2, &
ariavg1, ariavg2
double precision, dimension(0: Nstepsmax):: s1, s2
beta = 1.0D0
s1(0) = Log(spot)
s2(0) = Log(spot)
delta_t = tau / days / Nsteps
my = (rd - rf - 0.5D0 * vol * vol) * delta_t
volrt = vol * Sqrt(delta_t)
f1 = 0.D0
f2 = 0.D0
fc1 = 0.D0
fc2 = 0.D0
call random_seed() !initialize random number generator
do j = 0, Nsim
do i = 1, Nsteps
call random_number(zeta)
zeta = wy_ncinv(zeta) !convert uniform zeta into standard normal zeta
! simulate path and antithetic path according to Euler scheme
s1(i) = s1(i - 1) + my + volrt * zeta
s2(i) = s2(i - 1) + my - volrt * zeta
end do
select case (payofftype)
case (0) !European put/call
f1 = f1 + wy_payoff(Exp(s1(Nsteps)), strike, phi)
f2 = f2 + wy_payoff(Exp(s2(Nsteps)), strike, phi)
case (1) !Goemetric fixed strike Asian put/call
geoavg1 = 0.0D0
geoavg2 = 0.0D0
do i = 0, Nsteps
geoavg1 = geoavg1 + s1(i)
geoavg2 = geoavg2 + s2(i)
end do
select case (greek)
case (0) !value
f1 = f1 + wy_payoff(Exp(geoavg1/(Nsteps+1)), strike, phi)
f2 = f2 + wy_payoff(Exp(geoavg2/(Nsteps+1)), strike, phi)
case (1) !delta
if (phi*geoavg1/(Nsteps+1) .gt. phi*log(strike)) then
f1 = f1 + Exp(geoavg1/(Nsteps+1))/spot
endif
if (phi*geoavg2/(Nsteps+1) .gt. phi*log(strike)) then
f2 = f2 + Exp(geoavg2/(Nsteps+1))/spot
endif
end select
case (2) !Arithmetic fixed strike Asian put/call
geoavg1 = 0.0D0
geoavg2 = 0.0D0
ariavg1 = 0.0D0
ariavg2 = 0.0D0
do i = 0, Nsteps
geoavg1 = geoavg1 + s1(i)
geoavg2 = geoavg2 + s2(i)
ariavg1 = ariavg1 + exp(s1(i))
ariavg2 = ariavg2 + exp(s2(i))
end do
fc1 = fc1 + wy_payoff(Exp(geoavg1/(Nsteps+1)), strike, phi)
fc2 = fc2 + wy_payoff(Exp(geoavg2/(Nsteps+1)), strike, phi)
f1 = f1 + wy_payoff(ariavg1/(Nsteps+1), strike, phi)
f2 = f2 + wy_payoff(ariavg2/(Nsteps+1), strike, phi)
end select
end do
! get value of option with antithetic technique
f = 0.5D0 * (f1 + f2) / Nsim * Exp(-rd * tau / days)
if (payofftype .eq. 2) then
f = f + beta * (wy_Asiangeo(spot, strike, vol, rd, rf, tau / days, 0.0D0, spot, phi, &
convention, quotation, 0) - 0.5D0 * (fc1 + fc2) / Nsim * Exp(-rd * tau / days))
end if
End function
!------------------------------------------------------------------------------------------------------------
! Asiangeo computes value and greek of fixed strike European puts and calls on the continuously sampled
! geometric average taken between -s and T (-backtau and tau), check Wilmott:Derivatives
!------------------------------------------------------------------------------------------------------------
recursive double precision Function wy_Asiangeo(spot, strike, vol, rd, rf, tau, backtau, currentavg, phi, &
convention, quotation, greek) result(f)
!DEC$ ATTRIBUTES DLLEXPORT :: WY_ASIANGEO
double precision :: spot, strike, vol, rd, rf, tau, backtau, currentavg
integer :: phi, convention, quotation, greek
double precision :: d1, d2, forward, disc, convrd, convrf, dconvrd, dconvrf, valuefactor, geofactor, alpha
alpha = tau / (tau + backtau)
Select Case (convention)
Case (0) !linear for for mat less than one year, annually compounded for higher maturities
If (tau .le. 1.D0) Then
forward = spot * (1.0D0 + rd * tau * 365.0D0 / 360.0D0) / (1.0D0 + rf * tau * 365.0D0 / 360.0D0)
disc = 1.0D0 / (1.0D0 + rd * tau * 365 / 360)
convrd = (rd * 365.0D0 / 360.0D0) / (1.0D0 + rd * tau * 365.0D0 / 360.0D0)
convrf = (rf * 365.0D0 / 360.0D0) / (1.0D0 + rf * tau * 365.0D0 / 360.0D0)
dconvrd = (365.0D0 / 360.0D0) / (1.0D0 + rd * tau * 365.0D0 / 360.0D0)
dconvrf = (365.0D0 / 360.0D0) / (1.0D0 + rf * tau * 365.0D0 / 360.0D0)
geofactor = Exp(-0.5D0 * tau * alpha * (Log(1.0D0 + rd) - Log(1.0D0 + rf) + vol * vol / 2.0D0 * (1.0D0 - 2.0D0 * alpha / 3.0D0)))
Else
forward = spot * Exp((Log(1.0D0 + rd) - Log(1.0D0 + rf)) * tau)
disc = Exp(-Log(1.0D0 + rd) * tau)
convrd = Log(1.0D0 + rd)
convrf = Log(1.0D0 + rf)
dconvrd = 1.0D0 / (1.0D0 + rd)
dconvrf = 1.0D0 / (1.0D0 + rf)
geofactor = Exp(-0.5D0 * tau * alpha * (Log(1.0D0 + rd) - Log(1.0D0 + rf) + vol * vol / 2.0D0 * (1.0D0 - 2.0D0 * alpha / 3.0D0)))
End If
Case (1) !annually compounded
forward = spot * Exp((Log(1.0D0 + rd) - Log(1.0D0 + rf)) * tau)
disc = Exp(-Log(1.0D0 + rd) * tau)
convrd = Log(1.0D0 + rd)
convrf = Log(1.0D0 + rf)
dconvrd = 1.0D0 / (1.0D0 + rd)
dconvrf = 1.0D0 / (1.0D0 + rf)
geofactor = Exp(-0.5D0 * tau * alpha * (Log(1.0D0 + rd) - Log(1.0D0 + rf) + vol * vol / 2.0D0 * (1.0D0 - 2.0D0 * alpha / 3.0D0)))
Case (2) !continuously compounded
forward = spot * Exp((rd - rf) * tau)
disc = Exp(-rd * tau)
convrd = rd
convrf = rf
dconvrd = 1.0D0
dconvrf = 1.0D0
geofactor = Exp(-0.5D0 * tau * alpha * (rd - rf + vol * vol / 2.0D0 * (1.0D0 - 2.0D0 * alpha / 3.0D0)))
End Select
Select Case (greek)
Case (0) !value
f = disc ** (1.0D0 - alpha) * wy_vanilla(geofactor * spot ** alpha * currentavg ** (1.0D0 - alpha), strike, &
alpha * vol / Sqrt(3.0D0), alpha * rd, alpha * rf, tau, phi, &
convention, 0, 0)
Case (1) !spotdelta
f= disc ** (1.0D0 - alpha) * wy_vanilla(geofactor * spot ** alpha * currentavg ** (1.0D0 - alpha), strike, &
alpha * vol / Sqrt(3.0D0), alpha * rd, alpha * rf, tau, phi, &
convention, 0, 1) * geofactor * alpha * (currentavg / spot) ** (1.0D0 - alpha)
End Select
End Function
!----------------------------------------------------------------------------------------------------------
