! Greeks computation
      use feynman_kac_int
      use hqsval_int
      use mp_types
      use anordf_int
      use const_int
      use umach_int

      implicit none
      
      real(kind(1d0)), external :: fkcoef, fkinitcond 
      external fkbc

! The set of strike prices
      real(kind(1d0)) :: ks(3) = (/15.0d0,20.0d0,25.0d0/)
! The set of sigma values
      real(kind(1d0)) :: sigma(3) = (/0.2d0, 0.3d0, 0.4d0/)
! The set of model diffusion powers: alpha = 2.0 <==> Black Scholes
      real(kind(1d0)) :: alpha(3) = (/2.0d0, 1.0d0, 0.0d0/)
! Time values for the options
      integer, parameter :: nt = 3
      real(kind(1d0)) :: time(nt)=(/1.d0/12., 4.d0/12., 7.d0/12./)
! Values of the min and max underlying values modeled
      real(kind(1d0)) :: x_min = 0.0d0, x_max = 60.0d0
! Value of the interest rate and continuous dividend
      real(kind(1d0)) :: r = 0.05d0, dividend = 0.0d0
! Values of the underlying where evaluations are made
      integer, parameter :: nv = 3
      real(kind(1d0)) :: eval_points(nv) = &
                                  (/19.0d0, 20.0d0, 21.0d0/)
! Define parameters for the integration step.
      integer, parameter :: nx = 121, nint = nx-1, n = 3*nx
      real(kind(1d0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt)
      type(d_fcn_data) fcn_data
! Number of left/right boundary conditions
      integer, parameter :: nlbc = 3, nrbc = 3
! Further parameters for the integration step
      real(kind(1d0)) :: dx, dx2, pi, sqrt2pi
! used to calc derivatives
      real(kind(1d0)) :: epsfrac = .001d0
      character(len=6) :: greek_name(12) = (/&
             " Value", " Delta", " Gamma", " Theta",&
             " Charm", " Color", "  Vega", " Volga",&
             " Vanna", "   Rho", " Speed", "Speed2"/)
! Time values for the options
      real(kind(1d0)) :: rex(12), reavg(12)
      integer :: irect(12)
      integer :: i, i2, i3, j, ig, iquart, nout

      real(kind(1d0)) ::&
                  spline_values(nv,nt,12), spline_values1(nv,nt),&
                  spline_valuesp(nv,nt), spline_valuesm(nv,nt),&
                  spline_valuespp(nv,nt), spline_valuesmm(nv,nt),&
                  xgridp(nx), xgridm(nx),xgridpp(nx), xgridmm(nx),&
                  eval_pointsp(nv), eval_pointspp(nv),&
                  eval_pointsm(nv), eval_pointsmm(nv),&
                  BS_values(nv,nt,12), sVo_array(nt)

      call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
      allocate(fcn_data % rdata (6))

      pi = const('pi')
      sqrt2pi = sqrt(2.0 * pi)
      dx2 = epsfrac * 0.5d0 * (x_min + x_max)

! Compute Constant Elasticity of Variance Model for Vanilla Call
! Define 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

      write(nout,'(T05,A)') "Constant Elasticity of Variance"//&
                            " Model for Vanilla Call Option"
      write(nout,'(T10,"Interest Rate: ", F7.3, T38,'//&
                   '"Continuous Dividend: ", F7.3  )') r, dividend
      write(nout,'(T10,"Minimum and Maximum Prices of '//&
                   ' Underlying: ", 2F7.2)') x_min, x_max
      write(nout,'(T10,"Number of equally spaced spline knots:",'//&
                   'I4)') nx - 1
      write(nout,'(T10,"Number of unknowns: ",I4)') n
      write(nout,*)
      write(nout,'(/T10,"Time in Years Prior to Expiration: ",'//&
                   '2X,3F7.4)') time
      write(nout,'(T10,"Option valued at Underlying Prices:  ",'//&
                   '3F7.2)') eval_points
      write(nout,*)

!
! iquart = 0 : derivatives estimated with 3-point fitted parabola
! iquart = 1 : derivatives estimated with 5-point fitted quartic
!              polynomial
!

      iquart = 0
      if (iquart == 0) then
        write(nout,'(T10,"3 point (parabolic) estimate of '//&
                    'parameter derivatives")')
      else
        write(nout,'(T10,"5 point (quartic) estimate of parameter'//&
                  ' derivatives")')
      end if
      write(nout, '(T10,"epsfrac = ", F11.8)') epsfrac

      irect = 0
      reavg = 0.0d0
      rex = 0.0d0
    
! alpha:  Black-Scholes
      do i2 = 1, 3
! Loop over volatility
        do i3 = 1, 3
! Loop over strike price
          call calc_Greeks(i2, i3, iquart)
        end do
      end do

      write(nout,*)
      do ig = 1, 12
        reavg(ig) = reavg(ig)/irect(ig)
        write (nout, '(" Greek: ", A6, ";  avg rel err: ",'//&
                     'F15.12, ";  max rel err: ", F15.12)')&
                     greek_name(ig), reavg(ig), rex(ig)
      end do

      contains

      subroutine calc_Greeks(volatility, strike_price, iquart)
      implicit none
      integer, intent(in) :: volatility, strike_price, iquart
    ! Local variables
      integer :: i1 = 1, j, iSderiv, gNi, l, k
      integer :: nt = 3
      real(kind(1d0)) :: x_maxp, x_maxm, x_maxpp, x_maxmm
      real(kind(1d0)) :: eps, tau, sigsqtau, sqrt_sigsqtau, sigsq
      real(kind(1d0)) :: d1, d2, n01pdf_d1, nu, relerr, relerrmax
      real(kind(1d0)) :: sVo, BSVo, S

      if ((volatility == 1) .and. (strike_price == 1)) then
        write(nout,*)
        write(nout,'(/T10,"Strike = ",F5.2,", Sigma =", F5.2,'//&
                    '", Alpha =",F5.2,":")') ks(strike_price),&
                    sigma(volatility), alpha(i1)
        write(nout,*)
        write(nout,'(T10,"years to expiration: ", 3F7.4)')&
                   (time(j), j=1,3)
        write(nout,*)
      end if

      fcn_data % rdata = (/ks(strike_price), x_max,&
              sigma(volatility), alpha(i1), r, dividend/)
                           
      call feynman_kac(xgrid, time, nlbc, nrbc, fkcoef,&
                       fkinitcond, fkbc, y, yprime,&
                       FCN_DATA = fcn_data)
       
! Compute Value, Delta, Gamma, Theta, Charm and Color
      do l = 0,2
        do i=1,NT
          spline_values(:,i,l+1) = hqsval(eval_points, xgrid,&
                                          y(:,i), IDERIV=l)
          spline_values(:,i,l+4) = hqsval(eval_points, xgrid,&
                                          yprime(:,i), IDERIV=l)
        end do
      end do
! Signs for Charm and Color must be inverted because FEYNMAN_KAC
! computes -d/dt instead of d/dt
      spline_values(:,:,5:6) = -spline_values(:,:,5:6)
              
! Speed2
      do i=1,nt
        spline_values(:,i,12) = hqsval(eval_points, xgrid, Y(:,i),&
                                       IDERIV=3)
      end do
       
! Compute Vega, Volga, Vanna, Rho, Speed
!     l =   7     8      9     10    11
      do l = 7,11
        xgridp = xgrid
        xgridm = xgrid
        eval_pointsp = eval_points
        eval_pointsm = eval_points
        x_maxp = x_max
        x_maxm = x_max
        fcn_data % rdata(3) = sigma(volatility)
        fcn_data % rdata(5) = r
        iSderiv = 0
        if (l == 9) iSderiv = 1 ! Vanna
        if (l == 11) iSderiv = 2 ! Speed
        if (l == 10) then
          fcn_data % rdata(5) = r * (1.0 + epsfrac) ! Rho
        else if (l < 10) then
          fcn_data % rdata(3) = sigma(volatility) * (1.0 + epsfrac)
        end if
          if (l == 11) then
             xgridp = xgrid + dx2
             xgridm = xgrid - dx2
             eval_pointsp = eval_points + dx2
             eval_pointsm = eval_points - dx2
             x_maxp = x_max + dx2
             x_maxm = x_max - dx2
          end if
          fcn_data % rdata(2) = x_maxp
          call feynman_kac(xgridp, time, nlbc, nrbc, fkcoef,&
                      fkinitcond, fkbc, y, yprime,&
                      FCN_DATA = fcn_data)
          do i=1,nt
            spline_valuesp(:,i) = hqsval(eval_pointsp, xgridp,&
                                         y(:,i), IDERIV=iSderiv)
          end do
          
          if (l == 10) then
             fcn_data % rdata(5) = r * (1.0 - epsfrac) ! Rho
          else if (l < 10) then
             fcn_data % rdata(3) = sigma(volatility) *&
                                   (1.0 - epsfrac)
          end if
          fcn_data % rdata(2) = x_maxm
! calculate spline values for sigmaM = sigmai2-1*(1. - epsfrac) or
! rM = r*(1. - epsfrac):
          call feynman_kac(xgridm, time, nlbc, nrbc, fkcoef,&
                           fkinitcond, fkbc, y, yprime,&
                           FCN_DATA = fcn_data)
          do i=1,NT
            spline_valuesm(:,i) = hqsval(eval_pointsm, xgridm,&
                                         y(:,i), IDERIV=iSderiv)
          end do

          if (iquart == 1) then
            xgridpp = xgrid
            xgridmm = xgrid
            eval_pointspp = eval_points
            eval_pointsmm = eval_points
            x_maxpp = x_max
            x_maxmm = x_max
              
            if (l == 11) then  ! Speed
              xgridpp = xgrid + 2.0 * dx2
              xgridmm = xgrid - 2.0 * dx2
              eval_pointspp = eval_points + 2.0 * dx2
              eval_pointsmm = eval_points - 2.0 * dx2
              x_maxpp = x_max + 2.0 * dx2
              x_maxmm = x_max - 2.0 * dx2
            end if      

            fcn_data % rdata(2) = x_maxpp
            if (l == 10) then
! calculate spline values for rPP = r*(1. + 2.*epsfrac):
              fcn_data % rdata(5) = r * (1.0 + 2.0 * epsfrac)
            else if (l < 10) then
! calculate spline values for sigmaPP = sigma(i2)*(1. + 2.*epsfrac):
              fcn_data % rdata(3) = sigma(volatility) *& 
                                     (1.0 + 2.0 * epsfrac)
            end if
            call feynman_kac (xgridpp, time, nlbc, nrbc, fkcoef,&
                              fkinitcond, fkbc, y, yprime,&
                              FCN_DATA = fcn_data)

            do i=1,nt
              spline_valuespp(:,i) = hqsval(eval_pointspp, xgridpp,&
                                            Y(:,i), IDERIV=iSderiv)
            end do

            fcn_data % rdata(2) = x_maxmm
! calculate spline values for sigmaMM = sigmai2-1*(1. - 2.*epsfrac)
! or rMM = r*(1. - 2.*epsfrac)
            if (l == 10) then
              fcn_data % rdata(5) = r * (1.0 - 2.0 * epsfrac)
            else if (l < 10) then
              fcn_data % rdata(3) = sigma(volatility) *& 
                                    (1.0 - 2.0 * epsfrac)
            end if
            call feynman_kac (xgridmm, time, nlbc, nrbc, fkcoef,&
                              fkinitcond, fkbc, y, yprime,&
                              FCN_DATA = fcn_data)

            do i=1,nt
              spline_valuesmm(:,i) = hqsval(eval_pointsmm, xgridmm,&
                                            y(:,i), IDERIV=iSderiv)
            end do
          end if ! if (iquart == 1)
          
          if (l /= 8) then
            eps = sigma(volatility) * epsfrac
            if (l == 10) eps = r * epsfrac ! Rho
            if (l == 11) eps = dx2 ! Speed

            spline_values(:,:,l) = &
               (spline_valuesp - spline_valuesm) / (2.0 * eps)
            if (iquart /= 0)  then
              spline_values1 =&
                (spline_valuespp - spline_valuesmm) /(4.0 * eps)
              spline_values(:,:,l)  =&
                (4.0 * spline_values(:,:,l) - spline_values1) / 3.0
            end if
          end if

          if (l == 8) then ! Volga
            eps = sigma(volatility) * epsfrac                        
            spline_values(:,:,l) =&
                 (spline_valuesp + spline_valuesm - 2.0 * &
                  spline_values(:,:,1)) / (eps * eps)
                        
             if (iquart /= 0) then
               spline_values1 =&
                 (spline_valuespp + spline_valuesmm - 2.0 * & 
                   spline_values(:,:,1)) / (4.0 * eps * eps)
               spline_values(:,:,l) = &
                (4.0 * spline_values(:,:,l) - spline_values1) / 3.0
          end if
        end if
      end do
! Evaluate BS solution at vector eval_points, 
! at each time value prior to expiration.

      do i = 1,nt
!
!  Black-Scholes (BS) European call option
!  value = ValBSEC(S,t) = exp(-q*tau)*S*N01CDF(d1) -
!                         exp(-r*tau)*K*N01CDF(d2),
!  where:
!    tau = time to maturity = T - t;
!    q = annual dividend yield;
!    r = risk free rate;
!    K = strike price;
!    S = stock price;
!    N01CDF(x) = N(0,1)_CDF(x);
!    d1 = ( log( S/K ) +
!         ( r - q + 0.5*sigma**2 )*tau ) /
!         ( sigma*sqrt(tau) );
!    d2 = d1 - sigma*sqrt(tau)
!
! BS option values for tau = time(i):
        tau = time(i)
        sigsqtau = (sigma(volatility)**2) * tau
        sqrt_sigsqtau = sqrt(sigsqtau)
        sigsq = sigma(volatility) * sigma(volatility)
        do j = 1, nv
! Values of the underlying price where evaluations are made:
          S = eval_points(j)
          d1 = (log(S / ks(strike_price)) + (r - dividend)& 
                * tau + 0.5 * sigsqtau) / sqrt_sigsqtau
          n01pdf_d1 = exp((-0.5) * d1 * d1) / sqrt2pi
          nu = exp((-dividend) * tau) * S * n01pdf_d1 * sqrt(tau)
                          
          d2 = d1 - sqrt_sigsqtau
          BS_values(j,i,1) = exp((-dividend) * tau) * S *& 
                                  anordf(d1) - exp((-r) * tau) *&
                                  ks(strike_price) * anordf(d2)
! greek = Rho
          BS_values(j,i,10) = exp((-r) * tau) * ks(strike_price) *&
                                    tau * anordf(d2)
! greek = Vega
          BS_values(j,i,7) = nu
! greek = Volga
          BS_values(j,i,8) = nu * d1 * d2 / sigma(volatility)
! greek = delta
          BS_values(j,i,2) = exp((-dividend) * tau) * anordf(d1)
! greek = Vanna
          BS_values(j,i,9) = (nu / S) * (1.0 - d1 / sqrt_sigsqtau)
! greek = dgamma
          BS_values(j,i,3) = exp((-dividend) * tau) *& 
                                 n01pdf_d1 /(S * sqrt_sigsqtau)
! greek = speed
          BS_values(j,i,11) = (-exp((-dividend) * tau)) *&
                         n01pdf_d1 * (1.0 + d1 / sqrt_sigsqtau)&
                         / (S * S * sqrt_sigsqtau)
! greek = speed
          BS_values(j,i,12) = (-exp((-dividend) * tau)) * & 
                       n01pdf_d1 * (1.0 + d1 / sqrt_sigsqtau) / & 
                       (S * S * sqrt_sigsqtau)
          d2 = d1 - sqrt_sigsqtau
! greek = theta
          BS_values(j,i,4) = exp((-dividend) * tau) * S * & 
                          (dividend * anordf(d1) - 0.5 * sigsq * &
                           n01pdf_d1 / sqrt_sigsqtau) - r * &
                           exp((-r) * tau) * ks(strike_price) * & 
                           anordf(d2)
! greek = charm
          BS_values(j,i,5) = exp((-dividend) * tau) * ((-dividend)&
                           * anordf(d1) + n01pdf_d1 *&
                          ((r - dividend) * tau - 0.5 * d2 *&
                           sqrt_sigsqtau) / (tau * sqrt_sigsqtau))
! greek = color
          BS_values(j,i,6) = & 
                          (-exp((-dividend) * tau)) * n01pdf_d1 *&
                          (2.0 * dividend * tau + 1.0 + d1 *&
                          (2.0 * (r - dividend) * tau - d2 *&
                           sqrt_sigsqtau) / sqrt_sigsqtau) / &
                          (2.0 * S * tau * sqrt_sigsqtau)
        end do
      end do

      do l=1,12
        relerrmax = 0.0
        do i = 1,nv
          do j = 1,nt
            sVo = spline_values(i,j,l) 
            BSVo = BS_values(i,j,l)
            relerr = abs((sVo - BSVo) / BSVo)
            if (relerr > relerrmax) relerrmax = relerr
            reavg(l) = reavg(l) + relerr
            irect(l) = irect(l) + 1
          end do
        end do
        if (relerrmax > rex(l)) rex(l) = relerrmax

        if ((volatility == 1) .and. (strike_price == 1)) then
          do i=1,nv
            sVo_array(1:nt) = spline_values(i,1:nt,l)
            write(nout,'("underlying price: ", F4.1,"; FK ",'//&
                  'A6,": ", 3(F13.10,TR1))') eval_points(i),&
                  greek_name(l),&
                  (sVo_array(k), k=1,nt)
            write(nout, '(T25, "BS ", A6,": ", 3(F13.10,TR1))')&
                  greek_name(l), (BS_values(i,k,l), k=1,nt)
           end do
        end if
      end do
      end subroutine calc_Greeks
      end

! These functions and routines  define the coefficients, payoff and boundary conditions.
      function fkcoef (x, tx, iflag, fcn_data)
      use mp_types
      implicit none
      real(kind(1d0)), intent(in) :: x, tx
      integer, intent(inout) :: iflag
      type(d_fcn_data), optional :: fcn_data
      real(kind(1d0)) :: fkcoef

      real(kind(1d0)) :: sigma, interest_rate, alpha, dividend,&
                         half = 0.5d0
      sigma = fcn_data % rdata(3)
      alpha = fcn_data % rdata(4)
      interest_rate = fcn_data % rdata(5)
      dividend = fcn_data % rdata(6)
      select case (iflag)
        case (1)
! The coefficient derivative d(sigma)/dx
          fkcoef = half*alpha*sigma*x**(alpha*half-1.0d0)
! The coefficient sigma(x)	
        case (2)
          fkcoef = sigma*x**(alpha*half) 
        case (3)
! The coefficient mu(x)
          fkcoef = (interest_rate - dividend) * x
        case (4)
! The coefficient kappa(x)
          fkcoef = interest_rate
      end select
! Note that there is no time dependence
      iflag = 0
      return 
      end function fkcoef

      function fkinitcond(x, fcn_data)
      use mp_types
      implicit none
      real(kind(1d0)), intent(in) :: x
      type (d_fcn_data), optional :: fcn_data
      real(kind(1d0)) :: fkinitcond
      real(kind(1d0)) :: zero = 0.0d0
      real(kind(1d0)) :: strike_price

      strike_price = fcn_data % rdata(1)
! The payoff function
      fkinitcond = max(x - strike_price, zero)
      return
      end function fkinitcond

      subroutine fkbc (tx, iflag, bccoefs, fcn_data)
      use mp_types
      implicit none
      real(kind(1d0)), intent(in) :: tx
      integer, intent(inout) :: iflag
      real(kind(1d0)), dimension(:,:), intent(out) :: bccoefs
      type (d_fcn_data), optional :: fcn_data
      real(kind(1d0)) :: x_max, df, interest_rate, strike_price

      strike_price = fcn_data % rdata(1)
      x_max = fcn_data % rdata(2)
      interest_rate = fcn_data % rdata(5)
      select case (iflag)
        case (1)
          bccoefs(1,1:4) = (/1.0d0, 0.0d0, 0.0d0, 0.0d0/)
          bccoefs(2,1:4) = (/0.0d0, 1.0d0, 0.0d0, 0.0d0/)
          bccoefs(3,1:4) = (/0.0d0, 0.0d0, 1.0d0, 0.0d0/)
! Note no time dependence at left end
          iflag = 0
        case (2)
          df = exp(interest_rate * tx)
          bccoefs(1,1:4) = (/1.0d0, 0.0d0, 0.0d0, x_max -&
                             df*strike_price/)   
          bccoefs(2,1:4) = (/0.0d0, 1.0d0, 0.0d0, 1.0d0/)
          bccoefs(3,1:4) = (/0.0d0, 0.0d0, 1.0d0, 0.0d0/)
      end select
      end subroutine fkbc
