几个简单方法实现欧洲电话的毕苏期权定价方程模式

!********************************************************************
!********************************************************************
! 几个简单方法实现欧洲电话的毕苏期权定价方程模式
!********************************************************************
! 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

几个简单方法实现欧洲电话的毕苏期权定价方程模式.

广州