SDEClass.f90 Source File

module SDEClass
   use RandomClass
   use ArrayClass
   use MathClass
   implicit none

   integer(int32) :: PF_SDE_GeometricBrown = 1
   type :: SDE_
      real(real64), allocatable :: process(:)
      real(real64) :: c = 0.0d0
      real(real64) :: sigma = 0.0d0
      integer(int32) :: SDEType = 0
   contains
      procedure, public :: init => initSDE
      procedure, public :: solve => solveSDE
   end type
contains

   subroutine initSDE(obj, SDEType, c, sigma)
      class(SDE_), intent(inout) :: obj
      integer(int32), intent(in) :: SDEType
      real(real64), optional, intent(in) :: c, sigma

      if (SDEType == PF_SDE_GeometricBrown) then
         ! Geometric Brown Motion
         if (.not. present(c)) then
            print *, "ERROR :: initSDE >> arg c is necessary."
            stop
         end if
         if (.not. present(sigma)) then
            print *, "ERROR :: initSDE >> arg sigma is necessary."
            stop
         end if
         obj%c = c
         obj%sigma = sigma
         obj%SDEType = PF_SDE_GeometricBrown

      end if

   end subroutine

   function solveSDE(obj, X0, dt, step) result(ret)
      class(SDE_), intent(inout) :: obj
      real(real64), intent(in) :: X0, dt
      integer(int32), intent(in) :: step
      real(real64) ::t, Bt
      real(real64), allocatable :: ret(:)
      integer(int32) :: n, i
      type(Random_) :: random

      if (obj%SDEType == PF_SDE_GeometricBrown) then
         ! Geometric Brown Motion

         n = step
         ret = zeros(n)

         Bt = 0.0d0
         t = 0.0d0
         ret(1) = X0
         do i = 2, n
            t = t + dt
            Bt = Bt + random%gauss(mu=0.0d0, sigma=1.0d0)
            ret(i) = X0*exp((obj%c - 0.50d0*obj%sigma*obj%sigma)*t + obj%sigma*Bt)
         end do
      end if

   end function

end module SDEClass