!********************************************************************
!********************************************************************
! 几个简单方法实现欧洲电话的毕苏期权定价方程模式
!********************************************************************
! BLACK_SCHOLES Simple Approaches to the Black-Scholes Equation,
! which demonstrates several approaches to
! the valuation of a European call.
!
! SINOMACH
!
! PotsyYZ, LilyZ, ShuaiZ, YibingW, JHA, XingtongC, XiaolanG，QihaoZ，
! 周勇，WeichuH，Dr.Mr.Guo,JinhongL,JiananF,ShaolanZ，WenzheL，GuoxinH，
! GuofZ,MingsenL，XiminL，XingH，DetaoW，XinyongL，Pengrino，YongjianL,
! HuxionC, YideZ,Angle
! Sept. 2022
!*********************************************************************
!c $$ Run, Output:
!*********************************************************************
!*********************************************************************
program main
!*********************************************************************72
! MAIN is the main program for BLACK_SCHOLES_PRB.
!
! Discussion:
!
! BLACK_SCHOLES_PRB tests the BLACK_SCHOLES library.
!
implicit none
call timestamp ( );
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BLACK_SCHOLES_PRB:'
write ( *, '(a)' ) ' FORTRAN95,2013&18 version'
write ( *, '(a)' ) ' Test the BLACK_SCHOLES library.'
call asset_path_test ( )
call binomial_test ( )
call bsf_test ( )
call forward_test ( )
call mc_test ( )
! Terminate.
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BLACK_SCHOLES_PRB:'
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ' '
call timestamp ( )
return
end
subroutine asset_path_test ( )
!*********************************************************************72
! ASSET_PATH_TEST tests ASSET_PATH.
implicit none
integer n
parameter ( n = 100 )
double precision mu
character * ( 100 ) output_filename
double precision s(0:n)
double precision s0
integer seed
double precision sigma
double precision t1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'ASSET_PATH_TEST:'
write ( *, '(a)' ) ' Demonstrate the simulated of an asset price path.'
s0 = 2.0D+00
mu = 0.1D+00
sigma = 0.3D+00
t1 = 1.0D+00
seed = 123456789
write ( *, '(a,g14.6)' ) ' '
write ( *, '(a,g14.6)' ) ' The asset price at time 0 S0 = ', s0
write ( *, '(a,g14.6)' ) ' The asset expected growth rate MU = ', mu
write ( *, '(a,g14.6)' ) ' The asset volatility SIGMA = ', sigma
write ( *, '(a,g14.6)' ) ' The expiry date T1 = ', t1
write ( *, '(a,i6)' ) ' The number of time steps N = ', n
write ( *, '(a,i12)' ) ' The random number seed was SEED = ', seed
call asset_path ( s0, mu, sigma, t1, n, seed, s )
call r8vec_print_part ( n + 1, s, 10, ' Partial results:' )
output_filename = 'asset_path.txt'
call r8vec_write ( output_filename, n + 1, s )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Full results written to "' // trim ( output_filename ) // '".'
return
end
subroutine binomial_test ( )
!*********************************************************************72
! BINOMIAL_TEST tests BINOMIAL.
implicit none
double precision c
double precision e
integer m
double precision r
double precision s0
double precision sigma
double precision t1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'BINOMIAL_TEST:'
write ( *, '(a)' ) ' A demonstration of the binomial method'
write ( *, '(a)' ) ' for option valuation.'
s0 = 2.0D+00
e = 1.0D+00
r = 0.05D+00
sigma = 0.25D+00
t1 = 3.0D+00
m = 256
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' The asset price at time 0 S0 = ', s0
write ( *, '(a,g14.6)' ) ' The exercise price E = ', e
write ( *, '(a,g14.6)' ) ' The interest rate R = ', r
write ( *, '(a,g14.6)' ) ' The asset volatility SIGMA = ', sigma
write ( *, '(a,g14.6)' ) ' The expiry date T1 = ', t1
write ( *, '(a,i8)' ) ' The number of intervals M = ', m
call binomial ( s0, e, r, sigma, t1, m, c )
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' The option value is ', c
return
end
subroutine bsf_test ( )
!*********************************************************************72
! BSF_TEST tests BSF.
!
implicit none
double precision c
double precision e
double precision r
double precision s0
double precision sigma
double precision t0
double precision t1
write ( *, '(a)' )
write ( *, '(a)' ) 'BSF_TEST:'
write ( *, '(a)' ) ' A demonstration of the Black-Scholes formula'
write ( *, '(a)' ) ' for option valuation.'
s0 = 2.0D+00
t0 = 0.0D+00
e = 1.0D+00
r = 0.05D+00
sigma = 0.25D+00
t1 = 3.0D+00
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' The asset price at time T0 S0 = ', s0
write ( *, '(a,g14.6)' ) ' The time T0 = ', t0
write ( *, '(a,g14.6)' ) ' The exercise price E = ', e
write ( *, '(a,g14.6)' ) ' The interest rate R = ', r
write ( *, '(a,g14.6)' ) ' The asset volatility SIGMA = ', sigma
write ( *, '(a,g14.6)' ) ' The expiry date T1 = ', t1
call bsf ( s0, t0, e, r, sigma, t1, c )
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' The option value C = ', c
return
end
subroutine forward_test ( )
!*********************************************************************72
! FORWARD_TEST tests FORWARD.
implicit none
integer nt
parameter ( nt = 29 )
integer nx
parameter ( nx = 11 )
double precision e
integer i
double precision r
double precision s
double precision sigma
double precision smax
double precision smin
double precision t1
double precision u(nx-1,nt+1)
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'FORWARD_TEST:'
write ( *, '(a)' ) ' A demonstration of the forward difference method'
write ( *, '(a)' ) ' for option valuation.'
e = 4.0D+00
r = 0.03D+00
sigma = 0.50D+00
t1 = 1.0D+00
smax = 10.0D+00
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6)' ) ' The exercise price E = ', e
write ( *, '(a,g14.6)' ) ' The interest rate R = ', r
write ( *, '(a,g14.6)' ) ' The asset volatility SIGMA = ', sigma;
write ( *, '(a,g14.6)' ) ' The expiry date T1 = ', t1
write ( *, '(a,i8)' ) ' The number of space steps NX = ', nx
write ( *, '(a,i8)' ) ' The number of time steps NT = ', nt
write ( *, '(a,g14.6)' ) ' The value of SMAX = ', smax
call forward ( e, r, sigma, t1, nx, nt, smax, u )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Initial Option'
write ( *, '(a)' ) ' Value Value'
write ( *, '(a)' ) ' '
smin = 0.0D+00
do i = 1, nx - 1
s = ( ( nx - i - 1 ) * smin + i * smax ) / dble ( nx - 1 )
write ( *, '(2x,g14.6,2x,g14.6)' ) s, u(i,nt+1)
end do
return
end
subroutine mc_test ( )
!*********************************************************************72
! MC_TEST tests MC.
implicit none
double precision conf(2)
double precision e
integer m
double precision r
double precision s0
integer seed
double precision sigma
double precision t1
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'MC_TEST:'
write ( *, '(a)' ) ' A demonstration of the Monte Carlo method'
write ( *, '(a)' ) ' for option valuation.'
s0 = 2.0D+00
e = 1.0D+00
r = 0.05D+00
sigma = 0.25D+00
t1 = 3.0D+00
m = 1000000
seed = 123456789
write ( *, '(a)' ) ' '
write ( *, '(a, g14.6)' ) ' The asset price at time 0, S0 = ', s0
write ( *, '(a, g14.6)' ) ' The exercise price E = ', e
write ( *, '(a, g14.6)' ) ' The interest rate R = ', r
write ( *, '(a, g14.6)' ) ' The asset volatility SIGMA = ', sigma
write ( *, '(a, g14.6)' ) ' The expiry date T1 = ', t1
write ( *, '(a, i8)' ) ' The number of simulations M = ', m
call mc ( s0, e, r, sigma, t1, m, seed, conf )
write ( *, '(a)' ) ' '
write ( *, '(a,g14.6,a,g14.6,a)' ) &
' The confidence interval is [', conf(1), ',', conf(2), '].'
return
end
subroutine asset_path ( s0, mu, sigma, t1, n, seed, s )
!*********************************************************************72
!
! ASSET_PATH simulates the behavior of an asset price over time.
!
! Reference:
!
! Desmond Higham,
! Black-Scholes for Scientific Computing Students,
! Computing in Science and Engineering,
! November/December 2004, Volume 6, Number 6, pages 72-79.
!
! Parameters:
!
! Input, double precision S0, the asset price at time 0.
!
! Input, double precision MU, the expected growth rate.
!
! Input, double precision SIGMA, the volatility of the asset.
!
! Input, double precision T1, the expiry date.
!
! Input, integer N, the number of steps to take between 0 and T1.
!
! Input/output, integer SEED, a seed for the random number generator.
!
! Output, double precision S(0:N), the option values from time 0 to T1
! in equal steps.
!
implicit none
integer n
double precision dt
integer i
double precision mu
double precision p
double precision r(n)
double precision s(0:n)
double precision s0
integer seed
double precision sigma
double precision t1
dt = t1 / dble ( n )
call r8vec_normal_01 ( n, seed, r )
s(0) = s0
p = s0
do i = 1, n
p = p * exp ( ( mu - sigma * sigma ) * dt + sigma * sqrt ( dt ) * r(i) )
s(i) = p
end do
return
end
subroutine binomial ( s0, e, r, sigma, t1, m, c )
!*********************************************************************72
!
! BINOMIAL uses the binomial method for a European call.
!
! Reference:
!
! Desmond Higham,
! Black-Scholes for Scientific Computing Students,
! Computing in Science and Engineering,
! November/December 2004, Volume 6, Number 6, pages 72-79.
!
! Parameters:
!
! Input, double precision S0, the asset price at time 0.
!
! Input, double precision E, the exercise price.
!
! Input, double precision R, the interest rate.
!
! Input, double precision SIGMA, the volatility of the asset.
!
! Input, double precision T1, the expiry date.
!
! Input, integer M, the number of steps to take
! between 0 and T1.
!
! Output, double precision C, the option value at time 0.
implicit none
double precision a
double precision b
double precision c
double precision d
double precision dt
double precision e
integer i
integer m
integer n
double precision p
double precision r
double precision s0
double precision sigma
double precision t1
double precision u
double precision w(1:m+1)
! Time stepsize.
dt = t1 / dble ( m )
a = 0.5D+00 * ( exp ( - r * dt ) + exp ( ( r + sigma**2 ) * dt ) )
d = a - sqrt ( a * a - 1.0D+00 )
u = a + sqrt ( a * a - 1.0D+00 )
p = ( exp ( r * dt ) - d ) / ( u - d )
do i = 1, m + 1
w(i) = max ( s0 * d**(m+1-i) * u**(i-1) - e, 0.0D+00 )
end do
! Trace backwards to get the option value at time 0.
do n = m, 1, -1
do i = 1, n
w(i) = ( 1.0D+00 - p ) * w(i) + p * w(i+1)
end do
end do
do i = 1, m + 1
w(i) = exp ( - r * t1 ) * w(i)
end do
c = w(1)
return
end
subroutine bsf ( s0, t0, e, r, sigma, t1, c )
!*********************************************************************72
! BSF evaluates the Black-Scholes formula for a European call.
! Reference:
! Desmond Higham,
! Black-Scholes for Scientific Computing Students,
! Computing in Science and Engineering,
! November/December 2004, Volume 6, Number 6, pages 72-79.
! Parameters:
! Input, double precision S0, the asset price at time T0.
! Input, double precision T0, the time at which the asset price is known.
! Input, double precision E, the exercise price.
! Input, double precision R, the interest rate.
! Input, double precision SIGMA, the volatility of the asset.
! Input, double precision T1, the expiry date.
! Output, double precision C, the value of the call option.
implicit none
double precision c
double precision d1
double precision d2
double precision e
double precision n1
double precision n2
double precision r
double precision s0
double precision sigma
double precision t0
double precision t1
double precision tau
tau = t1 - t0
if ( 0.0D+00 .lt. tau ) then
d1 = ( log ( s0 / e ) + ( r + 0.5D+00 * sigma * sigma ) * tau )&
/ ( sigma * sqrt ( tau ) )
d2 = d1 - sigma * sqrt ( tau )
n1 = 0.5D+00 * ( 1.0D+00 + erf ( d1 / sqrt ( 2.0D+00 ) ) )
n2 = 0.5D+00 * ( 1.0D+00 + erf ( d2 / sqrt ( 2.0D+00 ) ) )
c = s0 * n1 - e * exp ( - r * tau ) * n2
else
c = max ( s0 - e, 0.0D+00 )
end if
return
end
subroutine forward ( e, r, sigma, t1, nx, nt, smax, u )
!*********************************************************************72
!
! FORWARD uses the forward difference method to value a European call option.
! Reference:
! Desmond Higham,
! Black-Scholes for Scientific Computing Students,
! Computing in Science and Engineering,
! November/December 2004, Volume 6, Number 6, pages 72-79.
! Parameters:
! Input, double precision E, the exercise price.
! Input, double precision R, the interest rate.
! Input, double precision SIGMA, the volatility of the asset.
! Input, double precision T1, the expiry date.
! Input, integer NX, the number of "space" steps used to
! divide the interval [0,L].
! Input, integer NT, the number of time steps.
! Input, double precision SMAX, the maximum value of S to consider.
! Output, double precision U(NX-1,NT+1), the value of the European
! call option.
implicit none
integer nt
integer nx
double precision a(2:nx-1)
double precision b(1:nx-1)
double precision c(1:nx-2)
double precision dt
double precision dx
double precision e
integer i
integer j
double precision p
double precision r
double precision sigma
double precision smax
double precision t
double precision t1
double precision u(nx-1,nt+1)
double precision u0
dt = t1 / dble ( nt )
dx = smax / dble ( nx )
do i = 1, nx - 1
b(i) = 1.0D+00 - r * dt - dt * ( sigma * i )**2
end do
do i = 1, nx - 2
c(i) = 0.5D+00 * dt * ( sigma * i )**2 + 0.5D+00 * dt * r * i
end do
do i = 2, nx - 1
a(i) = 0.5D+00 * dt * ( sigma * i )**2 - 0.5D+00 * dt * r * i
end do
u0 = 0.0D+00
do i = 1, nx - 1
u0 = u0 + dx
u(i,1) = max ( u0 - e, 0.0D+00 )
end do
do j = 1, nt
t = dble ( j - 1 ) * t1 / dble ( nt )
p = 0.5D+00 * dt * ( nx - 1 ) * ( sigma * sigma * ( nx - 1 ) + r ) &
* ( smax - e * exp ( - r * t ) )
do i = 1, nx - 1
u(i,j+1) = b(i) * u(i,j)
end do
do i = 1, nx - 2
u(i,j+1) = u(i,j+1) + c(i) * u(i+1,j)
end do
do i = 2, nx - 1
u(i,j+1) = u(i,j+1) + a(i) * u(i-1,j)
end do
u(nx-1,j+1) = u(nx-1,j+1) + p
end do
return
end
subroutine get_unit ( iunit )
!*********************************************************************72
! GET_UNIT returns a free FORTRAN unit number.
! Discussion:
! A "free" FORTRAN unit number is a value between 1 and 99 which
! is not currently associated with an I/O device. A free FORTRAN unit
! number is needed in order to open a file with the OPEN command.
! If IUNIT = 0, then no free FORTRAN unit could be found, although
! all 99 units were checked (except for units 5, 6 and 9, which
! are commonly reserved for console I/O).
! Otherwise, IUNIT is a value between 1 and 99, representing a
! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6
! are special, and will never return those values.
! Parameters:
! Output, integer IUNIT, the free unit number.
implicit none
integer i
integer iunit
logical value
iunit = 0
do i = 1, 99
if ( i .ne. 5 .and. i .ne. 6 .and. i .ne. 9 ) then
inquire ( unit = i, opened = value, err = 10 )
if ( .not. value ) then
iunit = i
return
end if
end if
10 continue
end do
return
end
subroutine mc ( s0, e, r, sigma, t1, m, seed, conf )
!*********************************************************************72
! MC uses Monte Carlo valuation on a European call.
! Reference:
! Desmond Higham,
! Black-Scholes for Scientific Computing Students,
! Computing in Science and Engineering,
! November/December 2004, Volume 6, Number 6, pages 72-79.
! Parameters:
! Input, double precision S0, the asset price at time 0.
! Input, double precision E, the exercise price.
! Input, double precision R, the interest rate.
! Input, double precision SIGMA, the volatility of the asset.
! Input, double precision T1, the expiry date.
! Input, integer M, the number of simulations.
! Input/output, integer SEED, a seed for the random
! number generator.
! Output, double precision CONF(2), the estimated range of the valuation.
implicit none
integer m
double precision conf(2)
double precision e
integer i
double precision pmean
double precision pvals(m)
double precision r
double precision s0
integer seed
double precision sigma
double precision std
double precision svals(m)
double precision t1
double precision u(m)
double precision width
call r8vec_normal_01 ( m, seed, u )
do i = 1, m
svals(i) = s0 * exp ( ( r - 0.5D+00 * sigma * sigma ) * t1 +&
sigma * sqrt ( t1 ) * u(i) )
end do
do i = 1, m
pvals(i) = exp ( - r * t1 ) * max ( svals(i) - e, 0.0D+00 )
end do
pmean = 0.0D+00
do i = 1, m
pmean = pmean + pvals(i)
end do
pmean = pmean / dble ( m )
std = 0.0D+00
do i = 1, m
std = std + ( pvals(i) - pmean ) ** 2
end do
std = sqrt ( std / dble ( m - 1 ) )
width = 1.96D+00 * std / sqrt ( dble ( m ) )
conf(1) = pmean - width
conf(2) = pmean + width
return
end
function r8_normal_01 ( seed )
!*********************************************************************72
! R8_NORMAL_01 returns a unit pseudonormal R8.
! Discussion:
! Because this routine uses the Box Muller method, it requires pairs
! of uniform random values to generate a pair of normal random values.
! This means that on every other call, the code can use the second
! value that it calculated.
! However, if the user has changed the SEED value between calls,
! the routine automatically resets itself and discards the saved data.
! Parameters:
! Input/output, integer SEED, a seed for the random number generator.
! Output, double precision R8_NORMAL_01, a sample of the standard normal PDF.
implicit none
double precision pi
parameter ( pi = 3.141592653589793D+00 )
double precision r1
double precision r2
double precision r8_normal_01
double precision r8_uniform_01
integer seed
integer seed1
integer seed2
integer seed3
integer used
double precision v1
double precision v2
save seed1
save seed2
save seed3
save used
save v2
data seed2 / 0 /
data used / 0 /
data v2 / 0.0D+00 /
!
! If USED is odd, but the input SEED does not match
! the output SEED on the previous call, then the user has changed
! the seed. Wipe out internal memory.
if ( mod ( used, 2 ) .eq. 1 ) then
if ( seed .ne. seed2 ) then
used = 0
seed1 = 0
seed2 = 0
seed3 = 0
v2 = 0.0D+00
end if
end if
!
! If USED is even, generate two uniforms, create two normals,
! return the first normal and its corresponding seed.
!
if ( mod ( used, 2 ) .eq. 0 ) then
seed1 = seed
r1 = r8_uniform_01 ( seed )
if ( r1 .eq. 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8_NORMAL_01 - Fatal error!'
write ( *, '(a)' ) ' R8_UNIFORM_01 returned a value of 0.'
stop
end if
seed2 = seed
r2 = r8_uniform_01 ( seed )
seed3 = seed
v1 = sqrt ( -2.0D+00 * log ( r1 ) ) * cos ( 2.0D+00 * pi * r2 )
v2 = sqrt ( -2.0D+00 * log ( r1 ) ) * sin ( 2.0D+00 * pi * r2 )
r8_normal_01 = v1
seed = seed2
!
! If USED is odd (and the input SEED matched the output value from
! the previous call), return the second normal and its corresponding seed.
!
else
r8_normal_01 = v2
seed = seed3
end if
used = used + 1
return
end
function r8_uniform_01 ( seed )
!*********************************************************************72
! R8_UNIFORM_01 returns a unit pseudorandom R8.
! Discussion:
! This routine implements the recursion
! seed = 16807 * seed mod ( 231 - 1 )
! r8_uniform_01 = seed / ( 231 - 1 )
! The integer arithmetic never requires more than 32 bits,
! including a sign bit.
! If the initial seed is 12345, then the first three computations are
! Input Output R8_UNIFORM_01
! SEED SEED
! 12345 207482415 0.096616
! 207482415 1790989824 0.833995
! 1790989824 2035175616 0.947702
! Reference:
! Paul Bratley, Bennett Fox, Linus Schrage,
! A Guide to Simulation,
! Springer Verlag, pages 201-202, 1983.
! Pierre L'Ecuyer,
! Random Number Generation,
! in Handbook of Simulation,
! edited by Jerry Banks,
! Wiley Interscience, page 95, 1998.
! Bennett Fox,
! Algorithm 647:
! Implementation and Relative Efficiency of Quasirandom
! Sequence Generators,
! ACM Transactions on Mathematical Software,
! Volume 12, Number 4, pages 362-376, 1986.
! Peter Lewis, Allen Goodman, James Miller,
! A Pseudo-Random Number Generator for the System/360,
! IBM Systems Journal,
! Volume 8, pages 136-143, 1969.
! Parameters:
! Input/output, integer SEED, the "seed" value, which should NOT be 0.
! On output, SEED has been updated.
! Output, double precision R8_UNIFORM_01, a new pseudorandom variate,
! strictly between 0 and 1.
implicit none
double precision r8_uniform_01
integer k
integer seed
if ( seed .eq. 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!'
write ( *, '(a)' ) ' Input value of SEED = 0.'
stop
end if
k = seed / 127773
seed = 16807 * ( seed - k * 127773 ) - k * 2836
if ( seed .lt. 0 ) then
seed = seed + 2147483647
end if
r8_uniform_01 = dble ( seed ) * 4.656612875D-10
return
end
subroutine r8vec_normal_01 ( n, seed, x )
!*********************************************************************72
! R8VEC_NORMAL_01 returns a unit pseudonormal R8VEC.
! Discussion:
! An R8VEC is a vector of R8's.
! The standard normal probability distribution function (PDF) has
! mean 0 and standard deviation 1.
! This routine can generate a vector of values on one call. It
! has the feature that it should provide the same results
! in the same order no matter how we break up the task.
! The Box-Muller method is used, which is efficient, but
! generates an even number of values each time. On any call
! to this routine, an even number of new values are generated.
! Depending on the situation, one value may be left over.
! In that case, it is saved for the next call.
! Parameters:
! Input, integer N, the number of values desired. If N is negative,
! then the code will flush its internal memory; in particular,
! if there is a saved value to be used on the next call, it is
! instead discarded. This is useful if the user has reset the
! random number seed, for instance.
! Input/output, integer SEED, a seed for the random number generator.
! Output, double precision X(N), a sample of the standard normal PDF.
! Local parameters:
! Local, integer MADE, records the number of values that have
! been computed. On input with negative N, this value overwrites
! the return value of N, so the user can get an accounting of
! how much work has been done.
! Local, integer SAVED, is 0 or 1 depending on whether there is a
! single saved value left over from the previous call.
! Local, integer X_LO_INDEX, X_HI_INDEX, records the range of entries of
! X that we need to compute. This starts off as 1:N, but is adjusted
! if we have a saved value that can be immediately stored in X(1),
! and so on.
! Local, double precision Y, the value saved from the previous call, if
! SAVED is 1.
implicit none
integer n
integer i
integer m
integer made
double precision pi
parameter ( pi = 3.141592653589793D+00 )
double precision r(2)
double precision r8_uniform_01
integer saved
integer seed
double precision x(n)
integer x_hi_index
integer x_lo_index
double precision y
save made
save saved
save y
data made / 0 /
data saved / 0 /
data y / 0.0D+00 /
!
! I'd like to allow the user to reset the internal data.
! But this won't work properly if we have a saved value Y.
! I'm making a crock option that allows the user to signal
! explicitly that any internal memory should be flushed,
! by passing in a negative value for N.
!
if ( n .lt. 0 ) then
n = made
made = 0
saved = 0
y = 0.0D+00
return
else if ( n .eq. 0 ) then
return
end if
! Record the range of X we need to fill in.
x_lo_index = 1
x_hi_index = n
! Use up the old value, if we have it.
if ( saved .eq. 1 ) then
x(1) = y
saved = 0
x_lo_index = 2
end if
! Maybe we don't need any more values.
if ( x_hi_index - x_lo_index + 1 .eq. 0 ) then
! If we need just one new value, do that here to avoid null arrays.
else if ( x_hi_index - x_lo_index + 1 .eq. 1 ) then
r(1) = r8_uniform_01 ( seed )
if ( r(1) .eq. 0.0D+00 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8VEC_NORMAL_01 - Fatal errorc'
write ( *, '(a)' ) ' R8_UNIFORM_01 returned a value of 0.'
stop
end if
r(2) = r8_uniform_01 ( seed )
x(x_hi_index) =sqrt ( -2.0D+00 * log ( r(1) ) )&
* cos ( 2.0D+00 * pi * r(2) )
y = sqrt ( -2.0D+00 * log ( r(1) ) )* sin ( 2.0D+00 * pi * r(2) )
saved = 1
made = made + 2
! If we require an even number of values, that's easy.
else if ( mod ( x_hi_index - x_lo_index + 1, 2 ) .eq. 0 ) then
do i = x_lo_index, x_hi_index, 2
call r8vec_uniform_01 ( 2, seed, r )
x(i) = sqrt ( -2.0D+00 * log ( r(1) ) )* cos ( 2.0D+00 * pi * r(2) )
x(i+1) = sqrt ( -2.0D+00 * log ( r(1) ) )* sin ( 2.0D+00 * pi * r(2) )
end do
made = made + x_hi_index - x_lo_index + 1
!
! If we require an odd number of values, we generate an even number,
! and handle the last pair specially, storing one in X(N), and
! saving the other for later.
!
else
do i = x_lo_index, x_hi_index - 1, 2
call r8vec_uniform_01 ( 2, seed, r )
x(i) = sqrt ( -2.0D+00 * log ( r(1) ) )* cos ( 2.0D+00 * pi * r(2) )
x(i+1) = sqrt ( -2.0D+00 * log ( r(1) ) ) * sin ( 2.0D+00 * pi * r(2) )
end do
call r8vec_uniform_01 ( 2, seed, r )
x(n) = sqrt ( -2.0D+00 * log ( r(1) ) ) * cos ( 2.0D+00 * pi * r(1) )
y = sqrt ( -2.0D+00 * log ( r(2) ) ) * sin ( 2.0D+00 * pi * r(2) )
saved = 1
made = made + x_hi_index - x_lo_index + 2
end if
return
end
subroutine r8vec_print_part ( n, a, max_print, title )
!*********************************************************************72
!
! R8VEC_PRINT_PART prints "part" of an R8VEC.
! Discussion:
! The user specifies MAX_PRINT, the maximum number of lines to print.
! If N, the size of the vector, is no more than MAX_PRINT, then
! the entire vector is printed, one entry per line.
! Otherwise, if possible, the first MAX_PRINT-2 entries are printed,
! followed by a line of periods suggesting an omission,
! and the last entry.
! Parameters:
! Input, integer N, the number of entries of the vector.
! Input, double precision A(N), the vector to be printed.
! Input, integer MAX_PRINT, the maximum number of lines to print.
! Input, character*(*) TITLE, a title.
!
implicit none
integer n
double precision a(n)
integer i
integer max_print
character*(*) title
if ( max_print .le. 0 ) then
return
end if
if ( n .le. 0 ) then
return
end if
write ( *, '(a)' ) ' '
write ( *, '(a)' ) title
write ( *, '(a)' ) ' '
if ( n .le. max_print ) then
do i = 1, n
write ( *, '(2x,i8,a1,1x,g14.6)' ) i, ':', a(i)
end do
else if ( 3 .le. max_print ) then
do i = 1, max_print - 2
write ( *, '(2x,i8,a1,1x,g14.6)' ) i, ':', a(i)
end do
write ( *, '(a)' ) ' ........ ..............'
i = n
write ( *, '(2x,i8,a1,1x,g14.6)' ) i, ':', a(i)
else
do i = 1, max_print - 1
write ( *, '(2x,i8,a1,1x,g14.6)' ) i, ':', a(i)
end do
i = max_print
write ( *, '(2x,i8,a1,1x,g14.6,a)' ) i, ':', a(i), '...more entries...'
end if
return
end
subroutine r8vec_uniform_01 ( n, seed, r )
!*********************************************************************72
!
! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
! Discussion:
! An R8VEC is a vector of R8's.
! Reference:
! Paul Bratley, Bennett Fox, Linus Schrage,
! A Guide to Simulation,
! Springer Verlag, pages 201-202, 1983.
! Bennett Fox,
! Algorithm 647:
! Implementation and Relative Efficiency of Quasirandom
! Sequence Generators,
! ACM Transactions on Mathematical Software,
! Volume 12, Number 4, pages 362-376, 1986.
! Peter Lewis, Allen Goodman, James Miller,
! A Pseudo-Random Number Generator for the System/360,
! IBM Systems Journal,
! Volume 8, pages 136-143, 1969.
! Parameters:
! Input, integer N, the number of entries in the vector.
! Input/output, integer SEED, the "seed" value, which should NOT be 0.
! On output, SEED has been updated.
! Output, double precision R(N), the vector of pseudorandom values.
!
implicit none
integer n
integer i
integer k
integer seed
double precision r(n)
do i = 1, n
k = seed / 127773
seed = 16807 * ( seed - k * 127773 ) - k * 2836
if ( seed .lt. 0 ) then
seed = seed + 2147483647
end if
r(i) = dble ( seed ) * 4.656612875D-10
end do
return
end
subroutine r8vec_write ( output_filename, n, x )
!*********************************************************************72
!
! R8VEC_WRITE writes an R8VEC file.
! Discussion:
! An R8VEC is a vector of R8's.
! Parameters:
! Input, character * ( * ) OUTPUT_FILENAME, the output file name.
! Input, integer N, the number of points.
! Input, double precision X(N), the data.
implicit none
integer m
integer n
integer j
character * ( * ) output_filename
integer output_unit
double precision x(n)
! Open the file.
call get_unit ( output_unit )
open ( unit = output_unit, file = output_filename, status = 'replace' )
! Create the format string.
if ( 0 .lt. n ) then
! Write the data.
do j = 1, n
write ( output_unit, '(2x,g24.16)' ) x(j)
end do
end if
! Close the file.
close ( unit = output_unit )
return
end
subroutine timestamp ( )
!*********************************************************************72
! TIMESTAMP prints out the current YMDHMS date as a timestamp.
! Parameters:
! None
implicit none
character * ( 8 ) ampm
integer d
character * ( 8 ) date
integer h
integer m
integer mm
character * ( 9 ) month(12)
integer n
integer s
character * ( 10 ) time
integer y
save month
data month /'January ', 'February ', 'March ', 'April ',&
'May ', 'June ', 'July ', 'August ',&
'September', 'October ', 'November ', 'December ' /
call date_and_time ( date, time )
read ( date, '(i4,i2,i2)' ) y, m, d
read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm
if ( h .lt. 12 ) then
ampm = 'AM'
else if ( h .eq. 12 ) then
if ( n .eq. 0 .and. s .eq. 0 ) then
ampm = 'Noon'
else
ampm = 'PM'
end if
else
h = h - 12
if ( h .lt. 12 ) then
ampm = 'PM'
else if ( h .eq. 12 ) then
if ( n .eq. 0 .and. s .eq. 0 ) then
ampm = 'Midnight'
else
ampm = 'AM'
end if
end if
end if
write ( *,'(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) d,&
month(m), y, h, ':', n, ':', s, '.', mm, ampm
return
end
!c$$ asset_path.txt, 图形数据文档:
11 September 2022 3:30:47.075 PM BLACK_SCHOLES_PRB: FORTRAN95,2013&18 version Test the BLACK_SCHOLES library.
ASSET_PATH_TEST: Demonstrate the simulated of an asset price path.
The asset price at time 0 S0 = 2.00000 The asset expected growth rate MU = 0.100000 The asset volatility SIGMA = 0.300000
The expiry date T1 = 1.00000 The number of time steps N
= 100 The random number seed was SEED = 123456789 Partial results:
1: 2.00000 2: 2.10353 3: 2.07412 4: 2.03940 5: 2.02551 6: 2.10078 7: 2.13498 8: 2.21808 ........
.............. 101: 2.44083 Full results written to "asset_path.txt". BINOMIAL_TEST:
A demonstration of the binomial method for option valuation. The asset price at time 0 S0 =
2.00000 The exercise price E = 1.00000 The interest rate R = 0.500000E-01 The asset volatility SIGMA = 0.250000
The expiry date T1 = 3.00000 The number of intervals M = 256 The
option value is 1.14476 BSF_TEST: A demonstration of the Black-Scholes formula for option valuation.
The asset price at time T0 S0 = 2.00000 The time T0 = 0.00000 The exercise
price E = 1.00000 The interest rate R = 0.500000E-01 The asset volatility SIGMA = 0.250000
The expiry date T1 = 3.00000 The option value C = 1.14474 FORWARD_TEST: A
demonstration of the forward difference method for option valuation. The exercise price E = 4.00000
The interest rate R = 0.300000E-01 The asset volatility SIGMA = 0.500000 The
expiry date T1 = 1.00000 The number of space steps NX = 11 The number of time steps NT = 29 The value of SMAX = 10.0000
Initial Option Value Value 1.00000 0.139363E-02 2.00000
0.373367E-01 3.00000 0.223638 4.00000 0.627210 5.00000 1.20992 6.00000 1.91439 7.00000 2.69543 8.00000
3.52261 9.00000 4.37638 10.0000 5.24428 MC_TEST: A demonstration of the
Monte Carlo method for option valuation. The asset price at time 0, S0 = 2.00000 The exercise price E = 1.00000
The interest rate R = 0.500000E-01 The asset volatility SIGMA =
0.250000 The expiry date T1 = 3.00000 The number of simulations M = 1000000 The confidence interval is [ 1.14311 , 1.14663 ].
BLACK_SCHOLES_PRB: Normal end of execution. 11
September 2022 3:30:47.161 PM
几个简单方法实现欧洲电话的毕苏期权定价方程模式.
广州
更多回帖
