|
|
|
This function solves the generalized Feynman-Kac PDE on a rectangular grid using a finite element Galerkin method. Initial and boundary conditions are satisfied. The solution is represented by a series of C2 Hermite quintic splines.
Synopsis
#include <imsl.h>
void imsl_f_feynman_kac (int nxgrid, int ntgrid, int nlbcd, int nrbcd, float xgrid[], float tgrid[], void fcn_fkcfiv(), void fcn_fkbcp(), float y[], float y_prime[],…, 0)
The type double function is imsl_d_feynman_kac.
Required Arguments
int nxgrid
(Input)
Number of grid lines in the x-direction. nxgrid must be
at least 2.
int ntgrid
(Input)
Number of time points where an approximate solution is returned. The
value ntgrid
must be at least 1.
int nlbcd
(Input)
Number of left boundary conditions. It is required that 1 ≤ nlbcd ≤ 3.
int nrbcd
(Input)
Number of right boundary conditions. It is required that 1 ≤ nrbcd ≤ 3.
float xgrid[]
(Input)
Array of length nxgrid containing the
breakpoints for the Hermite quintic splines used in the x discretization.
The points in xgrid must be
strictly increasing. The values xgrid[0] and xgrid[nxgrid-1] are the endpoints of
the interval.
float tgrid[]
(Input)
Array of length ntgrid containing the
set of time points (in time-remaining units) where an approximate solution is
returned. The points in tgrid must be positive
and strictly increasing.
void fcn_fkcfiv
(float x,
float t,
int *iflag,
float *value)
User-supplied function to compute the values of the coefficients 
for the Feynman-Kac PDE and the
initial data function 
.
float x
(Input)
Space
variable.
float t
(Input)
Time
variable.
int *iflag
(Input/Output)
On entry, iflag indicates which
coefficient or data function is to be computed. The following table shows which
value has to be returned by fcn_fkcfiv for all
possible values of iflag:
|
iflag |
Computed coefficient |
|
-1 |
|
|
0 |
|
|
1 |
|
|
2 |
|
|
3 |
|
For non-zero input values of iflag, note when a coefficient does not depend on t. This is done 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.
float *value
(Output)
Value of the coefficient or initial data function. Which value is
computed depends on the input value for iflag, see description
of iflag.
void fcn_fkbcp
(int nbc,
float t, int *iflag,
float values[])
User-supplied
function to define boundary values that the solution of the differential
equation must satisfy. There are nlbcd conditions
specified at the left end, 
, and nrbcd conditions at
the right end, 
. The boundary conditions are
int nbc
(Input)
Number of boundary conditions. At xmin, nbc=nlbcd, at
xmax,
nbc = nrbcd.
float t (Input)
Time point of the
boundary conditions.
int *iflag
(Input/Output)
On entry, iflag indicates
whether the coefficients for the left or right boundary conditions are to be
computed:
|
iflag |
Computed boundary conditions |
|
1 |
Left end, x = xmin |
|
0 |
Right end, x = xmax |
If there is no time dependence for one of the boundaries then set iflag = 0 after the array values is defined for either end point. This flag can avoid unneeded continued computation of the finite element matrices.
float values[] (Output)
Array of
length 
containing the values of
the boundary condition coefficients in its first 4*nbc locations. The
coefficients for xmin are stored row-wise according to the
following scheme:
.
The coefficients for xmax are stored similarly.
float y[]
(Output)
An array of size (ntgrid+1) by (3*nxgrid) containing
the coefficients of the Hermite representation of the approximate solution for
the Feynman-Kac PDE at time points 0, tgrid[0], …, tgrid[ntgrid-1]. The approximate
solution is given by
for
The representation for the initial data at t = 0 is
.
The 
by 
matrix
is stored row-wise in array y.
After the integration,
use row i of
array y as input
argument coef in function
imsl_f_feynman_kac_evaluate
to obtain an array of values for f(x, t) or
its partials 
at time point t=0,tgrid[i-1], i=1,…,ntgrid.
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”.
float y_prime[]
(Output)
An array of size (ntgrid + 1)
by (3*nxgrid)
containing the first derivatives of y at time points
0,
tgrid[0],…,tgrid[ntgrid - 1], i.e.
for
,
and
for 
.
The 
by 
matrix
is stored row-wise in array y_prime.
After the integration, use row i of array y_prime as input
argument coef in function
imsl_f_feynman_kac_evaluate
to obtain an array of values for the partials 
at time point t=tgrid[i-1], i=1,…,ntgrid, and row 0 for the partials at
t=0.
Synopsis with Optional Arguments
#include <imsl.h>
void imsl_f_feynman_kac (int nxgrid, int ntgrid, int nlbcd, int nrbcd, float xgrid[], float tgrid[], void fcn_fkcfiv(), void fcn_fkbcp(), float y[], float y_prime[],
IMSL_FCN_FKCFIV_W_DATA, void fcn_fkcfiv(), void *data,
IMSL_FCN_FKBCP_W_DATA, void fcn_fkbcp(), void *data,
IMSL_FCN_INIT, void fcn_fkinit(),
IMSL_FCN_INIT_W_DATA, void fcn_fkinit(), void *data,
IMSL_FCN_FORCE, void fcn_fkforce(),
IMSL_FCN_FORCE_W_DATA, void fcn_fkforce(), void *data,
IMSL_ATOL_RTOL_SCALARS, float atol, float rtol,
IMSL_ATOL_RTOL_ARRAYS, float atol[], float rtol[]
IMSL_NDEGREE, int ndeg,
IMSL_TDEPEND, int tdepend[],
IMSL_MAX_STEP, float max_stepsize,
IMSL_INITIAL_STEPSIZE, float init_stepsize,
IMSL_MAX_NUMBER_STEPS, int max_steps,
IMSL_STEP_CONTROL, int step_control,
IMSL_MAX_BDF_ORDER, int max_bdf_order,
IMSL_T_BARRIER, float t_barrier,
IMSL_ISTATE, int istate[],
IMSL_EVALS, int nval[],
0)
Optional Arguments
IMSL_FCN_FKCFIV_W_DATA,
void fcn_fkcfiv(float x,
float t,
int *iflag,
float *value,
void *data), void *data
(Input)
User-supplied function to compute the values of the coefficients
for the Feynman-Kac PDE and the
initial data function 
. This function also accepts a
pointer to the object data supplied by the
user.
IMSL_FCN_FKBCP_W_DATA,
void fcn_fkbcp(int nbc,
float t, int *iflag,
float values[],
void *data),
void *data
(Input)
User-supplied function to define boundary values that the solution of
the differential equation must satisfy. This function also accepts a pointer to
data supplied by
the user.
IMSL_FCN_INIT,
void fcn_fkinit(int nxgrid,
int ntgrid, float xgrid[],
float tgrid[],
float time,
float yprime[],
float y[],
float atol[], float rtol[])
(Input)
User-supplied function for adjustment of initial data or
as an opportunity for output during the integration steps. The solution values
of the model parameters are presented in the arrays y and yprime at the
integration points time = 0, tgrid[j], j=0,…,ntgrid-1. At the
initial point, integral least-squares estimates are made for representing the
initial data 
. If this is not satisfactory, fcn_fkinit can change
the contents of y[] to match data for
any reason.
int nxgrid
(Input)
Number of grid lines in the x-direction.
int ntgrid
(Input)
Number of time points where an approximate solution is returned.
float xgrid[]
(Input)
Vector of length nxgrid containing the
breakpoints for the Hermite quintic splines used in the x
discretization.
float tgrid[]
(Input)
Vector of length ntgrid containing the
set of time points (in time-remaining units) where an approximate solution is
returned.
float time
(Input)
Time variable.
float yprime[]
(Input)
Vector of length 3*nxgrid containing the
derivative of the coefficients of the Hermite quintic spline at time point time.
float y[]
(Input/Output)
Vector of length 3*nxgrid containing the
coefficients of the Hermite quintic spline at time point time.
float atol[]
(Input/Output)
Array of length
3*nxgrid
containing absolute error tolerances.
float rtol[]
(Input/Output)
Array of length 3*nxgrid containing
relative error tolerances.
IMSL_FCN_INIT_W_DATA,
void fcn_fkinit(int nxgrid,
int ntgrid,
float xgrid[], float tgrid[],
float time,
float yprime[], float y[],
float atol[],
float rtol[],
void *data),
void *data
(Input)
User-supplied function for adjustment of initial data or as an
opportunity for output during the integration steps which also accepts a pointer
to data supplied by the user. For an explanation of the other arguments of
function fcn_fkinit, see
optional argument IMSL_FCN_INIT.
IMSL_FCN_FORCE,
void fcn_fkforce(int interval,
int ndeg, int nxgrid,
float y[],
float time,
float width,
float xlocal[],
float qw[],
float u[],
float phi[],
float dphi[])
(Input)
Function fcn_fkforce is used in
case there is a non-zero term 
in the Feynman-Kac
differential equation. If function fcn_fkforce is not used, it is assumed that is identically zero.
int interval
(Input)
Index
indicating the bounds xgrid[interval-1] and xgrid[interval] of the
integration interval, 1≤interval ≤ nxgrid-1.
int ndeg
(Input)
The degree used for the Gauss-Legendre formulas.
int nxgrid
(Input)
Number of grid lines in the x-direction.
float y[]
(Input)
Vector of length 3*nxgrid containing
the coefficients of the Hermite quintic spline representing the solution of the
Feynman-Kac equation at time point time.
float time
(Input)
Time variable.
float width (Input)
The interval
width, width = xgrid[interval] - xgrid[interval-1].
float xlocal[]
(Input)
Array of length ndeg containing the
Gauss-Legendre points translated and normalized to the interval [xgrid[interval-1], xgrid[interval]].
float qw[]
(Input)
Vector of
length ndeg
containing the Gauss-Legendre weights.
float u[]
(Input)
Array of dimension 12 by ndeg containing the
basis function values that define 
at the Gauss-Legendre points
xlocal. Setting
u[i+k*ndeg] and 
xlocal[i], 
is defined as
float phi[]
(Output)
Vector of length 6 containing Gauss-Legendre approximations for the
local contributions
where 
time and
.
Vector phi contains elements
,
for i=0,…,5.
float dphi[]
(Output)
Array of dimension 6 by 6 containing a
Gauss-Legendre approximation for the Jacobian of the local contributions
at time,
The approximation to this symmetric matrix is stored row-wise, i.e.
for i, j = 0, …, 5.
IMSL_FCN_FORCE_W_DATA,
void fcn_fkforce(int interval,
int ndeg,
int nxgrid,
float y[],
float time,
float width,
floa xlocal[],
float qw[],
float u[],
float phi[],
float dphi[],
void *data),
void *data
(Input)
Function fcn_fkforce is used in
case there is a non-zero term in the Feynman-Kac
differential equation. This function also accepts a data pointer to data
supplied by the user. For an explanation of the other arguments of function
fcn_fkforce, see
optional argument IMSL_FCN_FORCE.
IMSL_ATOL_RTOL_SCALARS,
float atol,
float rtol
(Input)
Scalar values that apply to the error estimates of all components of
the solution y
in the differential equation solver SDASLX. See optional
argument IMSL_ATOL_RTOL_ARRAYS
if separate tolerances are to be applied to each component of y.
Default: atol and rtol are set to 
in single precision and 
in double precision.
IMSL_ATOL_RTOL_ARRAYS,
float atol[],
float rtol[],
(Input)
Componentwise tolerances are used for the computation of solution
y in the
differential equation solver SDASLX. Arguments
atol and rtol are pointers to arrays of length 3*nxgrid to be used for
the absolute and relative tolerance and to be applied to each component of the
solution, y. See
optional argument IMSL_ATOL_RTOL_SCALARS
if scalar values of absolute and relative tolerances are to be
applied to all components.
Default: All elements of atol and rtol are set to
10-3 in
single precision and 10-5 in double precision.
IMSL_NDEGREE,
int ndeg
(Input)
The degree used for the Gauss-Legendre formulas for constructing the
finite element matrices. It is required that ndeg ≥ 6.
Default: ndeg = 6.
IMSL_TDEPEND,
int tdepend[]
(Output)
Vector of length 7 indicating time dependence of the coefficients,
boundary conditions and function 
in the Feynman-Kac equation.
If tdepend[i] = 0
then argument i
is not time dependent, if tdepend[i]=1 then
argument i is
time dependent.
|
i |
Computed coefficient |
|
0 |
|
|
1 |
|
|
2 |
|
|
3 |
|
|
4 |
Left boundary conditions |
|
5 |
Right boundary conditions |
|
6 |
|
IMSL_MAX_STEP,
float max_stepsize
(Input)
This is the maximum step size the integrator may take.
Default:
max_stepsize = imsl_f_machine(2), the
largest possible machine number.
IMSL_INITIAL_STEPSIZE,
float init_stepsize
(Input)
The starting step size for the integration. Must be less than zero
since the integration is internally done from t=0 to t=tgrid[ntgrid-1] in a
negative direction.
Default: init_stepsize = 
, where
is the machine precision.
IMSL_MAX_NUMBER_STEPS,
int max_steps
(Input)
The maximum number of steps between each output point of the
integration.
Default: maxsteps = 500000.
IMSL_STEP_CONTROL, int step_control
(Input)
Indicates which step control algorithm is used for the
integration. If step_control = 0, then the step control method of
Söderlind is used. If step_control = 1, then the method used by the
original Petzold code SASSL is
used.
Default: step_control = 0.
IMSL_MAX_BDF_ORDER,
int max_bdf_order
(Input)
Maximum order of the backward differentiation formulas used in the
integrator. It is required that 1≤ max_bdf_order ≤ 5.
Default:
max_bdf_order = 5.
IMSL_T_BARRIER, float t_barrier
(Input)
This optional argument controls whether the code should integrate
past a special point, t_barrier, and then
interpolate to get 
and 
at the points in tgrid[]. If this
optional argument is present, the integrator assumes the equations either change
on the alternate sides of t_barrier or they are undefined there. In
this case, the code creeps up to t_barrier in the
direction of integration. It is required that t_barrier ≥tgrid[ntgrid-1].
Default:
t_barrier = tgrid[ntgrid-1].
IMSL_ISTATE, int istate[]
(Output)
An
array of length 7 whose entries flag the state of computation for the matrices
and vectors required in the integration. For each entry, a zero indicates
that no computation has been done or 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 |
Entry in istate |
|
0 |
Mass Matrix, M |
|
1 |
Stiffness matrix, N |
|
2 |
Bending matrix, R |
|
3 |
Weighted mass, K |
|
4 |
Left boundary conditions |
|
5 |
Right- boundary conditions |
|
6 |
Forcing term |
Default: istate[i] = 0 for I = 0,…,6.
IMSL_EVALS,
int nval[]
(Output)
Array of length 3 summarizing the number of evaluations required
during the integration.
|
i |
nval[i] |
|
0 |
Number of residual function evaluations of the DAE used in the model. |
|
1 |
Number of factorizations of the differential matrix associated with solving the DAE. |
|
2 |
Number of linear system solve steps using the differential matrix. |
Description
The generalized Feynman-Kac differential equation has the form
where the initial data satisfies
. The derivatives are
, etc
The function imsl_f_feynman_kac uses a finite element Galerkin method over the rectangle
in 
to compute the approximate
solution. The interval 
is decomposed with a grid
On each subinterval the solution is represented by
The values
are time-dependent coefficients associated with each
interval. The basis functions 
are given for
by
The Galerkin principle is then applied. Using the provided initial and boundary conditions leads to an Index 1 differential-algebraic equation for the time-dependent coefficients
.
This system is integrated using the variable order, variable step algorithm DDASLX/SDASLX noted in Hanson and Krogh, R. (2008) Solving Constrained Differential-Algebraic Systems Using Projections. Solution values and their time derivatives are returned at a grid preceding time T, expressed in units of time remaining. For further details of deriving and solving the system see Hanson, R. (2008) Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements. To evaluate f or its partial derivatives at any time point in the grid, use the function imsl_f_feynman_kac_evaluate.
Examples
Example 1
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. The coefficients, payoff, boundary conditions and forcing term for American and European options are defined by the user functions fkcfiv_put, fkbcp_put and fkforce_put, respectively. 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”, 
.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define max(A,B) ((A) >= (B) ? (A) : (B))
#define NXGRID 61
#define NTGRID 2
#define NV 9
void fkcfiv_put(float, float, int *, float *);
void fkbcp_put(int, float, int *, float[]);
void fkforce_put(int, int, int, float[], float, float, float[],
float[], float[], float[],float[], void *);
int main()
{
/* Compute American Option Premium for Vanilla Put */
/* The strike price */
float KS = 10.0;
/* The sigma value */
float sigma = 0.4;
/* Time values for the options */
int nt = 2;
float time[] = { 0.25, 0.5};
/* Values of the underlying where evaluations are made */
float xs[] = { 0.0, 2.0, 4.0, 6.0, 8.0, 10.0,
12.0, 14.0, 16.0 };
/* Value of the interest rate */
float r = 0.1;
/* Values of the min and max underlying values modeled */
float x_min=0.0, x_max=30.0;
/* Define parameters for the integration step. */
int nint=NXGRID-1, n=3*NXGRID;
float xgrid[NXGRID];
float ye[(NTGRID+1)*3*NXGRID],yeprime[(NTGRID+1)*3*NXGRID];
float ya[(NTGRID+1)*3*NXGRID], yaprime[(NTGRID+1)*3*NXGRID];
float fe[NTGRID*NV], fa[NTGRID*NV];
float dx;
int i;
int nlbcd = 2, nrbcd = 3;
float atol = 0.2e-2;
/* Array for user-defined data */
float usr_data[3];
/* Define an equally-spaced grid of points for the
underlying price */
dx = (x_max-x_min)/((float) nint);
xgrid[0]=x_min;
xgrid[NXGRID-1]=x_max;
for (i=2; i<=NXGRID-1; i++) xgrid[i-1]=xgrid[i-2]+dx;
usr_data[0] = KS;
usr_data[1] = r;
usr_data[2] = atol;
imsl_f_feynman_kac(NXGRID, nt, nlbcd, nrbcd, xgrid, time,
fkcfiv_put, fkbcp_put, ye, yeprime,
IMSL_ATOL_RTOL_SCALARS, 0.5e-2, 0.5e-2, 0);
imsl_f_feynman_kac(NXGRID, nt, nlbcd, nrbcd, xgrid, time,
fkcfiv_put, fkbcp_put, ya, yaprime,
IMSL_FCN_FORCE_W_DATA, fkforce_put, usr_data,
IMSL_ATOL_RTOL_SCALARS, 0.5e-2, 0.5e-2, 0);
/* Evaluate solutions at vector of points XS(:), at each
time value prior to expiration. */
for (i=0; i<nt; i++)
{
imsl_f_feynman_kac_evaluate (NV, NXGRID, xgrid, xs, &ye[(i+1)*n],
IMSL_RETURN_USER, &fe[i*NV], 0);
imsl_f_feynman_kac_evaluate (NV, NXGRID, xgrid, xs, &ya[(i+1)*n],
IMSL_RETURN_USER, &fa[i*NV], 0);
}
printf("\nAmerican Option Premium for Vanilla Put, "
"3 and 6 Months Prior to Expiry\n");
printf("%7sNumber of equally spaced spline knots:%4d\n", " ",
NXGRID);
printf("%7sNumber of unknowns:%4d\n", " ", n);
printf("%7sStrike=%6.2f, sigma=%5.2f, Interest Rate=%5.2f\n\n",
" ",KS,sigma,r);
printf("%7s%10s%20s%20s\n", " ","Underlying","European","American");
for (i=0; i<NV; i++)
printf("%7s%10.4f%10.4f%10.4f%10.4f%10.4f\n", " ", xs[i], fe[i],
fe[i+NV], fa[i], fa[i+NV]);
}
/* These functions define the coefficients, payoff, boundary conditions
and forcing term for American and European Options. */
void fkcfiv_put(float x, float tx, int *iflag, float *value)
{
/* The sigma value */
float sigma=0.4;
/* Value of the interest rate and continuous dividend */
float strike_price=10.0, interest_rate=0.1, dividend=0.0;
float zero=0.0;
switch (*iflag)
{
case 0:
/* The payoff function */
*value = max(strike_price - x, zero);
break;
case -1:
/* The coefficient derivative d sigma/ dx */
*value = sigma;
break;
case 1:
/* The coefficient sigma(x) */
*value = sigma*x;
break;
case 2:
/* The coefficient mu(x) */
*value = (interest_rate - dividend) * x;
break;
case 3:
/* The coefficient kappa(x) */
*value = interest_rate;
break;
}
/* Note that there is no time dependence */
*iflag = 0;
return;
}s
void fkbcp_put(int nbc, float tx, int *iflag, float val[])
{
switch (*iflag)
{
case 1:
val[0] = 0.0; val[1] = 1.0; val[2] = 0.0;
val[3] = -1.0; val[4] = 0.0; val[5] = 0.0;
val[6] = 1.0; val[7] = 0.0;
break;
case 2:
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = 0.0; val[4] = 0.0; val[5] = 1.0;
val[6] = 0.0; val[7] = 0.0; val[8] = 0.0;
val[9] = 0.0; val[10] = 1.0; val[11] = 0.0;
break;
}
/* Note no time dependence */
*iflag = 0;
return;
}
void fkforce_put(int interval, int ndeg, int nxgrid,
float y[], float time, float width, float xlocal[],
float qw[], float u[], float phi[],float dphi[],
void *data_ptr)
{
int i, j, k, l;
const int local=6;
float yl[6], bf[6];
float value, strike_price, interest_rate, zero=0.0, one=1.0;
float rt, mu;
float *data = NULL;
data = (float *) data_ptr;
for (i=0; i<local; i++)
{
yl[i] = y[3*interval-3+i];
phi[i] = zero;
}
strike_price = data[0];
interest_rate = data[1];
value = data[2];
mu=2.0;
/* This is the local definition of the forcing term */
for (j=1; j<=local; j++)
for (l=1; l<=ndeg; l++)
{
bf[0] = u[(l-1)];
bf[1] = u[(l-1)+ndeg];
bf[2] = u[(l-1)+2*ndeg];
bf[3] = u[(l-1)+6*ndeg];
bf[4] = u[(l-1)+7*ndeg];
bf[5] = u[(l-1)+8*ndeg];
rt = 0.0;
for (k=0; k<local; k++)
rt += yl[k]*bf[k];
rt = value/(rt + value - (strike_price-xlocal[l-1]));
phi[j-1] += qw[l-1] * bf[j-1] * pow(rt,mu);
}
for (i=0; i<local; i++)
phi[i] = -phi[i]*width*interest_rate*strike_price;
/* This is the local derivative matrix for the forcing term */
for (i=0; i<local*local; i++)
dphi[i] = zero;
for (j=1; j<=local; j++)
for (i=1; i<=local; i++)
for (l=1; l<=ndeg; l++)
{
bf[0] = u[(l-1)];
bf[1] = u[(l-1)+ndeg];
bf[2] = u[(l-1)+2*ndeg];
bf[3] = u[(l-1)+6*ndeg];
bf[4] = u[(l-1)+7*ndeg];
bf[5] = u[(l-1)+8*ndeg];
rt = 0.0;
for (k=0; k<local; k++)
rt += yl[k]*bf[k];
rt = one/(rt + value-(strike_price-xlocal[l-1]));
dphi[i-1+(j-1)*local] += qw[l-1] * bf[i-1] * bf[j-1]
* pow(rt, mu+1.0);
}
for (i=0; i<local*local; i++)
dphi[i] = mu * dphi[i] * width * pow(value, mu) *
interest_rate * strike_price;
return;
}
Output
American Option Premium for Vanilla Put, 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
Underlying European American
0.0000 9.7534 9.5136 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.7615 3.5683 4.0000 4.0000
8.0000 1.9070 1.9168 2.0170 2.0865
10.0000 0.6529 0.8548 0.6771 0.9030
12.0000 0.1632 0.3371 0.1680 0.3519
14.0000 0.0372 0.1270 0.0373 0.1321
16.0000 0.0088 0.0483 0.0084 0.0501
Example 2
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”, 
.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define max(A,B) ((A) >= (B) ? (A) : (B))
#define NXGRID 121
#define NTGRID 3
#define NV 3
void fcn_fkcfiv(float, float, int *, float *, void *);
void fcn_fkbcp(int, float, int *, float[], void *);
int main()
{
/* Compute Constant Elasticity of Variance Model for Vanilla Call */
/* The set of strike prices */
float KS[] = { 15.0, 20.0, 25.0};
/* The set of sigma values */
float sigma[] = { 0.2, 0.3, 0.4};
/* The set of model diffusion powers */
float alpha[] = { 2.0, 1.0, 0.0};
/* Time values for the options */
int nt = 3;
float time[] = { 1.0/12.0, 4.0/12.0, 7.0/12.0 };
/* Values of the underlying where evaluations are made */
float xs[] = { 19.0, 20.0, 21.0};
/* Value of the interest rate and continuous dividend */
float r=0.05, dividend=0.0;
/* Values of the min and max underlying values modeled */
float x_min=0.0, x_max=60.0;
/* Define parameters for the integration step. */
int nint = NXGRID-1, n=3*NXGRID;
float xgrid[NXGRID], y[(NTGRID+1)*3*NXGRID];
float yprime[(NTGRID+1)*3*NXGRID], f[NTGRID*NV];
float dx;
/* Number of left/right boundary conditions */
int nlbcd = 3, nrbcd = 3;
float usr_data[6];
int i,i1,i2,i3,j;
/* Define equally-spaced grid of points for the underlying price */
dx = (x_max-x_min)/((float) nint);
xgrid[0] = x_min;
xgrid[NXGRID-1] = x_max;
for (i=2; i<=NXGRID-1; i++)
xgrid[i-1] = xgrid[i-2]+dx;
printf("%2sConstant Elasticity of Variance Model for Vanilla Call\n",
" ");
printf("%7sInterest Rate:%7.3f\tContinuous Dividend:%7.3f\n",
" ", r, dividend);
printf("%7sMinimum and Maximum Prices of Underlying:%7.2f%7.2f\n",
" ",x_min, x_max);
printf("%7sNumber of equally spaced spline knots:%4d \n", " ",
NXGRID-1);
printf("%7sNumber of unknowns:%4d\n\n", " ",n);
printf("%7sTime in Years Prior to Expiration:%7.4f%7.4f%7.4f\n",
" ",time[0], time[1], time[2]);
printf("%7sOption valued at Underlying Prices:%7.2f%7.2f%7.2f\n\n",
" ", xs[0], xs[1], xs[2]);
for (i1=1; i1<=3; i1++) /* Loop over power */
for (i2=1; i2<=3; i2++) /* Loop over volatility */
for (i3=1; i3<=3; i3++) /* Loop over strike price */
{
/* Pass data through into evaluation functions. */
usr_data[0] = KS[i3-1];
usr_data[1] = x_max;
usr_data[2] = sigma[i2-1];
usr_data[3] = alpha[i1-1];
usr_data[4] = r;
usr_data[5] = dividend;
imsl_f_feynman_kac(NXGRID, nt, nlbcd, nrbcd, xgrid, time,
NULL, NULL, y, yprime,
IMSL_FCN_FKCFIV_W_DATA, fcn_fkcfiv, usr_data,
IMSL_FCN_FKBCP_W_DATA, fcn_fkbcp, usr_data, 0);
/* Evaluate solution at vector of points xs, at each time
Value prior to expiration. */
for (i=0; i<NTGRID; i++)
imsl_f_feynman_kac_evaluate (NV, NXGRID, xgrid, xs,
&y[(i+1)*n], IMSL_RETURN_USER, &f[i*NV], 0);
printf("%2sStrike=%5.2f, Sigma=%5.2f, Alpha=%5.2f\n",
" ", KS[i3-1], sigma[i2-1], alpha[i1-1]);
for (i=0; i<NV; i++)
{
printf("%23sCall Option Values%2s", " ", " ");
for (j=0; j<nt; j++) printf("%7.4f ", f[j*NV+i]);
printf("\n");
}
printf("\n");
}
}
void fcn_fkcfiv(float x, float tx, int *iflag, float *value,
void *data_ptr)
{
float sigma, strike_price, interest_rate;
float alpha, dividend, zero=0.0, half=0.5;
float *data = NULL;
data = (float *)data_ptr;
strike_price = data[0];
sigma = data[2];
alpha = data[3];
interest_rate = data[4];
dividend = data[5];
switch (*iflag)
{
case 0:
/* The payoff function */
*value = max(x - strike_price, zero);
break;
case -1:
/* The coefficient derivative d sigma/ dx */
*value = half * alpha * sigma * pow(x, alpha*half-1.0);
break;
case 1:
/* The coefficient sigma(x) */
*value = sigma * pow(x, alpha*half);
break;
case 2:
/* The coefficient mu(x) */
*value = (interest_rate - dividend) * x;
break;
case 3:
/* The coefficient kappa(x) */
*value = interest_rate;
break;
}
data = NULL;
/* Note that there is no time dependence */
*iflag = 0;
return;
}
void fcn_fkbcp(int nbc, float tx, int *iflag, float val[],
void *data_ptr)
{
float x_max, df, interest_rate, strike_price;
float *data = NULL;
data = (float *)data_ptr;
strike_price = data[0];
x_max = data[1];
interest_rate = data[4];
switch (*iflag)
{
case 1:
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = 0.0; val[4] = 0.0; val[5] = 1.0;
val[6] = 0.0; val[7] = 0.0; val[8] = 0.0;
val[9] = 0.0; val[10] = 1.0; val[11] = 0.0;
/* Note no time dependence at left end */
*iflag = 0;
break;
case 2:
df = exp(interest_rate*tx);
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = x_max - df*strike_price; val[4] = 0.0;
val[5] = 1.0; val[6] = 0.0; val[7] = 1.0;
val[8] = 0.0; val[9] = 0.0; val[10] = 1.0;
val[11] = 0.0;
break;
}
data = NULL;
return;
}
Output
Constant Elasticity of Variance Model for Vanilla Call
Interest Rate: 0.050 Continuous Dividend: 0.000
Minimum and Maximum Prices of Underlying: 0.00 60.00
Number of equally spaced spline knots: 120
Number of unknowns: 363
Time in Years Prior to Expiration: 0.0833 0.3333 0.5833
Option valued at Underlying Prices: 19.00 20.00 21.00
Strike=15.00, Sigma= 0.20, Alpha= 2.00
Call Option Values 4.0624 4.2577 4.4730
Call Option Values 5.0624 5.2507 5.4491
Call Option Values 6.0624 6.2487 6.4386
Strike=20.00, Sigma= 0.20, Alpha= 2.00
Call Option Values 0.1310 0.5956 0.9699
Call Option Values 0.5028 1.0889 1.5100
Call Option Values 1.1979 1.7485 2.1751
Strike=25.00, Sigma= 0.20, Alpha= 2.00
Call Option Values 0.0000 0.0113 0.0742
Call Option Values 0.0000 0.0372 0.1614
Call Option Values 0.0007 0.1025 0.3132
Strike=15.00, Sigma= 0.30, Alpha= 2.00
Call Option Values 4.0637 4.3399 4.6623
Call Option Values 5.0626 5.2945 5.5788
Call Option Values 6.0624 6.2709 6.5241
Strike=20.00, Sigma= 0.30, Alpha= 2.00
Call Option Values 0.3103 1.0274 1.5500
Call Option Values 0.7315 1.5422 2.1024
Call Option Values 1.3758 2.1689 2.7385
Strike=25.00, Sigma= 0.30, Alpha= 2.00
Call Option Values 0.0006 0.1112 0.3547
Call Option Values 0.0038 0.2170 0.5552
Call Option Values 0.0184 0.3857 0.8225
Strike=15.00, Sigma= 0.40, Alpha= 2.00
Call Option Values 4.0759 4.5136 4.9673
Call Option Values 5.0664 5.4199 5.8321
Call Option Values 6.0635 6.3577 6.7294
Strike=20.00, Sigma= 0.40, Alpha= 2.00
Call Option Values 0.5116 1.4645 2.1286
Call Option Values 0.9623 1.9957 2.6944
Call Option Values 1.5815 2.6110 3.3230
Strike=25.00, Sigma= 0.40, Alpha= 2.00
Call Option Values 0.0083 0.3288 0.7790
Call Option Values 0.0285 0.5169 1.0657
Call Option Values 0.0813 0.7688 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.0241 0.1061
Call Option Values 0.1498 0.4102 0.6483
Call Option Values 1.0832 1.3313 1.5772
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.4310
Call Option Values 5.0624 5.2477 5.4310
Call Option Values 6.0624 6.2477 6.4310
Strike=20.00, Sigma= 0.30, Alpha= 1.00
Call Option Values 0.0016 0.0812 0.2214
Call Option Values 0.1981 0.4981 0.7535
Call Option Values 1.0836 1.3441 1.6018
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.0072 0.1556 0.3445
Call Option Values 0.2501 0.5919 0.8720
Call Option Values 1.0867 1.3783 1.6577
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.0007
Call Option Values 0.0000 0.0002 0.0059
Strike=15.00, Sigma= 0.20, Alpha= 0.00
Call Option Values 4.0625 4.2479 4.4311
Call Option Values 5.0625 5.2479 5.4312
Call Option Values 6.0623 6.2479 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.0624 4.2479 4.4311
Call Option Values 5.0625 5.2479 5.4311
Call Option Values 6.0623 6.2479 6.4311
Strike=20.00, Sigma= 0.30, Alpha= 0.00
Call Option Values 0.0000 -0.0000 0.0029
Call Option Values 0.0895 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.0624 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.0001 0.0111
Call Option Values 0.0985 0.3383 0.5781
Call Option Values 1.0830 1.3306 1.5749
Strike=25.00, Sigma= 0.40, Alpha= 0.00
Call Option Values 0.0000 0.0000 -0.0000
Call Option Values 0.0000 -0.0000 -0.0000
Call Option Values 0.0000 -0.0000 -0.0000
Example 3
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 =61, n = 3 × nx = 183.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define max(A,B) ((A) >= (B) ? (A) : (B))
#define NXGRID 61
#define NTGRID 2
#define NV 12
void fkcfiv_call(float, float, int *, float *, void *);
void fkbcp_call(int, float, int *, float[], void *);
int main()
{
int i;
/* The strike price */
float KS1 = 10.0;
/* The spread value */
float KS2 = 15.0;
/* The Bet for the Cash-or-Nothing Call */
float bet = 2.0;
/* The sigma value */
float sigma = 0.4;
/* Time values for the options */
int nt = 2;
float time[] = {0.25, 0.5};
/* Values of the underlying where evaluations are made */
float xs[NV];
/* Value of the interest rate and continuous dividend */
float r = 0.1, dividend=0.0;
/* Values of the min and max underlying values modeled */
float x_min=0.0, x_max=30.0;
/* Define parameters for the integration step. */
int nint=NXGRID-1, n=3*NXGRID;
float xgrid[NXGRID];
float yb[(NTGRID+1)*3*NXGRID], ybprime[(NTGRID+1)*3*NXGRID];
float yv[(NTGRID+1)*3*NXGRID], yvprime[(NTGRID+1)*3*NXGRID];
float fb[NTGRID*NV], fv[NTGRID*NV];
float dx;
/* Number of left/right boundary conditions. */
int nlbcd = 3, nrbcd = 3;
/* Structure for the evaluation functions. */
struct {
int idope[1];
float rdope[7];
} usr_data;
/* Define an equally-spaced grid of points for the underlying
price */
dx = (x_max-x_min)/((float)(nint));
xgrid[0]=x_min;
xgrid[NXGRID-1]=x_max;
for (i=2; i<=NXGRID-1; i++)
xgrid[i-1] = xgrid[i-2]+dx;
for (i=1; i<=NV; i++)
xs[i-1] = 2.0+(i-1)*2.0;
usr_data.rdope[0] = KS1;
usr_data.rdope[1] = bet;
usr_data.rdope[2] = KS2;
usr_data.rdope[3] = x_max;
usr_data.rdope[4] = sigma;
usr_data.rdope[5] = r;
usr_data.rdope[6] = dividend;
/* Flag the difference in payoff functions */
/* 1 for the Bet, 2 for the Vertical Spread */
usr_data.idope[0] = 1;
imsl_f_feynman_kac(NXGRID, NTGRID, nlbcd, nrbcd, xgrid, time,
NULL, NULL, yb, ybprime,
IMSL_FCN_FKCFIV_W_DATA, fkcfiv_call, &usr_data,
IMSL_FCN_FKBCP_W_DATA, fkbcp_call, &usr_data,
0);
usr_data.idope[0] = 2;
imsl_f_feynman_kac(NXGRID, NTGRID, nlbcd, nrbcd, xgrid, time,
NULL, NULL, yv, yvprime,
IMSL_FCN_FKCFIV_W_DATA, fkcfiv_call, &usr_data,
IMSL_FCN_FKBCP_W_DATA, fkbcp_call, &usr_data,
0);
/* Evaluate solutions at vector of points XS(:), at each time value
prior to expiration. */
for (i=0; i<NTGRID; i++)
{
imsl_f_feynman_kac_evaluate (NV, NXGRID, xgrid, xs, &yb[(i+1)*n],
IMSL_RETURN_USER, &fb[i*NV], 0);
imsl_f_feynman_kac_evaluate (NV, NXGRID, xgrid, xs, &yv[(i+1)*n],
IMSL_RETURN_USER, &fv[i*NV], 0);
}
printf("%2sEuropean Option Value for A Bet\n", " ");
printf("%3sand a Vertical Spread, 3 and 6 Months Prior to Expiry\n",
" ");
printf("%5sNumber of equally spaced spline knots:%4d\n", " ",
NXGRID);
printf("%5sNumber of unknowns:%4d\n", " ", n);
printf("%5sStrike=%5.2f, Sigma=%5.2f, Interest Rate=%5.2f\n",
" ", KS1, sigma, r);
printf("%5sBet=%5.2f, Spread Value=%5.2f\n\n", " ", bet, KS2);
printf("%17s%18s%18s\n", "Underlying", "A Bet", "Vertical Spread");
for (i=0; i<NV; i++)
printf("%8s%9.4f%9.4f%9.4f%9.4f%9.4f\n", " ", xs[i], fb[i],
fb[i+NV], fv[i], fv[i+NV]);
}
/* These functions define the coefficients, payoff, boundary conditions
and forcing term for American and European Options. */
void fkcfiv_call(float x, float tx, int *iflag, float *value,
void *data_ptr)
{
float sigma, strike_price, interest_rate;
float spread, bet, dividend, zero=0.0;
float *data_real = NULL;
int *data_int = NULL;
struct struct_data {
int idope[1];
float rdope[7];
};
struct struct_data *data = NULL;
data = data_ptr;
data_real = data->rdope;
data_int = data->idope;
strike_price = data_real[0];
bet = data_real[1];
spread = data_real[2];
sigma = data_real[4];
interest_rate = data_real[5];
dividend = data_real[6];
switch (*iflag)
{
case 0:
/* The payoff function - Use flag passed to decide which */
switch (data_int[0])
{
case 1:
/* After reaching the strike price the payoff jumps
from zero to the bet value. */
*value = zero;
if (x > strike_price) *value = bet;
break;
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);
break;
}
break;
case -1:
/* The coefficient derivative d sigma/ dx */
*value = sigma;
break;
case 1:
/* The coefficient sigma(x) */
*value = sigma*x;
break;
case 2:
/* The coefficient mu(x) */
*value = (interest_rate - dividend)*x;
break;
case 3:
/* The coefficient kappa(x) */
*value = interest_rate;
break;
}
/* Note that there is no time dependence */
*iflag = 0;
data_real = NULL;
data_int = NULL;
data = NULL;
return;
}
void fkbcp_call(int nbc, float tx, int *iflag, float
val[],
void *data_ptr)
{
float strike_price, spread, bet, interest_rate, df;
int *data_int = NULL;
float *data_real = NULL;
struct struct_data {
int idope[1];
float rdope[7];
};
struct struct_data *data = NULL;
data = data_ptr;
data_int = data->idope;
data_real = data->rdope;
strike_price = data_real[0];
bet = data_real[1];
spread = data_real[2];
interest_rate = data_real[5];
switch (*iflag)
{
case 1:
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = 0.0; val[4] = 0.0; val[5] = 1.0;
val[6] = 0.0; val[7] = 0.0; val[8] = 0.0;
val[9] = 0.0; val[10] = 1.0; val[11] = 0.0;
/* Note no time dependence in case (1) for IFLAG */
*iflag = 0;
break;
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 */
switch (data_int[0])
{
case 1:
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = bet*df;
break;
case 2:
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = (spread-strike_price)*df;
break;
}
val[4] = 0.0; val[5] = 1.0; val[6] = 0.0;
val[7] = 0.0; val[8] = 0.0; val[9] = 0.0;
val[10] = 1.0; val[11] = 0.0;
break;
}
data_real = NULL;
data_int = NULL;
data = NULL;
return;
}
Output
European Option Value for A Bet
and a Vertical Spread, 3 and 6 Months Prior to Expiry
Number of equally spaced spline knots: 61
Number of unknowns: 183
Strike=10.00, Sigma= 0.40, Interest Rate= 0.10
Bet= 2.00, Spread Value=15.00
Underlying A Bet Vertical Spread
2.0000 0.0000 0.0000 0.0000 0.0000
4.0000 0.0000 0.0014 0.0000 0.0006
6.0000 0.0110 0.0722 0.0039 0.0446
8.0000 0.2690 0.4304 0.1479 0.3831
10.0000 0.9948 0.9781 0.8909 1.1927
12.0000 1.6095 1.4287 2.1911 2.2274
14.0000 1.8654 1.6924 3.4255 3.1551
16.0000 1.9337 1.8177 4.2264 3.8263
18.0000 1.9476 1.8700 4.6264 4.2492
20.0000 1.9501 1.8903 4.7911 4.4922
22.0000 1.9505 1.8979 4.8497 4.6232
24.0000 1.9506 1.9007 4.8685 4.6909
Example 4
This example evaluates the price of a convertible bond. Here,
convertibility means that the bond may, at any time of the holder’s choosing, be
converted to a multiple of the specified asset. Thus a convertible bond
with price
returns an amountat 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 function fkforce_cbond.
The terminal conditions are provided with the user function fkinit_cbond.
#include <imsl.h>
#include <stdio.h>
#include <math.h>
#define max(A,B) ((A) >= (B) ? (A) : (B))
#define NXGRID 61
#define NTGRID 2
#define NV 13
void fkcfiv_cbond(float, float, int *, float *, void *);
void fkbcp_cbond(int, float, int *, float[], void *);
void fkforce_cbond(int, int, int, float[], float, float,
float[], float[], float[], float[], float[], void *);
void fkinit_cbond(int, int, float[], float[], float,
float[], float[], float[], float[], void *);
int main()
{
int i;
/* Compute value of a Convertible Bond */
/* The face value */
float KS = 1.0e0;
/* The sigma or volatility value */
float sigma = 0.25e0;
/* Time values for the options */
float time[] = { 0.5, 1.0};
/* Values of the underlying where evaluation are made */
float xs[NV];
/* Value of the interest rate, continuous dividend and factor */
float r = 0.1, dividend=0.02, factor = 1.125;
/* Values of the min and max underlying values modeled */
float x_min = 0.0, x_max = 4.0;
/* Define parameters for the integration step. */
int nint = NXGRID-1, n=3*NXGRID;
float xgrid[NXGRID];
float y[(NTGRID+1)*3*NXGRID], yprime[(NTGRID+1)*3*NXGRID];
float f[(NTGRID+1)*NV], dx;
/* Array for user-defined data */
float usr_data[8];
float atol;
/* Number of left/right boundary conditions. */
int nlbcd = 3, nrbcd = 3;
/*
* Define an equally-spaced grid of points for the
* underlying price
*/
dx = (x_max-x_min)/((float) nint);
xgrid[0] = x_min;
xgrid[NXGRID-1] = x_max;
for (i=2; i<=NXGRID-1; i++) xgrid[i-1] = xgrid[i-2] + dx;
for (i=1; i<=NV; i++) xs[i-1] = (i-1)*0.25;
/* Pass the data for evaluation */
usr_data[0] = KS;
usr_data[1] = x_max;
usr_data[2] = sigma;
usr_data[3] = r;
usr_data[4] = dividend;
usr_data[5] = factor;
/* Use a pure absolute error tolerance for the integration */
atol = 1.0e-3;
usr_data[6] = atol;
/* Compute value of convertible bond */
imsl_f_feynman_kac(NXGRID, NTGRID, nlbcd, nrbcd, xgrid, time,
NULL, NULL, y, yprime,
IMSL_FCN_FKCFIV_W_DATA, fkcfiv_cbond, usr_data,
IMSL_FCN_FKBCP_W_DATA, fkbcp_cbond, usr_data,
IMSL_FCN_INIT_W_DATA, fkinit_cbond, usr_data,
IMSL_FCN_FORCE_W_DATA, fkforce_cbond, usr_data,
IMSL_ATOL_RTOL_SCALARS, 1.0e-3, 0.0e0,
0);
/*
* Evaluate and display solutions at vector of points XS(:),
* at each time value prior to expiration.
*/
for (i=0; i<=NTGRID; i++)
imsl_f_feynman_kac_evaluate (NV, NXGRID, xgrid, xs, &y[i*n],
IMSL_RETURN_USER, &f[i*NV], 0);
printf("%2sConvertible Bond Value, 0+, 6 and 12 Months Prior "
"to Expiry\n", " ");
printf("%5sNumber of equally spaced spline knots:%4d\n", " ",
NXGRID);
printf("%5sNumber of unknowns:%4d\n", " ",n);
printf("%5sStrike=%5.2f, Sigma=%5.2f\n", " ", KS, sigma);
printf("%5sInterest Rate=%5.2f, Dividend=%5.2f, Factor=%6.3f\n\n",
" ", r, dividend, factor);
printf("%15s%18s\n", "Underlying", "Bond Value");
for (i=0; i<NV; i++)
printf("%7s%8.4f%8.4f%8.4f%8.4f\n",
" ", xs[i], f[i], f[i+NV], f[i+2*NV]);
}
/*
* These functions define the coefficients, payoff, boundary conditions
* and forcing term.
*/
void fkcfiv_cbond(float x, float tx, int *iflag, float *value,
void *data_ptr)
{
float sigma, strike_price, interest_rate;
float dividend, factor, zero=0.0;
float *data = NULL;
data = (float *) data_ptr;
strike_price = data[0];
sigma = data[2];
interest_rate = data[3];
dividend = data[4];
factor = data[5];
switch(*iflag)
{
case 0:
/* The payoff function - */
*value = max(factor * x, strike_price);
break;
case -1:
/* The coefficient derivative d sigma/ dx */
*value = sigma;
break;
case 1:
/* The coefficient sigma(x) */
*value = sigma*x;
break;
case 2:
/* The coefficient mu(x) */
*value = (interest_rate - dividend) * x;
break;
case 3:
/* The coefficient kappa(x) */
*value = interest_rate;
break;
}
/* Note that there is no time dependence */
*iflag = 0;
}
void fkbcp_cbond(int nbc, float tx, int *iflag, float val[],
void *data_ptr)
{
float interest_rate, strike_price, dp, factor, x_max;
float *data = NULL;
data = (float *) data_ptr;
switch (*iflag)
{
case 1:
strike_price = data[0];
interest_rate = data[3];
dp = strike_price * exp(tx*interest_rate);
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = dp; val[4] = 0.0; val[5] = 1.0;
val[6] = 0.0; val[7] = 0.0; val[8] = 0.0;
val[9] = 0.0; val[10] = 1.0; val[11] = 0.0;
break;
case 2:
x_max = data[1];
factor = data[5];
val[0] = 1.0; val[1] = 0.0; val[2] = 0.0;
val[3] = factor*x_max; val[4] = 0.0; val[5] = 1.0;
val[6] = 0.0; val[7] = factor; val[8] = 0.0;
val[9] = 0.0; val[10] = 1.0; val[11] = 0.0;
/* Note no time dependence */
*iflag = 0;
break;
}
return;
}
void fkforce_cbond(int interval, int ndeg, int nxgrid, float
y[],
float time, float width,
float xlocal[], float
qw[],
float u[], float phi[], float dphi[], void *data_ptr)
{
int i, j, k, l;
const int local=6;
float yl[6], bf[6];
float value, strike_price, interest_rate, zero=0.e0;
float one=1.0e0, rt, mu, factor;
float *data = NULL;
data = (float *) data_ptr;
for (i=0; i<local; i++)
{
yl[i] = y[3*interval-3+i];
phi[i] = zero;
}
for (i=0; i<local*local; i++)
dphi[i] = zero;
value = data[6];
strike_price = data[0];
interest_rate = data[3];
factor = data[5];
mu = 2.0;
/*
* This is the local definition of the forcing term -
* It "forces" the constraint f >= factor*x.
*/
for (j=1; j<=local; j++)
for (l=1; l<=ndeg; l++)
{
bf[0] = u[(l-1)];
bf[1] = u[(l-1)+ndeg];
bf[2] = u[(l-1)+2*ndeg];
bf[3] = u[(l-1)+6*ndeg];
bf[4] = u[(l-1)+7*ndeg];
bf[5] = u[(l-1)+8*ndeg];
rt = 0.0;
for (k=0; k<local; k++)
rt += yl[k]*bf[k];
rt = value/(rt + value - factor * xlocal[l-1]);
phi[j-1] += qw[l-1] * bf[j-1] * pow(rt,mu);
}
for (i=0; i<local; i++)
phi[i] = -phi[i]*width*factor*strike_price;
/*
* This is the local derivative matrix for the forcing term
*/
for (j=1; j<=local; j++)
for (i=1; i<=local; i++)
for (l=1; l<=ndeg; l++)
{
bf[0] = u[(l-1)];
bf[1] = u[(l-1)+ndeg];
bf[2] = u[(l-1)+2*ndeg];
bf[3] = u[(l-1)+6*ndeg];
bf[4] = u[(l-1)+7*ndeg];
bf[5] = u[(l-1)+8*ndeg];
rt = 0.0;
for (k=0; k<local; k++) rt += yl[k]*bf[k];
rt = one/(rt + value - factor * xlocal[l-1]);
dphi[i-1+(j-1)*local] += qw[l-1] * bf[i-1] *
bf[j-1] * pow(value*rt, mu) * rt;
}
for (i=0; i<local*local; i++)
dphi[i] = -mu * dphi[i] * width * factor * strike_price;
return;
}
void fkinit_cbond(int nxgrid, int ntgrid, float xgrid[], float
tgrid[],
float time, float yprime[], float y[], float
atol[],
float rtol[], void *data_ptr)
{
int i;
float *data = NULL;
data = (float *) data_ptr;
if (time == 0.0)
{
/* Set initial data precisely. */
for (i=1; i<=nxgrid; i++)
{
if (xgrid[i-1] * data[5] < data[0])
{
y[3*i-3] = data[0];
y[3*i-2] = 0.0;
y[3*i-1] = 0.0;
}
else
{
y[3*i-3] = xgrid[i-1] * data[5];
y[3*i-2] = data[5];
y[3*i-1] = 0.0;
}
}
}
return;
}
Output
Convertible Bond Value, 0+, 6 and 12 Months Prior to Expiry
Number of equally spaced spline knots: 61
Number of unknowns: 183
Strike= 1.00, Sigma= 0.25
Interest Rate= 0.10, Dividend= 0.02, Factor= 1.125
Underlying Bond Value
0.0000 1.0000 0.9512 0.9048
0.2500 1.0000 0.9512 0.9049
0.5000 1.0000 0.9513 0.9065
0.7500 1.0000 0.9737 0.9605
1.0000 1.1250 1.1416 1.1464
1.2500 1.4063 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
