Program g05xdfe ! G05XDF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. 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 ! .. 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(:,:) Real (Kind=nag_wp) :: c(d) = (/1.0_nag_wp/) Real (Kind=nag_wp) :: diff(d) = (/0.0_nag_wp/) Integer, Allocatable :: move(:) ! .. Intrinsic Procedures .. Intrinsic :: exp, real, size ! .. 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*(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,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,state) ! .. Use Statements .. Use nag_library, Only: g05kff ! .. Scalar Arguments .. Integer, Intent (In) :: genid, subid ! .. Array Arguments .. Integer, Intent (In) :: seed(:) Integer, Allocatable, Intent (Out) :: state(:) ! .. Local Scalars .. Integer :: ifail, lseed, lstate ! .. Executable Statements .. lseed = size(seed,1) ! 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