! Compute European Option Premium for a Cash-or-Nothing
! and a Vertical Spread Call.
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use umach_int
      implicit none
! The strike price
      real(kind(1e0)) :: ks1 = 10.0e0
! The spread value
      real(kind(1e0)) :: ks2 = 15.0e0
! The Bet for the Cash-or-Nothing Call
      real(kind(1e0)) :: bet = 2.0e0
! The sigma value
      real(kind(1e0)) :: sigma = 0.4e0

! Time values for the options
      integer, parameter :: nt = 2
      real(kind(1e0)) :: time(nt)=(/0.25e0, 0.5e0/)
! Values of the underlying where evaluation are made
      integer, parameter :: nv = 12, nlbc = 3, nrbc = 3
      integer :: i
      real(kind(1e0)) :: xs(nv) = (/(2+(I-1)*2.0e0,I=1,NV)/)
! Value of the interest rate and continuous dividend -
      real(kind(1e0)) :: r = 0.1e0, dividend = 0.0e0
! Values of the min and max underlying values modeled -
      real(kind(1e0)) :: x_min = 0.0e0, x_max = 30.0e0
 
! Define parameters for the integration step.  
      integer, parameter :: nx = 61, nint = nx-1, n=3*nx
      real(kind(1e0)) :: xgrid(nx), yb(n,0:nt), ybprime(n,0:nt),&
                         yv(n,0:nt), yvprime(n,0:nt),&
                         dx, fb(nv,nt), fv(nv,nt)
      type(s_fcn_data) fcn_data
      integer :: nout
      real(kind(1e0)), external :: fkcoef_call, fkinitcond_call
      external fkbc_call
      
      call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
      allocate(fcn_data % rdata (7), fcn_data % idata (1))
! Define an equally-spaced grid of points for the underlying price
      dx = (x_max-x_min)/real(nint)
      xgrid(1) = x_min
      xgrid(nx) = x_max
      do i = 2,nx-1
        xgrid(i) = xgrid(i-1) + dx
      end do

      fcn_data % rdata = (/ks1,bet,ks2,x_max,sigma,r,dividend/)

! Flag the difference in payoff functions - 
! 1 for the Bet, 2 for the Vertical Spread
      fcn_data % idata(1) = 1
      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_call,&
                        fkinitcond_call, fkbc_call, yb, ybprime,&
                        FCN_DATA = fcn_data)

      fcn_data % idata(1) = 2 
      call feynman_kac (Xgrid, time, nlbc, nrbc, fkcoef_call,&
                        fkinitcond_call, fkbc_call, yv, yvprime,&
                        FCN_DATA = fcn_data)

! Evaluate solutions at vector of points XS(:), at each time value
! prior to expiration.                        
      do i=1,nt
        fb(:,i) = hqsval (xs, xgrid, yb(:,I))
        fv(:,i) = hqsval (xs, xgrid, yv(:,I))
      end do 
      write(nout,'(T05,A)') "European Option Value for A Bet",&
                       " and a Vertical Spread, 3 and 6 Months "//&
                       "Prior to Expiry"
      write(nout,'(T08,"Number of equally spaced spline knots "'//&
                  ',I4,/T08,"Number of unknowns ",I4)') NX,N
      write(nout,'(T08,"Strike = ",F5.2,", Sigma =", F5.2,'//&
           '", Interest Rate =",F5.2,'//&
           '/T08,"Bet = ",F5.2,", Spread Value = ", F5.2/'//&
           '/T10,"Underlying", T32,"A Bet", T40,"Vertical Spread",'//&
           '/(T10,5F9.4))') ks1, sigma, r, bet, ks2, &
             (xs(i), fb(i,1:nt), fv(i,1:nt),i=1,nv)
      end

! These routines define the coefficients, payoff, boundary
! conditions and forcing term for American and European Options.
      function fkcoef_call (x, tx, iflag, fcn_data) result(value)
      use mp_types
      implicit none

      integer, intent(inout) :: iflag
      real(kind(1e0)), intent(in) :: x, tx
      type(s_fcn_data), optional :: fcn_data
      real(kind(1e0)) :: value

      real(kind(1e0)) :: sigma, interest_rate, dividend
! Data passed through using allocated components
! of the derived type s_fcn_data
      sigma = fcn_data % rdata(5)
      interest_rate = fcn_data % rdata(6)
      dividend = fcn_data % rdata(7)
      select case (iflag)
        case (1)
! The coefficient derivative d(sigma)/dx		
          value = sigma
! The coefficient sigma(x)				
        case (2)
          value = sigma * x
        case (3)
! The coefficient mu(x)			
          value = (interest_rate - dividend) * x
        case (4)
! The coefficient kappa(x)
          value = interest_rate
      end select
! Note that there is no time dependence
      iflag = 0
      return 
      end function fkcoef_call

      function fkinitcond_call(x, fcn_data) result(value)
      use mp_types
      implicit none

      real(kind(1e0)), intent(in) :: x
      type(s_fcn_data), optional :: fcn_data
      real(kind(1e0)) :: value

      real(kind(1e0)) :: strike_price, spread, bet
      real(kind(1e0)), parameter :: zero = 0.0e0

      strike_price = fcn_data % rdata(1)
      bet = fcn_data % rdata(2)
      spread = fcn_data % rdata(3)
! The payoff function - Use flag passed to decide which
      select case (fcn_data % idata(1))
        case(1)
! After reaching the strike price the payoff jumps
! from zero to the bet value.                     
        value = zero
        if (x > strike_price) value = bet
        case(2)
! Function is zero up to strike price.
! Then linear between strike price and spread.
! Then has constant value Spread-Strike Price after
! the value Spread.                    		 
          value = max(x-strike_price, zero) - max(x-spread, zero)
      end select
      return
      end function fkinitcond_call

      subroutine fkbc_call (TX, iflag, bccoefs, fcn_data)
      use mp_types
      implicit none

      real(kind(1e0)), intent(in) :: tx
      integer, intent(inout) :: iflag
      real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
      type(s_fcn_data), optional :: fcn_data

      real(kind(1e0)) :: strike_price, spread, bet,&
                         interest_rate, df

      strike_price = fcn_data % rdata(1)
      bet = fcn_data % rdata(2)
      spread = fcn_data % rdata(3)
      interest_rate = fcn_data % rdata(6)
      select case (iflag)
        case (1)
          bccoefs(1,1:4) = ((/1.0e0, 0.0e0, 0.0e0, 0.0e0/))
          bccoefs(2,1:4) = ((/0.0e0, 1.0e0, 0.0e0, 0.0e0/))
          bccoefs(3,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
        case (2)
! This is the discount factor using the risk-free
! interest rate
          df = exp(interest_rate * tx)
! Use flag passed to decide on boundary condition -
          select case (fcn_data % idata(1))
            case(1) 
              bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, bet*df/) 
            case(2)			 
              bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,&
                               (spread-strike_price)*df/) 
          end select
          bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 0.0e0/)
          bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
          return        			 
      end select
! Note no time dependence in case (1) for iflag
      iflag = 0
      end subroutine fkbc_call
