Example 4: Convertible bonds

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 \nu x,\, \nu \ge 1 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:

  1. Bond face value K = {1} , conversion factor  \nu = 1.125
  2. Volatility  \sigma = {0.25}
  3. Times until expiration = {1/2, 1}
  4. Interest rate r = 0.1, dividend D = 0.02
  5. x_{\min} = 0,\, x_{\max} = 4
  6. nxGrid = 61,\, n = 3 \times nxGrid = 183
  7. Boundary conditions f(0,t) = K \exp(r-(T-t)),\, f(x_{\max},t)=\nu x_{\max}
  8. Terminal data f(x,T)=\max(K,\nu x)
  9. Constraint for bond holder f(x,t) \ge \nu x
Note that the error tolerance is set to a pure absolute error of value 10^{-3}. The free boundary constraint f(x,t) \ge \nu x 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;
      }
   }
}

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