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.0461
Link to C# source.