Chapter 5: Differential Equations

FEYNMAN_KAC

TextonSpedometerHPClogo.gif



Solves the generalized Feynman-Kac PDE on a rectangular grid using a finite element Galerkin method. Initial and boundary conditions are provided. The solution is represented by a series of C2 Hermite quintic splines.

Required Arguments

XGRID — Rank-1 array containing the set of breakpoints that define the end points for the Hermite quintic splines.   (Input)
Let m = size(XGRID). The points in XGRID must be in strictly increasing order, and m  2.

TGRID — Rank-1 array containing the set of time points (in time-remaining units) at which an approximate solution is computed.   (Input)
Let n = size(TGRID). The points in TGRID must be strictly positive and in strictly increasing order and n 1.

NLBC — The number of left boundary conditions.   (Input)
1≤  NLBC ≤ 3.

NRBC — The number of right boundary conditions.   (Input)
1≤ NRBC ≤ 3.

FKCOEF — User-supplied FUNCTION to evaluate the coefficients and of the Feynman-Kac PDE. The usage is FKCOEF (X, TX, IFLAG[]), where

Function Return Value

FKCOEF — Value of the coefficient function. Which value is computed depends on the input value for IFLAG, see description of IFLAG.

Required Arguments

X — Point in the x-space at which the coefficient is to be evaluated.   (Input)

TX — Time point at which the coefficient is to be evaluated.   (Input)

IFLAG — Flag related to the coefficient that has to be computed   (Input/Output).
On entry, IFLAG indicates which coefficient is to be computed. The following table shows which value has to be returned by FKCOEF for all possible values of IFLAG:

 

IFLAG

Computed coefficient

1



2



3



4



 

                One indicates when a coefficient does not depend on t by setting IFLAG = 0 after the coefficient is defined. If there is time dependence, the value of IFLAG should not be changed. This action will usually yield a more efficient algorithm because some finite element matrices do not have to be reassembled for each t value.

Optional Arguments

FCN_DATA —  A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied function. For a detailed description of this argument see FCN_DATA below.   (Input/Output)

     FKCOEF must be declared EXTERNAL in the calling program.

FKINITCOND — User-supplied FUNCTION to evaluate the initial condition function  in the Feynman-Kac PDE. The usage is FKINITCOND (X[,…]), where

Function Return Value

FKINITCOND — Value of the initial condition function .

Required Arguments

X — Point in the x-space at which the initial condition is to be evaluated.   (Input)

Optional Arguments

FCN_DATA —  A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied function. For a detailed description of this argument see FCN_DATA below.   (Input/Output)

     FKINITCOND must be declared EXTERNAL in the calling program.

FKBC — User-supplied subroutine to evaluate the coefficients for the left and right boundary conditions the Feynman-Kac PDE must satisfy. There are NLBC conditions specified at the left end, , and NRBC conditions at the right end, . The boundary conditions can be vectors of dimension 1, 2 or 3 and are defined by

            The usage is FKBC (TX, IFLAG, BCCOEFS[]) where

Required Arguments

TX — Time point at which the coefficients are to be evaluated. (Input)

IFLAG — Flag related to the boundary conditions that have to be computed (Input/Output).
On input, IFLAG indicates whether the coefficients for the left or right boundary conditions have to be computed:

 

IFLAG

Computed boundary conditions

1

Left end, x = xmin

2

Right end, x = xmax

 

If there is no time dependence for one of the boundaries then set IFLAG = 0 after the array BCCOEFS is defined for either end point.  This can avoid unneeded continued computation of the finite element matrices.

BCCOEFS — Array of size 3 Χ 4 containing the coefficients of the left or right boundary conditions in its first NLBC or NRBC rows, respectively.   (Output)
The coefficients for  are stored row-wise according to the following matrix-scheme:

The coefficients for are stored similarly.

Optional Arguments

FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a detailed description of this argument see FCN_DATA below.   (Input/Output)

     FKBC must be declared EXTERNAL in the calling program. 

Y — Array of size (3*m) by (n+1) containing the coefficients of the Hermite representation of the approximate solution for the Feynman-Kac PDE at time points (in time-remaining units) 0, TGRID(1), …, TGRID(n). (Output)
For , the coefficients are stored in columns 1,…,n of array Y and  the  approximate solution is given by

.

            The coefficients of the representation for the initial data are given in column 0 of array Y and are defined by

.

            The starting coefficients  are estimated using least-squares.

            After the integrations, use Y(:,0) and Y(:,j) as input argument COEFFS to function HQSVAL to obtain an array of values for f(x, t) or its partials  at time points t=0 and t=TGRID(j), j=1,…,n, respectively.

            The expressions for the basis functions are represented piece-wise and can be found in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”.

YPRIME — Array of size (3*m) by (n + 1) containing the first derivatives of the coefficients of the Hermite representation of the approximate solution for the Feynman-Kac PDE at time points (in time-remaining units) 0, TGRID(1), …, TGRID(n).  (Output)
For  and , the derivatives of the coefficients are stored in column 0 and columns 1 to n of array YPRIME, respectively. The  columns in YPRIME represent

  for ,

            and

  for .

            After the integrations, use YPRIME(:,j) as input argument COEFFS to function HQSVAL to obtain an array of values for the partials  at time points t = TGRID(j), j = 1,…,n, and YPRIME(:,0) for the partials at t = 0.

Optional Arguments

FKINIT — User-supplied subroutine that allows for adjustment of initial data or as an opportunity for output during the integration steps.

            The usage is CALL FKINIT (XGRID, TGRID, TX, YPRIME, Y, ATOL, RTOL[,…]) where

Required Arguments

XGRID — Array of size m containing the set of breakpoints that define the end points for the Hermite quintic splines.   (Input)

TGRID — Array of size n containing the set of time points (in time-remaining units) at which an approximate solution is computed.   (Input)

TX — Time point for the evaluation.   (Input)
Possible values are 0 (the initial or “terminal” time point) and all values in array TGRID.

YPRIME — Array of length 3*m containing the derivatives of the Hermite quintic spline coefficients at time point TX.   (Input)

Y — Array of length 3* m containing the Hermite quintic spline coefficients at time point TX.   (Input/Output)
For the initial time point TX=0 this array can be used to reset the Hermite quintic spline coefficients to user defined values. For all other values of TX array Y is an input array.

ATOL — Array of length 3* m containing absolute error tolerances used in the integration routine that determines the Hermite quintic spline coefficients and its derivatives.   (Input/Output)

RTOL — Array of length 3*m containing relative error tolerances used in the integration routine that determines the Hermite quintic spline coefficients and its derivatives.   (Input/Output)

Optional Arguments

FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied function. For a detailed description of this argument see FCN_DATA below.   (Input/Output)

     FKINIT must be declared EXTERNAL in the calling program.

FKFORCE — User-supplied subroutine that computes local contributions

 and .

            The usage is CALL FKFORCE (I, T, WIDTH, Y, XLOCAL, QW, U, PHI, DPHI[,]) where

Required Arguments

I — Index related to the integration interval (XGRID(I), XGRID(I+1)).   (Input)

T — Time point at which the local contributions are computed.   (Input)

WIDTH — Width of the integration interval I, WIDTH= XGRID(I+1)- XGRID(I).   (Input)

Y — Array of length 3*m containing the coefficients of the Hermite quintic spline representing the solution of the Feynman-Kac PDE at time point T.   (Input)
For each

,

the approximate solution is locally defined by

Here, the functions  are basis polynomials of order 5 and

.

The values

,

are stored as successive triplets in array Y.

XLOCAL — Array containing the Gauss-Legendre points translated and normalized to the interval [XGRID(I),XGRID(I+1)].   (Input)
The size of the array is equal to the degree of the Gauss-Legendre polynomials used for constructing the finite element matrices.

QW — Array containing the Gauss-Legendre weights.   (Input)
The size of the array is equal to the degree of the Gauss-Legendre polynomials used for constructing the finite element matrices.

U — Array of size size(XLOCAL) Χ 12 containing the basis function values that define at the Gauss-Legendre points XLOCAL.   (Input)
Let

.

Using the local approximation in the I-th interval, defined by

,

and setting

,  and ,

vector  is defined as

PHI — Array of size 6 containing a Gauss-Legendre approximation for the local contribution, where T and .   (Output)

Setting NDEG:=SIZE(XLOCAL) and , array PHI contains elements

PHI(i) = WIDTH

for i= 1,…, 6.

DPHI — Array of size 6Χ6, a Gauss-Legendre approximation for the Jacobian of the local contribution  at time point = T,

.   (Output)

The approximation to this symmetric matrix is stored row-wise, i.e.

DPHI(i,j) = WIDTH

or i, j = 1,…,6.

Optional Arguments

FCN_DATA —  A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a detailed description of this argument see FCN_DATA below.   (Input/Output)

     FKFORCE must be declared EXTERNAL in the calling program.

            If subroutine FKFORCE is not used as an optional argument then it is assumed that the forcing term  in the Feynman-Kac equation is identically zero.

ATOL — Array of non-negative values containing absolute error tolerances used in the computation of each column of solution array Y via integration routine DASPH.   (Input)
The size of array ATOL can be 1 or 3Χm. In the first case, ATOL(1:1) is applied to all solution components, in the latter each component of ATOL is assigned to the corresponding solution component allowing for individual control of the error tolerances. At least one entry in arrays ATOL or RTOL must be greater than 0.
Default: ATOL(1:1) = 1.0e-3 for single and 1.0d-5 for double precision.

RTOL — Array of non-negative values containing relative error tolerances used in the computation of each column of solution array Y via integration routine DASPH.   (Input)
The size of array RTOL can be 1 or 3Χm. In the first case, RTOL(1:1) is applied to all solution components, in the latter each component of RTOL is assigned to the corresponding solution component allowing for individual control of the error tolerances. At least one entry in arrays ATOL or RTOL must be greater than 0.
Default: RTOL(1:1) = 1.0e-3 for single and 1.0d-5 for double precision.

NDEG — Degree of the Gauss-Legendre formulas used for constructing the finite element matrices.   (Input)
NDEG ≥ 6.
Default: NDEG = 6.

RINITSTEPSIZE — Starting step size for the integration.   (Input)
RINITSTEPSIZE must be strictly negative since the integration is internally done from T = 0 to T = TGRID(n) in a negative direction.
Default: Program defined initial stepsize.

MAXBDFORDER — Maximum order of the backward differentiation formulas (BDF) used in the integrator DASPH.   (Input)
1 MAXBDFORDER 5.
Default: MAXBDFORDER = 5.

RMAXSTEPSIZE — Maximum step size the integrator may take.   (Input)
RMAXSTEPSIZE must be strictly positive.
Default: RMAXSTEPSIZE = AMACH(2), the largest possible machine number.

MAXIT — Maximum number of internal integration steps between two consecutive time points in TGRID.   (Input)
MAXIT must be strictly positive.
Default: MAXIT = 500000.

IMETHSTEPCTRL — Indicates which step control algorithm is used in the integration.   (Input)
If IMETHSTEPCTRL = 0, then the step control method of Sφderlind is used. If  IMETHSTEPCTRL =1, then the method used by the original Petzold code SASSL is used.

 

IMETHSTEPCTRL

Method used

0

Method of Sφderlind..

1

Method from Petzold code SASSL.

 

Default: IMETHSTEPCTRL = 0.

TBARRIER — Time barrier past which the integration routine DASPH will not go during integration.   (Input)
TBARRIER ≥ TGRID(n).
Default: TBARRIER = TGRID(n).

ISTATE — Array of size 5 whose entries flag the state of computation for the matrices and vectors required in the integration.  (Output)
For each entry, a zero indicates that no computation has been done or that there is a time dependence. A one indicates that the entry has been computed and there is no time dependence.
The ISTATE entries are as follows:

 

I

ISTATE(I)

1

State of computation of Mass matrix, M.

2

State of computation of Stiffness matrix, N.

3

State of computation of Bending matrix, R.

4

State of computation of Weighted mass matrix, K.

5

State of computation of initial data.

 

NVAL — Array of size 3 summarizing the number of evaluations required during the integration. (Output)

 

I

NVAL(I)

1

Number of residual function evaluations
of the DAE used in the model.

2

Number of factorizations of the differential
matrix associated with solving the DAE.

3

Number of linear system solve steps using
the differential matrix.

 

ITDEPEND — Logical array of size 7 indicating time dependence of the coefficients, boundary conditions and forcing term in the Feynman-Kac equation.   (Output)
If ITDEPEND(I)=.FALSE. then argument I is not time dependent, if ITDEPEND(I) =.TRUE. then argument I is time dependent.

 

I

ITDEPEND(I)

1

Time dependence of .

2

Time dependence of .

3

Time dependence of .

4

Time dependence of .

5

Time dependence of left boundary conditions.

6

Time dependence of right boundary conditions.

7

Time dependence of .

 

FCN_DATA — A derived type, s_fcn_data, which may be used to pass additional information to/from the user-supplied function.   (Input/Output)
The derived type, s_fcn_data, is defined as:

type s_fcn_data

   real(kind(1e0)), pointer, dimension(:) :: rdata

   integer, pointer, dimension(:) :: idata

end type

            in module mp_types. The double precision counterpart to s_fcn_data is named d_fcn_data. The user must include a use mp_types statement in the calling program to define this derived type.

            Note that if user-supplied data are required in one of the user-defined functions or subroutines available for routine FEYNMAN_KAC then these data must be defined via FCN_DATA.

FORTRAN 90 Interface

Generic:          CALL FEYNMAN_KAC (XGRID, TGRID, NLBC, NRBC, FKCOEF, FKINITCOND, FKBC, Y, YPRIME [,…])

Specific:         The specific interface names are S_FEYNMAN_KAC and D_FEYNMAN_KAC.

Description

The generalized Feynman-Kac differential equation has the form

,

where the initial data satisfies

.

The derivatives are  etc.

FEYNMAN_KAC uses a finite element Galerkin method over the rectangle[3]

in  to compute the approximate solution. The interval  is decomposed with a grid

.

On each subinterval the solution is represented by

The values  are time-dependent coefficients associated with each interval. The basis functions are given for ,

by

The Galerkin principle is then applied. Using the provided initial and boundary conditions leads to an index 1 differential-algebraic equation (DAE) for the time-dependent coefficients

.

This system is integrated using the variable order, variable step algorithm DASPH. Solution values and their time derivatives are returned at a grid preceding time T, expressed in units of time remaining.

More mathematical details are found in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”.

Examples

Example 1 – A Diffusion Model For Call Options

In Beckers (1980) there is a model for a Stochastic Differential Equation of option pricing.  The idea is a “constant elasticity of variance diffusion (or CEV) class”

The Black-Scholes model is the limiting case.  A numerical solution of this diffusion model yields the price of a call option.  Various values of the strike price, time values, and power coefficientare used to evaluate the option price at values of the underlying price.  The sets of parameters in the computation are:

1.   power

2.   strike price

3.   volatility

4.   times until expiration

5.   underlying prices

6.   interest rate

7.  

8.  

With this model the Feynman-Kac differential equation is defined by identifying:

The payoff function is the “vanilla option”, .

Link to example source (feynman_kac_ex1.f90)

 

! Compute Constant Elasticity of Variance Model for Vanilla Call

      use feynman_kac_int

      use hqsval_int

      use mp_types

      use umach_int

 

      implicit none

 

! The set of strike prices

      real(kind(1e0)) :: ks(3)=(/15.0e0,20.0e0,25.0e0/)

! The set of sigma values

      real(kind(1e0)) :: sigma(3) = (/0.2e0, 0.3e0, 0.4e0/)

! The set of model diffusion powers

      real(kind(1e0)) :: alpha(3) = (/2.0e0,1.0e0,0.0e0/)

! Time values for the options

      integer, parameter :: nt = 3

      real(kind(1e0)) :: time(nt)=(/1.e0/12., 4.e0/12., 7.e0/12./)

! Values of the underlying where evaluation are made

      integer, parameter :: nv = 3, nlbc = 3, nrbc = 3

      real(kind(1e0)) :: xs(nv) = (/19.0e0,20.0e0,21.0e0/)

! Value of the interest rate and continuous dividend

      real(kind(1e0)) :: r = 0.05e0, dividend = 0.0e0

! Values of the min and max underlying values modeled

      real(kind(1e0)) :: x_min = 0.0e0, x_max = 60.0e0

 

! Define parameters for the integration step.

      integer, parameter :: nx = 121, nint = nx-1, n = 3*nx

      real(kind(1e0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt),&

                         dx, f(nv,nt)

      type(s_fcn_data) fcn_data

      integer :: nout

      real(kind(1e0)), external :: fkcoef, fkinitcond              

      external fkbc     

 

      integer :: i,i1,i2,i3,j

! Allocate space inside the derived type for holding

! data values. These are for the evaluation routines.

      allocate(fcn_data % rdata (6))

! 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

 

      call umach(2, nout)

      write(nout,'(T05,A)') "Constant Elasticity of Variance Model "//&

                            "for Vanilla Call"

      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,'//&

                   '/T10,"Number of unknowns ",I4)')&

                   nx-1,n

      write(nout,'(/T10,"Time in Years Prior to Expiration ",2X,'//&

                   '3F7.4)') time

      write(nout,'(T10,"Option valued at Underlying Prices  ",'//&

                   '3F7.2)') xs

 

      do i1 = 1,3       ! Loop over power

        do i2=1,3     ! Loop over volatility

          do i3=1,3  ! Loop over strike price

! Pass data through into evaluation routines.      

            fcn_data % rdata =&

               (/ks(i3),x_max,sigma(i2),alpha(i1),r,dividend/)

            call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef,&

                              fkinitcond, fkbc, y, yprime,&

                              FCN_DATA = fcn_data)

! Evaluate solution at vector of points XS(:), at each time value

! prior to expiration.                       

            do i=1,nt

              f(:,i) = hqsval (xs, xgrid, y(:,i))

            end do

            write(nout,'(/T05,"Strike=",F5.2," Sigma=", F5.2,'//&

              '" Alpha=", F5.2,/(T25," Call Option Values ",'//&

              'X,3F7.4))') ks(I3),sigma(I2),&

              alpha(i1),(f(i,:),i=1,nv)

          end do !i3 - Strike price loop

        end do !i2 - Sigma loop

      end do !i1 - Alpha loop

      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(1e0)), intent(in) :: X, TX

        integer, intent(inout) :: iflag

         type(s_fcn_data), optional :: fcn_data

         real(kind(1e0)) :: fkcoef

         real(kind(1e0)) :: sigma, interest_rate, alpha, dividend,&

                            zero = 0.0e0, half = 0.5e0

                  

        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.0e0)

! 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(1e0)), intent(in) :: x

        type (s_fcn_data), optional :: fcn_data

        real(kind(1e0)) :: fkinitcond

        real(kind(1e0)) :: zero = 0.0e0

        real(kind(1e0)) :: 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(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)) :: 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.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/)

! Note no time dependence at left end 

            iflag = 0

          case (2)

            df = exp(interest_rate * tx)

            bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,&

                             x_max - df*strike_price/)  

            bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 1.0e0/)

            bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)

        end select

      end subroutine fkbc

Output

    Constant Elasticity of Variance Model for Vanilla Call

         Interest Rate   0.050       Continuous Dividend   0.000

         Minimum and Maximum Prices of Underlying    0.00  60.00

         Number of equally spaced spline knots  120

         Number of unknowns  363

 

         Time in Years Prior to Expiration    0.0833 0.3333 0.5833

         Option valued at Underlying Prices    19.00  20.00  21.00

 

    Strike=15.00 Sigma= 0.20 Alpha= 2.00

                         Call Option Values   4.0624 4.2575 4.4730

                         Call Option Values   5.0624 5.2506 5.4490

                         Call Option Values   6.0624 6.2486 6.4385

 

    Strike=20.00 Sigma= 0.20 Alpha= 2.00

                         Call Option Values   0.1310 0.5955 0.9699

                         Call Option Values   0.5018 1.0887 1.5101

                         Call Option Values   1.1977 1.7483 2.1752

 

    Strike=25.00 Sigma= 0.20 Alpha= 2.00

                         Call Option Values   0.0000 0.0112 0.0745

                         Call Option Values   0.0000 0.0372 0.1621

                         Call Option Values   0.0007 0.1027 0.3141

 

    Strike=15.00 Sigma= 0.30 Alpha= 2.00

                         Call Option Values   4.0637 4.3398 4.6622

                         Call Option Values   5.0626 5.2944 5.5786

                         Call Option Values   6.0624 6.2708 6.5240

 

    Strike=20.00 Sigma= 0.30 Alpha= 2.00

                         Call Option Values   0.3109 1.0276 1.5494

                         Call Option Values   0.7326 1.5424 2.1017

                         Call Option Values   1.3765 2.1690 2.7379

 

    Strike=25.00 Sigma= 0.30 Alpha= 2.00

                         Call Option Values   0.0006 0.1112 0.3543

                         Call Option Values   0.0038 0.2169 0.5548

                         Call Option Values   0.0184 0.3857 0.8222

 

    Strike=15.00 Sigma= 0.40 Alpha= 2.00

                         Call Option Values   4.0755 4.5138 4.9675

                         Call Option Values   5.0662 5.4201 5.8326

                         Call Option Values   6.0634 6.3579 6.7301

 

    Strike=20.00 Sigma= 0.40 Alpha= 2.00

                         Call Option Values   0.5115 1.4640 2.1273

                         Call Option Values   0.9621 1.9951 2.6929

                         Call Option Values   1.5814 2.6105 3.3216

 

    Strike=25.00 Sigma= 0.40 Alpha= 2.00

                         Call Option Values   0.0083 0.3286 0.7790

                         Call Option Values   0.0285 0.5167 1.0657

                         Call Option Values   0.0813 0.7687 1.4103

 

    Strike=15.00 Sigma= 0.20 Alpha= 1.00

                         Call Option Values   4.0624 4.2479 4.4311

                         Call Option Values   5.0624 5.2479 5.4311

                         Call Option Values   6.0624 6.2479 6.4311

 

    Strike=20.00 Sigma= 0.20 Alpha= 1.00

                         Call Option Values   0.0000 0.0218 0.1045

                         Call Option Values   0.1498 0.4109 0.6485

                         Call Option Values   1.0832 1.3314 1.5773

 

    Strike=25.00 Sigma= 0.20 Alpha= 1.00

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

 

    Strike=15.00 Sigma= 0.30 Alpha= 1.00

                         Call Option Values   4.0624 4.2477 4.4309

                         Call Option Values   5.0624 5.2477 5.4309

                         Call Option Values   6.0624 6.2477 6.4309

 

    Strike=20.00 Sigma= 0.30 Alpha= 1.00

                         Call Option Values   0.0011 0.0781 0.2201

                         Call Option Values   0.1994 0.5000 0.7543

                         Call Option Values   1.0835 1.3443 1.6023

 

    Strike=25.00 Sigma= 0.30 Alpha= 1.00

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0005

 

    Strike=15.00 Sigma= 0.40 Alpha= 1.00

                         Call Option Values   4.0624 4.2479 4.4312

                         Call Option Values   5.0624 5.2479 5.4312

                         Call Option Values   6.0624 6.2479 6.4312

 

    Strike=20.00 Sigma= 0.40 Alpha= 1.00

                         Call Option Values   0.0076 0.1563 0.3452

                         Call Option Values   0.2495 0.5907 0.8706

                         Call Option Values   1.0868 1.3779 1.6571

 

    Strike=25.00 Sigma= 0.40 Alpha= 1.00

                         Call Option Values   0.0000 0.0000 0.0001

                         Call Option Values   0.0000 0.0000 0.0008

                         Call Option Values   0.0000 0.0003 0.0063

 

    Strike=15.00 Sigma= 0.20 Alpha= 0.00

                         Call Option Values   4.0626 4.2479 4.4311

                         Call Option Values   5.0623 5.2480 5.4311

                         Call Option Values   6.0624 6.2480 6.4312

 

    Strike=20.00 Sigma= 0.20 Alpha= 0.00

                         Call Option Values   0.0001 0.0001 0.0002

                         Call Option Values   0.0816 0.3316 0.5748

                         Call Option Values   1.0818 1.3308 1.5748

 

    Strike=25.00 Sigma= 0.20 Alpha= 0.00

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

 

    Strike=15.00 Sigma= 0.30 Alpha= 0.00

                         Call Option Values   4.0625 4.2479 4.4312

                         Call Option Values   5.0623 5.2479 5.4312

                         Call Option Values   6.0624 6.2479 6.4312

 

    Strike=20.00 Sigma= 0.30 Alpha= 0.00

                         Call Option Values   0.0000 0.0000 0.0029

                         Call Option Values   0.0894 0.3326 0.5753

                         Call Option Values   1.0826 1.3306 1.5749

 

    Strike=25.00 Sigma= 0.30 Alpha= 0.00

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

 

    Strike=15.00 Sigma= 0.40 Alpha= 0.00

                         Call Option Values   4.0624 4.2479 4.4312

                         Call Option Values   5.0623 5.2479 5.4312

                         Call Option Values   6.0624 6.2479 6.4312

 

    Strike=20.00 Sigma= 0.40 Alpha= 0.00

                         Call Option Values   0.0000 0.0002 0.0113

                         Call Option Values   0.0985 0.3383 0.5781

                         Call Option Values   1.0830 1.3306 1.5749

 

    Strike=25.00 Sigma= 0.40 Alpha= 0.00

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

                         Call Option Values   0.0000 0.0000 0.0000

 

Example 2 – American Option vs. European Option On a Vanilla Put

The value of the American Option on a Vanilla Put can be no smaller than its European counterpart.  That is due to the American Option providing the opportunity to exercise at any time prior to expiration.  This example compares this difference – or premium value of the American Option – at two time values using the Black-Scholes model. The example is based on Wilmott et al. (1996, p. 176), and uses the non-linear forcing or weighting term described in Hanson, R. (2008), “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”, for evaluating the price of the American Option. A call to the subroutine fkinit_put sets the initial conditions. One breakpoint is set exactly at the strike price.

The sets of parameters in the computation are:

1.     Strike price

2.     Volatility

3.     Times until expiration

4.     Interest rate

5.    

6.    

The payoff function is the “vanilla option”, .

Link to example source (feynman_kac_ex2.f90)

 

! Compute American Option Premium for Vanilla Put

      use feynman_kac_int

      use hqsval_int

      use mp_types

      use umach_int

      implicit none

! The strike price

      real(kind(1e0)) :: ks = 10.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 evaluations are made

      integer, parameter :: nv = 9

      integer, parameter :: nlbc = 2, nrbc = 3, ndeg = 6

 

      integer :: i

      real(kind(1e0)) :: xs(nv) = (/((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

      real(kind(1e0)) :: atol(1), rtol(1)

 

! Define parameters for the integration step. 

      integer, parameter :: nx = 121, nint = nx-1, n = 3*nx

      real(kind(1e0)) :: xgrid(nx), ye(n,0:nt), yeprime(n,0:nt),&

                         ya(n,0:nt), yaprime(n,0:nt),&

                         dx, fe(nv,nt), fa(nv,nt)

      type(s_fcn_data) fcn_data

      integer :: nout

      real(kind(1e0)), external :: fkcoef_put, fkinitcond_put

      external fkbc_put, fkinit_put, fkforce_put

     

      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), 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

 

! Place a breakpoint at the strike price.

      do i = 1,nx

        if (xgrid(i) > ks) then

          xgrid(i-1) = ks

          exit

        end if

      end do

! Request less accuracy than the default values provide.

      atol(1) = 0.5e-2

      rtol(1) = 0.5e-2

      fcn_data % rdata = (/ks,x_max,sigma,r,dividend,atol(1)/)

      fcn_data % idata = (/ndeg/)

 

! Compute European then American Put Option Values.

      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&

                        fkinitcond_put, fkbc_put, ye, yeprime,&

                        FKINIT=fkinit_put, ATOL=atol,RTOL=rtol,&

                        FCN_DATA = fcn_data)

      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&

                        fkinitcond_put, fkbc_put, ya, yaprime,&

                        FKINIT=fkinit_put, ATOL=atol, RTOL=rtol,&

                        FKFORCE=fkforce_put, FCN_DATA = fcn_data)

 

! Evaluate solutions at vector of points XS(:), at each time value

! prior to expiration.                        

      do i=1,nt

        fe(:,i) = hqsval (xs, xgrid, ye(:,I))

        fa(:,I) = hqsval (xs, xgrid, ya(:,I))

      end do

      write(nout,'(T05,A,/,T05,A)')&

         "American Option Premium for Vanilla Put, 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,"Underlying", T26,"European",'//&

               'T42,"American",/(T10,5F8.4))') ks,sigma,r,&

               (xs(i), fe(i,1:nt), fa(i,1:nt),i=1,nv)

      end

! These routines define the coefficients, payoff, boundary

! conditions, forcing term and initial conditions for American and

! European Options.

      function fkcoef_put(x, tx, iflag, fcn_data)

        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)) :: fkcoef_put

 

        real(kind(1e0)) :: sigma, strike_price, interest_rate, &

                           dividend, zero=0.e0

        sigma = fcn_data % rdata(3)

        interest_rate = fcn_data % rdata(4)

        dividend = fcn_data % rdata(5)

        select case (iflag)

          case (1)

! The coefficient derivative d(sigma)/dx

            fkcoef_put = sigma

! The coefficient sigma(x)

          case (2)

            fkcoef_put = sigma*x

          case (3)

! The coefficient mu(x)

            fkcoef_put = (interest_rate - dividend)*x

          case (4)

! The coefficient kappa(x)

            fkcoef_put = interest_rate

        end select

! Note that there is no time dependence           

        iflag = 0

        return

      end function fkcoef_put

 

      function fkinitcond_put(x, fcn_data)

        use mp_types

        implicit none

        real(kind(1e0)), intent(in) :: x

        type (s_fcn_data), optional :: fcn_data

        real(kind(1e0)) :: fkinitcond_put

        real(kind(1e0)) :: zero = 0.0e0

        real(kind(1e0)) :: strike_price

 

        strike_price = fcn_data % rdata(1)

! The payoff function                

        fkinitcond_put = max(strike_price - x, zero)

        return

      end function fkinitcond_put

 

      subroutine fkbc_put (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

 

        select case (iflag)

          case (1)      

            bccoefs(1,1:4) = ((/0.0e0, 1.0e0, 0.0e0, -1.0e0/))

            bccoefs(2,1:4) = ((/0.0e0, 0.0e0, 1.0e0,  0.0e0/))

          case (2)

            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/))

        end select

! Note no time dependence

        iflag = 0             

      end subroutine fkbc_put

 

      subroutine fkforce_put (interval, t, hx, y, xlocal, qw, u,&

                              phi, dphi, fcn_data)

        use mp_types

        implicit none

        integer, parameter :: local = 6

        integer :: i, j, l, ndeg

        integer, intent(in) :: interval

        real(kind(1e0)), intent(in) :: y(:), t, hx, qw(:),&

                                       xlocal(:), u(:,:)

        real(kind(1e0)), intent(out) :: phi(:), dphi(:,:)

        type (s_fcn_data), optional :: fcn_data

 

        real(kind(1e0)) :: yl(local), bf(local)

        real(kind(1e0)) :: value, strike_price, interest_rate,&

                           zero=0.0e0, one=1.0e0, rt, mu

 

        yl = y(3*interval-2:3*interval+3)

        phi = zero

 

        value = fcn_data % rdata(6)

        strike_price = fcn_data % rdata(1)

        interest_rate = fcn_data % rdata(4)

        ndeg = fcn_data % idata(1)

 

        mu = 2

! This is the local definition of the forcing term

        do j=1,local

          do l=1,ndeg

            bf(1:3) = u(l,1:3)

            bf(4:6) = u(l,7:9)

            rt = dot_product(yl,bf)

            rt = value/(rt + value - (strike_price - xlocal(l)))

            phi(j) = phi(j) + qw(l) * bf(j) * rt**mu

          end do

        end do

        phi = -phi*hx*interest_rate*strike_price

 

! This is the local derivative matrix for the forcing term

 

        dphi = zero

        do j =1,local

          do i = 1,local

            do l=1,ndeg

              bf(1:3) = u(l,1:3)

              bf(4:6) = u(l,7:9)

              rt = dot_product(yl,bf)

              rt = one/(rt + value - (strike_price - xlocal(l)))

              dphi(i,j) = dphi(i,j) + qw(l) * bf(I) * bf(j) *&

                          rt**(mu+1)

            end do

          end do

        end do

        dphi = mu*dphi*hx*value**mu*interest_rate*strike_price

        return

      end subroutine fkforce_put

 

      subroutine fkinit_put(xgrid,tgrid,t,yprime,y,atol,rtol,&

                            fcn_data)

        use mp_types

        implicit none

        real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&

                                       yprime(:)

        real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)

        type (s_fcn_data), optional :: fcn_data

        integer :: i

 

        if (t == 0.0e0) then

! Set initial data precisely.  The strike price is a breakpoint.

! Average the derivative limit values from either side.

          do i=1,size(xgrid)

            if (xgrid(i) < fcn_data % rdata(1)) then

              y(3*i-2) = fcn_data % rdata(1) - xgrid(i)

              y(3*i-1) = -1.0e0

              y(3*i)= 0.0e0

            else if (xgrid(i) ==  fcn_data % rdata(1)) then

              y(3*i-2) = 0.0e0

              y(3*i-1) = -0.5e0

              y(3*i) = 0.0e0

            else

              y(3*i-2) = 0.0e0

              y(3*i-1) = 0.0e0

              y(3*i) = 0.0e0

            end if

          end do

        end if

      end subroutine fkinit_put

Output

    American Option Premium for Vanilla Put, 3 and 6 Months Prior to

    Expiry

       Number of equally spaced spline knots  121

       Number of unknowns  363

       Strike= 10.00, Sigma= 0.40, Interest Rate= 0.10

       Underlying        European        American

           0.0000  9.7536  9.5137 10.0000 10.0000

           2.0000  7.7536  7.5138  8.0000  8.0000

           4.0000  5.7537  5.5156  6.0000  6.0000

           6.0000  3.7614  3.5680  4.0000  4.0000

           8.0000  1.9064  1.9162  2.0214  2.0909

          10.0000  0.6516  0.8540  0.6767  0.9034

          12.0000  0.1625  0.3365  0.1675  0.3515

          14.0000  0.0369  0.1266  0.0374  0.1322

          16.0000  0.0088  0.0481  0.0086  0.0504

 

Example 3 – European Option With Two Payoff Strategies

This example evaluates the price of a European Option using two payoff strategies: Cash-or-Nothing and Vertical Spread.  In the first case the payoff function is

.

The value B is regarded as the bet on the asset price, see Wilmott et al. (1995, p. 39-40). The second case has the payoff function

Both problems use the same boundary conditions.  Each case requires a separate integration of the Black-Scholes differential equation, but only the payoff function evaluation differs in each case.  The sets of parameters in the computation are:

1.     Strike and bet prices K1={10.0}, K2 = {15.0}, and B = {2.0}

2.     Volatility σ = {0.4}.

3.     Times until expiration = {1/4, 1/2}.

4.     Interest rate r = 0.1.

5.     xmin = 0, xmax = 30.

6.     nx =121, n = 3 Χ nx = 363.

Link to example source (feynman_kac_ex3.f90)

 

! 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

Output

    European Option Value for A Bet

     and a Vertical Spread, 3 and 6 Months Prior to Expiry

       Number of equally spaced spline knots   61

       Number of unknowns  183

       Strike= 10.00, Sigma= 0.40, Interest Rate= 0.10

       Bet =  2.00, Spread Value = 15.00

 

         Underlying            A Bet   Vertical Spread

            2.0000   0.0000   0.0000   0.0000   0.0000

            4.0000   0.0000   0.0014   0.0000   0.0006

            6.0000   0.0110   0.0723   0.0039   0.0447

            8.0000   0.2691   0.4302   0.1478   0.3832

           10.0000   0.9948   0.9781   0.8909   1.1926

           12.0000   1.6094   1.4290   2.1911   2.2273

           14.0000   1.8655   1.6922   3.4254   3.1553

           16.0000   1.9338   1.8175   4.2263   3.8264

           18.0000   1.9476   1.8700   4.6264   4.2492

           20.0000   1.9501   1.8904   4.7911   4.4921

           22.0000   1.9505   1.8979   4.8497   4.6231

           24.0000   1.9506   1.9007   4.8685   4.6909

 

Example 4 – Convertible Bonds

This example evaluates the price of a convertible bond. Here, convertibility means that the bond may, at any time of the holder’s choosing, be converted to a multiple of the specified asset.  Thus a convertible bond with pricereturns an amountat timeunless the owner has converted the bond to units of the asset at some time prior to. This definition, the differential equation and boundary conditions are given in Chapter 18 of  Wilmott et al. (1996). Using a constant interest rate and volatility factor, the parameters and boundary conditions are:

1.     Bond face value , conversion factor

2.     Volatility

3.     Times until expiration

4.     Interest rate , dividend

5.    

6.    

7.     Boundary conditions

8.     Terminal data

9.     Constraint for bond holder

Note that the error tolerance is set to a pure absolute error of value.  The free boundary constraint  is achieved by use of a non-linear forcing term in the subroutine fkforce_cbond.  The terminal conditions are provided with the user subroutine fkinit_cbond.

Link to example source (feynman_kac_ex4.f90)

 

! Compute value of a Convertible Bond

      use feynman_kac_int

      use hqsval_int

      use mp_types

      use umach_int

 

      implicit none

! The face value

      real(kind(1e0)) :: ks = 1.0e0

! The sigma  or volatility value

      real(kind(1e0)) :: sigma = 0.25e0

 

! Time values for the options

      integer, parameter :: nt = 2

      real(kind(1e0)) :: time(nt)=(/0.5e0, 1.0e0/)

! Values of the underlying where evaluation are made

      integer, parameter :: nv = 13

      integer, parameter :: nlbc = 3, nrbc = 3, ndeg = 6

      integer :: i

      real(kind(1e0)) :: xs(nv) = (/((i-1)*0.25e0,i=1,nv)/)

! Value of the interest rate, continuous dividend and factor

      real(kind(1e0)) :: r = 0.1e0, dividend = 0.02e0,&

                         factor =1.125e0

! Values of the min and max underlying values modeled

      real(kind(1e0)) :: x_min = 0.0e0, x_max = 4.0e0

! Define parameters for the integration step. 

      integer, parameter :: nx = 61, nint = nx-1, n = 3*nx

      real(kind(1e0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt),&

                         dx, f(nv,0:nt)

! Relative and absolute error tolerances

      real(kind(1e0)) :: atol(1), rtol(1)

      type(s_fcn_data) fcn_data

      real(kind(1e0)), external :: fkcoef_cbond, fkinitcond_cbond

      external fkbc_cbond, fkforce_cbond, fkinit_cbond

 

      integer :: nout

 

      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

 

! Use a pure absolute error tolerance for the integration

! The default values require too much integation time.

      atol(1) = 1.0e-3

      rtol(1) = 0.0e0

 

! Pass the data for evaluation 

      fcn_data % rdata = (/ks,x_max,sigma,r,dividend,factor,&

                           atol(1)/)

      fcn_data % idata = (/ndeg/)

         

! Compute value of convertible bond                

      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_cbond,&

                     fkinitcond_cbond, fkbc_cbond, y, yprime,&

                     ATOL=atol, RTOL=rtol, FKINIT = fkinit_cbond,&

                     FKFORCE = fkforce_cbond, FCN_DATA = fcn_data)

  

! Evaluate and display solutions at vector of points XS(:), at each

! time value prior to expiration.                       

      do i=0,nt

        f(:,i) = hqsval (xs, xgrid, y(:,i))

      end do

 

      write(nout,'(T05,A)')&

      "Convertible Bond Value, 0+, 6 and 12 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,/'//&

         'T08,"Interest Rate =",F5.2,", Dividend =",F5.2,'//&

         '", Factor = ",F5.3,//T08,"Underlying", T26,"Bond Value",'//&

         '/(T10,4F8.4))') ks,sigma,r,dividend,factor,&

           (xs(i), f(i,0:nt),i=1,nv)

      end

 

! These routines define the coefficients, payoff, boundary

! conditions and forcing term.

      function fkcoef_cbond(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, zero = 0.e0

 

      sigma = fcn_data % rdata(3)

      interest_rate = fcn_data % rdata(4)

      dividend = fcn_data % rdata(5)

 

      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_cbond

 

      function fkinitcond_cbond(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, factor

 

      strike_price = fcn_data % rdata(1)

      factor = fcn_data % rdata(6)

      value = max(factor * x, strike_price)

      return

      end function fkinitcond_cbond

 

      subroutine fkbc_cbond (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)) :: interest_rate, strike_price, dp,&

                         factor, x_max

 

      select case (iflag)

        case (1)

          strike_price = fcn_data % rdata(1)

          interest_rate = fcn_data % rdata(4)

          dp = strike_price * exp(tx*interest_rate)

          bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,    dp/)

          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

        case (2)

          x_max = fcn_data % rdata(2)

          factor = fcn_data % rdata(6)

          bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, factor * x_max/)

          bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, factor/)

          bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)

      end select

! Note no time dependence                 

      iflag = 0

      return

      end subroutine fkbc_cbond

 

      subroutine fkforce_cbond (interval, t, hx, y, xlocal, qw, u,&

                                phi, dphi, fcn_data)

      use mp_types

      implicit none

      integer :: i, j, l

      integer, parameter :: local = 6

      integer, intent(in) :: interval

      real(kind(1.e0)), intent(in) :: y(:), t, hx, qw(:),xlocal(:),&

                                      u(:,:)

      real(kind(1.e0)), intent(out) :: phi(:), dphi(:,:)

 

      integer :: ndeg

      real(kind(1.e0)) :: yl(local), bf(local)

      real(kind(1.e0)) :: value, strike_price, interest_rate,&

                          zero = 0.0e0, one = 1.0e0, rt, mu, factor

      type(s_fcn_data), optional :: fcn_data

 

      yl = y(3*interval-2:3*interval+3)

      phi = zero

      dphi = zero

      value = fcn_data % rdata(7)

      strike_price = fcn_data % rdata(1)

      interest_rate = fcn_data % rdata(4)

      factor = fcn_data % rdata(6)

      ndeg = fcn_data % idata(1)

      mu = 2

! This is the local definition of the forcing term

! It "forces" the constraint f >= factor*x.       

      do j=1,local

        do l = 1,ndeg

          bf(1:3) = u(l,1:3)

          bf(4:6) = u(l,7:9)

          rt = dot_product(yl,bf)

          rt = value/(rt + value - factor * xlocal(l))

          phi(j) = phi(j) + qw(l) * bf(j) * rt**mu

        end do

      end do

      phi = -phi * hx * factor * strike_price

! This is the local derivative matrix for the forcing term -       

      do j=1,local

        do i = 1,local

          do l=1,ndeg

            bf(1:3) = u(L,1:3)

            bf(4:6) = u(L,7:9)

            rt = dot_product(yl,bf)

            rt = one/(rt + value - factor * xlocal(l))

            dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j)&

                     * (value * rt)**mu * rt

          end do

        end do

      end do

      dphi = -mu * dphi * hx * factor * strike_price

      return                

      end subroutine fkforce_cbond

         

      subroutine fkinit_cbond(xgrid,tgrid,t,yprime,y,atol,rtol,&

                              fcn_data)

      use mp_types

      implicit none

      real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)

      real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), yprime(:),&

                                     t

      type(s_fcn_data), optional :: fcn_data

 

      integer :: i

      if (t == 0.0e0) then

! Set initial data precisely.

        do i=1,size(Xgrid)

          if (xgrid(i)*fcn_data % rdata(6) <&

              fcn_data % rdata(1)) then

            y(3*i-2) = fcn_data % rdata(1)

            y(3*i-1) = 0.0e0

            y(3*i  ) = 0.0e0

          else

            y(3*i-2) = xgrid(i) * fcn_data % rdata(6)

            y(3*i-1) = fcn_data % rdata(6)

            y(3*i  ) = 0.0e0

          end if

        end do

      end if

      end subroutine fkinit_cbond

Output

    Convertible Bond Value, 0+, 6 and 12 Months Prior to Expiry

       Number of equally spaced spline knots   61

       Number of unknowns  183

       Strike=  1.00, Sigma= 0.25

       Interest Rate= 0.10, Dividend= 0.02, Factor= 1.125

           

       Underlying        Bond Value

           0.0000  1.0000  0.9512  0.9048

           0.2500  1.0000  0.9512  0.9049

           0.5000  1.0000  0.9513  0.9065

           0.7500  1.0000  0.9737  0.9605

           1.0000  1.1250  1.1416  1.1464

           1.2500  1.4062  1.4117  1.4121

           1.5000  1.6875  1.6922  1.6922

           1.7500  1.9688  1.9731  1.9731

           2.0000  2.2500  2.2540  2.2540

           2.2500  2.5312  2.5349  2.5349

           2.5000  2.8125  2.8160  2.8160

           2.7500  3.0938  3.0970  3.0970

           3.0000  3.3750  3.3781  3.3781

Example 5 – A Non-Standard American Option

This example illustrates a method for evaluating a certain “Bermudan Style” or non-standard American option.  These options are American Style options restricted to certain dates where the option may be exercised.  Since this agreement gives the holder more opportunity than a European option, it is worth more.  But since the holder can only exercise at certain times it is worth no more than the American style option value that can be exercised at any time.  Our solution method uses the same model and data as in Example 2, but allows exercise at weekly intervals.  Thus we integrate, for half a year, over each weekly interval using a European style Black-Scholes model, but with initial data at each new week taken from the corresponding values of the American style option.

Link to example source (feynman_kac_ex5.f90)

 

! Compute Bermudan-Style Option Premium for Vanilla Put

      use feynman_kac_int

      use hqsval_int

      use mp_types

      use umach_int

 

      implicit none

      integer :: nout

! The strike price

      real(kind(1e0)) :: ks = 10.0e0

! The sigma value

      real(kind(1e0)) :: sigma = 0.4e0

! Program working stores

      real(kind(1e0)) :: week

! Time values for the options

      integer, parameter :: nt = 26

      integer, parameter :: ndeg = 6

      real(kind(1e0)) :: time(nt), time_end = 0.5e0

! Values of the underlying where evaluation are made

      integer, parameter :: nv = 9, nlbc = 2, nrbc = 3

      integer :: i

      real(kind(1e0)) :: xs(nv) = (/((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),&

                         ya(n,0:nt), yaprime(n,0:nt),&

                         ytemp(n,0:1), ytempprime(n,0:1),&

                         dx, fb(nv,nt), fa(nv,nt)

      real(kind(1e0)) :: atol

      type(s_fcn_data) fcn_data_amer, fcn_data_berm

      real(kind(1e0)), external :: fkcoef_put, fkinitcond_put

      external fkbc_put, fkforce_put, fkinit_amer_put, fkinit_berm_put

 

      call umach(2, nout)

! Allocate space inside the derived type for holding

! data values. These are for the evaluation routines.

      allocate(fcn_data_amer % rdata (6), fcn_data_amer % 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

 

! Place a breakpoint at the strike price.

      do i=1,nx

        if (xgrid(i) > ks) then

          xgrid(i-1) = ks

          exit

        end if

      end do

 

! Compute time values where American option is computed

      week = time_end/real(nt,kind(week))

      time(1) = week

      do i=2,nt-1

        time(i) = time(i-1) + week

      end do

      time(nt) = time_end

 

      atol = 1.0e-3

      fcn_data_amer % rdata = (/ks,x_max,sigma,r,dividend,atol/)

      fcn_data_amer % idata = (/ndeg/)

 

! Compute American Put Option Values at the weekly grid.

      call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&

                        fkinitcond_put, fkbc_put, ya, yaprime,&

                        FKINIT = fkinit_amer_put,&

                        FKFORCE = fkforce_put,&

                        FCN_DATA = fcn_data_amer)

! Integrate once again over the weekly grid, using the American

! Option values as initial data for a piece-wise European option

!integration.

 

! Allocate space to hold coefficient data and initial values.

      allocate(fcn_data_berm % rdata(5+n))

      fcn_data_berm % rdata(1:5) = fcn_data_amer % rdata(1:5)

! Copy initial data so the payoff value is the same for

! American and Bermudan option values.       

      yb(1:n,0) = ya(1:n,0)

      ybprime(1:n,0) = ya(1:n,0)

 

      do i=0,nt-1

! Move American Option values into place as initial conditions,

! but now integrating with European style over each period of

! the weekly grid.      

        fcn_data_berm % rdata(6:) = ya(1:n,i)

        if (i .eq. 0) then

          call feynman_kac (xgrid, (/time(1)/), nlbc, nrbc,&

                            fkcoef_put, fkinitcond_put, fkbc_put,&

                            ytemp(:,0:1), ytempprime(:,0:1),&

                            FKINIT = fkinit_berm_put,&

                            FCN_DATA = fcn_data_berm)

        else

          call feynman_kac (xgrid, (/time(i+1)-time(i)/),&

                            nlbc, nrbc, fkcoef_put,&

                            fkinitcond_put, fkbc_put,&

                            ytemp(:,0:1), ytempprime(:,0:1),&

                            FKINIT = fkinit_berm_put,&

                            FCN_DATA = fcn_data_berm)

        end if

! Record values of the Bermudan option at the end of each integration.

        yb(1:n,i+1) = ytemp(1:n,1)

        ybprime(1:n,i+1) = ytempprime(1:n,1)

      end do                                           

! Evaluate solutions at vector of points XS(:), at each time value

! prior to expiration.                       

      do i=1,nt

        fa(:,i) = hqsval (xs, xgrid, ya(:,i))

        fb(:,i) = hqsval (xs, xgrid, yb(:,i))

      end do

 

      write(nout,'(T05,A)')&

        "American Option Premium for Vanilla Put, 6 Months "//&

        "Prior to Expiry"

      write(nout,'(T05,A)')&

        "Exercise Opportunities At Weekly Intervals"

      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,"Underlying",'//&

                   'T20,"Bermudan Style", T42,"American",'//&

                   '/(T10,F8.4, T26, F8.4, T42, F8.4))')&

                   KS,SIGMA,R,&

                  (xs(i), fb(i,nt:nt), fa(i,nt:nt),i=1,nv)

      end

 

! These subprograms set the coefficients, payoff, boundary

! conditions and forcing term for American and European Options.

      function fkcoef_put(x, tx, iflag, fcn_data_amer)&

                          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_amer

      real(kind(1e0)) :: value

 

      real(kind(1e0)) :: sigma, interest_rate, dividend, zero=0.0e0

 

      sigma = fcn_data_amer % rdata(3)

      interest_rate = fcn_data_amer % rdata(4)

      dividend = fcn_data_amer % rdata(5)

      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_put

 

      function fkinitcond_put(x, fcn_data_amer) result(value)

      use mp_types

      implicit none

 

      real(kind(1e0)), intent(in) :: x

      type (s_fcn_data), optional :: fcn_data_amer

 

      real(kind(1e0)) :: value

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

 

      strike_price = fcn_data_amer % rdata(1)

! The payoff function

      value = max(strike_price - x, zero)

      return

      end function fkinitcond_put

 

      subroutine fkbc_put (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

 

      select case (iflag)

        case (1)

          bccoefs(1,1:4) = ((/0.0e0, 1.0e0, 0.0e0, -1.0e0/))

          bccoefs(2,1:4) = ((/0.0e0, 0.0e0, 1.0e0,  0.0e0/))

        case (2)

          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/))

      end select

! Note no time dependence

      iflag = 0

      end subroutine fkbc_put

 

      subroutine fkforce_put (interval, t, hx, y, xlocal, qw, u,&

                              phi, dphi, fcn_data_amer)

      use mp_types

      implicit none

      integer, parameter :: local = 6

      integer :: i, j, l, ndeg

      integer, intent(in) :: interval

      real(kind(1.e0)), intent(in) :: y(:), t, hx, qw(:),&

                                      xlocal(:), u(:,:)

      real(kind(1.e0)), intent(out) :: phi(:), dphi(:,:)

      type (s_fcn_data), optional :: fcn_data_amer

 

      real(kind(1.e0)) :: yl(local), bf(local)

      real(kind(1.e0)) :: value, strike_price, interest_rate,&

                          zero = 0.e0, one = 1.e0, rt, mu

 

      yl = y(3*interval-2:3*interval+3)

      phi = zero

      value = fcn_data_amer % rdata(6)

      strike_price = fcn_data_amer % rdata(1)

      interest_rate = fcn_data_amer % rdata(4)

      ndeg = fcn_data_amer % idata(1)

 

      mu = 2

! This is the local definition of the forcing term

      do j=1,local

        do l=1,ndeg

          bf(1:3) = U(L,1:3)

          bf(4:6) = U(L,7:9)

          rt = dot_product(YL,BF)

          rt = value/(rt + value-(strike_price-xlocal(l)))

          phi(j) = phi(j) + qw(l) * bf(j) * rt**mu

        end do

      end do

      phi = -phi * hx * interest_rate * strike_price

! This is the local derivative matrix for the forcing term

      dphi = zero

      do j=1,local

        do i = 1,local

          do l=1,ndeg

            bf(1:3) = u(L,1:3)

            bf(4:6) = u(L,7:9)

            rt = dot_product(yl,bf)

            rt = one/(rt + value - (strike_price - xlocal(l)))

            dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j) *&

                        rt**(mu+1)

          end do

        end do

      end do

      dphi = mu * dphi * hx * value**mu * interest_rate *&

             strike_price

      end subroutine fkforce_put

 

      subroutine fkinit_amer_put(xgrid,tgrid,t,yprime,y,atol,rtol,&

                                 fcn_data_amer)

      use mp_types

      implicit none

 

      real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&

                                     yprime(:)

      real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)

      type(s_fcn_data), optional :: fcn_data_amer

      integer :: i

 

      if (t == 0.0e0) then

! Set initial data precisely.  The strike price is a breakpoint.

! Average the derivative limit values from either side.

        do i=1,size(xgrid)

          if (xgrid(i) < fcn_data_amer % rdata(1)) then

            y(3*i-2) = fcn_data_amer % rdata(1) - xgrid(i)

            y(3*i-1) = -1.0e0

            y(3*i) = 0.0e0

          else if (xgrid(i) == fcn_data_amer % rdata(1)) then

            y(3*i-2) = 0.0e0

            y(3*i-1) = -0.5e0

            y(3*i) = 0.0e0

          else

            y(3*i-2) = 0.0e0

            y(3*i-1) = 0.0e0

            y(3*i) = 0.0e0

          end if

        end do

      end if

      end subroutine fkinit_amer_put

 

      subroutine fkinit_berm_put(xgrid,tgrid,t,yprime,y,atol,rtol,&

                                 fcn_data_berm)

      use mp_types

      implicit none

 

      real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&

                                     yprime(:)

      real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)

      type(s_fcn_data), optional :: fcn_data_berm

      integer :: i

 

      if (t == 0.0e0) then

! Set initial data for each week at the previously computed

! American Option values.  These coefficients are passed

! in the derived type fcn_data_berm.

        do i=1,size(xgrid)

          y(3*i-2) = fcn_data_berm % rdata(3+3*i)

          y(3*i-1) = fcn_data_berm % rdata(4+3*i)

          y(3*i  ) = fcn_data_berm % rdata(5+3*i)

        end do

      end if

      end subroutine fkinit_berm_put

Output

    American Option Premium for Vanilla Put, 6 Months Prior to Expiry

    Exercise Opportunities At Weekly Intervals

       Number of equally spaced spline knots   61

       Number of unknowns  183

       Strike= 10.00, Sigma= 0.40, Interest Rate= 0.10

 

       Underlying  Bermudan Style        American

           0.0000          9.9808         10.0000

           2.0000          7.9808          8.0000

           4.0000          5.9808          6.0000

           6.0000          3.9808          4.0000

           8.0000          2.0924          2.0926

          10.0000          0.9139          0.9138

          12.0000          0.3570          0.3569

          14.0000          0.1309          0.1309

          16.0000          0.0468          0.0469

Example 6 – Oxygen Diffusion Problem

Our previous examples are from the field of financial engineering.   A final example is a physical model.  The Oxygen Diffusion Problem is summarized in Crank [4], p. 10-20, 261-262.  We present the numerical treatment of the transformed one-dimensional system

A slight difference from the Crank development is that we have reflected the time variable to match our form of the Feynman-Kac equation. We have a free boundary problem because the interfaceis implicit.  This interface is implicitly defined by solving the variational relation.  The first factor is zero for and the second factor is zero for.  We list the Feynman-Kac equation coefficients, forcing term and boundary conditions, followed by comments.

The  forcing term has the property of being almost the value 1 when the solution is larger than the factor.  As the solution, the forcing term  is almost the value zero. These properties combine to approximately achieve the variational relation that defines the free boundary.  Note that the arc of the free boundary is not explicitly available with this numerical method.  We have usedthe requested absolute error tolerance for the numerical integration.

The boundary condition is discontinuous as , since the initial data yields.  For the numerical integration we have chosen a boundary value function that starts with the value  at  and rapidly blends to the value zero as the integration proceeds in the negative direction.  It is necessary to give the integrator the opportunity to severely limit the step size near  to achieve continuity during the blending.

In the example code, values of  are checked against published values for certain values of .  Also checked are values of at published values of the free boundary, for the same values of .

Link to example source (feynman_kac_ex6.f90)

 

! Integrate Oxygen Diffusion Model found in the book

! Crank, John. Free and Moving Boundary Problems,

! Oxford Univ. Press, (1984), p. 19-20 and p. 261-262.

      use feynman_kac_int

      use hqsval_int

      use mp_types

      use norm_int

      use umach_int

      implicit none

 

      integer :: nout

      real(kind(1e0)), allocatable :: xgrid(:), tgrid(:), y(:,:),&

                                      yprime(:,:), f(:,:), s(:)

      real(kind(1e0)) :: dx, x_min, x_max, zero=0.0e0, one=1.0e0

      real(kind(1e0)) :: atol(1), rtol(1)

      type(s_fcn_data) :: fk_ox2

      integer :: i, nint, n, nunit, ntimes = 8

      integer, parameter :: ndeg = 6, nlbc = 1, nrbc = 2

      real(kind(1e0)), external :: fkcoef_ox2, fkinitcond_ox2

      external fkbc_ox2, fkforce_ox2

     

      call umach(2,nout)

 

! Define number of equally spaced intervals for elements

      nint = 100

! Allocate the space needed for the integration process

      n = 3*(nint+1)

      allocate(xgrid(nint+1), y(n,0:ntimes), yprime(n,0:ntimes),&

               tgrid(ntimes), f(1,ntimes), s(ntimes))

 

! Allocate space inside the derived type for holding

! data values. These are for the evaluation routines.

      allocate(fk_ox2 % rdata (1), fk_ox2 % idata (1))

 

      atol(1) = 0.5e-2

      rtol(1) = 0.5e-2

      fk_ox2 % rdata(1) = atol(1)

      fk_ox2 % idata(1) = ndeg

 

! Define interval endpoints

      x_min = zero

      x_max = one

! Define interval widths

      dx = (x_max-x_min)/real(nint)

      xgrid(1) = x_min

      xgrid(nint+1) = x_max

! Define grid points of interval

      do i=2,nint

        xgrid(i) = xgrid(i-1) + dx

      end do

! Define time integration output points

! These correspond to published values in Crank's book, p. 261-262

      tgrid = (/0.04e0,0.06e0,0.10e0,0.12e0,0.14e0,&

                0.16e0,0.18e0,0.185e0/)

      call feynman_kac (xgrid, tgrid, nlbc, nrbc, fkcoef_ox2,&

                        fkinitcond_ox2, fkbc_ox2, y, yprime,&

                        ATOL = atol, RTOL = rtol,&

                        FKFORCE = fkforce_ox2, FCN_DATA = fk_ox2)

! Summarize output at the left end

      do i=1,ntimes

        f(:,i)= hqsval ((/zero/), xgrid, y(:,i))

      end do

! Check differences of evaluation and published left end values

      f(1,:) = f(1,:) - (/2.743e-01, 2.236e-01, 1.432e-01,&

           1.091e-01, 7.786e-02, 4.883e-02, 2.179e-02, 1.534e-02/)

      write(nout,*) "Oxygen Depletion Model, from Crank's "//&

                    "Book, p. 261-262,"

      write(nout,*) "'Free and Moving Boundary Value Problems'"

      if (norm(f(1,:)) < ntimes * atol(1)) then

        write(nout,*) "FEYNMAN_KAC Example 6 - Fixed Sealed "//&

                      "Surface Values are correct"

      else

        write(nout,*) "FEYNMAN_KAC Example 6 - does not agree with"//&

                      " published left end values"

      end if

! Define known position of free boundary at the time points

      s = (/0.9992e0,0.9918e0,0.9350e0,0.8792e0,&

            0.7989e0,0.6834e0,0.5011e0,0.4334e0/)

                      

! Evaluate and verify solution is small near free boundary -

      do i=1,ntimes

        f(:,i) = hqsval ((/s(i)/), xgrid, y(:,i))

      end do

 

      if (norm(f(1,:)) < ntimes * atol(1)) then

        write(nout,*) "FEYNMAN_KAC Example 6 - Free Boundary "//&

                   "Position Values are correct"

      else

        write(nout,*) "FEYNMAN_KAC Example 6 - does not agree "//&

                   "with published free boundary values"

      end if

      end

 

      function fkcoef_ox2 (x, tx, iflag, fk_ox2) result(value)

      use mp_types

      implicit none

! Coefficient valuation routine for the Oxygen Diffusion

! Model found in Crank's book, p. 20

! Input/Ouput variables

      integer, intent(inout) :: iflag

      real(kind(1e0)), intent(in) :: x, tx

      type(s_fcn_data), optional :: fk_ox2

      real(kind(1e0)) :: value

 

! Local variables

      real(kind(1e0)) :: zero = 0.0e0, two = 2.0e0        

 

      select case (iflag)

        case (1) ! Factor DSigma/Dx(x,t)

          value = zero

        case (2) !  Factor Sigma(x,t)

          value = sqrt(two)

        case (3) !  Factor Mu (x,t)

          value = zero

        case (4) !  Factor Kappa (x,t)

          value = zero

      end select

! Signal no dependence on tx=t=time for any coefficient.

      iflag = 0

      return

      end function fkcoef_ox2

 

      function fkinitcond_ox2(x, fk_ox2) result(value)

      use mp_types

      implicit none

      real(kind(1e0)), intent(in) :: x

      type (s_fcn_data), optional :: fk_ox2

      real(kind(1e0)) :: value

 

      real(kind(1e0)) :: half = 0.5e0, one = 1.0e0

 

      value = half * (one - x)**2

      return          

      end function fkinitcond_ox2

 

      subroutine fkbc_ox2 (tx, iflag, bccoefs, fk_ox2)

      use mp_types

      implicit none

! Evaluation routine for Oxygen Diffusion Model

! boundary conditions. 

! Input/Ouput variables

      real(kind(1e0)), intent(in) :: tx

      integer, intent(inout) :: iflag

      real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs

      type (s_fcn_data), optional :: fk_ox2

 

! Local variables

      real(kind(1e0)) :: zero = 0.0e0, one = 1.0e0, atol

 

      atol = fk_ox2 % rdata(1)

      select case (iflag)

        case (1) ! Left Boundary Condition,  at X_min=0

! There is a rapid blending of the boundary condition to achieve

! a zero derivative value at the left end.

! The initial data has the derivative with value one.

! This boundary condition essentially abruptly changes that

! derivative value to zero.

! Returning iflag=1 signals time dependence.  This is important

! for this problem.

          bccoefs(1,1:4) = (/0.0e0, one, 0.0e0, exp(tx/atol**2)/)

          return                    

        case (2) ! Right Boundary Condition, at X_max=1

          bccoefs(1,1:4) = (/one, 0.0e0, 0.0e0, 0.0e0/)

          bccoefs(2,1:4) = (/0.0e0, one, 0.0e0, 0.0e0/)

      end select

      iflag = 0 ! Signal no dependence on tx=time.

      end subroutine fkbc_ox2

 

      subroutine fkforce_ox2 (interval, t, hx, y, xlocal, qw, u,&

                              phi, dphi, fk_ox2)

! Evaluation routine for Oxygen Diffusion model forcing function.

      use mp_types

      implicit none

      integer, parameter :: local = 6

      integer :: i, j, l, mu, ndeg

      integer, intent(in) :: interval

 

      real(kind(1e0)), intent(in) :: y(:), t, hx, qw(:),&

                                     xlocal(:), u(:,:)

      real(kind(1e0)), intent(out) :: phi(:), dphi(:,:)

      type (s_fcn_data), optional :: fk_ox2

 

      real(kind(1e0)) :: yl(local), bf(local)

      real(kind(1e0)) :: value, zero = 0.0e0, one = 1.0e0, rt

 

      yl = y(3*interval-2:3*interval+3)

      phi = zero

 

      value = fk_ox2 % rdata(1)

      ndeg = fk_ox2 % idata(1)

 

      mu = 2

      do j=1,local

        do l=1,ndeg

          bf(1:3) = u(l,1:3)

          bf(4:6) = u(l,7:9)

          rt = dot_product(yl,bf)

          rt = one - (value/(rt + value))**mu             

          phi(j) = phi(j) + qw(l) * bf(j) * RT

        end do

      end do

      phi = phi * hx

! This is the local derivative matrix for the forcing term -       

      dphi = zero

      do j=1,local

        do i = 1,local

          do l=1,ndeg

            bf(1:3) = u(l,1:3)

            bf(4:6) = u(l,7:9)

            rt = dot_product(yl,bf)

            rt = one/(rt + value)

            dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j) *&

                        rt**(mu+1)

          end do

        end do

      end do

      dphi = mu * dphi * hx * value**mu

      return                               

      end subroutine fkforce_ox2

Output

Oxygen Depletion Model, from Crank's Book, p. 261-262,

'Free and Moving Boundary Value Problems'

FEYNMAN_KAC Example 6 - Fixed Sealed Surface Values are correct

FEYNMAN_KAC Example 6 - Free Boundary Position Values are correct

Example 7 – Calculating the Greeks

In this example, routine FEYNMAN_KAC is used to solve for the Greeks, i.e. various derivatives of Feynman-Kac (FK) solutions applicable to the pricing of options and related financial derivatives. In order to illustrate and verify these calculations, the Greeks are calculated by two methods. The first method involves the FK solution to the diffusion model for call options given in Example 1 for the Black-Scholes (BS) case, i.e. . The second method calculates the Greeks using the closed-form BS evaluations which can be found at http://en.wikipedia.org/wiki/The_Greeks.

This example calculates FK and BS solutions  to the BS problem and the following Greeks:

•      Delta =

 is the first derivative of the Value, , of a portfolio of derivative security derived from underlying instrument with respect to the underlying instrument’s price S;

•      Gamma = ;

•      Theta =  is the negative first derivative of V with respect to time t;

•      Charm = ;

•      Color = ;

•      Rho =  is the first derivative of V with respect to the risk-free rate r;

•      Vega =  measures sensivity to volatility parameter  of the underlying S;

•      Volga = ;

•      Vanna = ;

•      Speed = .

Intrinsic Greeks include derivatives involving only S and t, the intrinsic FK arguments. In the above list, Value, Delta, Gamma, Theta, Charm, Color and Speed are all intrinsic Greeks. As is discussed in Hanson, R. (2008) “Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements”, the expansion of the FK solution function  in terms of quintic polynomial functions defined on S-grid subintervals and subject to continuity constraints in derivatives 0, 1 and 2 across the boundaries of these subintervals allows Value, Delta, Gamma, Theta, Charm and Color to be calculated directly by routines FEYNMAN_KAC and HQSVAL.

Non-intrinsic Greeks are derivatives of V involving FK parameters other than the intrinsic arguments S and t, such as r and . Non-intrinsic Greeks in the above list include Rho, Vega, Volga and Vanna. In order to calculate non-intrinsic Greek (parameter) derivatives or intrinsic Greek S-derivatives beyond the second (such as Speed) or t-derivatives beyond the first, the entire FK solution must be calculated 3 times (for a parabolic fit) or five times (for a quartic fit), at the point where the derivative is to be evaluated and at nearby points located symmetrically on either side.

Using a Taylor series expansion of   truncated to  terms (to allow an m-degree polynomial fit of  m+1 data points),

,

we are able to derive the following parabolic (3 point) estimation of first and second derivatives  and   in terms of the three values ,  and , where  and 0 < ε frac <<1:

 ,

.

Similarly, the quartic (5 point) estimation of  and  in terms of , ,,  and  is:

.

For our example, the quartic estimate does not appear to be significantly better than the parabolic estimate, so we have produced only parabolic estimates by setting variable iquart to 0. The user may try the example with the quartic estimate simply by setting iquart to 1.

As is pointed out in Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements, the quintic polynomial expansion function used by FEYNMAN_KAC only allows for continuous derivatives through the second derivative. While up to fifth derivatives can be calculated from the quintic expansion (indeed  function HQSVAL will allow the third derivative to be calculated by setting optional argument IDERIV to 3, as is done in this example), the accuracy is compromised by the inherent lack of continuity across grid points (i.e. subinterval boundaries).

The accurate second derivatives in S returned by function HQSVAL can be leveraged into a third derivative estimate by calculating three FK second derivative solutions, the first solution for grid and evaluation point  set  and the second and third solutions for solution grid and evaluation point sets   and , where the solution grid and evaluation point sets are shifted up and down by . In this example, is set to , where  is the average value of S over the range of grid values and . The third derivative solution can then be obtained using the parabolic estimate

This procedure is implemented in the current example to calculate the Greek Speed. (For comparison purposes, Speed  is also calculated directly by setting the optional argument IDERIV to 3 in the call to function HQSVAL. The output from this direct calculation is called “Speed2”.)

To reach better accuracy results, all computations are done in double precision.

The average and maximum relative errors (defined as the absolute value of the difference between the BS and FK values divided by the BS value) for each of the Greeks is given at the end of the output. (These relative error statistics are given for nine combinations of Strike Price and volatility, but only one of the nine combinations is actually printed in the output.) Both intrinsic and non-intrinsic Greeks have good accuracy (average relative error is in the range 0.01 – 0.00001) except for Volga, which has an average relative error of about 0.05. This is probably a result of the fact that Volga involves differences of differences, which will more quickly erode accuracy than calculations using only one difference to approximate a derivative. Possible ways to improve upon the 2 to 4 significant digits of accuracy achieved in this example include increasing FK integration accuracy by reducing the initial stepsize (via optional argument RINITSTEPSIZE), by choosing more closely spaced S and t grid points (by adjusting FEYNMAN_KAC’s input arrays XGRID and TGRID) and by adjusting   so that the central differences used to calculate the derivatives are not too small to compromise accuracy.

Link to example source (feynman_kac_ex7.f90)

 

! 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

Output

    Constant Elasticity of Variance Model for Vanilla Call Option

         Interest Rate:   0.050      Continuous Dividend:   0.000

         Minimum and Maximum Prices of Underlying:    0.00  60.00

         Number of equally spaced spline knots:  120

         Number of unknowns:  363

 

 

         Time in Years Prior to Expiration:    0.0833 0.3333 0.5833

         Option valued at Underlying Prices:    19.00  20.00  21.00

 

         3 point (parabolic) estimate of parameter derivatives

         epsfrac =  0.00100000

 

 

         Strike =15.00 Sigma = 0.20 Alpha = 2.00:

 

         years to expiration:  0.0833 0.3333 0.5833

 

underlying price: 19.0; FK  Value:  4.0623732450  4.2575924184  4.4733805278

                        BS  Value:  4.0623732509  4.2575929678  4.4733814062

underlying price: 20.0; FK  Value:  5.0623700127  5.2505145764  5.4492418798

                        BS  Value:  5.0623700120  5.2505143129  5.4492428547

underlying price: 21.0; FK  Value:  6.0623699727  6.2485587059  6.4385718831

                        BS  Value:  6.0623699726  6.2485585270  6.4385720688

underlying price: 19.0; FK  Delta:  0.9999864098  0.9877532309  0.9652249945

                        BS  Delta:  0.9999863811  0.9877520034  0.9652261127

underlying price: 20.0; FK  Delta:  0.9999998142  0.9964646548  0.9842482622

                        BS  Delta:  0.9999998151  0.9964644003  0.9842476147

underlying price: 21.0; FK  Delta:  0.9999999983  0.9990831687  0.9932459040

                        BS  Delta:  0.9999999985  0.9990834124  0.9932451927

underlying price: 19.0; FK  Gamma:  0.0000543456  0.0144908955  0.0264849216

                        BS  Gamma:  0.0000547782  0.0144911447  0.0264824761

underlying price: 20.0; FK  Gamma:  0.0000008315  0.0045912854  0.0129288434

                        BS  Gamma:  0.0000008437  0.0045925328  0.0129280372

underlying price: 21.0; FK  Gamma:  0.0000000080  0.0012817012  0.0058860348

                        BS  Gamma:  0.0000000077  0.0012818272  0.0058865489

underlying price: 19.0; FK  Theta: -0.7472631891 -0.8301000450 -0.8845209253

                        BS  Theta: -0.7472638978 -0.8301108199 -0.8844992143

underlying price: 20.0; FK  Theta: -0.7468881086 -0.7706770630 -0.8152217385

                        BS  Theta: -0.7468880640 -0.7706789470 -0.8152097697

underlying price: 21.0; FK  Theta: -0.7468815742 -0.7479185416 -0.7728950748

                        BS  Theta: -0.7468815673 -0.7479153725 -0.7728982104

underlying price: 19.0; FK  Charm: -0.0014382828 -0.0879903285 -0.0843323992

                        BS  Charm: -0.0014397520 -0.0879913927 -0.0843403333

underlying price: 20.0; FK  Charm: -0.0000284881 -0.0364107814 -0.0547260337

                        BS  Charm: -0.0000285354 -0.0364209077 -0.0547074804

underlying price: 21.0; FK  Charm: -0.0000003396 -0.0126436426 -0.0313343015

                        BS  Charm: -0.0000003190 -0.0126437838 -0.0313252716

underlying price: 19.0; FK  Color:  0.0051622176  0.0685064195  0.0299871130

                        BS  Color:  0.0051777484  0.0684737183  0.0300398444

underlying price: 20.0; FK  Color:  0.0001188761  0.0355826975  0.0274292189

                        BS  Color:  0.0001205713  0.0355891884  0.0274307898

underlying price: 21.0; FK  Color:  0.0000015432  0.0143174420  0.0190897159

                        BS  Color:  0.0000015141  0.0143247729  0.0190752019

underlying price: 19.0; FK   Vega:  0.0003289870  0.3487168323  1.1153520921

                        BS   Vega:  0.0003295819  0.3487535501  1.1153536190

underlying price: 20.0; FK   Vega:  0.0000056652  0.1224632724  0.6032458218

                        BS   Vega:  0.0000056246  0.1224675413  0.6033084039

underlying price: 21.0; FK   Vega:  0.0000000623  0.0376974472  0.3028275297

                        BS   Vega:  0.0000000563  0.0376857196  0.3028629419

underlying price: 19.0; FK  Volga:  0.0286254576  8.3705173459 16.7944554708

                        BS  Volga:  0.0286064650  8.3691191978 16.8219823169

underlying price: 20.0; FK  Volga:  0.0007137402  4.2505025277 12.9315441466

                        BS  Volga:  0.0007186004  4.2519372748 12.9612638820

underlying price: 21.0; FK  Volga:  0.0000100364  1.7613083436  8.6626161799

                        BS  Volga:  0.0000097963  1.7617504949  8.6676581034

underlying price: 19.0; FK  Vanna: -0.0012418872 -0.3391850563 -0.6388552010

                        BS  Vanna: -0.0012431594 -0.3391932673 -0.6387423326

underlying price: 20.0; FK  Vanna: -0.0000244490 -0.1366771953 -0.3945466661

                        BS  Vanna: -0.0000244825 -0.1367114682 -0.3945405194

underlying price: 21.0; FK  Vanna: -0.0000002904 -0.0466333335 -0.2187406645

                        BS  Vanna: -0.0000002726 -0.0466323413 -0.2187858632

underlying price: 19.0; FK    Rho:  1.2447807022  4.8365676561  8.0884594648

                        BS    Rho:  1.2447806658  4.8365650322  8.0884502627

underlying price: 20.0; FK    Rho:  1.2448021850  4.8929216544  8.3041708173

                        BS    Rho:  1.2448021908  4.8929245641  8.3041638392

underlying price: 21.0; FK    Rho:  1.2448024992  4.9107294560  8.4114197621

                        BS    Rho:  1.2448024996  4.9107310444  8.4114199038

underlying price: 19.0; FK  Speed: -0.0002124684 -0.0156265453 -0.0179534748

                        BS  Speed: -0.0002123854 -0.0156192867 -0.0179536520

underlying price: 20.0; FK  Speed: -0.0000037247 -0.0055877024 -0.0097502607

                        BS  Speed: -0.0000037568 -0.0055859333 -0.0097472434

underlying price: 21.0; FK  Speed: -0.0000000385 -0.0017085830 -0.0048143174

                        BS  Speed: -0.0000000378 -0.0017082128 -0.0048130214

underlying price: 19.0; FK Speed2: -0.0002310655 -0.0156276977 -0.0179516855

                        BS Speed2: -0.0002123854 -0.0156192867 -0.0179536520

underlying price: 20.0; FK Speed2: -0.0000043215 -0.0055923924 -0.0097502997

                        BS Speed2: -0.0000037568 -0.0055859333 -0.0097472434

underlying price: 21.0; FK Speed2: -0.0000000475 -0.0017117661 -0.0048153107

                        BS Speed2: -0.0000000378 -0.0017082128 -0.0048130214

 

 Greek:  Value;  avg rel err:  0.000146171196;  max rel err:  0.009030737566

 Greek:  Delta;  avg rel err:  0.000035817272;  max rel err:  0.001158483076

 Greek:  Gamma;  avg rel err:  0.001088392379;  max rel err:  0.044845800289

 Greek:  Theta;  avg rel err:  0.000054196359;  max rel err:  0.001412847300

 Greek:  Charm;  avg rel err:  0.001213347059;  max rel err:  0.064576457415

 Greek:  Color;  avg rel err:  0.003323954467;  max rel err:  0.136355681544

 Greek:   Vega;  avg rel err:  0.001514753397;  max rel err:  0.106255126885

 Greek:  Volga;  avg rel err:  0.058531380389;  max rel err:  1.639564208085

 Greek:  Vanna;  avg rel err:  0.001061525805;  max rel err:  0.065629483069

 Greek:    Rho;  avg rel err:  0.000146868262;  max rel err:  0.009438788128

 Greek:  Speed;  avg rel err:  0.002065441607;  max rel err:  0.073086615101

 Greek: Speed2;  avg rel err:  0.008429883935;  max rel err:  0.255746328973



http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260