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:
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.3781Link to C# source.