C.hr stng1.f
C....*...1.........2.........3.........4.........5.........6.........7.*
C    
C     stng1.f
C
C     Copyright (C) 1995, 2003.
C
C     A. Ronald Gallant, P.O. Box 659, Chapel Hill NC 27514-0659, USA
C
C     Permission to use, copy, modify, and distribute this software and 
C     its documentation for any purpose and without fee is hereby 
C     granted, provided that the above copyright notice appear in all
C     copies and that both that copyright notice and this permission 
C     notice appear in supporting documentation.
C     
C     This software is provided "as is" without any expressed or implied 
C     warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*
C.hr stng1
C@
*....*...1.........2.........3.........4.........5.........6.........7.*
*
*     stng1     12/23/94
*
*     purpose
*     compute a realization of length n from an explicit order 1 strong
*     scheme for the nonautonomous Ito stochastic differential equation
*
*       dX = A(t,X)dt + B(t,X)dW
*
*     where X has dimension d and W has dimension m.
*
*     usage
*     call stng1(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y)
*
*     libf calls
*     unsk
*
*     arguments
*     drift  - subroutine with usage call drift(t,X,d,A) that is
*              declared external in the calling program.
*              input, external
*     difuse - subroutine with usage call difuse(t,X,d,m,B) that is
*              declared external in the calling program.
*              input, external
*     t      - input to drift and difuse, on return t=t0+n*dt.
*              workspace, real*8
*     X      - input to drift and difuse, a vector of length d.  on
*              return X(i)=Y(i,n) for i=1,...,d.
*              workspace, real*8
*     A      - output of drift, a vector of length d.
*              workspace, real*8
*     B      - output of difuse, a matrix of order d by m.
*              workspace, real*8
*     d      - dimension of X and A, row dimension of B and Y.  d must
*              be less than parameter md which is currently set to 15.
*              input, integer*4
*     m      - column dimension of B.  m must be less than parameter mm
*              which is currently set to 5.
*              input, integer*4
*     t0     - time of initial condition.
*              input, real*8
*     X0     - initial condition at time t=t0, a vector of length d.
*              input, real*8
*     dt     - time increment for the simulations.
*              input, real*8
*     iseed  - seed for the simulations.
*              input and output, integer*4
*     n      - number of simulations, column dimension of Y.
*              input, integer*4
*     Y      - computed simulations, a d by n matrix.
*              output, real*8
*
*     remarks
*     To reduce memory requirements for a computation such as
*     (1/T)(integral from 0 to T of func[X(t)]) where T = n*dt,
*     dimension as Y(d,n0) and use
*       do k=1,n,n0
*         call stng1(drift,difuse,t0,X0,A,B,d,m,t0,X0,dt,iseed,n0,Y)
*         do it=1,n0
*           do i=1,d
*             ave(i) = ave(i) + func(Y(i,it))/n
*           end do
*         end do
*       end do
*     instead of dimensioning as Y(d,n) and using
*       call stng1(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y)
*       do it=1,n
*         do i=1,d
*           ave(i) = ave(i) + func(Y(i,it))/n
*         end do
*       end do
*     
*     reference
*     Kloeden, Peter E., and Eckhard Platen (1992), Numerical Solution
*     of Stochastic Differential Equations, Berlin, Springer-Verlag,
*     p. 347 (integral approximation) and p. 376 (explicit scheme).
*
      subroutine stng1(drift,difuse,t,X,A,B,d,m,t0,X0,dt,iseed,n,Y)
      implicit none

*     arguments

      external drift, difuse

      real*8 X, A, B, X0, Y
      real*8 t,t0,dt

      integer*4 d, m, iseed, n

      dimension X(d), A(d), B(d,m), X0(d), Y(d,n)

*     parameters

      real*8 pi
      parameter (pi=3.141592653589793d0)

      integer*4 p
      parameter (p=50)
*     p should be greater than 0.05/dt, 50=0.05/1.d-3

      integer*4 md, mm
      integer*4 nout
      parameter (md=15, mm=5)
      parameter (nout=3)

*     local variables

      intrinsic DSQRT

      real*4 z,unsk
      real*8 z0,z1,z2

      real*8 dW,Ip
      real*8 Tn,BTn
      real*8 sum

      real*8 tmp,r,left,right
      real*8 rdt,rho,rrho,rtwo

      integer*4 i,j,k,it,j1,j2

      dimension z0(mm),z1(mm),z2(mm)
      dimension dW(mm),Ip(mm,mm)
      dimension Tn(md,mm),BTn(md*mm)
      dimension sum(md)

*     error traps

      if (m.gt.mm) then
        write(nout,'(a,i5)') 'Error, stng1, mm < m =', m
        stop
      end if

      if (d.gt.md) then
        write(nout,'(a,i5)') 'Error, stng1, md < d =', d
        stop
      end if

*     explicit order 1 strong scheme

      rho = 0.d0
      do k=1,p
        r = k
        rho =  rho + 1.d0/(r**2)
      end do
      rho = 1.d0/12.d0 - rho/(2.d0*(pi**2))

      rdt  = DSQRT(dt)
      rrho = DSQRT(rho)
      rtwo = DSQRT(2.d0)

      do i=1,d
        X(i)=X0(i)
      end do

      t = t0

      do it=1,n

        call drift(t,X,d,A)
        call difuse(t,X,d,m,B)

        do j=1,m
          do i=1,d
            Tn(i,j) = X(i) + A(i)*dt + B(i,j)*rdt
          end do
        end do

        do j=1,m
          z = unsk(iseed)
          z0(j) = z
          dW(j) = rdt*z
        end do

        if (m.eq.1) then

          Ip(1,1) = 0.5d0*(dW(1)**2 - dt)

        else

          do j=1,m
            Ip(j,j) = 0.5d0*(dW(j)**2 - dt)
            z = unsk(iseed)
            z1(j) = z
          end do

          do j1=1,m
          do j2=1,j1-1
            left  = 0.5d0*z0(j1)*z0(j2)
            right = rrho*(z1(j1)*z0(j2) - z1(j2)*z0(j1))
            Ip(j1,j2) = dt*(left + right)
            Ip(j2,j1) = dt*(left - right)
          end do
          end do

          do k=1, p
            r = k

            do j=1,m
              z = unsk(iseed)
              z1(j) = z
              z = unsk(iseed)
              z2(j) = z
            end do

            do j1=1,m
            do j2=1,j1-1
              left  = z2(j1)*(rtwo*z0(j2) + z1(j2))
              right = z2(j2)*(rtwo*z0(j1) + z1(j1))
              tmp = (0.5d0/pi)*(dt/r)*(left - right)
              Ip(j1,j2) = Ip(j1,j2) + tmp
              Ip(j2,j1) = Ip(j2,j1) - tmp
            end do
            end do

          end do

        end if

        do i=1,d
          sum(i) = 0.d0
        end do

        do j=1,m
          do i=1,d
            sum(i) = sum(i) + B(i,j)*dW(j)
          end do
        end do

        do j1=1,m
          call difuse(t,Tn(1,j1),d,m,BTn)
        do j2=1,m
          do i=1,d
            sum(i) = sum(i)
     &       + (Ip(j1,j2)/rdt) * ( BTn(-d+d*j2+i) - B(i,j2) )
          end do
        end do
        end do

        do i=1,d
          Y(i,it) = X(i) + A(i)*dt + sum(i)
        end do

*       For a stiff system the above loop is replaced by an equation
*       solver to solve Y = X + A(it,Y)dt + sum for Y.  See p. 407 of
*       Kloeden and Platen (1992).  Subroutine drift should be revised
*       to return first derivatives with respect to its argument X to
*       allow use of a solver that uses analytic first derivatives.
 
        do i=1,d
          X(i) = Y(i,it)
        end do

        t = t + dt

      end do
 
      return
      end