FEYNMAN_KAC
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)
NLBC  3.
NRBC — The number of right boundary conditions. (Input)
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 p(x) in the Feynman-Kac PDE. The usage is FKINITCOND (X []), where
Function Return Value
FKINITCOND — Value of the initial condition function p(x).
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, xmin, and NRBC conditions at the right end, xmax. 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 xmin are stored row-wise according to the following matrix-scheme:
The coefficients for xmax 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) 0TGRID(1)TGRID(n). (Output)
For t = TGRID(j), j = 1, n, the coefficients are stored in columns 1n 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 Y(i,0), i = 1, m 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=1n, respectively.
The expressions for the basis functions βi(x) 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) 0TGRID(1)TGRID(n). (Output)
For t = 0 and t = TGRID(j)j=1n , the derivatives of the coefficients are stored in column 0 and columns 1 to n of array YPRIME, respectively. The columns in YPRIME represent
and
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 = 1n, 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
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
vector is defined as
PHI — Array of size 6 containing a Gauss-Legendre approximation for the local contribution , where t = T and . (Output)
Setting NDEG:=SIZE(XLOCAL) and xj = XLOCAL(j), array PHI contains elements
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 = T, (Output)
The approximation to this symmetric matrix is stored row-wise, i.e.
DPHI(i,j) = WIDTH
for ij = 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 evaluationsof 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
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 α  2. A numerical solution of this diffusion model yields the price of a call option. Various values of the strike price K, time values σ, and power coefficient α are used to evaluate the option price at values of the underlying price. The sets of parameters in the computation are:
1. power α = {2.0,1.0,0.0}
2. strike price K = {15.0,20.0,25.0}
3. volatility σ = {0.2,0.3,0.4}
4. times until expiration = {1/12,4/12,7/12}
5. underlying prices = {19.0,20.0,21.0}
6. interest rate r = 0.05
7. xmin = 0, xmax = 60
8. nx = 121, n = 3 × nx = 363
With this model the Feynman-Kac differential equation is defined by identifying:
The payoff function is the “vanilla option”, .
(Example 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 K = {10.0}
2. Volatility σ = {0.4}
3. Times until expiration = {1/4,1/2}
4. Interest rate r = 0.1
5. xmin = 0.0, xmax = 30.0
6. nx = 121, n = 3 × nx = 363
The payoff function is the “vanilla option”, .
(Example 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.
(Example 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
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,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)), parameter :: zero = 0.0e0

strike_price = fcn_data % rdata(1)
bet = fcn_data % rdata(2)
! 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
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

interest_rate, df

strike_price = fcn_data % rdata(1)
bet = fcn_data % rdata(2)
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,&
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

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 price x returns an amount K at time T unless the owner has converted the bond to νx, ν  1, units of the asset at some time prior to T. 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 K = {1}, conversion factor ν = 1.125
2. Volatility σ = {0.25}
3. Times until expiration = {1/2,1}
4. Interest rate r = 0.1, dividend D = 0.02
5. xmin = 0, xmax = 4
6. nx = 61, n = 3 × nx = 183
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 10-3. The free boundary constraint f(x,t) = νx 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.
(Example 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.
(Example 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 t  -t to match our form of the Feynman-Kac equation. We have a free boundary problem because the interface s(t) is implicit. This interface is implicitly defined by solving the variational relation (ft + fxx 1)f = 0, f  0. The first factor is zero for 0  x < s(t) and the second factor is zero for s(t x  1. 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 used ɛ = ATOL, the requested absolute error tolerance for the numerical integration.
The boundary condition fx(t,0) = 0, t < 0 is discontinuous as , since the initial data yields fx(0,0) = 1. For the numerical integration we have chosen a boundary value function that starts with the value 1 at t = 0 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 t = 0 to achieve continuity during the blending.
In the example code, values of f(0,t) are checked against published values for certain values of t. Also checked are values of f(0,s(t)) = 0 at published values of the free boundary, for the same values of t.
(Example 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 V(S,t) to the BS problem and the following Greeks:
Delta =
is the first derivative of the Value, V(S,t), 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 V(S,t) 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 f(σ + ɛ) truncated to m + 1 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 f (1)(σ) and f (2)(σ) in terms of the three values f(σ  ɛ ), f(σ) and f(σ + ɛ), where ɛ = ɛfracσ and 0 < ɛfrac << 1:
Similarly, the quartic (5 point) estimation of f (1)(σ) and f (2)(σ) in terms of f(σ  2ɛ), f(σ  ɛ ), f(σ), f(σ + ɛ), and f(σ  + 2ɛ) 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 {Sf (2)(S)} and the second and third solutions for solution grid and evaluation point sets {S + ɛf (2)(S + ɛ)} and {S + ɛ,f (2)(S  ɛ)} , 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 0 < ɛfrac << 1. 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 ɛfrac so that the central differences used to calculate the derivatives are not too small to compromise accuracy.
(Example 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