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:
.
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]);
}
}
}
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.