This example evaluates the price of a convertible bond. Here, convertibility means that the bond may, at any time of the holder's choosing, be converted to a multiple of the specified asset. Thus a convertible bond with price x
returns an amount K
at time T
unless
the owner has converted the bond to 
units of the asset at some time prior to T
. This definition, the differential equation and boundary conditions are given in Chapter 18 of Wilmott et al. (1996). Using a constant interest rate and volatility factor, the parameters and boundary conditions are:
, conversion factor 
. The free boundary constraint 
is achieved by use of a non-linear forcing term in interface IForcingTerm
. The coefficient values of the Hermite quintic spline representing the approximate solution of the differential equation at the initial time point are provided with the interface IInitialData
.
using System;
using Imsl.Math;
public class FeynmanKacEx4
{
public static void Main(String[] args)
{
int i, j;
int nxGrid = 61;
int ntGrid = 2;
int nv = 13;
// Compute value of a Convertible Bond
// The face value
double KS = 1.0;
// The sigma or volatility value
double sigma = 0.25;
// Time values for the options
double[] tGrid = { 0.5, 1.0 };
// Values of the underlying where evaluation are made
double[] evaluationPoints = new double[nv];
// Value of the interest rate, continuous dividend and factor
double r = 0.1, dividend = 0.02, factor = 1.125;
// Values of the min and max underlying values modeled
double xMin = 0.0, xMax = 4.0;
// Define parameters for the integration step.
int nint = nxGrid - 1, n = 3 * nxGrid;
double[] xGrid = new double[nxGrid];
// Array for user-defined data
double[] rData = new double[8];
double dx, atol;
// Number of left/right boundary conditions.
int nlbcd = 3, nrbcd = 3;
bool[] timeDependence = new bool[7];
/*
* Define an equally-spaced grid of points for the
* underlying price
*/
dx = (xMax - xMin) / ((double)nint);
xGrid[0] = xMin;
xGrid[nxGrid - 1] = xMax;
for (i = 2; i <= nxGrid - 1; i++)
{
xGrid[i - 1] = xGrid[i - 2] + dx;
}
for (i = 1; i <= nv; i++)
{
evaluationPoints[i - 1] = (i - 1) * 0.25;
}
// Pass the data for evaluation
rData[0] = KS;
rData[1] = xMax;
rData[2] = sigma;
rData[3] = r;
rData[4] = dividend;
rData[5] = factor;
// Use a pure absolute error tolerance for the integration
atol = 1.0e-3;
rData[6] = atol;
// Compute value of convertible bond
FeynmanKac convertibleBond = new FeynmanKac(new MyCoefficients(rData));
MyForcingTerm forceTerm = new MyForcingTerm(rData);
MyInitialData initialData = new MyInitialData(rData);
convertibleBond.SetForcingTerm(forceTerm);
convertibleBond.SetInitialData(initialData);
convertibleBond.SetAbsoluteErrorTolerances(1.0e-3);
convertibleBond.SetRelativeErrorTolerances(0.0);
// Only the left boundary conditions are time dependent
timeDependence[4] = true;
convertibleBond.SetTimeDependence(timeDependence);
convertibleBond.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(rData), xGrid, tGrid);
double[,] bondCoefficients = convertibleBond.GetSplineCoefficients();
double[] bondTimeCoeffs = new double[n];
double[][] bondSplineValues = new double[ntGrid + 1][];
/*
* Evaluate and display solutions at vector of points XS(:),
* at each time value prior to expiration.
*/
for (i = 0; i <= ntGrid; i++)
{
for (j = 0; j < n; j++)
bondTimeCoeffs[j] = bondCoefficients[i, j];
bondSplineValues[i] = convertibleBond.GetSplineValue(
evaluationPoints, bondTimeCoeffs, 0);
}
Console.Out.WriteLine(" Convertible Bond Value, 0+, 6 and 12 Months "
+ "Prior to Expiry", " ");
Console.Out.WriteLine(
" Number of equally spaced spline knots:{0,4:d}", nxGrid);
Console.Out.WriteLine(" Number of unknowns:{0,4:d}", n);
Console.Out.WriteLine(" Strike={0,5:f2}, Sigma={1,5:f2}",
KS, sigma);
Console.Out.WriteLine(" Interest Rate={0,5:f2}, " +
"Dividend={1,5:f2}, Factor={0,6:f3}", r, dividend, factor);
Console.Out.WriteLine(" Underlying Bond Value");
for (i = 0; i < nv; i++)
{
Console.Out.WriteLine(" {0,8:f4}{1,8:f4}{2,8:f4}{3,8:f4}",
evaluationPoints[i],
bondSplineValues[0][i],
bondSplineValues[1][i],
bondSplineValues[2][i]);
}
}
/*
* These classes define the coefficients, payoff, boundary conditions
* and forcing term.
*/
class MyCoefficients : FeynmanKac.IPdeCoefficients
{
double sigma, strikePrice, interestRate;
double dividend, factor;
public MyCoefficients(double[] rData)
{
this.strikePrice = rData[0];
this.sigma = rData[2];
this.interestRate = rData[3];
this.dividend = rData[4];
this.factor = rData[5];
}
// The coefficient sigma(x)
public double Sigma(double x, double t)
{
return (sigma * x);
}
// The coefficient derivative d(sigma) / dx
public double SigmaPrime(double x, double t)
{
return sigma;
}
// The coefficient mu(x)
public double Mu(double x, double t)
{
return ((interestRate - dividend) * x);
}
// The coefficient kappa(x)
public double Kappa(double x, double t)
{
return interestRate;
}
}
class MyBoundaries : FeynmanKac.IBoundaries
{
double interestRate, strikePrice, dp, factor, xMax;
public MyBoundaries(double[] rData)
{
this.strikePrice = rData[0];
this.xMax = rData[1];
this.interestRate = rData[3];
this.factor = rData[5];
}
public void LeftBoundaries(double time, double[,] bndCoeffs)
{
dp = strikePrice * Math.Exp(time * interestRate);
bndCoeffs[0, 0] = 1.0;
bndCoeffs[0, 1] = 0.0;
bndCoeffs[0, 2] = 0.0;
bndCoeffs[0, 3] = dp;
bndCoeffs[1, 0] = 0.0;
bndCoeffs[1, 1] = 1.0;
bndCoeffs[1, 2] = 0.0;
bndCoeffs[1, 3] = 0.0;
bndCoeffs[2, 0] = 0.0;
bndCoeffs[2, 1] = 0.0;
bndCoeffs[2, 2] = 1.0;
bndCoeffs[2, 3] = 0.0;
return;
}
public void RightBoundaries(double time, double[,] bndCoeffs)
{
bndCoeffs[0, 0] = 1.0;
bndCoeffs[0, 1] = 0.0;
bndCoeffs[0, 2] = 0.0;
bndCoeffs[0, 3] = factor * xMax;
bndCoeffs[1, 0] = 0.0;
bndCoeffs[1, 1] = 1.0;
bndCoeffs[1, 2] = 0.0;
bndCoeffs[1, 3] = factor;
bndCoeffs[2, 0] = 0.0;
bndCoeffs[2, 1] = 0.0;
bndCoeffs[2, 2] = 1.0;
bndCoeffs[2, 3] = 0.0;
return;
}
public double Terminal(double x)
{
// The payoff function
double val = Math.Max(factor * x, strikePrice);
return val;
}
}
class MyForcingTerm : FeynmanKac.IForcingTerm
{
const double zero = 0.0;
const double one = 1.0;
double val, strikePrice, interestRate;
double rt, mu, factor;
public MyForcingTerm(double[] rData)
{
this.val = rData[6];
this.strikePrice = rData[0];
this.interestRate = rData[3];
this.factor = rData[5];
}
public void Force(int interval, double[] y, double time, double width,
double[] xlocal, double[] qw, double[,] u,
double[] phi, double[,] dphi)
{
const int local = 6;
int ndeg = xlocal.Length;
double[] yl = new double[6];
double[] bf = new double[6];
for (int i = 0; i < local; i++)
{
yl[i] = y[3 * interval - 3 + i];
phi[i] = zero;
}
for (int i = 0; i < local; i++)
{
for (int j = 0; j < local; j++)
{
dphi[i, j] = zero;
}
}
mu = 2.0;
/*
* This is the local definition of the forcing term -
* It "forces" the constraint f >= factor*x.
*/
for (int j = 1; j <= local; j++)
{
for (int l = 1; l <= ndeg; l++)
{
bf[0] = u[0, l - 1];
bf[1] = u[1, l - 1];
bf[2] = u[2, l - 1];
bf[3] = u[6, l - 1];
bf[4] = u[7, l - 1];
bf[5] = u[8, l - 1];
rt = 0.0;
for (int k = 0; k < local; k++)
{
rt += yl[k] * bf[k];
}
rt = val / (rt + val - factor * xlocal[l - 1]);
phi[j - 1] += qw[l - 1] * bf[j - 1] * Math.Pow(rt, mu);
}
}
for (int i = 0; i < local; i++)
{
phi[i] = (-phi[i]) * width * factor * strikePrice;
}
/*
* This is the local derivative matrix for the forcing term
*/
for (int j = 1; j <= local; j++)
{
for (int i = 1; i <= local; i++)
{
for (int l = 1; l <= ndeg; l++)
{
bf[0] = u[0, l - 1];
bf[1] = u[1, l - 1];
bf[2] = u[2, l - 1];
bf[3] = u[6, l - 1];
bf[4] = u[7, l - 1];
bf[5] = u[8, l - 1];
rt = 0.0;
for (int k = 0; k < local; k++)
{
rt += yl[k] * bf[k];
}
rt = one / (rt + val - factor * xlocal[l - 1]);
dphi[j - 1, i - 1] += qw[l - 1] * bf[i - 1] * bf[j - 1] *
Math.Pow(val * rt, mu) * rt;
}
}
}
for (int i = 0; i < local; i++)
{
for (int j = 0; j < local; j++)
{
dphi[i, j] = (-mu) * dphi[i, j] * width * factor *
strikePrice;
}
}
return;
}
}
class MyInitialData : FeynmanKac.IInitialData
{
private double[] data;
public MyInitialData(double[] rData)
{
data = new double[rData.Length];
for (int i = 0; i < rData.Length; i++)
{
data[i] = rData[i];
}
}
public void Init(double[] xGrid, double[] tGrid, double tp,
double[] yprime, double[] y, double[] atol, double[] rtol)
{
int nxGrid = xGrid.Length;
if (tp == 0.0)
// Set initial data precisely.
{
for (int 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;
}
}
}
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= 0.100
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.5313 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
Link to C# source.