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:
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]); } } }
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.0461Link to C# source.