```    Program g05xdfe

!     G05XDF Example Program Text

!     Mark 27.0 Release. NAG Copyright 2019.

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