MathFinance - Financial Functions Library in Fortran 90

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
!----------------------------------------------------------------------------------------------------------
 
 






























MathFinance logo Footer Pic 1 Footer Pic 2 Footer Pic 3 Footer Pic 4 © MathFinance AG
Privacy Policy  |  Disclaimer  |  AGB / Terms and Conditions