feynman

CNLMath : Differential Equations : feynman_kac

feynman_kac

Synopsis

Required Arguments

Synopsis with Optional Arguments

more...

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 typedouble 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 4 * max (nlbcd, nrbcd) containing the values of the boundary condition coefficients in its first 4*nbc locations. The coefficients for xminare 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 (ntgrid + 1) by (3 * nxgrid) 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

The (ntgrid + 1) by (3 *nxgrid) 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[],

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

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, float 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 10-3 in single precision and 10-5 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

The function imsl_f_feynman_kac uses a finite element Galerkin method over the rectangle

to compute the approximate solution. The interval

is decomposed with a grid

On each subinterval the solution is represented by

The values

are time-dependent coefficients associated with each interval. The basis functions

are given for

The Galerkin principle is then applied. Using the provided initial and boundary conditions leads to an Index 1 differential-algebraic equation 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

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;

}

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

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

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

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

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

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

Fatal Errors

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm. User flag = "#".