This example evaluates the price of a European Option using two payoff strategies: Cash-or-Nothing and Vertical Spread. In the first case the payoff function is
using System;
using Imsl.Math;
public class FeynmanKacEx3
{
public static void Main(String[] args)
{
int i, j;
int nxGrid = 61;
int ntGrid = 2;
int nv = 12;
// The strike price
double strikePrice = 10.0;
// The spread value
double spreadValue = 15.0;
// The Bet for the Cash-or-Nothing Call
double bet = 2.0;
// The sigma value
double sigma = 0.4;
// Time values for the options
double[] tGrid = { 0.25, 0.5 };
// Values of the underlying where evaluations are made
double[] evaluationPoints = new double[nv];
// Value of the interest rate and continuous dividend
double r = 0.1, dividend = 0.0;
// Values of the min and max underlying values modeled
double xMin = 0.0, xMax = 30.0;
// Define parameters for the integration step.
int nint = nxGrid - 1, n = 3 * nxGrid;
double[] xGrid = new double[nxGrid];
double dx;
// Number of left/right boundary conditions.
int nlbcd = 3, nrbcd = 3;
// Structure for the evaluation routines.
int iData;
double[] rData = new double[7];
bool[] timeDependence = new bool[7];
// Define an equally-spaced grid of points for the
// underlying price
dx = (xMax - xMin) / ((double)(nint));
xGrid[0] = xMin;
xGrid[nxGrid - 1] = xMax;
for (i = 2; i <= nxGrid - 1; i++)
{
xGrid[i - 1] = xGrid[i - 2] + dx;
}
for (i = 1; i <= nv; i++)
{
evaluationPoints[i - 1] = 2.0 + (i - 1) * 2.0;
}
rData[0] = strikePrice;
rData[1] = bet;
rData[2] = spreadValue;
rData[3] = xMax;
rData[4] = sigma;
rData[5] = r;
rData[6] = dividend;
// Flag the difference in payoff functions
// 1 for the Bet, 2 for the Vertical Spread
// In both cases, no time dependencies for
// the coefficients and the left boundary
// conditions
timeDependence[5] = true;
iData = 1;
FeynmanKac betOption = new FeynmanKac(new MyCoefficients(rData));
betOption.SetTimeDependence(timeDependence);
betOption.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(rData, iData), xGrid, tGrid);
double[,] betOptionCoefficients = betOption.GetSplineCoefficients();
double[][] splineValuesBetOption = new double[ntGrid][];
double[] betOptionTimeCoeffs = new double[n];
iData = 2;
FeynmanKac spreadOption = new FeynmanKac(new MyCoefficients(rData));
spreadOption.SetTimeDependence(timeDependence);
spreadOption.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(rData, iData), xGrid, tGrid);
double[,] spreadOptionCoefficients =
spreadOption.GetSplineCoefficients();
double[][] splineValuesSpreadOption = new double[ntGrid][];
double[] spreadOptionTimeCoeffs = new double[n];
// Evaluate solutions at vector evaluationPoints, at each time value
// prior to expiration.
for (i = 0; i < ntGrid; i++)
{
for (j = 0; j < n; j++)
{
betOptionTimeCoeffs[j] = betOptionCoefficients[i + 1, j];
spreadOptionTimeCoeffs[j] = spreadOptionCoefficients[i + 1, j];
}
splineValuesBetOption[i] = betOption.GetSplineValue(
evaluationPoints, betOptionTimeCoeffs, 0);
splineValuesSpreadOption[i] = spreadOption.GetSplineValue(
evaluationPoints, spreadOptionTimeCoeffs, 0);
}
Console.Out.WriteLine(" European Option Value for A Bet");
Console.Out.WriteLine(
" and a Vertical Spread, 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,5:f2}, Sigma={1,5:f2}, Interest Rate={2,5:f2}",
strikePrice, sigma, r);
Console.Out.WriteLine(" Bet={0,5:f2}, Spread Value={1,5:f2}",
bet, spreadValue);
Console.Out.WriteLine();
Console.Out.WriteLine(
" Underlying A Bet Vertical Spread");
for (i = 0; i < nv; i++)
{
Console.Out.WriteLine(
" {0,9:f4}{1,9:f4}{2,9:f4}{3,9:f4}{4,9:f4}",
evaluationPoints[i], splineValuesBetOption[0][i],
splineValuesBetOption[1][i], splineValuesSpreadOption[0][i],
splineValuesSpreadOption[1][i]);
}
}
// These routines define the coefficients, payoff, boundary conditions
// and forcing term for American and European Options.
class MyCoefficients : FeynmanKac.IPdeCoefficients
{
const double zero = 0.0;
double sigma, strikePrice, interestRate;
double spread, bet, dividend;
public MyCoefficients(double[] rData)
{
this.strikePrice = rData[0];
this.bet = rData[1];
this.spread = rData[2];
this.sigma = rData[4];
this.interestRate = rData[5];
this.dividend = rData[6];
}
// The coefficient sigma(x)
public double Sigma(double x, double t)
{
return (sigma * x);
}
// The coefficient derivative d(sigma) / dx
public double SigmaPrime(double x, double t)
{
return sigma;
}
// The coefficient mu(x)
public double Mu(double x, double t)
{
return ((interestRate - dividend) * x);
}
// The coefficient kappa(x)
public double Kappa(double x, double t)
{
return interestRate;
}
}
class MyBoundaries : FeynmanKac.IBoundaries
{
double strikePrice, spread, bet, interestRate, df;
int dataInt;
public MyBoundaries(double[] rData, int iData)
{
this.strikePrice = rData[0];
this.bet = rData[1];
this.spread = rData[2];
this.interestRate = rData[5];
this.dataInt = iData;
}
public void LeftBoundaries(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;
return;
}
public void RightBoundaries(double time, double[,] bndCoeffs)
{
// This is the discount factor using the risk-free
// interest rate
df = Math.Exp(interestRate * time);
// Use flag passed to decide on boundary condition
switch (dataInt)
{
case 1:
bndCoeffs[0, 0] = 1.0;
bndCoeffs[0, 1] = 0.0;
bndCoeffs[0, 2] = 0.0;
bndCoeffs[0, 3] = bet * df;
break;
case 2:
bndCoeffs[0, 0] = 1.0;
bndCoeffs[0, 1] = 0.0;
bndCoeffs[0, 2] = 0.0;
bndCoeffs[0, 3] = (spread - strikePrice) * df;
break;
}
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;
return;
}
public double Terminal(double x)
{
const double zero = 0.0;
double val = 0.0;
switch (dataInt)
{
// The payoff function - Use flag passed to decide which
case 1:
// After reaching the strike price the payoff jumps
// from zero to the bet value.
val = zero;
if (x > strikePrice)
{
val = bet;
}
break;
case 2:
/* Function is zero up to strike price.
Then linear between strike price and spread.
Then has constant value Spread-Strike Price after
the value Spread. */
val = Math.Max(x - strikePrice, zero) -
Math.Max(x - spread, zero);
break;
}
return val;
}
}
}
European Option Value for A Bet
and a Vertical Spread, 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
Bet= 2.00, Spread Value=15.00
Underlying A Bet Vertical Spread
2.0000 0.0000 0.0000 0.0000 0.0000
4.0000 0.0000 0.0013 0.0000 0.0005
6.0000 0.0112 0.0729 0.0038 0.0452
8.0000 0.2686 0.4291 0.1486 0.3833
10.0000 0.9948 0.9781 0.8898 1.1907
12.0000 1.6103 1.4301 2.1904 2.2267
14.0000 1.8650 1.6926 3.4267 3.1567
16.0000 1.9335 1.8171 4.2274 3.8282
18.0000 1.9477 1.8696 4.6261 4.2499
20.0000 1.9501 1.8902 4.7903 4.4918
22.0000 1.9505 1.8979 4.8493 4.6222
24.0000 1.9506 1.9008 4.8685 4.6901
Link to C# source.