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 IPdeCoefficients , IBoundaries and IForcingTerm , 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. nxGrid = 61, \, n=3 \times nxGrid = 183
The payoff function is the "vanilla option", p(x) = \max(K-x,0) .
using System;
using Imsl.Math;

public class FeynmanKacEx1 : FeynmanKac.IPdeCoefficients, 
   FeynmanKac.IBoundaries
{
   // 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;
   }

   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 val = Math.Max(strikePrice - x, zero);

      return val;
   }

   class MyForcingTerm : FeynmanKac.IForcingTerm
   {
      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;
         const double zero = 0.0;
         const double one = 1.0;
         double[] yl = new double[local];
         double[] bf = new double[local];
         double val, 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;
         val = 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 = val / (rt + val - (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 + val - (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(val, mu)
                             * interestRate * strikePrice;
            }
         }
         return;
      }
   }

   public static void Main(String[] args)
   {
      // 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;
      }

      FeynmanKacEx1 ex1 = new FeynmanKacEx1();

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

      MyForcingTerm forceTerm = new MyForcingTerm();

      american.SetForcingTerm(forceTerm);

      european.ComputeCoefficients(nlbcd, nrbcd, ex1, xGrid, tGrid);
      american.ComputeCoefficients(nlbcd, nrbcd, ex1, 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[] europeanTimeCoeffs = new double[n];
      double[] americanTimeCoeffs = new double[n];

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

      for (int i = 0; i < nt; i++)
      {
         // Exctract spline coefficients for individual times
         for (int j = 0; j < n; j++)
         {
            europeanTimeCoeffs[j] = europeanCoefficients[i + 1, j];
            americanTimeCoeffs[j] = americanCoefficients[i + 1, j];
         }

         splineValuesEuropean[i] = 
            european.GetSplineValue(evaluationPoints, europeanTimeCoeffs, 0);
         splineValuesAmerican[i] = 
            american.GetSplineValue(evaluationPoints, americanTimeCoeffs, 0);
      }

      Console.Out.WriteLine();
      Console.Out.WriteLine("American Option Premium for Vanilla Put, " +
                            "3 and 6 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,6:f2}, sigma={1,5:f2}, " +
                            "Interest Rate={2,5:f2}", KS, sigma, r);
      Console.Out.WriteLine();
      Console.Out.WriteLine("       Underlying            European" +
                            "            American");
      for (int i = 0; i < nv; i++)
      {
         Console.Out.WriteLine("       {0,10:f4}{1,10:f4}{2,10:f4}" +
                               "{3,10:f4}{4,10:f4}",
                                  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 C# source.