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 ForcingTerm
. 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 InitialData
.
import com.imsl.math.*;
import java.util.*;
public class FeynmanKacEx4 {
public static void main(String args[]) throws Exception {
int i;
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;
boolean[] timeDependence = new boolean[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.setErrorTolerances(1.0e-3,0.0);
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[][] 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++) {
bondSplineValues[i] = convertibleBond.getSplineValue(evaluationPoints,
bondCoefficients[i], 0);
}
System.out.printf("%2sConvertible Bond Value, 0+, 6 and 12 Months " +
"Prior to Expiry%n", " ");
System.out.printf("%5sNumber of equally spaced spline knots:%4d%n",
" ", nxGrid);
System.out.printf("%5sNumber of unknowns:%4d%n", " ", n);
System.out.printf(Locale.ENGLISH, "%5sStrike=%5.2f, Sigma=%5.2f%n",
" ", KS, sigma);
System.out.printf(Locale.ENGLISH, "%5sInterest Rate=%5.2f, " +
"Dividend=%5.2f, Factor=%6.3f%n%n",
" ", r, dividend, factor);
System.out.printf("%15s%18s%n", "Underlying", "Bond Value");
for (i = 0; i < nv; i++) {
System.out.printf(Locale.ENGLISH, "%7s%8.4f%8.4f%8.4f%8.4f%n",
" ", evaluationPoints[i],
bondSplineValues[0][i],
bondSplineValues[1][i],
bondSplineValues[2][i]);
}
}
/*
* These classes define the coefficients, payoff, boundary conditions
* and forcing term.
*/
static class MyCoefficients implements FeynmanKac.PdeCoefficients {
private double sigma, strikePrice, interestRate;
private double dividend, factor, zero = 0.0;
private double value = 0.0;
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;
}
}
static class MyBoundaries implements FeynmanKac.Boundaries {
private 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 value = Math.max(factor * x, strikePrice);
return value;
}
}
static class MyForcingTerm implements FeynmanKac.ForcingTerm {
private final double zero = 0.0, one = 1.0;
private double value, strikePrice, interestRate;
private double rt, mu, factor;
public MyForcingTerm(double[] rData) {
this.value = 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) {
final 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 = value / (rt + value - 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 + value - factor * xlocal[l - 1]);
dphi[j - 1][i - 1] += qw[l - 1] * bf[i - 1] *
bf[j - 1] * Math.pow(value * 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;
}
}
static class MyInitialData implements FeynmanKac.InitialData {
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= 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.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 Java source.