MOLCH

Deprecated Routine: MOLCH is a deprecated routine and has been replaced with MMOLCH. To view the deprecated documentation, see molch.pdf on the IMSL website. You can also access a local copy in your IMSL documentation directory at pdf\deprecated_routines\math\molch.pdf.

FEYNMAN_KAC

more...

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

Required Arguments

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

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

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

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

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

Function Return Value

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

Required Arguments

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

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

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

IFLAG	Computed coefficient
1
2	σ
3	μ
4	κ

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

Optional Arguments

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

FKCOEF must be declared EXTERNAL in the calling program.

FKINITCOND — User-supplied FUNCTION to evaluate the initial condition function 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

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) 0, TGRID(1), …, TGRID(n). (Output)
For t = TGRID(j), j = 1, …, n, the coefficients are stored in columns 1, …, n of array Y and the approximate solution is given by

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

The starting coefficients 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=1, …, n, 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) 0, TGRID(1), …, TGRID(n). (Output)
For t = 0 and t = TGRID(j), j=1, …, n , 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 = 1, …, n, and YPRIME(:,0) for the partials at t = 0.

Optional Arguments

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

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

Required Arguments

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

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

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

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

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

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

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

Optional Arguments

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 i, j = 1, …, 6.

Optional Arguments

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 ,

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:

power α = {2.0,1.0,0.0}

strike price K = {15.0,20.0,25.0}

volatility σ = {0.2,0.3,0.4}

times until expiration = {1/12,4/12,7/12}

underlying prices = {19.0,20.0,21.0}

interest rate r = 0.05

xmin = 0, xmax = 60

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

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.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

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

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

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:

Strike price K = {10.0}

Volatility σ = {0.4}

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

Interest rate r = 0.1

xmin = 0.0, xmax = 30.0

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

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

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:

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

Volatility σ = {0.4}.

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

Interest rate r = 0.1.

xmin = 0, xmax = 30.

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

! The spread value

real(kind(1e0)) :: ks2 = 15.0e0

! The Bet for the Cash-or-Nothing Call

real(kind(1e0)) :: bet = 2.0e0

! The sigma value

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

! Time values for the options

integer, parameter :: nt = 2

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

! Values of the underlying where evaluation are made

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

integer :: i

real(kind(1e0)) :: xs(nv) = (/(2+(I-1)*2.0e0,I=1,NV)/)

! Value of the interest rate and continuous dividend -

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

! Values of the min and max underlying values modeled -

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

! Define parameters for the integration step.

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

real(kind(1e0)) :: xgrid(nx), yb(n,0:nt), ybprime(n,0:nt),&

yv(n,0:nt), yvprime(n,0:nt),&

dx, fb(nv,nt), fv(nv,nt)

type(s_fcn_data) fcn_data

integer :: nout

real(kind(1e0)), external :: fkcoef_call, fkinitcond_call

external fkbc_call

call umach(2, nout)

! Allocate space inside the derived type for holding

! data values. These are for the evaluation routines.

allocate(fcn_data % rdata (7), fcn_data % idata (1))

! Define an equally-spaced grid of points for the underlying price

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

xgrid(1) = x_min

xgrid(nx) = x_max

do i = 2,nx-1

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

end do

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

! Flag the difference in payoff functions -

! 1 for the Bet, 2 for the Vertical Spread

fcn_data % idata(1) = 1

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

fkinitcond_call, fkbc_call, yb, ybprime,&

FCN_DATA = fcn_data)

fcn_data % idata(1) = 2

call feynman_kac (Xgrid, time, nlbc, nrbc, fkcoef_call,&

fkinitcond_call, fkbc_call, yv, yvprime,&

FCN_DATA = fcn_data)

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

! prior to expiration.

do i=1,nt

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

fv(:,i) = hqsval (xs, xgrid, yv(:,I))

end do

write(nout,'(T05,A)') "European Option Value for A Bet",&

" and a Vertical Spread, 3 and 6 Months "//&

"Prior to Expiry"

write(nout,'(T08,"Number of equally spaced spline knots "'//&

',I4,/T08,"Number of unknowns ",I4)') NX,N

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

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

'/T08,"Bet = ",F5.2,", Spread Value = ", F5.2/'//&

'/T10,"Underlying", T32,"A Bet", T40,"Vertical Spread",'//&

'/(T10,5F9.4))') ks1, sigma, r, bet, ks2, &

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

end

! These routines define the coefficients, payoff, boundary

! conditions and forcing term for American and European Options.

function fkcoef_call (x, tx, iflag, fcn_data) result(value)

use mp_types

implicit none

integer, intent(inout) :: iflag

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

type(s_fcn_data), optional :: fcn_data

real(kind(1e0)) :: value

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

! Data passed through using allocated components

! of the derived type s_fcn_data

sigma = fcn_data % rdata(5)

interest_rate = fcn_data % rdata(6)

dividend = fcn_data % rdata(7)

select case (iflag)

case (1)

! The coefficient derivative d(sigma)/dx

value = sigma

! The coefficient sigma(x)

case (2)

value = sigma * x

case (3)

! The coefficient mu(x)

value = (interest_rate - dividend) * x

case (4)

! The coefficient kappa(x)

value = interest_rate

end select

! Note that there is no time dependence

iflag = 0

return

end function fkcoef_call

function fkinitcond_call(x, fcn_data) result(value)

use mp_types

implicit none

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

type(s_fcn_data), optional :: fcn_data

real(kind(1e0)) :: value

real(kind(1e0)) :: strike_price, spread, bet

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

strike_price = fcn_data % rdata(1)

bet = fcn_data % rdata(2)

spread = fcn_data % rdata(3)

! The payoff function - Use flag passed to decide which

select case (fcn_data % idata(1))

case(1)

! After reaching the strike price the payoff jumps

! from zero to the bet value.

value = zero

if (x > strike_price) value = bet

case(2)

! Function is zero up to strike price.

! Then linear between strike price and spread.

! Then has constant value Spread-Strike Price after

! the value Spread.

value = max(x-strike_price, zero) - max(x-spread, zero)

end select

return

end function fkinitcond_call

subroutine fkbc_call (TX, iflag, bccoefs, fcn_data)

use mp_types

implicit none

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

integer, intent(inout) :: iflag

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

type(s_fcn_data), optional :: fcn_data

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

interest_rate, df

strike_price = fcn_data % rdata(1)

bet = fcn_data % rdata(2)

spread = fcn_data % rdata(3)

interest_rate = fcn_data % rdata(6)

select case (iflag)

case (1)

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

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

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

case (2)

! This is the discount factor using the risk-free

! interest rate

df = exp(interest_rate * tx)

! Use flag passed to decide on boundary condition -

select case (fcn_data % idata(1))

case(1)

bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, bet*df/)

case(2)

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

(spread-strike_price)*df/)

end select

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

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

return

end select

! Note no time dependence in case (1) for iflag

iflag = 0

end subroutine fkbc_call

Output

European Option Value for A Bet

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

Number of equally spaced spline knots 61

Number of unknowns 183

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

Bet = 2.00, Spread Value = 15.00

Underlying A Bet Vertical Spread

2.0000 0.0000 0.0000 0.0000 0.0000

4.0000 0.0000 0.0014 0.0000 0.0006

6.0000 0.0110 0.0723 0.0039 0.0447

8.0000 0.2691 0.4302 0.1478 0.3832

10.0000 0.9948 0.9781 0.8909 1.1926

12.0000 1.6094 1.4290 2.1911 2.2273

14.0000 1.8655 1.6922 3.4254 3.1553

16.0000 1.9338 1.8175 4.2263 3.8264

18.0000 1.9476 1.8700 4.6264 4.2492

20.0000 1.9501 1.8904 4.7911 4.4921

22.0000 1.9505 1.8979 4.8497 4.6231

24.0000 1.9506 1.9007 4.8685 4.6909

Example 4– Convertible Bonds

This example evaluates the price of a convertible bond. Here, convertibility means that the bond may, at any time of the holder’s choosing, be converted to a multiple of the specified asset. Thus a convertible bond with 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:

Bond face value K = {1}, conversion factor ν = 1.125

Volatility σ = {0.25}

Times until expiration = {1/2,1}

Interest rate r = 0.1, dividend D = 0.02

xmin = 0, xmax = 4

nx = 61, n = 3 × nx = 183

Boundary conditions

Terminal data

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

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

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(:),&

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

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

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

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

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 {S, f (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

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

! 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

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

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

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