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. nx = 61,\, n = 3 \times nx = 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 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;
      }
   }
}

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