Example 1: American Option vs. European Option on a Vanilla Put

The value of the American Option on a Vanilla Put can be no smaller than its European counterpart. That is due to the American Option providing the opportunity to exercise at any time prior to expiration. This example compares this difference - or premium value of the American Option - at two time values using the Black-Scholes model. The example is based on Wilmott et al. (1996, p. 176), and uses the non-linear forcing or weighting term described in Hanson, R. (2008), Integrating Feynman-Kac Equations Using Hermite Quintic Finite Elements , for evaluating the price of the American Option. The coefficients, payoff, boundary conditions and forcing term for American and European options are defined through interfaces PdeCoefficients , Boundaries and ForcingTerm , respectively. One breakpoint is set exactly at the strike price. The sets of parameters in the computation are:

  1. Strike price K = {10.0}
  2. Volatility  \sigma = {0.4}
  3. Times until expiration = {1/ 4, 1/ 2}
  4. Interest rate r = 0.1
  5. x_{\min}=0.0,\, x_{\max}=30.0
  6. nx = 61, \, n=3 \times nx = 183
The payoff function is the "vanilla option", p(x) = \max(K-x,0) .
import com.imsl.math.*;
import java.util.*;

public class FeynmanKacEx1 {

   public static void main(String args[]) throws Exception {
      // Compute American Option Premium for Vanilla Put
      // The strike price
      double KS = 10.0;
      // The sigma value
      double sigma = 0.4;
      // Time values for the options
      int nt = 2;
      double[] tGrid = {0.25, 0.5};
      // Values of the underlying where evaluations are made
      double[] evaluationPoints = {0.0, 2.0, 4.0, 6.0, 8.0, 10.0,
         12.0, 14.0, 16.0};
      // Value of the interest rate
      double r = 0.1;
      // Values of the min and max underlying values modeled
      double xMin = 0.0, xMax = 30.0;
      // Define parameters for the integration step.
      int nxGrid = 61;
      int nv = 9;
      int nint = nxGrid - 1, n = 3 * nxGrid;
      double[] xGrid = new double[nxGrid];
      double dx;
      int nlbcd = 2, nrbcd = 3;
      // Define an equally-spaced grid of points for the
      // underlying price
      dx = (xMax - xMin) / ((double) nint);
      xGrid[0] = xMin;
      xGrid[nxGrid - 1] = xMax;
      for (int i = 2; i <= nxGrid - 1; i++) {
         xGrid[i - 1] = xGrid[i - 2] + dx;
      }

      FeynmanKac.Boundaries boundaries = new FeynmanKac.Boundaries() {

         public void leftBoundaries(double time, double[][] bndCoeffs) {
            bndCoeffs[0][0] = 0.0;
            bndCoeffs[0][1] = 1.0;
            bndCoeffs[0][2] = 0.0;
            bndCoeffs[0][3] = -1.0;
            bndCoeffs[1][0] = 0.0;
            bndCoeffs[1][1] = 0.0;
            bndCoeffs[1][2] = 1.0;
            bndCoeffs[1][3] = 0.0;
         }

         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] = 0.0;
            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;
         }

         public double terminal(double x) {
            double zero = 0.0;
            // Strike price
            double strikePrice = 10.0;
            // The payoff function
            double value = Math.max(strikePrice - x, zero);

            return value;
         }
      };

      FeynmanKac.PdeCoefficients pdeCoefficients =
            new FeynmanKac.PdeCoefficients() {
               // The coefficient sigma(x)

               public double sigma(double x, double t) {
                  double sigma = 0.4;
                  return (sigma * x);
               }

               // The coefficient derivative d(sigma) / dx
               public double sigmaPrime(double x, double t) {
                  double sigma = 0.4;
                  return sigma;
               }

               // The coefficient mu(x)
               public double mu(double x, double t) {
                  double interestRate = 0.1;
                  double dividend = 0.0;
                  return ((interestRate - dividend) * x);
               }

               // The coefficient kappa(x)
               public double kappa(double x, double t) {
                  double interestRate = 0.1;
                  return interestRate;
               }
            };

      FeynmanKac.ForcingTerm forceTerm = new FeynmanKac.ForcingTerm() {

         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;
            final double zero = 0.0,  one = 1.0;
            double[] yl = new double[local];
            double[] bf = new double[local];
            double value, strikePrice, interestRate;
            double rt, mu;
            int nxGrid = y.length / 3;
            int ndeg = xlocal.length;

            for (int i = 0; i < local; i++) {
               yl[i] = y[3 * interval - 3 + i];
               phi[i] = zero;
            }
            strikePrice = 10.0;
            interestRate = 0.1;
            value = 1.0e-5;
            mu = 2.0;
            // This is the local definition of the forcing term
            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 - (strikePrice - 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 * interestRate * strikePrice;
            }
            // This is the local derivative matrix for the forcing term
            for (int i = 0; i < local; i++) {
               for (int j = 0; j < local; j++) {
                  dphi[i][j] = zero;
               }
            }
            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 - (strikePrice - xlocal[l - 1]));
                     dphi[j - 1][i - 1] += qw[l - 1] * bf[i - 1] * bf[j - 1] *
                           Math.pow(rt, mu + 1.0);
                  }
               }
            }
            for (int i = 0; i < local; i++) {
               for (int j = 0; j < local; j++) {
                  dphi[i][j] = mu * dphi[i][j] * width * Math.pow(value, mu) *
                        interestRate * strikePrice;
               }
            }
            return;
         }
      };

      FeynmanKac european = new FeynmanKac(pdeCoefficients);
      FeynmanKac american = new FeynmanKac(pdeCoefficients);

      american.setForcingTerm(forceTerm);

      european.computeCoefficients(nlbcd, nrbcd, boundaries, xGrid, tGrid);
      american.computeCoefficients(nlbcd, nrbcd, boundaries, xGrid, tGrid);

      // Evaluate solutions at vector of points evaluationPoints, at each
      // time value prior to expiration.

      double[][] europeanCoefficients = european.getSplineCoefficients();
      double[][] americanCoefficients = american.getSplineCoefficients();

      double[][] splineValuesEuropean = new double[nt][];
      double[][] splineValuesAmerican = new double[nt][];

      for (int i = 0; i < nt; i++) {
         splineValuesEuropean[i] =
               european.getSplineValue(evaluationPoints,
               europeanCoefficients[i + 1], 0);
         splineValuesAmerican[i] =
               american.getSplineValue(evaluationPoints,
               americanCoefficients[i + 1], 0);
      }

      System.out.printf("%nAmerican Option Premium for Vanilla Put, " +
            "3 and 6 Months Prior to Expiry%n");
      System.out.printf("%7sNumber of equally spaced spline knots:" +
            "%4d%n", " ", nxGrid);
      System.out.printf("%7sNumber of unknowns:%4d%n", " ", n);
      System.out.printf(Locale.ENGLISH, "%7sStrike=%6.2f, sigma=%5.2f," +
            " Interest Rate=%5.2f%n%n", " ", KS, sigma, r);
      System.out.printf("%7s%10s%20s%20s%n", " ", "Underlying",
            "European", "American");
      for (int i = 0; i < nv; i++) {
         System.out.printf(Locale.ENGLISH, "%7s%10.4f%10.4f%10.4f" +
               "%10.4f%10.4f%n", " ",
               evaluationPoints[i],
               splineValuesEuropean[0][i],
               splineValuesEuropean[1][i],
               splineValuesAmerican[0][i],
               splineValuesAmerican[1][i]);
      }
   }
}

Output


American Option Premium for Vanilla Put, 3 and 6 Months Prior to Expiry
       Number of equally spaced spline knots:  61
       Number of unknowns: 183
       Strike= 10.00, sigma= 0.40, Interest Rate= 0.10

       Underlying            European            American
           0.0000    9.7531    9.5123   10.0000   10.0000
           2.0000    7.7531    7.5123    8.0000    8.0000
           4.0000    5.7531    5.5128    6.0000    6.0000
           6.0000    3.7569    3.5583    4.0000    4.0000
           8.0000    1.9024    1.9181    2.0202    2.0954
          10.0000    0.6694    0.8703    0.6923    0.9219
          12.0000    0.1675    0.3477    0.1712    0.3625
          14.0000    0.0326    0.1279    0.0332    0.1321
          16.0000    0.0054    0.0448    0.0055    0.0461
Link to Java source.