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