|
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.
XGRID Rank-1
array containing the set of breakpoints that define the end points for
the Hermite quintic splines. (Input)
Let m = size(XGRID). The points in
XGRID must be in
strictly increasing order, and m 2.
TGRID Rank-1
array containing the set of time points (in time-remaining units) at which an
approximate solution is computed. (Input)
Let n =
size(TGRID). The
points in TGRID
must be strictly positive and in strictly increasing order and n
≥
1.
NLBC The number
of left boundary conditions. (Input)
1≤ NLBC ≤ 3.
NRBC The number
of right boundary conditions. (Input)
1≤ NRBC ≤ 3.
FKCOEF
User-supplied FUNCTION to evaluate
the coefficients and
of the Feynman-Kac PDE. The usage is FKCOEF (X, TX, IFLAG[,
]), where
Function Return Value
FKCOEF Value of the coefficient function. Which value is computed depends on the input value for IFLAG, see description of IFLAG.
Required Arguments
X Point in the x-space at which the coefficient is to be evaluated. (Input)
TX Time point at which the coefficient is to be evaluated. (Input)
IFLAG Flag related to the
coefficient that has to be computed (Input/Output).
On entry,
IFLAG indicates
which coefficient is to be computed. The following table shows which value has
to be returned by FKCOEF for all
possible values of IFLAG:
IFLAG |
Computed coefficient |
1 |
|
2 |
|
3 |
|
4 |
|
One indicates when a coefficient does not depend on t by setting IFLAG = 0 after the coefficient is defined. If there is time dependence, the value of IFLAG should not be changed. This action will usually yield a more efficient algorithm because some finite element matrices do not have to be reassembled for each t value.
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied function. For a detailed description of this argument see FCN_DATA below. (Input/Output)
FKCOEF must be declared EXTERNAL in the calling program.
FKINITCOND
User-supplied FUNCTION to evaluate
the initial condition function in the Feynman-Kac
PDE. The usage is FKINITCOND (X[,
]), where
Function Return Value
FKINITCOND Value of the initial
condition function .
Required Arguments
X Point in the x-space at which the initial condition is to be evaluated. (Input)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied function. For a detailed description of this argument see FCN_DATA below. (Input/Output)
FKINITCOND must be declared EXTERNAL in the calling program.
FKBC
User-supplied subroutine to evaluate
the coefficients for the left and right boundary conditions the Feynman-Kac PDE
must satisfy. There are NLBC conditions
specified at the left end, , and NRBC conditions at the
right end,
. The boundary conditions
can be vectors of dimension 1, 2 or 3 and are defined by
The usage is FKBC (TX, IFLAG, BCCOEFS[, ]) where
Required Arguments
TX Time point at which the coefficients are to be evaluated. (Input)
IFLAG Flag related
to the boundary conditions that have to be computed (Input/Output).
On input,
IFLAG indicates
whether the coefficients for the left or right boundary conditions have to be
computed:
IFLAG |
Computed boundary conditions |
1 |
Left end, x = xmin |
2 |
Right end, x = xmax |
If there is no time dependence for one of the boundaries then set IFLAG = 0 after the array BCCOEFS is defined for either end point. This can avoid unneeded continued computation of the finite element matrices.
BCCOEFS Array
of size 3 Χ 4 containing the coefficients of the left or right boundary
conditions in its first NLBC or NRBC rows,
respectively. (Output)
The coefficients for are stored row-wise according to the following
matrix-scheme:
The coefficients
for are stored similarly.
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a detailed description of this argument see FCN_DATA below. (Input/Output)
FKBC must be declared EXTERNAL in the calling program.
Y Array of size
(3*m) by
(n+1) containing the coefficients of the Hermite representation of the
approximate solution for the Feynman-Kac PDE at time points (in time-remaining
units) 0,
TGRID(1),
, TGRID(n). (Output)
For
, the coefficients are
stored in columns 1,
,n of array Y and the
approximate solution is given by
.
The coefficients of the representation for the initial data are given in column 0 of array Y and are defined by
.
The starting coefficients are estimated using
least-squares.
After the integrations, use Y(:,0) and Y(:,j) as
input argument COEFFS to function
HQSVAL to obtain
an array of values for f(x, t) or its partials at time points t=0 and t=TGRID(j), j=1,
,n, respectively.
The expressions for the basis functions are represented piece-wise and can be found in Hanson, R. (2008)
Integrating
Feynman-Kac Equations Using Hermite Quintic Finite Elements.
YPRIME Array of
size (3*m) by (n + 1) containing the
first derivatives of the coefficients of the Hermite representation of the
approximate solution for the Feynman-Kac PDE at time points (in time-remaining
units) 0,
TGRID(1),
, TGRID(n).
(Output)
For and
, the derivatives of the coefficients are stored in column 0 and
columns 1 to n of array YPRIME, respectively.
The columns in YPRIME represent
for
,
and
for
.
After the integrations, use YPRIME(:,j) as input
argument COEFFS
to function HQSVAL
to obtain an array of values for the partials at time points t = TGRID(j), j
= 1,
,n, and YPRIME(:,0) for the partials at
t =
0.
FKINIT User-supplied subroutine that allows for adjustment of initial data or as an opportunity for output during the integration steps.
The usage is CALL FKINIT (XGRID, TGRID, TX, YPRIME, Y, ATOL, RTOL[, ]) where
Required Arguments
XGRID Array of size m containing the set of breakpoints that define the end points for the Hermite quintic splines. (Input)
TGRID Array of size n containing the set of time points (in time-remaining units) at which an approximate solution is computed. (Input)
TX Time point for
the evaluation. (Input)
Possible values are 0 (the initial or
terminal time point) and all values in array TGRID.
YPRIME Array of length 3*m containing the derivatives of the Hermite quintic spline coefficients at time point TX. (Input)
Y Array of length
3* m containing the Hermite quintic spline
coefficients at time point TX.
(Input/Output)
For the initial time point TX=0 this array can
be used to reset the Hermite quintic spline coefficients
to user defined values. For all other values of TX array Y is an input array.
ATOL Array of length 3* m containing absolute error tolerances used in the integration routine that determines the Hermite quintic spline coefficients and its derivatives. (Input/Output)
RTOL Array of length 3*m containing relative error tolerances used in the integration routine that determines the Hermite quintic spline coefficients and its derivatives. (Input/Output)
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied function. For a detailed description of this argument see FCN_DATA below. (Input/Output)
FKINIT must be declared EXTERNAL in the calling program.
FKFORCE User-supplied subroutine that computes local contributions
and
.
The usage is CALL FKFORCE (I, T, WIDTH, Y, XLOCAL, QW, U, PHI, DPHI[,, ]) where
Required Arguments
I Index related to the integration interval (XGRID(I), XGRID(I+1)). (Input)
T Time point at which the local contributions are computed. (Input)
WIDTH Width of the integration interval I, WIDTH= XGRID(I+1)- XGRID(I). (Input)
Y
Array of length 3*m containing the coefficients of the Hermite
quintic spline representing the solution of the Feynman-Kac PDE at time
point T.
(Input)
For each
,
the approximate solution is locally defined by
Here, the functions
are basis polynomials
of order 5 and
.
The values
,
are stored as successive triplets in array Y.
XLOCAL Array
containing the Gauss-Legendre points translated and normalized to the
interval
[XGRID(I),XGRID(I+1)]. (Input)
The size of the array is
equal to the degree of the Gauss-Legendre polynomials used for constructing the
finite element matrices.
QW Array containing
the Gauss-Legendre weights. (Input)
The size of the array is
equal to the degree of the Gauss-Legendre polynomials used for constructing the
finite element matrices.
U Array of size
size(XLOCAL) Χ
12 containing the basis function values that define at the Gauss-Legendre points XLOCAL.
(Input)
Let
.
Using the local approximation in the I-th interval, defined by
,
and setting
,
and
,
vector is defined as
PHI Array of size 6
containing a Gauss-Legendre approximation for the local contribution, where
T and
. (Output)
Setting NDEG:=SIZE(XLOCAL) and
, array PHI contains
elements
PHI(i) = WIDTH
for i= 1, , 6.
DPHI Array of size
6Χ6, a Gauss-Legendre approximation for the Jacobian of the local contribution
at time point
= T,
. (Output)
The approximation to this symmetric matrix is stored row-wise, i.e.
DPHI(i,j) = WIDTH
or i, j = 1, ,6.
Optional Arguments
FCN_DATA A derived type, s_fcn_data, which may be used to pass additional integer or floating point information to or from the user-supplied subroutine. For a detailed description of this argument see FCN_DATA below. (Input/Output)
FKFORCE must be declared EXTERNAL in the calling program.
If subroutine FKFORCE is not used as
an optional argument then it is assumed that the forcing term in the Feynman-Kac equation is identically zero.
ATOL Array of
non-negative values containing absolute error tolerances used in the computation
of each column of solution array Y via integration
routine DASPH.
(Input)
The size of array ATOL can be 1 or
3Χm. In the first
case, ATOL(1:1)
is applied to
all solution components, in the latter each component of ATOL is assigned to
the corresponding solution component allowing for individual control of the
error tolerances. At least one entry in arrays ATOL or RTOL must be greater
than 0.
Default: ATOL(1:1) = 1.0e-3 for
single and 1.0d-5 for double precision.
RTOL Array of
non-negative values containing relative error tolerances used in the computation
of each column of solution array Y via integration
routine DASPH.
(Input)
The size of array RTOL can be 1 or
3Χm. In the first case, RTOL(1:1) is applied to all
solution components, in the latter each component of RTOL is assigned to
the corresponding solution component allowing for individual control of the
error tolerances. At least one entry in arrays ATOL or RTOL must be greater
than 0.
Default: RTOL(1:1) = 1.0e-3 for
single and 1.0d-5 for double precision.
NDEG Degree of
the Gauss-Legendre formulas used for constructing the finite element
matrices. (Input)
NDEG ≥ 6.
Default:
NDEG = 6.
RINITSTEPSIZE
Starting step size for the integration. (Input)
RINITSTEPSIZE must be
strictly negative since the integration is internally done from T = 0 to T = TGRID(n) in a negative direction.
Default:
Program defined initial stepsize.
MAXBDFORDER
Maximum order of the backward differentiation formulas (BDF) used in the
integrator DASPH.
(Input)
1 ≤ MAXBDFORDER ≤ 5.
Default: MAXBDFORDER = 5.
RMAXSTEPSIZE
Maximum step size the integrator may take. (Input)
RMAXSTEPSIZE must be
strictly positive.
Default: RMAXSTEPSIZE = AMACH(2), the largest
possible machine number.
MAXIT Maximum
number of internal integration steps between two consecutive time points in
TGRID.
(Input)
MAXIT must be strictly
positive.
Default: MAXIT = 500000.
IMETHSTEPCTRL Indicates
which step control algorithm is used in the integration.
(Input)
If IMETHSTEPCTRL = 0,
then the step control method of Sφderlind is used. If IMETHSTEPCTRL =1, then
the method used by the original Petzold code SASSL is used.
IMETHSTEPCTRL |
Method used |
0 |
Method of Sφderlind.. |
1 |
Method from Petzold code SASSL. |
Default: IMETHSTEPCTRL = 0.
TBARRIER Time
barrier past which the integration routine DASPH will not go
during integration. (Input)
TBARRIER ≥ TGRID(n).
Default: TBARRIER = TGRID(n).
ISTATE Array of
size 5 whose entries flag the state of computation for the matrices and vectors
required in the integration. (Output)
For each entry, a zero indicates
that no computation has been done or that there is a time dependence. A one
indicates that the entry has been computed and there is no time
dependence.
The ISTATE entries are as
follows:
I |
ISTATE(I) |
1 |
State of computation of Mass matrix, M. |
2 |
State of computation of Stiffness matrix, N. |
3 |
State of computation of Bending matrix, R. |
4 |
State of computation of Weighted mass matrix, K. |
5 |
State of computation of initial data. |
NVAL Array of size 3 summarizing the number of evaluations required during the integration. (Output)
I |
NVAL(I) |
1 |
Number of residual function evaluations |
2 |
Number of factorizations of the differential |
3 |
Number of linear system solve steps using |
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.
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.
The generalized Feynman-Kac differential equation has the form
,
where the initial data satisfies
.
The derivatives are etc.
FEYNMAN_KAC uses a finite element Galerkin method over the rectangle[3]
in to compute the
approximate solution. The interval
is decomposed with a grid
.
On each subinterval the solution is represented by
The values are time-dependent
coefficients associated with each interval. The basis functions
are given for
,
by
The Galerkin principle is then applied. Using the provided initial and boundary conditions leads to an index 1 differential-algebraic equation (DAE) for the time-dependent coefficients
.
This system is integrated using the variable order, variable step algorithm DASPH. Solution values and their time derivatives are returned at a grid preceding time T, expressed in units of time remaining.
More mathematical details are found in Hanson, R. (2008) Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements.
In Beckers (1980) there is a model for a Stochastic Differential Equation of option pricing. The idea is a constant elasticity of variance diffusion (or CEV) class
The Black-Scholes model is the limiting case. A numerical solution of this diffusion model yields the
price of a call option. Various values of the strike price
, time values,
and power coefficient
are used to evaluate the option price at values of the
underlying price. The sets of parameters in the computation are:
1.
power
2.
strike price
3.
volatility
4.
times until expiration
5.
underlying prices
6.
interest rate
7.
8.
With this model the Feynman-Kac differential equation is defined by identifying:
The payoff function is the vanilla option, .
Link to example source (feynman_kac_ex1.f90)
! Compute Constant Elasticity of Variance Model for Vanilla Call
use feynman_kac_int
use hqsval_int
use mp_types
use umach_int
implicit none
! The set of strike prices
real(kind(1e0)) :: ks(3)=(/15.0e0,20.0e0,25.0e0/)
! The set of sigma values
real(kind(1e0)) :: sigma(3) = (/0.2e0, 0.3e0, 0.4e0/)
! The set of model diffusion powers
real(kind(1e0)) :: alpha(3) = (/2.0e0,1.0e0,0.0e0/)
! Time values for the options
integer, parameter :: nt = 3
real(kind(1e0)) :: time(nt)=(/1.e0/12., 4.e0/12., 7.e0/12./)
! Values of the underlying where evaluation are made
integer, parameter :: nv = 3, nlbc = 3, nrbc = 3
real(kind(1e0)) :: xs(nv) = (/19.0e0,20.0e0,21.0e0/)
! Value of the interest rate and continuous dividend
real(kind(1e0)) :: r = 0.05e0, dividend = 0.0e0
! Values of the min and max underlying values modeled
real(kind(1e0)) :: x_min = 0.0e0, x_max = 60.0e0
! Define parameters for the integration step.
integer, parameter :: nx = 121, nint = nx-1, n = 3*nx
real(kind(1e0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt),&
dx, f(nv,nt)
type(s_fcn_data) fcn_data
integer :: nout
real(kind(1e0)), external :: fkcoef, fkinitcond
external fkbc
integer :: i,i1,i2,i3,j
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fcn_data % rdata (6))
! Define an equally-spaced grid of points for the underlying price
dx = (x_max-x_min)/real(nint)
xgrid(1) = x_min
xgrid(nx) = x_max
do i = 2,nx-1
xgrid(i) = xgrid(i-1) + dx
end do
call umach(2, nout)
write(nout,'(T05,A)') "Constant Elasticity of Variance Model "//&
"for Vanilla Call"
write(nout,'(T10,"Interest Rate ", F7.3, T38,"Continuous '//&
'Dividend ", F7.3 )') r, dividend
write(nout,'(T10,"Minimum and Maximum Prices of Underlying ",'//&
'2F7.2)') x_min, x_max
write(nout,'(T10,"Number of equally spaced spline knots ",I4,'//&
'/T10,"Number of unknowns ",I4)')&
nx-1,n
write(nout,'(/T10,"Time in Years Prior to Expiration ",2X,'//&
'3F7.4)') time
write(nout,'(T10,"Option valued at Underlying Prices ",'//&
'3F7.2)') xs
do i1 = 1,3 ! Loop over power
do i2=1,3 ! Loop over volatility
do i3=1,3 ! Loop over strike price
! Pass data through into evaluation routines.
fcn_data % rdata =&
(/ks(i3),x_max,sigma(i2),alpha(i1),r,dividend/)
call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef,&
fkinitcond, fkbc, y, yprime,&
FCN_DATA = fcn_data)
! Evaluate solution at vector of points XS(:), at each time value
! prior to expiration.
do i=1,nt
f(:,i) = hqsval (xs, xgrid, y(:,i))
end do
write(nout,'(/T05,"Strike=",F5.2," Sigma=", F5.2,'//&
'" Alpha=", F5.2,/(T25," Call Option Values ",'//&
'X,3F7.4))') ks(I3),sigma(I2),&
alpha(i1),(f(i,:),i=1,nv)
end do !i3 - Strike price loop
end do !i2 - Sigma loop
end do !i1 - Alpha loop
end
! These functions and routines define the coefficients, payoff
! and boundary conditions.
function fkcoef (X, TX, iflag, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: X, TX
integer, intent(inout) :: iflag
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: fkcoef
real(kind(1e0)) :: sigma, interest_rate, alpha, dividend,&
zero = 0.0e0, half = 0.5e0
sigma = fcn_data % rdata(3)
alpha = fcn_data % rdata(4)
interest_rate = fcn_data % rdata(5)
dividend = fcn_data % rdata(6)
select case (iflag)
case (1)
! The coefficient derivative d(sigma)/dx
fkcoef = half*alpha*sigma*x**(alpha*half-1.0e0)
! The coefficient sigma(x)
case (2)
fkcoef = sigma*x**(alpha*half)
case (3)
! The coefficient mu(x)
fkcoef = (interest_rate - dividend) * x
case (4)
! The coefficient kappa(x)
fkcoef = interest_rate
end select
! Note that there is no time dependence
iflag = 0
return
end function fkcoef
function fkinitcond(x, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: x
type (s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: fkinitcond
real(kind(1e0)) :: zero = 0.0e0
real(kind(1e0)) :: strike_price
strike_price = fcn_data % rdata(1)
! The payoff function
fkinitcond = max(x - strike_price, zero)
return
end function fkinitcond
subroutine fkbc (tx, iflag, bccoefs, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
type (s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: x_max, df, interest_rate, strike_price
strike_price = fcn_data % rdata(1)
x_max = fcn_data % rdata(2)
interest_rate = fcn_data % rdata(5)
select case (iflag)
case (1)
bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, 0.0e0/)
bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 0.0e0/)
bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
! Note no time dependence at left end
iflag = 0
case (2)
df = exp(interest_rate * tx)
bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,&
x_max - df*strike_price/)
bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 1.0e0/)
bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
end select
end subroutine fkbc
Constant Elasticity of Variance Model for Vanilla Call
Interest Rate 0.050 Continuous Dividend 0.000
Minimum and Maximum Prices of Underlying 0.00 60.00
Number of equally spaced spline knots 120
Number of unknowns 363
Time in Years Prior to Expiration 0.0833 0.3333 0.5833
Option valued at Underlying Prices 19.00 20.00 21.00
Strike=15.00 Sigma= 0.20 Alpha= 2.00
Call Option Values 4.0624 4.2575 4.4730
Call Option Values 5.0624 5.2506 5.4490
Call Option Values 6.0624 6.2486 6.4385
Strike=20.00 Sigma= 0.20 Alpha= 2.00
Call Option Values 0.1310 0.5955 0.9699
Call Option Values 0.5018 1.0887 1.5101
Call Option Values 1.1977 1.7483 2.1752
Strike=25.00 Sigma= 0.20 Alpha= 2.00
Call Option Values 0.0000 0.0112 0.0745
Call Option Values 0.0000 0.0372 0.1621
Call Option Values 0.0007 0.1027 0.3141
Strike=15.00 Sigma= 0.30 Alpha= 2.00
Call Option Values 4.0637 4.3398 4.6622
Call Option Values 5.0626 5.2944 5.5786
Call Option Values 6.0624 6.2708 6.5240
Strike=20.00 Sigma= 0.30 Alpha= 2.00
Call Option Values 0.3109 1.0276 1.5494
Call Option Values 0.7326 1.5424 2.1017
Call Option Values 1.3765 2.1690 2.7379
Strike=25.00 Sigma= 0.30 Alpha= 2.00
Call Option Values 0.0006 0.1112 0.3543
Call Option Values 0.0038 0.2169 0.5548
Call Option Values 0.0184 0.3857 0.8222
Strike=15.00 Sigma= 0.40 Alpha= 2.00
Call Option Values 4.0755 4.5138 4.9675
Call Option Values 5.0662 5.4201 5.8326
Call Option Values 6.0634 6.3579 6.7301
Strike=20.00 Sigma= 0.40 Alpha= 2.00
Call Option Values 0.5115 1.4640 2.1273
Call Option Values 0.9621 1.9951 2.6929
Call Option Values 1.5814 2.6105 3.3216
Strike=25.00 Sigma= 0.40 Alpha= 2.00
Call Option Values 0.0083 0.3286 0.7790
Call Option Values 0.0285 0.5167 1.0657
Call Option Values 0.0813 0.7687 1.4103
Strike=15.00 Sigma= 0.20 Alpha= 1.00
Call Option Values 4.0624 4.2479 4.4311
Call Option Values 5.0624 5.2479 5.4311
Call Option Values 6.0624 6.2479 6.4311
Strike=20.00 Sigma= 0.20 Alpha= 1.00
Call Option Values 0.0000 0.0218 0.1045
Call Option Values 0.1498 0.4109 0.6485
Call Option Values 1.0832 1.3314 1.5773
Strike=25.00 Sigma= 0.20 Alpha= 1.00
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Strike=15.00 Sigma= 0.30 Alpha= 1.00
Call Option Values 4.0624 4.2477 4.4309
Call Option Values 5.0624 5.2477 5.4309
Call Option Values 6.0624 6.2477 6.4309
Strike=20.00 Sigma= 0.30 Alpha= 1.00
Call Option Values 0.0011 0.0781 0.2201
Call Option Values 0.1994 0.5000 0.7543
Call Option Values 1.0835 1.3443 1.6023
Strike=25.00 Sigma= 0.30 Alpha= 1.00
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0005
Strike=15.00 Sigma= 0.40 Alpha= 1.00
Call Option Values 4.0624 4.2479 4.4312
Call Option Values 5.0624 5.2479 5.4312
Call Option Values 6.0624 6.2479 6.4312
Strike=20.00 Sigma= 0.40 Alpha= 1.00
Call Option Values 0.0076 0.1563 0.3452
Call Option Values 0.2495 0.5907 0.8706
Call Option Values 1.0868 1.3779 1.6571
Strike=25.00 Sigma= 0.40 Alpha= 1.00
Call Option Values 0.0000 0.0000 0.0001
Call Option Values 0.0000 0.0000 0.0008
Call Option Values 0.0000 0.0003 0.0063
Strike=15.00 Sigma= 0.20 Alpha= 0.00
Call Option Values 4.0626 4.2479 4.4311
Call Option Values 5.0623 5.2480 5.4311
Call Option Values 6.0624 6.2480 6.4312
Strike=20.00 Sigma= 0.20 Alpha= 0.00
Call Option Values 0.0001 0.0001 0.0002
Call Option Values 0.0816 0.3316 0.5748
Call Option Values 1.0818 1.3308 1.5748
Strike=25.00 Sigma= 0.20 Alpha= 0.00
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Strike=15.00 Sigma= 0.30 Alpha= 0.00
Call Option Values 4.0625 4.2479 4.4312
Call Option Values 5.0623 5.2479 5.4312
Call Option Values 6.0624 6.2479 6.4312
Strike=20.00 Sigma= 0.30 Alpha= 0.00
Call Option Values 0.0000 0.0000 0.0029
Call Option Values 0.0894 0.3326 0.5753
Call Option Values 1.0826 1.3306 1.5749
Strike=25.00 Sigma= 0.30 Alpha= 0.00
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Strike=15.00 Sigma= 0.40 Alpha= 0.00
Call Option Values 4.0624 4.2479 4.4312
Call Option Values 5.0623 5.2479 5.4312
Call Option Values 6.0624 6.2479 6.4312
Strike=20.00 Sigma= 0.40 Alpha= 0.00
Call Option Values 0.0000 0.0002 0.0113
Call Option Values 0.0985 0.3383 0.5781
Call Option Values 1.0830 1.3306 1.5749
Strike=25.00 Sigma= 0.40 Alpha= 0.00
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
Call Option Values 0.0000 0.0000 0.0000
The value of the American Option on a Vanilla Put can be no smaller than its European counterpart. That is due to the American Option providing the opportunity to exercise at any time prior to expiration. This example compares this difference or premium value of the American Option at two time values using the Black-Scholes model. The example is based on Wilmott et al. (1996, p. 176), and uses the non-linear forcing or weighting term described in Hanson, R. (2008), Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements, for evaluating the price of the American Option. A call to the subroutine fkinit_put sets the initial conditions. One breakpoint is set exactly at the strike price.
The sets of parameters in the computation are:
1.
Strike price
2.
Volatility
3.
Times until expiration
4.
Interest rate
5.
6.
The payoff function is the vanilla option, .
Link to example source (feynman_kac_ex2.f90)
! Compute American Option Premium for Vanilla Put
use feynman_kac_int
use hqsval_int
use mp_types
use umach_int
implicit none
! The strike price
real(kind(1e0)) :: ks = 10.0e0
! The sigma value
real(kind(1e0)) :: sigma = 0.4e0
! Time values for the options
integer, parameter :: nt = 2
real(kind(1e0)) :: time(nt)=(/0.25e0, 0.5e0/)
! Values of the underlying where evaluations are made
integer, parameter :: nv = 9
integer, parameter :: nlbc = 2, nrbc = 3, ndeg = 6
integer :: i
real(kind(1e0)) :: xs(nv) = (/((i-1)*2.0e0,i=1,nv)/)
! Value of the interest rate and continuous dividend
real(kind(1e0)) :: r = 0.1e0, dividend = 0.0e0
! Values of the min and max underlying values modeled
real(kind(1e0)) :: x_min = 0.0e0, x_max = 30.0e0
real(kind(1e0)) :: atol(1), rtol(1)
! Define parameters for the integration step.
integer, parameter :: nx = 121, nint = nx-1, n = 3*nx
real(kind(1e0)) :: xgrid(nx), ye(n,0:nt), yeprime(n,0:nt),&
ya(n,0:nt), yaprime(n,0:nt),&
dx, fe(nv,nt), fa(nv,nt)
type(s_fcn_data) fcn_data
integer :: nout
real(kind(1e0)), external :: fkcoef_put, fkinitcond_put
external fkbc_put, fkinit_put, fkforce_put
call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fcn_data % rdata (6), fcn_data % idata (1))
! Define an equally-spaced grid of points for the underlying price
dx = (x_max-x_min)/real(nint)
xgrid(1) = x_min
xgrid(nx) = x_max
do i=2,nx-1
xgrid(i) = xgrid(i-1) + dx
end do
! Place a breakpoint at the strike price.
do i = 1,nx
if (xgrid(i) > ks) then
xgrid(i-1) = ks
exit
end if
end do
! Request less accuracy than the default values provide.
atol(1) = 0.5e-2
rtol(1) = 0.5e-2
fcn_data % rdata = (/ks,x_max,sigma,r,dividend,atol(1)/)
fcn_data % idata = (/ndeg/)
! Compute European then American Put Option Values.
call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&
fkinitcond_put, fkbc_put, ye, yeprime,&
FKINIT=fkinit_put, ATOL=atol,RTOL=rtol,&
FCN_DATA = fcn_data)
call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&
fkinitcond_put, fkbc_put, ya, yaprime,&
FKINIT=fkinit_put, ATOL=atol, RTOL=rtol,&
FKFORCE=fkforce_put, FCN_DATA = fcn_data)
! Evaluate solutions at vector of points XS(:), at each time value
! prior to expiration.
do i=1,nt
fe(:,i) = hqsval (xs, xgrid, ye(:,I))
fa(:,I) = hqsval (xs, xgrid, ya(:,I))
end do
write(nout,'(T05,A,/,T05,A)')&
"American Option Premium for Vanilla Put, 3 and 6 Months "//&
"Prior to", "Expiry"
write(nout,'(T08,"Number of equally spaced spline knots ",I4,'//&
'/T08,"Number of unknowns ",I4)') nx,n
write(nout,'(T08,"Strike= ",F5.2,", Sigma=", F5.2,", Interest'//&
' Rate=",F5.2,/T08,"Underlying", T26,"European",'//&
'T42,"American",/(T10,5F8.4))') ks,sigma,r,&
(xs(i), fe(i,1:nt), fa(i,1:nt),i=1,nv)
end
! These routines define the coefficients, payoff, boundary
! conditions, forcing term and initial conditions for American and
! European Options.
function fkcoef_put(x, tx, iflag, fcn_data)
use mp_types
implicit none
integer, intent(inout) :: iflag
real(kind(1e0)), intent(in) :: x, tx
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: fkcoef_put
real(kind(1e0)) :: sigma, strike_price, interest_rate, &
dividend, zero=0.e0
sigma = fcn_data % rdata(3)
interest_rate = fcn_data % rdata(4)
dividend = fcn_data % rdata(5)
select case (iflag)
case (1)
! The coefficient derivative d(sigma)/dx
fkcoef_put = sigma
! The coefficient sigma(x)
case (2)
fkcoef_put = sigma*x
case (3)
! The coefficient mu(x)
fkcoef_put = (interest_rate - dividend)*x
case (4)
! The coefficient kappa(x)
fkcoef_put = interest_rate
end select
! Note that there is no time dependence
iflag = 0
return
end function fkcoef_put
function fkinitcond_put(x, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: x
type (s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: fkinitcond_put
real(kind(1e0)) :: zero = 0.0e0
real(kind(1e0)) :: strike_price
strike_price = fcn_data % rdata(1)
! The payoff function
fkinitcond_put = max(strike_price - x, zero)
return
end function fkinitcond_put
subroutine fkbc_put (tx, iflag, bccoefs, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
type (s_fcn_data), optional :: fcn_data
select case (iflag)
case (1)
bccoefs(1,1:4) = ((/0.0e0, 1.0e0, 0.0e0, -1.0e0/))
bccoefs(2,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
case (2)
bccoefs(1,1:4) = ((/1.0e0, 0.0e0, 0.0e0, 0.0e0/))
bccoefs(2,1:4) = ((/0.0e0, 1.0e0, 0.0e0, 0.0e0/))
bccoefs(3,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
end select
! Note no time dependence
iflag = 0
end subroutine fkbc_put
subroutine fkforce_put (interval, t, hx, y, xlocal, qw, u,&
phi, dphi, fcn_data)
use mp_types
implicit none
integer, parameter :: local = 6
integer :: i, j, l, ndeg
integer, intent(in) :: interval
real(kind(1e0)), intent(in) :: y(:), t, hx, qw(:),&
xlocal(:), u(:,:)
real(kind(1e0)), intent(out) :: phi(:), dphi(:,:)
type (s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: yl(local), bf(local)
real(kind(1e0)) :: value, strike_price, interest_rate,&
zero=0.0e0, one=1.0e0, rt, mu
yl = y(3*interval-2:3*interval+3)
phi = zero
value = fcn_data % rdata(6)
strike_price = fcn_data % rdata(1)
interest_rate = fcn_data % rdata(4)
ndeg = fcn_data % idata(1)
mu = 2
! This is the local definition of the forcing term
do j=1,local
do l=1,ndeg
bf(1:3) = u(l,1:3)
bf(4:6) = u(l,7:9)
rt = dot_product(yl,bf)
rt = value/(rt + value - (strike_price - xlocal(l)))
phi(j) = phi(j) + qw(l) * bf(j) * rt**mu
end do
end do
phi = -phi*hx*interest_rate*strike_price
! This is the local derivative matrix for the forcing term
dphi = zero
do j =1,local
do i = 1,local
do l=1,ndeg
bf(1:3) = u(l,1:3)
bf(4:6) = u(l,7:9)
rt = dot_product(yl,bf)
rt = one/(rt + value - (strike_price - xlocal(l)))
dphi(i,j) = dphi(i,j) + qw(l) * bf(I) * bf(j) *&
rt**(mu+1)
end do
end do
end do
dphi = mu*dphi*hx*value**mu*interest_rate*strike_price
return
end subroutine fkforce_put
subroutine fkinit_put(xgrid,tgrid,t,yprime,y,atol,rtol,&
fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&
yprime(:)
real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)
type (s_fcn_data), optional :: fcn_data
integer :: i
if (t == 0.0e0) then
! Set initial data precisely. The strike price is a breakpoint.
! Average the derivative limit values from either side.
do i=1,size(xgrid)
if (xgrid(i) < fcn_data % rdata(1)) then
y(3*i-2) = fcn_data % rdata(1) - xgrid(i)
y(3*i-1) = -1.0e0
y(3*i)= 0.0e0
else if (xgrid(i) == fcn_data % rdata(1)) then
y(3*i-2) = 0.0e0
y(3*i-1) = -0.5e0
y(3*i) = 0.0e0
else
y(3*i-2) = 0.0e0
y(3*i-1) = 0.0e0
y(3*i) = 0.0e0
end if
end do
end if
end subroutine fkinit_put
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
This example evaluates the price of a European Option using two payoff strategies: Cash-or-Nothing and Vertical Spread. In the first case the payoff function is
.
The value B is regarded as the bet on the asset price, see Wilmott et al. (1995, p. 39-40). The second case has the payoff function
Both problems use the same boundary conditions. Each case requires a separate integration of the Black-Scholes differential equation, but only the payoff function evaluation differs in each case. The sets of parameters in the computation are:
1. Strike and bet prices K1={10.0}, K2 = {15.0}, and B = {2.0}
2. Volatility σ = {0.4}.
3. Times until expiration = {1/4, 1/2}.
4. Interest rate r = 0.1.
5. xmin = 0, xmax = 30.
6. nx =121, n = 3 Χ nx = 363.
Link to example source (feynman_kac_ex3.f90)
! Compute European Option Premium for a Cash-or-Nothing
! and a Vertical Spread Call.
use feynman_kac_int
use hqsval_int
use mp_types
use umach_int
implicit none
! The strike price
real(kind(1e0)) :: ks1 = 10.0e0
! The spread value
real(kind(1e0)) :: ks2 = 15.0e0
! The Bet for the Cash-or-Nothing Call
real(kind(1e0)) :: bet = 2.0e0
! The sigma value
real(kind(1e0)) :: sigma = 0.4e0
! Time values for the options
integer, parameter :: nt = 2
real(kind(1e0)) :: time(nt)=(/0.25e0, 0.5e0/)
! Values of the underlying where evaluation are made
integer, parameter :: nv = 12, nlbc = 3, nrbc = 3
integer :: i
real(kind(1e0)) :: xs(nv) = (/(2+(I-1)*2.0e0,I=1,NV)/)
! Value of the interest rate and continuous dividend -
real(kind(1e0)) :: r = 0.1e0, dividend = 0.0e0
! Values of the min and max underlying values modeled -
real(kind(1e0)) :: x_min = 0.0e0, x_max = 30.0e0
! Define parameters for the integration step.
integer, parameter :: nx = 61, nint = nx-1, n=3*nx
real(kind(1e0)) :: xgrid(nx), yb(n,0:nt), ybprime(n,0:nt),&
yv(n,0:nt), yvprime(n,0:nt),&
dx, fb(nv,nt), fv(nv,nt)
type(s_fcn_data) fcn_data
integer :: nout
real(kind(1e0)), external :: fkcoef_call, fkinitcond_call
external fkbc_call
call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fcn_data % rdata (7), fcn_data % idata (1))
! Define an equally-spaced grid of points for the underlying price
dx = (x_max-x_min)/real(nint)
xgrid(1) = x_min
xgrid(nx) = x_max
do i = 2,nx-1
xgrid(i) = xgrid(i-1) + dx
end do
fcn_data % rdata = (/ks1,bet,ks2,x_max,sigma,r,dividend/)
! Flag the difference in payoff functions -
! 1 for the Bet, 2 for the Vertical Spread
fcn_data % idata(1) = 1
call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_call,&
fkinitcond_call, fkbc_call, yb, ybprime,&
FCN_DATA = fcn_data)
fcn_data % idata(1) = 2
call feynman_kac (Xgrid, time, nlbc, nrbc, fkcoef_call,&
fkinitcond_call, fkbc_call, yv, yvprime,&
FCN_DATA = fcn_data)
! Evaluate solutions at vector of points XS(:), at each time value
! prior to expiration.
do i=1,nt
fb(:,i) = hqsval (xs, xgrid, yb(:,I))
fv(:,i) = hqsval (xs, xgrid, yv(:,I))
end do
write(nout,'(T05,A)') "European Option Value for A Bet",&
" and a Vertical Spread, 3 and 6 Months "//&
"Prior to Expiry"
write(nout,'(T08,"Number of equally spaced spline knots "'//&
',I4,/T08,"Number of unknowns ",I4)') NX,N
write(nout,'(T08,"Strike = ",F5.2,", Sigma =", F5.2,'//&
'", Interest Rate =",F5.2,'//&
'/T08,"Bet = ",F5.2,", Spread Value = ", F5.2/'//&
'/T10,"Underlying", T32,"A Bet", T40,"Vertical Spread",'//&
'/(T10,5F9.4))') ks1, sigma, r, bet, ks2, &
(xs(i), fb(i,1:nt), fv(i,1:nt),i=1,nv)
end
! These routines define the coefficients, payoff, boundary
! conditions and forcing term for American and European Options.
function fkcoef_call (x, tx, iflag, fcn_data) result(value)
use mp_types
implicit none
integer, intent(inout) :: iflag
real(kind(1e0)), intent(in) :: x, tx
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: value
real(kind(1e0)) :: sigma, interest_rate, dividend
! Data passed through using allocated components
! of the derived type s_fcn_data
sigma = fcn_data % rdata(5)
interest_rate = fcn_data % rdata(6)
dividend = fcn_data % rdata(7)
select case (iflag)
case (1)
! The coefficient derivative d(sigma)/dx
value = sigma
! The coefficient sigma(x)
case (2)
value = sigma * x
case (3)
! The coefficient mu(x)
value = (interest_rate - dividend) * x
case (4)
! The coefficient kappa(x)
value = interest_rate
end select
! Note that there is no time dependence
iflag = 0
return
end function fkcoef_call
function fkinitcond_call(x, fcn_data) result(value)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: x
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: value
real(kind(1e0)) :: strike_price, spread, bet
real(kind(1e0)), parameter :: zero = 0.0e0
strike_price = fcn_data % rdata(1)
bet = fcn_data % rdata(2)
spread = fcn_data % rdata(3)
! The payoff function - Use flag passed to decide which
select case (fcn_data % idata(1))
case(1)
! After reaching the strike price the payoff jumps
! from zero to the bet value.
value = zero
if (x > strike_price) value = bet
case(2)
! Function is zero up to strike price.
! Then linear between strike price and spread.
! Then has constant value Spread-Strike Price after
! the value Spread.
value = max(x-strike_price, zero) - max(x-spread, zero)
end select
return
end function fkinitcond_call
subroutine fkbc_call (TX, iflag, bccoefs, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: strike_price, spread, bet,&
interest_rate, df
strike_price = fcn_data % rdata(1)
bet = fcn_data % rdata(2)
spread = fcn_data % rdata(3)
interest_rate = fcn_data % rdata(6)
select case (iflag)
case (1)
bccoefs(1,1:4) = ((/1.0e0, 0.0e0, 0.0e0, 0.0e0/))
bccoefs(2,1:4) = ((/0.0e0, 1.0e0, 0.0e0, 0.0e0/))
bccoefs(3,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
case (2)
! This is the discount factor using the risk-free
! interest rate
df = exp(interest_rate * tx)
! Use flag passed to decide on boundary condition -
select case (fcn_data % idata(1))
case(1)
bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, bet*df/)
case(2)
bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0,&
(spread-strike_price)*df/)
end select
bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 0.0e0/)
bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
return
end select
! Note no time dependence in case (1) for iflag
iflag = 0
end subroutine fkbc_call
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
This example evaluates the price of a convertible bond. Here,
convertibility means that the bond may, at any time of the holders choosing, be
converted to a multiple of the specified asset. Thus a convertible bond
with pricereturns an amount
at time
unless the owner has
converted the bond to
units of the asset at some
time prior to
. This definition, the
differential equation and boundary conditions are given in Chapter 18 of
Wilmott et al. (1996). Using a constant interest rate and volatility factor, the
parameters and boundary conditions are:
1.
Bond face value , conversion factor
2.
Volatility
3.
Times until expiration
4.
Interest rate , dividend
5.
6.
7.
Boundary conditions
8.
Terminal data
9.
Constraint for bond holder
Note that the error tolerance is set to a pure absolute error
of value. The free boundary
constraint
is achieved by use of
a non-linear forcing term in the subroutine fkforce_cbond.
The terminal conditions are provided with the user subroutine fkinit_cbond.
Link to example source (feynman_kac_ex4.f90)
! Compute value of a Convertible Bond
use feynman_kac_int
use hqsval_int
use mp_types
use umach_int
implicit none
! The face value
real(kind(1e0)) :: ks = 1.0e0
! The sigma or volatility value
real(kind(1e0)) :: sigma = 0.25e0
! Time values for the options
integer, parameter :: nt = 2
real(kind(1e0)) :: time(nt)=(/0.5e0, 1.0e0/)
! Values of the underlying where evaluation are made
integer, parameter :: nv = 13
integer, parameter :: nlbc = 3, nrbc = 3, ndeg = 6
integer :: i
real(kind(1e0)) :: xs(nv) = (/((i-1)*0.25e0,i=1,nv)/)
! Value of the interest rate, continuous dividend and factor
real(kind(1e0)) :: r = 0.1e0, dividend = 0.02e0,&
factor =1.125e0
! Values of the min and max underlying values modeled
real(kind(1e0)) :: x_min = 0.0e0, x_max = 4.0e0
! Define parameters for the integration step.
integer, parameter :: nx = 61, nint = nx-1, n = 3*nx
real(kind(1e0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt),&
dx, f(nv,0:nt)
! Relative and absolute error tolerances
real(kind(1e0)) :: atol(1), rtol(1)
type(s_fcn_data) fcn_data
real(kind(1e0)), external :: fkcoef_cbond, fkinitcond_cbond
external fkbc_cbond, fkforce_cbond, fkinit_cbond
integer :: nout
call umach(2,nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fcn_data % rdata (7), fcn_data % idata (1))
! Define an equally-spaced grid of points for the underlying price
dx = (x_max - x_min)/real(nint)
xgrid(1) = x_min
xgrid(nx) = x_max
do i=2,nx-1
xgrid(i) = xgrid(i-1) + dx
end do
! Use a pure absolute error tolerance for the integration
! The default values require too much integation time.
atol(1) = 1.0e-3
rtol(1) = 0.0e0
! Pass the data for evaluation
fcn_data % rdata = (/ks,x_max,sigma,r,dividend,factor,&
atol(1)/)
fcn_data % idata = (/ndeg/)
! Compute value of convertible bond
call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_cbond,&
fkinitcond_cbond, fkbc_cbond, y, yprime,&
ATOL=atol, RTOL=rtol, FKINIT = fkinit_cbond,&
FKFORCE = fkforce_cbond, FCN_DATA = fcn_data)
! Evaluate and display solutions at vector of points XS(:), at each
! time value prior to expiration.
do i=0,nt
f(:,i) = hqsval (xs, xgrid, y(:,i))
end do
write(nout,'(T05,A)')&
"Convertible Bond Value, 0+, 6 and 12 Months Prior to Expiry"
write(nout,'(T08,"Number of equally spaced spline knots ",I4,'//&
'/T08,"Number of unknowns ",I4)') NX,N
write(nout,'(T08,"Strike = ",F5.2,", Sigma =", F5.2,/'//&
'T08,"Interest Rate =",F5.2,", Dividend =",F5.2,'//&
'", Factor = ",F5.3,//T08,"Underlying", T26,"Bond Value",'//&
'/(T10,4F8.4))') ks,sigma,r,dividend,factor,&
(xs(i), f(i,0:nt),i=1,nv)
end
! These routines define the coefficients, payoff, boundary
! conditions and forcing term.
function fkcoef_cbond(x, tx, iflag, fcn_data) result(value)
use mp_types
implicit none
integer, intent(inout) :: iflag
real(kind(1e0)), intent(in) :: x, tx
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: value
real(kind(1e0)) :: sigma, interest_rate, &
dividend, zero = 0.e0
sigma = fcn_data % rdata(3)
interest_rate = fcn_data % rdata(4)
dividend = fcn_data % rdata(5)
select case (iflag)
case (1)
! The coefficient derivative d(sigma)/dx
value = sigma
! The coefficient sigma(x)
case (2)
value = sigma * x
case (3)
! The coefficient mu(x)
value = (interest_rate - dividend) * x
case (4)
! The coefficient kappa(x)
value = interest_rate
end select
! Note that there is no time dependence
iflag = 0
return
end function fkcoef_cbond
function fkinitcond_cbond(x, fcn_data) result(value)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: x
type (s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: value
real(kind(1e0)) :: strike_price, factor
strike_price = fcn_data % rdata(1)
factor = fcn_data % rdata(6)
value = max(factor * x, strike_price)
return
end function fkinitcond_cbond
subroutine fkbc_cbond (tx, iflag, bccoefs, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
type(s_fcn_data), optional :: fcn_data
real(kind(1e0)) :: interest_rate, strike_price, dp,&
factor, x_max
select case (iflag)
case (1)
strike_price = fcn_data % rdata(1)
interest_rate = fcn_data % rdata(4)
dp = strike_price * exp(tx*interest_rate)
bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, dp/)
bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, 0.0e0/)
bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
return
case (2)
x_max = fcn_data % rdata(2)
factor = fcn_data % rdata(6)
bccoefs(1,1:4) = (/1.0e0, 0.0e0, 0.0e0, factor * x_max/)
bccoefs(2,1:4) = (/0.0e0, 1.0e0, 0.0e0, factor/)
bccoefs(3,1:4) = (/0.0e0, 0.0e0, 1.0e0, 0.0e0/)
end select
! Note no time dependence
iflag = 0
return
end subroutine fkbc_cbond
subroutine fkforce_cbond (interval, t, hx, y, xlocal, qw, u,&
phi, dphi, fcn_data)
use mp_types
implicit none
integer :: i, j, l
integer, parameter :: local = 6
integer, intent(in) :: interval
real(kind(1.e0)), intent(in) :: y(:), t, hx, qw(:),xlocal(:),&
u(:,:)
real(kind(1.e0)), intent(out) :: phi(:), dphi(:,:)
integer :: ndeg
real(kind(1.e0)) :: yl(local), bf(local)
real(kind(1.e0)) :: value, strike_price, interest_rate,&
zero = 0.0e0, one = 1.0e0, rt, mu, factor
type(s_fcn_data), optional :: fcn_data
yl = y(3*interval-2:3*interval+3)
phi = zero
dphi = zero
value = fcn_data % rdata(7)
strike_price = fcn_data % rdata(1)
interest_rate = fcn_data % rdata(4)
factor = fcn_data % rdata(6)
ndeg = fcn_data % idata(1)
mu = 2
! This is the local definition of the forcing term
! It "forces" the constraint f >= factor*x.
do j=1,local
do l = 1,ndeg
bf(1:3) = u(l,1:3)
bf(4:6) = u(l,7:9)
rt = dot_product(yl,bf)
rt = value/(rt + value - factor * xlocal(l))
phi(j) = phi(j) + qw(l) * bf(j) * rt**mu
end do
end do
phi = -phi * hx * factor * strike_price
! This is the local derivative matrix for the forcing term -
do j=1,local
do i = 1,local
do l=1,ndeg
bf(1:3) = u(L,1:3)
bf(4:6) = u(L,7:9)
rt = dot_product(yl,bf)
rt = one/(rt + value - factor * xlocal(l))
dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j)&
* (value * rt)**mu * rt
end do
end do
end do
dphi = -mu * dphi * hx * factor * strike_price
return
end subroutine fkforce_cbond
subroutine fkinit_cbond(xgrid,tgrid,t,yprime,y,atol,rtol,&
fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)
real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), yprime(:),&
t
type(s_fcn_data), optional :: fcn_data
integer :: i
if (t == 0.0e0) then
! Set initial data precisely.
do i=1,size(Xgrid)
if (xgrid(i)*fcn_data % rdata(6) <&
fcn_data % rdata(1)) then
y(3*i-2) = fcn_data % rdata(1)
y(3*i-1) = 0.0e0
y(3*i ) = 0.0e0
else
y(3*i-2) = xgrid(i) * fcn_data % rdata(6)
y(3*i-1) = fcn_data % rdata(6)
y(3*i ) = 0.0e0
end if
end do
end if
end subroutine fkinit_cbond
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
This example illustrates a method for evaluating a certain Bermudan Style or non-standard American option. These options are American Style options restricted to certain dates where the option may be exercised. Since this agreement gives the holder more opportunity than a European option, it is worth more. But since the holder can only exercise at certain times it is worth no more than the American style option value that can be exercised at any time. Our solution method uses the same model and data as in Example 2, but allows exercise at weekly intervals. Thus we integrate, for half a year, over each weekly interval using a European style Black-Scholes model, but with initial data at each new week taken from the corresponding values of the American style option.
Link to example source (feynman_kac_ex5.f90)
! Compute Bermudan-Style Option Premium for Vanilla Put
use feynman_kac_int
use hqsval_int
use mp_types
use umach_int
implicit none
integer :: nout
! The strike price
real(kind(1e0)) :: ks = 10.0e0
! The sigma value
real(kind(1e0)) :: sigma = 0.4e0
! Program working stores
real(kind(1e0)) :: week
! Time values for the options
integer, parameter :: nt = 26
integer, parameter :: ndeg = 6
real(kind(1e0)) :: time(nt), time_end = 0.5e0
! Values of the underlying where evaluation are made
integer, parameter :: nv = 9, nlbc = 2, nrbc = 3
integer :: i
real(kind(1e0)) :: xs(nv) = (/((i-1)*2.0e0,i=1,nv)/)
! Value of the interest rate and continuous dividend
real(kind(1e0)) :: r = 0.1e0, dividend = 0.0e0
! Values of the min and max underlying values modeled
real(kind(1e0)) :: x_min = 0.0e0, x_max = 30.0e0
! Define parameters for the integration step.
integer, parameter :: nx = 61, nint = nx-1, n = 3*nx
real(kind(1e0)) :: xgrid(nx), yb(n,0:nt), ybprime(n,0:nt),&
ya(n,0:nt), yaprime(n,0:nt),&
ytemp(n,0:1), ytempprime(n,0:1),&
dx, fb(nv,nt), fa(nv,nt)
real(kind(1e0)) :: atol
type(s_fcn_data) fcn_data_amer, fcn_data_berm
real(kind(1e0)), external :: fkcoef_put, fkinitcond_put
external fkbc_put, fkforce_put, fkinit_amer_put, fkinit_berm_put
call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fcn_data_amer % rdata (6), fcn_data_amer % idata (1))
! Define an equally-spaced grid of points for the underlying price
dx = (x_max-x_min)/real(nint)
xgrid(1) = x_min
xgrid(nx)= x_max
do i=2,nx-1
xgrid(i) = xgrid(i-1) + dx
end do
! Place a breakpoint at the strike price.
do i=1,nx
if (xgrid(i) > ks) then
xgrid(i-1) = ks
exit
end if
end do
! Compute time values where American option is computed
week = time_end/real(nt,kind(week))
time(1) = week
do i=2,nt-1
time(i) = time(i-1) + week
end do
time(nt) = time_end
atol = 1.0e-3
fcn_data_amer % rdata = (/ks,x_max,sigma,r,dividend,atol/)
fcn_data_amer % idata = (/ndeg/)
! Compute American Put Option Values at the weekly grid.
call feynman_kac (xgrid, time, nlbc, nrbc, fkcoef_put,&
fkinitcond_put, fkbc_put, ya, yaprime,&
FKINIT = fkinit_amer_put,&
FKFORCE = fkforce_put,&
FCN_DATA = fcn_data_amer)
! Integrate once again over the weekly grid, using the American
! Option values as initial data for a piece-wise European option
!integration.
! Allocate space to hold coefficient data and initial values.
allocate(fcn_data_berm % rdata(5+n))
fcn_data_berm % rdata(1:5) = fcn_data_amer % rdata(1:5)
! Copy initial data so the payoff value is the same for
! American and Bermudan option values.
yb(1:n,0) = ya(1:n,0)
ybprime(1:n,0) = ya(1:n,0)
do i=0,nt-1
! Move American Option values into place as initial conditions,
! but now integrating with European style over each period of
! the weekly grid.
fcn_data_berm % rdata(6:) = ya(1:n,i)
if (i .eq. 0) then
call feynman_kac (xgrid, (/time(1)/), nlbc, nrbc,&
fkcoef_put, fkinitcond_put, fkbc_put,&
ytemp(:,0:1), ytempprime(:,0:1),&
FKINIT = fkinit_berm_put,&
FCN_DATA = fcn_data_berm)
else
call feynman_kac (xgrid, (/time(i+1)-time(i)/),&
nlbc, nrbc, fkcoef_put,&
fkinitcond_put, fkbc_put,&
ytemp(:,0:1), ytempprime(:,0:1),&
FKINIT = fkinit_berm_put,&
FCN_DATA = fcn_data_berm)
end if
! Record values of the Bermudan option at the end of each integration.
yb(1:n,i+1) = ytemp(1:n,1)
ybprime(1:n,i+1) = ytempprime(1:n,1)
end do
! Evaluate solutions at vector of points XS(:), at each time value
! prior to expiration.
do i=1,nt
fa(:,i) = hqsval (xs, xgrid, ya(:,i))
fb(:,i) = hqsval (xs, xgrid, yb(:,i))
end do
write(nout,'(T05,A)')&
"American Option Premium for Vanilla Put, 6 Months "//&
"Prior to Expiry"
write(nout,'(T05,A)')&
"Exercise Opportunities At Weekly Intervals"
write(nout,'(T08,"Number of equally spaced spline knots ",'//&
'I4,/T08,"Number of unknowns ",I4)') nx, n
write(nout,'(T08,"Strike = ",F5.2,", Sigma =", F5.2,'//&
'", Interest Rate =",F5.2,//T08,"Underlying",'//&
'T20,"Bermudan Style", T42,"American",'//&
'/(T10,F8.4, T26, F8.4, T42, F8.4))')&
KS,SIGMA,R,&
(xs(i), fb(i,nt:nt), fa(i,nt:nt),i=1,nv)
end
! These subprograms set the coefficients, payoff, boundary
! conditions and forcing term for American and European Options.
function fkcoef_put(x, tx, iflag, fcn_data_amer)&
result(value)
use mp_types
implicit none
integer, intent(inout) :: iflag
real(kind(1e0)), intent(in) :: x, tx
type(s_fcn_data), optional :: fcn_data_amer
real(kind(1e0)) :: value
real(kind(1e0)) :: sigma, interest_rate, dividend, zero=0.0e0
sigma = fcn_data_amer % rdata(3)
interest_rate = fcn_data_amer % rdata(4)
dividend = fcn_data_amer % rdata(5)
select case (iflag)
case (1)
! The coefficient derivative d(sigma)/dx
value = sigma
! The coefficient sigma(x)
case (2)
value = sigma * x
case (3)
! The coefficient mu(x)
value = (interest_rate - dividend) * x
case (4)
! The coefficient kappa(x)
value = interest_rate
end select
! Note that there is no time dependence
iflag = 0
return
end function fkcoef_put
function fkinitcond_put(x, fcn_data_amer) result(value)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: x
type (s_fcn_data), optional :: fcn_data_amer
real(kind(1e0)) :: value
real(kind(1e0)) :: strike_price, zero = 0.0e0
strike_price = fcn_data_amer % rdata(1)
! The payoff function
value = max(strike_price - x, zero)
return
end function fkinitcond_put
subroutine fkbc_put (tx, iflag, bccoefs, fcn_data)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
type (s_fcn_data), optional :: fcn_data
select case (iflag)
case (1)
bccoefs(1,1:4) = ((/0.0e0, 1.0e0, 0.0e0, -1.0e0/))
bccoefs(2,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
case (2)
bccoefs(1,1:4) = ((/1.0e0, 0.0e0, 0.0e0, 0.0e0/))
bccoefs(2,1:4) = ((/0.0e0, 1.0e0, 0.0e0, 0.0e0/))
bccoefs(3,1:4) = ((/0.0e0, 0.0e0, 1.0e0, 0.0e0/))
end select
! Note no time dependence
iflag = 0
end subroutine fkbc_put
subroutine fkforce_put (interval, t, hx, y, xlocal, qw, u,&
phi, dphi, fcn_data_amer)
use mp_types
implicit none
integer, parameter :: local = 6
integer :: i, j, l, ndeg
integer, intent(in) :: interval
real(kind(1.e0)), intent(in) :: y(:), t, hx, qw(:),&
xlocal(:), u(:,:)
real(kind(1.e0)), intent(out) :: phi(:), dphi(:,:)
type (s_fcn_data), optional :: fcn_data_amer
real(kind(1.e0)) :: yl(local), bf(local)
real(kind(1.e0)) :: value, strike_price, interest_rate,&
zero = 0.e0, one = 1.e0, rt, mu
yl = y(3*interval-2:3*interval+3)
phi = zero
value = fcn_data_amer % rdata(6)
strike_price = fcn_data_amer % rdata(1)
interest_rate = fcn_data_amer % rdata(4)
ndeg = fcn_data_amer % idata(1)
mu = 2
! This is the local definition of the forcing term
do j=1,local
do l=1,ndeg
bf(1:3) = U(L,1:3)
bf(4:6) = U(L,7:9)
rt = dot_product(YL,BF)
rt = value/(rt + value-(strike_price-xlocal(l)))
phi(j) = phi(j) + qw(l) * bf(j) * rt**mu
end do
end do
phi = -phi * hx * interest_rate * strike_price
! This is the local derivative matrix for the forcing term
dphi = zero
do j=1,local
do i = 1,local
do l=1,ndeg
bf(1:3) = u(L,1:3)
bf(4:6) = u(L,7:9)
rt = dot_product(yl,bf)
rt = one/(rt + value - (strike_price - xlocal(l)))
dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j) *&
rt**(mu+1)
end do
end do
end do
dphi = mu * dphi * hx * value**mu * interest_rate *&
strike_price
end subroutine fkforce_put
subroutine fkinit_amer_put(xgrid,tgrid,t,yprime,y,atol,rtol,&
fcn_data_amer)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&
yprime(:)
real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)
type(s_fcn_data), optional :: fcn_data_amer
integer :: i
if (t == 0.0e0) then
! Set initial data precisely. The strike price is a breakpoint.
! Average the derivative limit values from either side.
do i=1,size(xgrid)
if (xgrid(i) < fcn_data_amer % rdata(1)) then
y(3*i-2) = fcn_data_amer % rdata(1) - xgrid(i)
y(3*i-1) = -1.0e0
y(3*i) = 0.0e0
else if (xgrid(i) == fcn_data_amer % rdata(1)) then
y(3*i-2) = 0.0e0
y(3*i-1) = -0.5e0
y(3*i) = 0.0e0
else
y(3*i-2) = 0.0e0
y(3*i-1) = 0.0e0
y(3*i) = 0.0e0
end if
end do
end if
end subroutine fkinit_amer_put
subroutine fkinit_berm_put(xgrid,tgrid,t,yprime,y,atol,rtol,&
fcn_data_berm)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: xgrid(:), tgrid(:), t,&
yprime(:)
real(kind(1e0)), intent(inout) :: y(:), atol(:), rtol(:)
type(s_fcn_data), optional :: fcn_data_berm
integer :: i
if (t == 0.0e0) then
! Set initial data for each week at the previously computed
! American Option values. These coefficients are passed
! in the derived type fcn_data_berm.
do i=1,size(xgrid)
y(3*i-2) = fcn_data_berm % rdata(3+3*i)
y(3*i-1) = fcn_data_berm % rdata(4+3*i)
y(3*i ) = fcn_data_berm % rdata(5+3*i)
end do
end if
end subroutine fkinit_berm_put
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
Our previous examples are from the field of financial engineering. A final example is a physical model. The Oxygen Diffusion Problem is summarized in Crank [4], p. 10-20, 261-262. We present the numerical treatment of the transformed one-dimensional system
A slight difference from the Crank development is that we
have reflected the time variable to match our form of
the Feynman-Kac equation. We have a free boundary problem because the
interface
is implicit. This
interface is implicitly defined by solving the variational relation
. The first factor is zero for
and the second factor is zero for
. We list the Feynman-Kac equation coefficients, forcing
term and boundary conditions, followed by comments.
The forcing term has the
property of being almost the value 1 when the solution is larger than the
factor
. As the solution
, the forcing term
is almost the value zero. These properties combine to
approximately achieve the variational relation that defines the free
boundary. Note that the arc of the free boundary is not explicitly
available with this numerical method. We have used
the requested absolute error tolerance for the numerical
integration.
The boundary condition is discontinuous as
, since the initial data yields
. For the numerical integration we have chosen a boundary
value function that starts with the value
at
and rapidly blends to
the value zero as the integration proceeds in the negative direction. It
is necessary to give the integrator the opportunity to severely limit the step
size near
to achieve continuity
during the blending.
In the example code, values of are checked against published values for certain values of
. Also checked are
values of
at published values of the
free boundary, for the same values of
.
Link to example source (feynman_kac_ex6.f90)
! Integrate Oxygen Diffusion Model found in the book
! Crank, John. Free and Moving Boundary Problems,
! Oxford Univ. Press, (1984), p. 19-20 and p. 261-262.
use feynman_kac_int
use hqsval_int
use mp_types
use norm_int
use umach_int
implicit none
integer :: nout
real(kind(1e0)), allocatable :: xgrid(:), tgrid(:), y(:,:),&
yprime(:,:), f(:,:), s(:)
real(kind(1e0)) :: dx, x_min, x_max, zero=0.0e0, one=1.0e0
real(kind(1e0)) :: atol(1), rtol(1)
type(s_fcn_data) :: fk_ox2
integer :: i, nint, n, nunit, ntimes = 8
integer, parameter :: ndeg = 6, nlbc = 1, nrbc = 2
real(kind(1e0)), external :: fkcoef_ox2, fkinitcond_ox2
external fkbc_ox2, fkforce_ox2
call umach(2,nout)
! Define number of equally spaced intervals for elements
nint = 100
! Allocate the space needed for the integration process
n = 3*(nint+1)
allocate(xgrid(nint+1), y(n,0:ntimes), yprime(n,0:ntimes),&
tgrid(ntimes), f(1,ntimes), s(ntimes))
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fk_ox2 % rdata (1), fk_ox2 % idata (1))
atol(1) = 0.5e-2
rtol(1) = 0.5e-2
fk_ox2 % rdata(1) = atol(1)
fk_ox2 % idata(1) = ndeg
! Define interval endpoints
x_min = zero
x_max = one
! Define interval widths
dx = (x_max-x_min)/real(nint)
xgrid(1) = x_min
xgrid(nint+1) = x_max
! Define grid points of interval
do i=2,nint
xgrid(i) = xgrid(i-1) + dx
end do
! Define time integration output points
! These correspond to published values in Crank's book, p. 261-262
tgrid = (/0.04e0,0.06e0,0.10e0,0.12e0,0.14e0,&
0.16e0,0.18e0,0.185e0/)
call feynman_kac (xgrid, tgrid, nlbc, nrbc, fkcoef_ox2,&
fkinitcond_ox2, fkbc_ox2, y, yprime,&
ATOL = atol, RTOL = rtol,&
FKFORCE = fkforce_ox2, FCN_DATA = fk_ox2)
! Summarize output at the left end
do i=1,ntimes
f(:,i)= hqsval ((/zero/), xgrid, y(:,i))
end do
! Check differences of evaluation and published left end values
f(1,:) = f(1,:) - (/2.743e-01, 2.236e-01, 1.432e-01,&
1.091e-01, 7.786e-02, 4.883e-02, 2.179e-02, 1.534e-02/)
write(nout,*) "Oxygen Depletion Model, from Crank's "//&
"Book, p. 261-262,"
write(nout,*) "'Free and Moving Boundary Value Problems'"
if (norm(f(1,:)) < ntimes * atol(1)) then
write(nout,*) "FEYNMAN_KAC Example 6 - Fixed Sealed "//&
"Surface Values are correct"
else
write(nout,*) "FEYNMAN_KAC Example 6 - does not agree with"//&
" published left end values"
end if
! Define known position of free boundary at the time points
s = (/0.9992e0,0.9918e0,0.9350e0,0.8792e0,&
0.7989e0,0.6834e0,0.5011e0,0.4334e0/)
! Evaluate and verify solution is small near free boundary -
do i=1,ntimes
f(:,i) = hqsval ((/s(i)/), xgrid, y(:,i))
end do
if (norm(f(1,:)) < ntimes * atol(1)) then
write(nout,*) "FEYNMAN_KAC Example 6 - Free Boundary "//&
"Position Values are correct"
else
write(nout,*) "FEYNMAN_KAC Example 6 - does not agree "//&
"with published free boundary values"
end if
end
function fkcoef_ox2 (x, tx, iflag, fk_ox2) result(value)
use mp_types
implicit none
! Coefficient valuation routine for the Oxygen Diffusion
! Model found in Crank's book, p. 20
! Input/Ouput variables
integer, intent(inout) :: iflag
real(kind(1e0)), intent(in) :: x, tx
type(s_fcn_data), optional :: fk_ox2
real(kind(1e0)) :: value
! Local variables
real(kind(1e0)) :: zero = 0.0e0, two = 2.0e0
select case (iflag)
case (1) ! Factor DSigma/Dx(x,t)
value = zero
case (2) ! Factor Sigma(x,t)
value = sqrt(two)
case (3) ! Factor Mu (x,t)
value = zero
case (4) ! Factor Kappa (x,t)
value = zero
end select
! Signal no dependence on tx=t=time for any coefficient.
iflag = 0
return
end function fkcoef_ox2
function fkinitcond_ox2(x, fk_ox2) result(value)
use mp_types
implicit none
real(kind(1e0)), intent(in) :: x
type (s_fcn_data), optional :: fk_ox2
real(kind(1e0)) :: value
real(kind(1e0)) :: half = 0.5e0, one = 1.0e0
value = half * (one - x)**2
return
end function fkinitcond_ox2
subroutine fkbc_ox2 (tx, iflag, bccoefs, fk_ox2)
use mp_types
implicit none
! Evaluation routine for Oxygen Diffusion Model
! boundary conditions.
! Input/Ouput variables
real(kind(1e0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1e0)), dimension(:,:), intent(out) :: bccoefs
type (s_fcn_data), optional :: fk_ox2
! Local variables
real(kind(1e0)) :: zero = 0.0e0, one = 1.0e0, atol
atol = fk_ox2 % rdata(1)
select case (iflag)
case (1) ! Left Boundary Condition, at X_min=0
! There is a rapid blending of the boundary condition to achieve
! a zero derivative value at the left end.
! The initial data has the derivative with value one.
! This boundary condition essentially abruptly changes that
! derivative value to zero.
! Returning iflag=1 signals time dependence. This is important
! for this problem.
bccoefs(1,1:4) = (/0.0e0, one, 0.0e0, exp(tx/atol**2)/)
return
case (2) ! Right Boundary Condition, at X_max=1
bccoefs(1,1:4) = (/one, 0.0e0, 0.0e0, 0.0e0/)
bccoefs(2,1:4) = (/0.0e0, one, 0.0e0, 0.0e0/)
end select
iflag = 0 ! Signal no dependence on tx=time.
end subroutine fkbc_ox2
subroutine fkforce_ox2 (interval, t, hx, y, xlocal, qw, u,&
phi, dphi, fk_ox2)
! Evaluation routine for Oxygen Diffusion model forcing function.
use mp_types
implicit none
integer, parameter :: local = 6
integer :: i, j, l, mu, ndeg
integer, intent(in) :: interval
real(kind(1e0)), intent(in) :: y(:), t, hx, qw(:),&
xlocal(:), u(:,:)
real(kind(1e0)), intent(out) :: phi(:), dphi(:,:)
type (s_fcn_data), optional :: fk_ox2
real(kind(1e0)) :: yl(local), bf(local)
real(kind(1e0)) :: value, zero = 0.0e0, one = 1.0e0, rt
yl = y(3*interval-2:3*interval+3)
phi = zero
value = fk_ox2 % rdata(1)
ndeg = fk_ox2 % idata(1)
mu = 2
do j=1,local
do l=1,ndeg
bf(1:3) = u(l,1:3)
bf(4:6) = u(l,7:9)
rt = dot_product(yl,bf)
rt = one - (value/(rt + value))**mu
phi(j) = phi(j) + qw(l) * bf(j) * RT
end do
end do
phi = phi * hx
! This is the local derivative matrix for the forcing term -
dphi = zero
do j=1,local
do i = 1,local
do l=1,ndeg
bf(1:3) = u(l,1:3)
bf(4:6) = u(l,7:9)
rt = dot_product(yl,bf)
rt = one/(rt + value)
dphi(i,j) = dphi(i,j) + qw(l) * bf(i) * bf(j) *&
rt**(mu+1)
end do
end do
end do
dphi = mu * dphi * hx * value**mu
return
end subroutine fkforce_ox2
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
In this example, routine FEYNMAN_KAC
is used to solve for the Greeks, i.e. various derivatives of Feynman-Kac (FK)
solutions applicable to the pricing of options and related financial
derivatives. In order to illustrate and verify these calculations, the Greeks
are calculated by two methods. The first method involves the FK solution to the
diffusion model for call options given in Example 1 for the Black-Scholes (BS)
case, i.e. . The second method
calculates the Greeks using the closed-form BS evaluations which can be found at
http://en.wikipedia.org/wiki/The_Greeks.
This example calculates FK and BS solutions to the BS problem and the following Greeks:
Delta =
is the
first derivative of the Value, , of a portfolio of derivative security derived from underlying
instrument with respect to the underlying instruments price S;
Gamma = ;
Theta = is the negative first
derivative of V with respect to time t;
Charm = ;
Color = ;
Rho = is the first
derivative of V with respect to the risk-free rate r;
Vega = measures sensivity to
volatility parameter
of the underlying
S;
Volga = ;
Vanna = ;
Speed = .
Intrinsic Greeks include derivatives involving only
S and t, the intrinsic FK arguments. In the above list,
Value, Delta, Gamma, Theta, Charm,
Color and Speed are all intrinsic Greeks. As is discussed in
Hanson, R. (2008) Integrating
Feynman-Kac Equations Using Hermite Quintic Finite Elements, the
expansion of the FK solution function in terms of quintic polynomial functions defined on
S-grid subintervals and subject to continuity constraints in derivatives
0, 1 and 2 across the boundaries of these subintervals allows Value,
Delta, Gamma, Theta, Charm and Color to be
calculated directly by routines FEYNMAN_KAC
and HQSVAL.
Non-intrinsic Greeks are derivatives of V
involving FK parameters other than the intrinsic arguments S and
t, such as r and . Non-intrinsic Greeks in
the above list include Rho, Vega, Volga and Vanna.
In order to calculate non-intrinsic Greek (parameter) derivatives or intrinsic
Greek S-derivatives beyond the second (such as Speed) or
t-derivatives beyond the first, the entire FK solution must be calculated
3 times (for a parabolic fit) or five times (for a quartic fit), at the point
where the derivative is to be evaluated and at nearby points located
symmetrically on either side.
Using a Taylor series expansion of truncated to
terms (to allow an m-degree polynomial fit of
m+1 data points),
,
we are able to derive the following parabolic (3 point)
estimation of first and second derivatives and
in terms of the three
values
,
and
, where
and 0
< ε frac
<<1:
,
.
Similarly, the quartic (5 point) estimation of and
in terms of
,
,
,
and
is:
.
For our example, the quartic estimate does not appear to be significantly better than the parabolic estimate, so we have produced only parabolic estimates by setting variable iquart to 0. The user may try the example with the quartic estimate simply by setting iquart to 1.
As is pointed out in Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements, the quintic polynomial expansion function used by FEYNMAN_KAC only allows for continuous derivatives through the second derivative. While up to fifth derivatives can be calculated from the quintic expansion (indeed function HQSVAL will allow the third derivative to be calculated by setting optional argument IDERIV to 3, as is done in this example), the accuracy is compromised by the inherent lack of continuity across grid points (i.e. subinterval boundaries).
The accurate second derivatives in S returned by
function HQSVAL
can be leveraged into a third derivative estimate by calculating three FK second
derivative solutions, the first solution for grid and evaluation point set
and the second and
third solutions for solution grid and evaluation point sets
and
, where the solution grid
and evaluation point sets are shifted up and down by
. In this example,
is set to
, where
is the average value of S over the range of grid
values and
. The third derivative
solution can then be obtained using the parabolic estimate
This procedure is implemented in the current example to calculate the Greek Speed. (For comparison purposes, Speed is also calculated directly by setting the optional argument IDERIV to 3 in the call to function HQSVAL. The output from this direct calculation is called Speed2.)
To reach better accuracy results, all computations are done in double precision.
The average and maximum relative errors (defined as the
absolute value of the difference between the BS and FK values divided by the BS
value) for each of the Greeks is given at the end of the output. (These relative
error statistics are given for nine combinations of Strike Price and volatility,
but only one of the nine combinations is actually printed in the output.) Both
intrinsic and non-intrinsic Greeks have good accuracy (average relative error is
in the range 0.01 0.00001) except for Volga, which has an average
relative error of about 0.05. This is probably a result of the fact that
Volga involves differences of differences, which will more quickly erode
accuracy than calculations using only one difference to approximate a
derivative. Possible ways to improve upon the 2 to 4 significant digits of
accuracy achieved in this example include increasing FK integration accuracy by
reducing the initial stepsize (via optional argument RINITSTEPSIZE),
by choosing more closely spaced S and t grid points (by adjusting
FEYNMAN_KACs
input arrays XGRID
and TGRID)
and by adjusting so that the central
differences used to calculate the derivatives are not too small to compromise
accuracy.
Link to example source (feynman_kac_ex7.f90)
! Greeks computation
use feynman_kac_int
use hqsval_int
use mp_types
use anordf_int
use const_int
use umach_int
implicit none
real(kind(1d0)), external :: fkcoef, fkinitcond
external fkbc
! The set of strike prices
real(kind(1d0)) :: ks(3) = (/15.0d0,20.0d0,25.0d0/)
! The set of sigma values
real(kind(1d0)) :: sigma(3) = (/0.2d0, 0.3d0, 0.4d0/)
! The set of model diffusion powers: alpha = 2.0 <==> Black Scholes
real(kind(1d0)) :: alpha(3) = (/2.0d0, 1.0d0, 0.0d0/)
! Time values for the options
integer, parameter :: nt = 3
real(kind(1d0)) :: time(nt)=(/1.d0/12., 4.d0/12., 7.d0/12./)
! Values of the min and max underlying values modeled
real(kind(1d0)) :: x_min = 0.0d0, x_max = 60.0d0
! Value of the interest rate and continuous dividend
real(kind(1d0)) :: r = 0.05d0, dividend = 0.0d0
! Values of the underlying where evaluations are made
integer, parameter :: nv = 3
real(kind(1d0)) :: eval_points(nv) = &
(/19.0d0, 20.0d0, 21.0d0/)
! Define parameters for the integration step.
integer, parameter :: nx = 121, nint = nx-1, n = 3*nx
real(kind(1d0)) :: xgrid(nx), y(n,0:nt), yprime(n,0:nt)
type(d_fcn_data) fcn_data
! Number of left/right boundary conditions
integer, parameter :: nlbc = 3, nrbc = 3
! Further parameters for the integration step
real(kind(1d0)) :: dx, dx2, pi, sqrt2pi
! used to calc derivatives
real(kind(1d0)) :: epsfrac = .001d0
character(len=6) :: greek_name(12) = (/&
" Value", " Delta", " Gamma", " Theta",&
" Charm", " Color", " Vega", " Volga",&
" Vanna", " Rho", " Speed", "Speed2"/)
! Time values for the options
real(kind(1d0)) :: rex(12), reavg(12)
integer :: irect(12)
integer :: i, i2, i3, j, ig, iquart, nout
real(kind(1d0)) ::&
spline_values(nv,nt,12), spline_values1(nv,nt),&
spline_valuesp(nv,nt), spline_valuesm(nv,nt),&
spline_valuespp(nv,nt), spline_valuesmm(nv,nt),&
xgridp(nx), xgridm(nx),xgridpp(nx), xgridmm(nx),&
eval_pointsp(nv), eval_pointspp(nv),&
eval_pointsm(nv), eval_pointsmm(nv),&
BS_values(nv,nt,12), sVo_array(nt)
call umach(2, nout)
! Allocate space inside the derived type for holding
! data values. These are for the evaluation routines.
allocate(fcn_data % rdata (6))
pi = const('pi')
sqrt2pi = sqrt(2.0 * pi)
dx2 = epsfrac * 0.5d0 * (x_min + x_max)
! Compute Constant Elasticity of Variance Model for Vanilla Call
! Define equally-spaced grid of points for the underlying price
dx = (x_max - x_min)/real(nint)
xgrid(1) = x_min
xgrid(nx) = x_max
do i = 2,nx-1
xgrid(i) = xgrid(i-1) + dx
end do
write(nout,'(T05,A)') "Constant Elasticity of Variance"//&
" Model for Vanilla Call Option"
write(nout,'(T10,"Interest Rate: ", F7.3, T38,'//&
'"Continuous Dividend: ", F7.3 )') r, dividend
write(nout,'(T10,"Minimum and Maximum Prices of '//&
' Underlying: ", 2F7.2)') x_min, x_max
write(nout,'(T10,"Number of equally spaced spline knots:",'//&
'I4)') nx - 1
write(nout,'(T10,"Number of unknowns: ",I4)') n
write(nout,*)
write(nout,'(/T10,"Time in Years Prior to Expiration: ",'//&
'2X,3F7.4)') time
write(nout,'(T10,"Option valued at Underlying Prices: ",'//&
'3F7.2)') eval_points
write(nout,*)
!
! iquart = 0 : derivatives estimated with 3-point fitted parabola
! iquart = 1 : derivatives estimated with 5-point fitted quartic
! polynomial
!
iquart = 0
if (iquart == 0) then
write(nout,'(T10,"3 point (parabolic) estimate of '//&
'parameter derivatives")')
else
write(nout,'(T10,"5 point (quartic) estimate of parameter'//&
' derivatives")')
end if
write(nout, '(T10,"epsfrac = ", F11.8)') epsfrac
irect = 0
reavg = 0.0d0
rex = 0.0d0
! alpha: Black-Scholes
do i2 = 1, 3
! Loop over volatility
do i3 = 1, 3
! Loop over strike price
call calc_Greeks(i2, i3, iquart)
end do
end do
write(nout,*)
do ig = 1, 12
reavg(ig) = reavg(ig)/irect(ig)
write (nout, '(" Greek: ", A6, "; avg rel err: ",'//&
'F15.12, "; max rel err: ", F15.12)')&
greek_name(ig), reavg(ig), rex(ig)
end do
contains
subroutine calc_Greeks(volatility, strike_price, iquart)
implicit none
integer, intent(in) :: volatility, strike_price, iquart
! Local variables
integer :: i1 = 1, j, iSderiv, gNi, l, k
integer :: nt = 3
real(kind(1d0)) :: x_maxp, x_maxm, x_maxpp, x_maxmm
real(kind(1d0)) :: eps, tau, sigsqtau, sqrt_sigsqtau, sigsq
real(kind(1d0)) :: d1, d2, n01pdf_d1, nu, relerr, relerrmax
real(kind(1d0)) :: sVo, BSVo, S
if ((volatility == 1) .and. (strike_price == 1)) then
write(nout,*)
write(nout,'(/T10,"Strike = ",F5.2,", Sigma =", F5.2,'//&
'", Alpha =",F5.2,":")') ks(strike_price),&
sigma(volatility), alpha(i1)
write(nout,*)
write(nout,'(T10,"years to expiration: ", 3F7.4)')&
(time(j), j=1,3)
write(nout,*)
end if
fcn_data % rdata = (/ks(strike_price), x_max,&
sigma(volatility), alpha(i1), r, dividend/)
call feynman_kac(xgrid, time, nlbc, nrbc, fkcoef,&
fkinitcond, fkbc, y, yprime,&
FCN_DATA = fcn_data)
! Compute Value, Delta, Gamma, Theta, Charm and Color
do l = 0,2
do i=1,NT
spline_values(:,i,l+1) = hqsval(eval_points, xgrid,&
y(:,i), IDERIV=l)
spline_values(:,i,l+4) = hqsval(eval_points, xgrid,&
yprime(:,i), IDERIV=l)
end do
end do
! Signs for Charm and Color must be inverted because FEYNMAN_KAC
! computes -d/dt instead of d/dt
spline_values(:,:,5:6) = -spline_values(:,:,5:6)
! Speed2
do i=1,nt
spline_values(:,i,12) = hqsval(eval_points, xgrid, Y(:,i),&
IDERIV=3)
end do
! Compute Vega, Volga, Vanna, Rho, Speed
! l = 7 8 9 10 11
do l = 7,11
xgridp = xgrid
xgridm = xgrid
eval_pointsp = eval_points
eval_pointsm = eval_points
x_maxp = x_max
x_maxm = x_max
fcn_data % rdata(3) = sigma(volatility)
fcn_data % rdata(5) = r
iSderiv = 0
if (l == 9) iSderiv = 1 ! Vanna
if (l == 11) iSderiv = 2 ! Speed
if (l == 10) then
fcn_data % rdata(5) = r * (1.0 + epsfrac) ! Rho
else if (l < 10) then
fcn_data % rdata(3) = sigma(volatility) * (1.0 + epsfrac)
end if
if (l == 11) then
xgridp = xgrid + dx2
xgridm = xgrid - dx2
eval_pointsp = eval_points + dx2
eval_pointsm = eval_points - dx2
x_maxp = x_max + dx2
x_maxm = x_max - dx2
end if
fcn_data % rdata(2) = x_maxp
call feynman_kac(xgridp, time, nlbc, nrbc, fkcoef,&
fkinitcond, fkbc, y, yprime,&
FCN_DATA = fcn_data)
do i=1,nt
spline_valuesp(:,i) = hqsval(eval_pointsp, xgridp,&
y(:,i), IDERIV=iSderiv)
end do
if (l == 10) then
fcn_data % rdata(5) = r * (1.0 - epsfrac) ! Rho
else if (l < 10) then
fcn_data % rdata(3) = sigma(volatility) *&
(1.0 - epsfrac)
end if
fcn_data % rdata(2) = x_maxm
! calculate spline values for sigmaM = sigmai2-1*(1. - epsfrac) or
! rM = r*(1. - epsfrac):
call feynman_kac(xgridm, time, nlbc, nrbc, fkcoef,&
fkinitcond, fkbc, y, yprime,&
FCN_DATA = fcn_data)
do i=1,NT
spline_valuesm(:,i) = hqsval(eval_pointsm, xgridm,&
y(:,i), IDERIV=iSderiv)
end do
if (iquart == 1) then
xgridpp = xgrid
xgridmm = xgrid
eval_pointspp = eval_points
eval_pointsmm = eval_points
x_maxpp = x_max
x_maxmm = x_max
if (l == 11) then ! Speed
xgridpp = xgrid + 2.0 * dx2
xgridmm = xgrid - 2.0 * dx2
eval_pointspp = eval_points + 2.0 * dx2
eval_pointsmm = eval_points - 2.0 * dx2
x_maxpp = x_max + 2.0 * dx2
x_maxmm = x_max - 2.0 * dx2
end if
fcn_data % rdata(2) = x_maxpp
if (l == 10) then
! calculate spline values for rPP = r*(1. + 2.*epsfrac):
fcn_data % rdata(5) = r * (1.0 + 2.0 * epsfrac)
else if (l < 10) then
! calculate spline values for sigmaPP = sigma(i2)*(1. + 2.*epsfrac):
fcn_data % rdata(3) = sigma(volatility) *&
(1.0 + 2.0 * epsfrac)
end if
call feynman_kac (xgridpp, time, nlbc, nrbc, fkcoef,&
fkinitcond, fkbc, y, yprime,&
FCN_DATA = fcn_data)
do i=1,nt
spline_valuespp(:,i) = hqsval(eval_pointspp, xgridpp,&
Y(:,i), IDERIV=iSderiv)
end do
fcn_data % rdata(2) = x_maxmm
! calculate spline values for sigmaMM = sigmai2-1*(1. - 2.*epsfrac)
! or rMM = r*(1. - 2.*epsfrac)
if (l == 10) then
fcn_data % rdata(5) = r * (1.0 - 2.0 * epsfrac)
else if (l < 10) then
fcn_data % rdata(3) = sigma(volatility) *&
(1.0 - 2.0 * epsfrac)
end if
call feynman_kac (xgridmm, time, nlbc, nrbc, fkcoef,&
fkinitcond, fkbc, y, yprime,&
FCN_DATA = fcn_data)
do i=1,nt
spline_valuesmm(:,i) = hqsval(eval_pointsmm, xgridmm,&
y(:,i), IDERIV=iSderiv)
end do
end if ! if (iquart == 1)
if (l /= 8) then
eps = sigma(volatility) * epsfrac
if (l == 10) eps = r * epsfrac ! Rho
if (l == 11) eps = dx2 ! Speed
spline_values(:,:,l) = &
(spline_valuesp - spline_valuesm) / (2.0 * eps)
if (iquart /= 0) then
spline_values1 =&
(spline_valuespp - spline_valuesmm) /(4.0 * eps)
spline_values(:,:,l) =&
(4.0 * spline_values(:,:,l) - spline_values1) / 3.0
end if
end if
if (l == 8) then ! Volga
eps = sigma(volatility) * epsfrac
spline_values(:,:,l) =&
(spline_valuesp + spline_valuesm - 2.0 * &
spline_values(:,:,1)) / (eps * eps)
if (iquart /= 0) then
spline_values1 =&
(spline_valuespp + spline_valuesmm - 2.0 * &
spline_values(:,:,1)) / (4.0 * eps * eps)
spline_values(:,:,l) = &
(4.0 * spline_values(:,:,l) - spline_values1) / 3.0
end if
end if
end do
! Evaluate BS solution at vector eval_points,
! at each time value prior to expiration.
do i = 1,nt
!
! Black-Scholes (BS) European call option
! value = ValBSEC(S,t) = exp(-q*tau)*S*N01CDF(d1) -
! exp(-r*tau)*K*N01CDF(d2),
! where:
! tau = time to maturity = T - t;
! q = annual dividend yield;
! r = risk free rate;
! K = strike price;
! S = stock price;
! N01CDF(x) = N(0,1)_CDF(x);
! d1 = ( log( S/K ) +
! ( r - q + 0.5*sigma**2 )*tau ) /
! ( sigma*sqrt(tau) );
! d2 = d1 - sigma*sqrt(tau)
!
! BS option values for tau = time(i):
tau = time(i)
sigsqtau = (sigma(volatility)**2) * tau
sqrt_sigsqtau = sqrt(sigsqtau)
sigsq = sigma(volatility) * sigma(volatility)
do j = 1, nv
! Values of the underlying price where evaluations are made:
S = eval_points(j)
d1 = (log(S / ks(strike_price)) + (r - dividend)&
* tau + 0.5 * sigsqtau) / sqrt_sigsqtau
n01pdf_d1 = exp((-0.5) * d1 * d1) / sqrt2pi
nu = exp((-dividend) * tau) * S * n01pdf_d1 * sqrt(tau)
d2 = d1 - sqrt_sigsqtau
BS_values(j,i,1) = exp((-dividend) * tau) * S *&
anordf(d1) - exp((-r) * tau) *&
ks(strike_price) * anordf(d2)
! greek = Rho
BS_values(j,i,10) = exp((-r) * tau) * ks(strike_price) *&
tau * anordf(d2)
! greek = Vega
BS_values(j,i,7) = nu
! greek = Volga
BS_values(j,i,8) = nu * d1 * d2 / sigma(volatility)
! greek = delta
BS_values(j,i,2) = exp((-dividend) * tau) * anordf(d1)
! greek = Vanna
BS_values(j,i,9) = (nu / S) * (1.0 - d1 / sqrt_sigsqtau)
! greek = dgamma
BS_values(j,i,3) = exp((-dividend) * tau) *&
n01pdf_d1 /(S * sqrt_sigsqtau)
! greek = speed
BS_values(j,i,11) = (-exp((-dividend) * tau)) *&
n01pdf_d1 * (1.0 + d1 / sqrt_sigsqtau)&
/ (S * S * sqrt_sigsqtau)
! greek = speed
BS_values(j,i,12) = (-exp((-dividend) * tau)) * &
n01pdf_d1 * (1.0 + d1 / sqrt_sigsqtau) / &
(S * S * sqrt_sigsqtau)
d2 = d1 - sqrt_sigsqtau
! greek = theta
BS_values(j,i,4) = exp((-dividend) * tau) * S * &
(dividend * anordf(d1) - 0.5 * sigsq * &
n01pdf_d1 / sqrt_sigsqtau) - r * &
exp((-r) * tau) * ks(strike_price) * &
anordf(d2)
! greek = charm
BS_values(j,i,5) = exp((-dividend) * tau) * ((-dividend)&
* anordf(d1) + n01pdf_d1 *&
((r - dividend) * tau - 0.5 * d2 *&
sqrt_sigsqtau) / (tau * sqrt_sigsqtau))
! greek = color
BS_values(j,i,6) = &
(-exp((-dividend) * tau)) * n01pdf_d1 *&
(2.0 * dividend * tau + 1.0 + d1 *&
(2.0 * (r - dividend) * tau - d2 *&
sqrt_sigsqtau) / sqrt_sigsqtau) / &
(2.0 * S * tau * sqrt_sigsqtau)
end do
end do
do l=1,12
relerrmax = 0.0
do i = 1,nv
do j = 1,nt
sVo = spline_values(i,j,l)
BSVo = BS_values(i,j,l)
relerr = abs((sVo - BSVo) / BSVo)
if (relerr > relerrmax) relerrmax = relerr
reavg(l) = reavg(l) + relerr
irect(l) = irect(l) + 1
end do
end do
if (relerrmax > rex(l)) rex(l) = relerrmax
if ((volatility == 1) .and. (strike_price == 1)) then
do i=1,nv
sVo_array(1:nt) = spline_values(i,1:nt,l)
write(nout,'("underlying price: ", F4.1,"; FK ",'//&
'A6,": ", 3(F13.10,TR1))') eval_points(i),&
greek_name(l),&
(sVo_array(k), k=1,nt)
write(nout, '(T25, "BS ", A6,": ", 3(F13.10,TR1))')&
greek_name(l), (BS_values(i,k,l), k=1,nt)
end do
end if
end do
end subroutine calc_Greeks
end
! These functions and routines define the coefficients, payoff and boundary conditions.
function fkcoef (x, tx, iflag, fcn_data)
use mp_types
implicit none
real(kind(1d0)), intent(in) :: x, tx
integer, intent(inout) :: iflag
type(d_fcn_data), optional :: fcn_data
real(kind(1d0)) :: fkcoef
real(kind(1d0)) :: sigma, interest_rate, alpha, dividend,&
half = 0.5d0
sigma = fcn_data % rdata(3)
alpha = fcn_data % rdata(4)
interest_rate = fcn_data % rdata(5)
dividend = fcn_data % rdata(6)
select case (iflag)
case (1)
! The coefficient derivative d(sigma)/dx
fkcoef = half*alpha*sigma*x**(alpha*half-1.0d0)
! The coefficient sigma(x)
case (2)
fkcoef = sigma*x**(alpha*half)
case (3)
! The coefficient mu(x)
fkcoef = (interest_rate - dividend) * x
case (4)
! The coefficient kappa(x)
fkcoef = interest_rate
end select
! Note that there is no time dependence
iflag = 0
return
end function fkcoef
function fkinitcond(x, fcn_data)
use mp_types
implicit none
real(kind(1d0)), intent(in) :: x
type (d_fcn_data), optional :: fcn_data
real(kind(1d0)) :: fkinitcond
real(kind(1d0)) :: zero = 0.0d0
real(kind(1d0)) :: strike_price
strike_price = fcn_data % rdata(1)
! The payoff function
fkinitcond = max(x - strike_price, zero)
return
end function fkinitcond
subroutine fkbc (tx, iflag, bccoefs, fcn_data)
use mp_types
implicit none
real(kind(1d0)), intent(in) :: tx
integer, intent(inout) :: iflag
real(kind(1d0)), dimension(:,:), intent(out) :: bccoefs
type (d_fcn_data), optional :: fcn_data
real(kind(1d0)) :: x_max, df, interest_rate, strike_price
strike_price = fcn_data % rdata(1)
x_max = fcn_data % rdata(2)
interest_rate = fcn_data % rdata(5)
select case (iflag)
case (1)
bccoefs(1,1:4) = (/1.0d0, 0.0d0, 0.0d0, 0.0d0/)
bccoefs(2,1:4) = (/0.0d0, 1.0d0, 0.0d0, 0.0d0/)
bccoefs(3,1:4) = (/0.0d0, 0.0d0, 1.0d0, 0.0d0/)
! Note no time dependence at left end
iflag = 0
case (2)
df = exp(interest_rate * tx)
bccoefs(1,1:4) = (/1.0d0, 0.0d0, 0.0d0, x_max -&
df*strike_price/)
bccoefs(2,1:4) = (/0.0d0, 1.0d0, 0.0d0, 1.0d0/)
bccoefs(3,1:4) = (/0.0d0, 0.0d0, 1.0d0, 0.0d0/)
end select
end subroutine fkbc
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
PHONE: 713.784.3131 FAX:713.781.9260 |