Program g05xdfe
! G05XDF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
! .. Use Statements ..
Use nag_library, Only: g05xcf, g05xdf, g05xef, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: a = 0, d = 1, nout = 6, rcord = 2
Real (Kind=nag_wp), Parameter :: c(d) = (/1.0_nag_wp/)
Real (Kind=nag_wp), Parameter :: diff(d) = (/0.0_nag_wp/)
! .. Local Scalars ..
Real (Kind=nag_wp) :: r, s0, sigma, t0, tend
Integer :: bgord, i, ifail, ldb, ldz, nmove, &
npaths, ntimesteps, p
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: analytic(:), b(:,:), rcomm(:), &
st(:,:), t(:), times(:), z(:,:)
Integer, Allocatable :: move(:)
! .. Intrinsic Procedures ..
Intrinsic :: exp, real, size, sqrt
! .. Executable Statements ..
ifail = 0
! We wish to solve the stochastic differential equation (SDE)
! dSt = r*St*dt + sigma*St*dXt
! where X is a one dimensional Wiener process.
! This means we have
! A = 0
! D = 1
! C = 1
! We now set the other parameters of the SDE and the Euler-Maruyama scheme
! Initial value of the process
s0 = 1.0_nag_wp
r = 0.05_nag_wp
sigma = 0.12_nag_wp
! Number of paths to simulate
npaths = 3
! The time interval [t0,T] on which to solve the SDE
t0 = 0.0_nag_wp
tend = 1.0_nag_wp
! The time steps in the discretization of [t0,T]
ntimesteps = 20
Allocate (t(ntimesteps))
Do i = 1, ntimesteps
t(i) = t0 + i*(tend-t0)/real(ntimesteps+1,kind=nag_wp)
End Do
! Make the bridge construction order
nmove = 0
Allocate (times(ntimesteps),move(nmove))
bgord = 3
Call g05xef(bgord,t0,tend,ntimesteps,t,nmove,move,times,ifail)
! Generate the input Z values and allocate memory for b
Call get_z(rcord,npaths,d,a,ntimesteps,z,b)
! Leading dimensions for the various input arrays
ldz = size(z,1)
ldb = size(b,1)
! Initialize the generator
Allocate (rcomm(12*(ntimesteps+1)))
Call g05xcf(t0,tend,times,ntimesteps,rcomm,ifail)
! Get the scaled increments of the Wiener process
Call g05xdf(npaths,rcord,d,a,diff,z,ldz,c,d,b,ldb,rcomm,ifail)
! Do the Euler-Maruyama time stepping
Allocate (st(npaths,ntimesteps+1),analytic(npaths))
! Do first time step
st(:,1) = s0 + (t(1)-t0)*(r*s0+sigma*s0*b(:,1))
Do i = 2, ntimesteps
Do p = 1, npaths
st(p,i) = st(p,i-1) + (t(i)-t(i-1))*(r*st(p,i-1)+sigma*st(p,i-1)*b(p &
,i))
End Do
End Do
! Do last time step
st(:,i) = st(:,i-1) + (tend-t(i-1))*(r*st(:,i-1)+sigma*st(:,i-1)*b(:,i))
! Compute the analytic solution:
! ST = S0*exp( (r-sigma**2/2)T + sigma WT )
analytic(:) = s0*exp((r-0.5_nag_wp*sigma*sigma)*tend+sigma*sqrt(tend-t0) &
*z(:,1))
! Display the results
Call display_results(ntimesteps,npaths,st,analytic)
Contains
Subroutine get_z(rcord,npaths,d,a,ntimes,z,b)
! .. Use Statements ..
Use nag_library, Only: g05yjf
! .. Scalar Arguments ..
Integer, Intent (In) :: a, d, npaths, ntimes, rcord
! .. Array Arguments ..
Real (Kind=nag_wp), Allocatable, Intent (Out) :: b(:,:), z(:,:)
! .. Local Scalars ..
Integer :: idim, ifail
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: std(:), tz(:,:), xmean(:)
Integer, Allocatable :: iref(:), state(:)
Integer :: seed(1)
! .. Intrinsic Procedures ..
Intrinsic :: transpose
! .. Executable Statements ..
idim = d*(ntimes+1-a)
! Allocate Z
If (rcord==1) Then
Allocate (z(idim,npaths))
Allocate (b(d*(ntimes+1),npaths))
Else
Allocate (z(npaths,idim))
Allocate (b(npaths,d*(ntimes+1)))
End If
! We now need to generate the input quasi-random points
! First initialize the base pseudorandom number generator
seed(1) = 1023401
Call initialize_prng(6,0,seed,size(seed),state)
! Scrambled quasi-random sequences preserve the good discrepancy
! properties of quasi-random sequences while counteracting the bias
! some applications experience when using quasi-random sequences.
! Initialize the scrambled quasi-random generator.
Call initialize_scrambled_qrng(1,2,idim,state,iref)
! Generate the quasi-random points from N(0,1)
Allocate (xmean(idim),std(idim))
xmean(1:idim) = 0.0_nag_wp
std(1:idim) = 1.0_nag_wp
If (rcord==1) Then
Allocate (tz(npaths,idim))
ifail = 0
Call g05yjf(xmean,std,npaths,tz,iref,ifail)
z(:,:) = transpose(tz)
Else
ifail = 0
Call g05yjf(xmean,std,npaths,z,iref,ifail)
End If
End Subroutine get_z
Subroutine initialize_prng(genid,subid,seed,lseed,state)
! .. Use Statements ..
Use nag_library, Only: g05kff
! .. Scalar Arguments ..
Integer, Intent (In) :: genid, lseed, subid
! .. Array Arguments ..
Integer, Intent (In) :: seed(lseed)
Integer, Allocatable, Intent (Out) :: state(:)
! .. Local Scalars ..
Integer :: ifail, lstate
! .. Executable Statements ..
! Initial call to initializer to get size of STATE array
lstate = 0
Allocate (state(lstate))
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
! Reallocate STATE
Deallocate (state)
Allocate (state(lstate))
! Initialize the generator to a repeatable sequence
ifail = 0
Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)
End Subroutine initialize_prng
Subroutine initialize_scrambled_qrng(genid,stype,idim,state,iref)
! .. Use Statements ..
Use nag_library, Only: g05ynf
! .. Scalar Arguments ..
Integer, Intent (In) :: genid, idim, stype
! .. Array Arguments ..
Integer, Allocatable, Intent (Out) :: iref(:)
Integer, Intent (Inout) :: state(*)
! .. Local Scalars ..
Integer :: ifail, iskip, liref, nsdigits
! .. Executable Statements ..
liref = 32*idim + 7
iskip = 0
nsdigits = 32
Allocate (iref(liref))
ifail = 0
Call g05ynf(genid,stype,idim,iref,liref,iskip,nsdigits,state,ifail)
End Subroutine initialize_scrambled_qrng
Subroutine display_results(ntimesteps,npaths,st,analytic)
! .. Scalar Arguments ..
Integer, Intent (In) :: npaths, ntimesteps
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: analytic(:), st(:,:)
! .. Local Scalars ..
Integer :: i, p
! .. Executable Statements ..
Write (nout,*) 'G05XDF Example Program Results'
Write (nout,*)
Write (nout,*) 'Euler-Maruyama solution for Geometric Brownian motion'
Write (nout,99999)('Path',p,p=1,npaths)
Do i = 1, ntimesteps + 1
Write (nout,99998) i, st(:,i)
End Do
Write (nout,*)
Write (nout,'(A)') 'Analytic solution at final time step'
Write (nout,99999)('Path',p,p=1,npaths)
Write (nout,'(4X,20(1X,F10.4))') analytic(:)
99999 Format (4X,20(5X,A,I2))
99998 Format (1X,I2,1X,20(1X,F10.4))
End Subroutine display_results
End Program g05xdfe