. The second method calculates the Greeks using the closed-form BS evaluations which can be found at http://en.wikipedia.org/wiki/The_Greeks.
In this example, the Feynman-Kac (FK) class is used to solve for the Greeks, i.e. various derivatives of FK solutions applicable to the pricing of options and related financial derivatives. In order to illustrate and verify these calculations, the Greeks are calculated by two methods. The first method involves the FK solution to the diffusion model for call options given in Example 2 for the Black-Scholes (BS) case, i.e. . The second method calculates the Greeks using the closed-form BS evaluations which can be found at http://en.wikipedia.org/wiki/The_Greeks.
This example calculates FK and BS solutions 
to the BS problem and the following Greeks:
is the first derivative of the Value
, 
, of a portfolio of derivative security derived from underlying instrument with respect to the underlying instrument's price S
;
;
is the negative first derivative of V
with respect to time t
;
;
;
is the first derivative of V
with respect to risk free rate r
;
measures sensitivity to volatility parameter 
of the underlying S
;
;
;
.Intrinsic Greeks
include derivatives involving only S
and t
, the intrinsic FK arguments. In the above list, Value, Delta, Gamma, Theta, Charm, Color,
and Speed
are all intrinsic Greeks. As is discussed in Hanson, R. (2008) Integrating Feynman-Kac equations Using Hermite Quintic Finite Elements
, the expansion of the FK solution function 
in terms of quintic polynomial functions defined on S
-grid subintervals and subject to continuity constraints in derivatives 0, 1, and 2 across the boundaries of these subintervals allows Value, Delta, Gamma, Theta, Charm,
and Color
to be calculated directly by the FK class methods GetSplineCoefficients
, GetSplineCoefficientsPrime
, and GetSplineValue
.
Non-intrinsic Greeks
are derivatives of V
involving FK parameters other than the intrinsic arguments S
and t
, such as r
and 
. Non-intrinsic Greeks in the above list include Rho, Vega, Volga,
and Vanna
. In order to calculate non-intrinsic Greek (parameter) derivatives or intrinsic Greek S
-derivatives beyond the second (such as Speed
) or t
-derivatives beyond the first, the entire FK solution must be calculated 3 times (for a parabolic fit) or five times (for a quartic fit), at the point where the derivative is to be evaluated and at nearby points located symmetrically on either side.
Using a Taylor series expansion of 
truncated to m
+1 terms (to allow an m
-degree polynomial fit of m
+1 data points):
we are able to derive the following parabolic (3 point) estimation of first and second derivatives 
and 
in terms of the three values: 
, 
, and 
, where 
and 
:
![f^{(1)}(\sigma)\; \equiv \;\frac{\partial f(\sigma)}{\partial \sigma}\;\approx\; f^{[1]}(\sigma, \epsilon)\; \equiv \;\frac{f(\sigma + \epsilon)-f(\sigma - \epsilon)}{2 \epsilon}](eqn/eqn_0218.png){kind=small}
![f^{(2)}(\sigma)\; \equiv \;\frac{\partial^2 f(\sigma)}{\partial \sigma^2}\;\approx\; f^{[2]}(\sigma, \epsilon)\; \equiv \;\frac{f(\sigma + \epsilon)+f(\sigma - \epsilon)-2f(\sigma)}{\epsilon^2}](eqn/eqn_0219.png){kind=small}
Similarly, the quartic (5 point) estimation of 
and 
in terms of 
, 
, 
, 
, and 
is:
![f^{(1)}(\sigma)\;\approx\;\frac{4}{3} f^{[1]}(\sigma, \epsilon) \;-\; \frac{1}{3}f^{[1]}(\sigma, 2\epsilon)](eqn/eqn_0227.png){kind=small}
![f^{(2)}(\sigma)\;\approx\;\frac{4}{3} f^{[2]}(\sigma, \epsilon) \;-\; \frac{1}{3}f^{[2]}(\sigma, 2\epsilon)](eqn/eqn_0228.png){kind=small}
For our example, the quartic estimate does not appear to be significantly better than the parabolic estimate, so we have produced only parabolic estimates by setting variable iquart
to 0. The user may try the example with the quartic estimate simply by setting iquart
to 1.
As is pointed out in Integrating Feynman-Kac equations Using Hermite Quintic Finite Elements
, the quintic polynomial expansion function used by Feynman-Kac only allows for continuous derivatives through the second derivative. While up to fifth derivatives can be calculated from the quintic expansion (indeed class FeynmanKac
method GetSplineValue
will allow the third derivative to be calculated by setting parameter ideriv
to 3
, as is done in this example), the accuracy is compromised by the inherent lack of continuity across grid points (i.e. subinterval boundaries).
The accurate second derivatives in S
returned by FK method GetSplineValue
can be leveraged into a third derivative estimate by calculating three FK second derivative solutions, the first solution for grid and evaluation point sets 
and the second and third solutions for solution grid and evaluation point sets 
and 
, where the solution grid and evaluation point sets are shifted up and down by 
. In this example, 
is set to 
, where 
is the average value of S
over the range of grid values and 
. The third derivative solution can then be obtained using the parabolic estimate:
This procedure is implemented in the current example to calculate the Greek Speed
. (For comparison purposes, Speed
is also calculated directly
by setting the GetSplineValue
input S derivative parameter iSDeriv
to 3
. The output from this direct calculation is called "Speed2
".)
The average and maximum relative errors (defined as the absolute value of the difference between the BS and FK values divided by the BS value) for each of the Greeks is given at the end of the output. (These relative error statistics are given for nine combinations of Strike Price and Volatility, but only one of the nine combinations is actually printed in the output.) Both intrinsic and non-intrinsic Greeks have good accuracy (average relative error is in the range 0.01 -- 0.0001) except for Volga
, which has an average relative error of about 0.05. This is probably a result of the fact that Volga
involves differences of differences, which will more quickly erode accuracy than calculations using only one difference to approximate a derivative. Possible ways to improve upon the 2 to 4 significant digits of accuracy achieved in this example include increasing FK integration accuracy by reducing the initial step size (using property InitialStepsize
), by choosing more closely spaced S
and t
grid points (by adjusting method ComputeCoefficients
input parameter arrays xGrid
and tGrid
), and by adjusting 
so that the central differences used to calculate the derivatives are not too small to compromise accuracy.
using System;
using Imsl.Math;
using Imsl.Stat;
public class FeynmanKacEx5
{
private static double[] strikePrices = { 15.0, 20.0, 25.0 };
// The set of sigma values
private static double[] sigma = { 0.2, 0.3, 0.4 };
// The set of model diffusion powers: alpha = 2.0 <==> Black Scholes
private static double[] alpha = { 2.0, 1.0, 0.0 };
// Time values for the options
private static double[] tGrid = { 1.0 / 12.0, 4.0 / 12.0, 7.0 / 12.0 };
// Values of the min and max underlying values modeled
private static double xMin = 0.0;
private static double xMax = 60.0;
// Value of the interest rate and continuous dividend
private static double r = 0.05;
private static double dividend = 0.0;
// Define parameters for the integration step.
private static int nXGgrid = 121;
private static int nTGrid = 3;
private static int nv = 3;
private static int nint = nXGgrid - 1;
private static int n = 3 * nXGgrid;
private static double[] xGrid = new double[nXGgrid];
private static double dx;
// Time dependency
private static bool[] timeDependence = new bool[7];
// Number of left/right boundary conditions
private static int nlbcd = 3;
private static int nrbcd = 3;
// Values of the underlying price where evaluations are made
private static double[] evaluationPoints = { 19.0, 20.0, 21.0 };
private static double epsfrac = .001; //used to calc derivatives
private static double dx2 = epsfrac * 0.5 * (xMin + xMax);
private static double sqrt2pi = Math.Sqrt(2.0 * Math.PI);
private static String[] greekName = new String[]{
" Value", " Delta", " Gamma", " Theta",
" Charm", " Color", " Vega", " Volga",
" Vanna", " Rho", " Speed", "Speed2"
};
// Time values for the options
private static double[] rex = new double[greekName.Length];
private static double[] reavg = new double[greekName.Length];
private static int[] irect = new int[greekName.Length];
// Compute Constant Elasticity of Variance Model for Vanilla Call
public static void Main(String[] args)
{
// Define equally-spaced grid of points for the underlying price
dx = (xMax - xMin) / ((double)nint);
xGrid[0] = xMin;
xGrid[nXGgrid - 1] = xMax;
for (int i = 2; i <= nXGgrid - 1; i++)
{
xGrid[i - 1] = xGrid[i - 2] + dx;
}
Console.Out.WriteLine(" Constant Elasticity of Variance Model" +
" for Vanilla Call Option");
Console.Out.WriteLine(
" Interest Rate:{0,7:f3} Continuous Dividend:{1,7:f3}",
r, dividend);
Console.Out.WriteLine(
" Minimum and Maximum Prices of Underlying:{0,7:f2}{1,7:f2}",
xMin, xMax);
Console.Out.WriteLine(
" Number of equally spaced spline knots:{0,4:d}",
nXGgrid - 1);
Console.Out.WriteLine(" Number of unknowns:{0,4:d}", n);
Console.Out.WriteLine();
Console.Out.WriteLine(
" Time in Years Prior to Expiration: {0,7:f4}{1,7:f4}{2,7:f4}",
tGrid[0], tGrid[1], tGrid[2]); // tGrid[] = tau == T - t
Console.Out.WriteLine(
" Option valued at Underlying Prices:{0,7:f2}{1,7:f2}{2,7:f2}",
evaluationPoints[0], evaluationPoints[1], evaluationPoints[2]);
Console.Out.WriteLine();
/*
* iquart = 0 : derivatives estimated with 3-point fitted parabola
* iquart = 1 : derivatives estimated with 5-point fitted quartic
* polynomial
*/
int iquart = 0;
if (iquart == 0)
{
Console.Out.WriteLine(
" 3 point (parabolic) estimate of parameter derivatives;");
}
else
{
Console.Out.WriteLine(
" 5 point (quartic) estimate of parameter derivatives");
}
Console.Out.WriteLine(" epsfrac = {0,11:f8}", epsfrac);
//alpha: Black-Scholes
for (int i2 = 1; i2 <= 3; i2++)
/* Loop over volatility */
for (int i3 = 1; i3 <= 3; i3++)
/* Loop over strike price */
calculateGreeks(i2, i3, iquart);
Console.Out.WriteLine();
for (int ig = 0; ig < 12; ig++)
{
reavg[ig] /= irect[ig];
Console.Out.WriteLine();
Console.Out.Write(" Greek: " + greekName[ig] +
"; avg rel err: {0,15:f12}; max rel err: {1,15:f12}",
reavg[ig], rex[ig]);
}
Console.Out.WriteLine();
} // end main
private static void calculateGreeks(
int volatility, int strikePrice, int iquart)
{
int i1 = 1;
int nt = 3;
if ((volatility == 1) && (strikePrice == 1))
{
Console.Out.WriteLine();
Console.Out.WriteLine(" Strike={0,5:f2}, Sigma=" +
"{1,5:f2}, Alpha={2,5:f2}:", strikePrices[strikePrice - 1],
sigma[volatility - 1], alpha[i1 - 1]);
Console.Out.WriteLine();
Console.Out.WriteLine(" years to expiration: " +
" {0,7:f4} {1,7:f4} {2,7:f4}",
tGrid[0], tGrid[1], tGrid[2]);
Console.Out.WriteLine();
}
/* Loop over t derivative index */
for (int iTDeriv = 0; iTDeriv < 2; iTDeriv++)
{
int iSDerivMax = 4 - iTDeriv;
/* Loop over S derivative index */
for (int iSDeriv = 0; iSDeriv < iSDerivMax; iSDeriv++)
{
int iSigDerivMax = 1;
if (iTDeriv == 0)
{
if (iSDeriv == 0)
iSigDerivMax = 3;
if (iSDeriv == 1)
iSigDerivMax = 2;
}
//Loop over sigma deriv index
for (int iSigDeriv = 0; iSigDeriv < iSigDerivMax; iSigDeriv++)
{
int iRDerivMax = 1;
if (iTDeriv == 0 && iSDeriv == 0 && iSigDeriv == 0)
iRDerivMax = 2;
// Loop over r derivative index
for (int iRDeriv = 0; iRDeriv < iRDerivMax; iRDeriv++)
{
int ispeedmin = 0;
int ispeedmax = 1;
if (iTDeriv == 0 && iSDeriv == 2)
ispeedmax = 2;
if (iTDeriv == 0 && iSDeriv == 3)
{
ispeedmin = 2;
ispeedmax = 3;
}
// Loop over speed index
for (int ispeed = ispeedmin; ispeed < ispeedmax; ispeed++)
{
// Pass data through into evaluation routines.
double[] userData = new double[6];
userData[0] = strikePrices[strikePrice - 1];
userData[1] = xMax;
userData[2] = sigma[volatility - 1];
userData[3] = alpha[i1 - 1];
userData[4] = r;
userData[5] = dividend;
double[][] splineValuesOption = new double[nTGrid][];
for (int i = 0; i < nTGrid; i++)
{
splineValuesOption[i] = new double[nv];
}
double[][] splineValuesOptionP = new double[nTGrid][];
double[][] splineValuesOptionM = new double[nTGrid][];
double[][] splineValuesOptionPP = new double[nTGrid][];
double[][] splineValuesOptionMM = new double[nTGrid][];
double[] xGridP = new double[nXGgrid];
double[] xGridM = new double[nXGgrid];
double[] evaluationPointsP = new double[nv];
double[] evaluationPointsM = new double[nv];
double[] xGridPP = new double[nXGgrid];
double[] xGridMM = new double[nXGgrid];
double[] evaluationPointsPP = new double[nv];
double[] evaluationPointsMM = new double[nv];
double xMaxP = xMax;
double xMaxM = xMax;
double xMaxPP = xMax, xMaxMM = xMax;
// Evaluate FK solution at vector evaluationPoints, at
// each time value prior to expiration.
if ((iSigDeriv != 1 || iSDeriv == 1 || iRDeriv != 1) &&
(ispeed == 0 || ispeed == 2))
{
FeynmanKac callOption =
new FeynmanKac(new MyCoefficients(userData));
//Right boundary condition time-dependent
timeDependence[5] = true;
callOption.SetTimeDependence(timeDependence);
callOption.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(userData), xGrid, tGrid);
double[] tmpArray = new double[3 * xGrid.Length];
double[,] optionCoefficients;
if (iTDeriv == 0)
{
optionCoefficients =
callOption.GetSplineCoefficients();
}
else
{
optionCoefficients =
callOption.GetSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < 3 * xGrid.Length; j++)
tmpArray[j] = optionCoefficients[i + 1, j];
// FK option values for tau = tGrid[i]:
splineValuesOption[i] =
callOption.GetSplineValue(evaluationPoints,
tmpArray, iSDeriv);
}
}
if (iSigDeriv > 0 || iRDeriv > 0 || ispeed == 1)
{
Array.Copy(xGrid, 0, xGridP, 0, nXGgrid);
Array.Copy(xGrid, 0, xGridM, 0, nXGgrid);
Array.Copy(evaluationPoints, 0,
evaluationPointsP, 0, nv);
Array.Copy(evaluationPoints, 0,
evaluationPointsM, 0, nv);
if (ispeed == 1)
{
for (int i = 0; i < nXGgrid; i++)
{
xGridP[i] = xGrid[i] + dx2;
xGridM[i] = xGrid[i] - dx2;
}
for (int i = 0; i < nv; i++)
{
evaluationPointsP[i] =
evaluationPoints[i] + dx2;
evaluationPointsM[i] =
evaluationPoints[i] - dx2;
}
xMaxP = xMax + dx2;
xMaxM = xMax - dx2;
}
userData[1] = xMaxP;
// calculate spline values for
// sigmaP = sigma[i2-1]*(1. + epsfrac)
// or rP = r*(1. + epsfrac)
if (iSigDeriv > 0)
userData[2] = sigma[volatility - 1] *
(1.0 + epsfrac);
else if (iRDeriv > 0)
userData[4] = r * (1.0 + epsfrac);
FeynmanKac callOptionP =
new FeynmanKac(new MyCoefficients(userData));
//Right boundary condition time-dependent
timeDependence[5] = true;
callOptionP.SetTimeDependence(timeDependence);
callOptionP.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(userData), xGridP, tGrid);
double[,] optionCoefficientsP;
double[] tmpArray2 = new double[3 * xGrid.Length];
if (iTDeriv == 0)
{
optionCoefficientsP =
callOptionP.GetSplineCoefficients();
}
else
{
optionCoefficientsP =
callOptionP.GetSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < 3 * xGrid.Length; j++)
tmpArray2[j] = optionCoefficientsP[i + 1, j];
// FK option values for tau = tGrid[i]:
splineValuesOptionP[i] =
callOptionP.GetSplineValue(evaluationPointsP,
tmpArray2, iSDeriv);
}
userData[1] = xMaxM;
// calculate spline values for
// sigmaM = sigma[i2-1]*(1. - epsfrac) or
// rM = r*(1. - epsfrac):
if (iSigDeriv > 0)
userData[2] = sigma[volatility - 1] *
(1.0 - epsfrac);
else
userData[4] = r * (1.0 - epsfrac);
FeynmanKac callOptionM =
new FeynmanKac(new MyCoefficients(userData));
//Right boundary condition time-dependent
timeDependence[5] = true;
callOptionM.SetTimeDependence(timeDependence);
callOptionM.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(userData), xGridM, tGrid);
double[,] optionCoefficientsM;
double[] tmpArray3 = new double[3 * xGrid.Length];
if (iTDeriv == 0)
{
optionCoefficientsM =
callOptionM.GetSplineCoefficients();
}
else
{
optionCoefficientsM =
callOptionM.GetSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < 3 * xGrid.Length; j++)
tmpArray3[j] = optionCoefficientsM[i + 1, j];
// FK option values for tau = tGrid[i]:
splineValuesOptionM[i] =
callOptionM.GetSplineValue(evaluationPointsM,
tmpArray3, iSDeriv);
}
if (iquart == 1)
{
Array.Copy(xGrid, 0, xGridPP, 0, nXGgrid);
Array.Copy(xGrid, 0, xGridMM, 0, nXGgrid);
Array.Copy(evaluationPoints, 0,
evaluationPointsPP, 0, nv);
Array.Copy(evaluationPoints, 0,
evaluationPointsMM, 0, nv);
if (ispeed == 1)
{
for (int i = 0; i < nXGgrid; i++)
{
xGridPP[i] = xGrid[i] + 2.0 * dx2;
xGridMM[i] = xGrid[i] - 2.0 * dx2;
}
for (int i = 0; i < nv; i++)
{
evaluationPointsPP[i] =
evaluationPoints[i] + 2.0 * dx2;
evaluationPointsMM[i] =
evaluationPoints[i] - 2.0 * dx2;
}
xMaxPP = xMax + 2.0 * dx2;
xMaxMM = xMax - 2.0 * dx2;
}
userData[1] = xMaxPP;
if (iSigDeriv > 0)
{
// calculate spline values for
// sigmaPP = sigma[i2-1]*(1. + 2.*epsfrac):
userData[2] = sigma[volatility - 1] *
(1.0 + 2.0 * epsfrac);
}
else if (iRDeriv > 0)
{
// calculate spline values for
// rPP = r*(1. + 2.*epsfrac):
userData[4] = r * (1.0 + 2.0 * epsfrac);
}
FeynmanKac callOptionPP =
new FeynmanKac(new MyCoefficients(userData));
//Right boundary condition time-dependent
timeDependence[5] = true;
callOptionPP.SetTimeDependence(timeDependence);
callOptionPP.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(userData), xGridPP, tGrid);
double[] tmpArray4 = new double[3 * xGrid.Length];
double[,] optionCoefficientsPP;
if (iTDeriv == 0)
{
optionCoefficientsPP =
callOptionPP.GetSplineCoefficients();
}
else
{
optionCoefficientsPP =
callOptionPP.GetSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < 3 * xGrid.Length; j++)
tmpArray4[j] =
optionCoefficientsPP[i + 1, j];
// FK option values for tau = tGrid[i]:
splineValuesOptionPP[i] =
callOptionPP.GetSplineValue(
evaluationPointsPP, tmpArray4, iSDeriv);
}
userData[1] = xMaxMM;
// calculate spline values for
// sigmaMM = sigma[i2-1]*(1. - 2.*epsfrac) or
// rMM = r*(1. - 2.*epsfrac)
if (iSigDeriv > 0)
userData[2] = sigma[volatility - 1] *
(1.0 - 2.0 * epsfrac);
else if (iRDeriv > 0)
userData[4] = r * (1.0 - 2.0 * epsfrac);
FeynmanKac callOptionMM =
new FeynmanKac(new MyCoefficients(userData));
//Right boundary condition time-dependent
timeDependence[5] = true;
callOptionMM.SetTimeDependence(timeDependence);
callOptionMM.ComputeCoefficients(nlbcd, nrbcd,
new MyBoundaries(userData), xGridMM, tGrid);
double[] tmpArray5 = new double[3 * xGrid.Length];
double[,] optionCoefficientsMM;
if (iTDeriv == 0)
{
optionCoefficientsMM =
callOptionMM.GetSplineCoefficients();
}
else
{
optionCoefficientsMM =
callOptionMM.GetSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < 3 * xGrid.Length; j++)
tmpArray5[j] =
optionCoefficientsMM[i + 1, j];
// FK option values for tau = tGrid[i]:
splineValuesOptionMM[i] =
callOptionMM.GetSplineValue(
evaluationPointsMM,tmpArray5, iSDeriv);
}
}
}
if (iSigDeriv == 1 || iRDeriv == 1 || ispeed == 1)
{
double eps = 0.0, f11 = 0.0, f12 = 0.0;
if (iSigDeriv == 1)
eps = sigma[volatility - 1] * epsfrac;
if (iRDeriv == 1)
eps = r * epsfrac;
if (ispeed == 1)
eps = dx2;
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < nv; j++)
{
f11 = (splineValuesOptionP[i][j] -
splineValuesOptionM[i][j]) / (2.0 * eps);
if (iquart == 0)
splineValuesOption[i][j] = f11;
else
{
f12 = (splineValuesOptionPP[i][j] -
splineValuesOptionMM[i][j]) / (4.0 * eps);
splineValuesOption[i][j] =
(4.0 * f11 - f12) / 3.0;
}
}
}
}
if (iSigDeriv == 2)
{
double eps = sigma[volatility - 1] * epsfrac;
double f21 = 0.0;
double f22 = 0.0;
for (int i = 0; i < nTGrid; i++)
{
for (int j = 0; j < nv; j++)
{
f21 = (splineValuesOptionP[i][j] +
splineValuesOptionM[i][j] - 2.0 *
splineValuesOption[i][j]) / (eps * eps);
if (iquart == 0)
splineValuesOption[i][j] = f21;
else
{
f22 = (splineValuesOptionPP[i][j] +
splineValuesOptionMM[i][j] - 2.0 *
splineValuesOption[i][j]) /
(4.0 * eps * eps);
splineValuesOption[i][j] =
(4.0 * f21 - f22) / 3.0;
}
}
}
}
// Evaluate BS solution at vector evaluationPoints,
// at each time value prior to expiration.
double[][] BSValuesOption = new double[nTGrid][];
for (int i7 = 0; i7 < nTGrid; i7++)
{
BSValuesOption[i7] = new double[nv];
}
for (int i = 0; i < nTGrid; i++)
{
/*
* Black-Scholes (BS) European call option
* value = ValBSEC(S,t) = exp(-q*tau)*S*N01CDF(d1) -
* exp(-r*tau)*K*N01CDF(d2),
* where:
* tau = time to maturity = T - t;
* q = annual dividend yield;
* r = risk free rate;
* K = strike price;
* S = stock price;
* N01CDF(x) = N(0,1)_CDF(x);
* d1 = ( log( S/K ) +
* ( r - q + 0.5*sigma**2 )*tau ) /
* ( sigma*sqrt(tau) );
* d2 = d1 - sigma*sqrt(tau)
*/
// BS option values for tau = tGrid[i]:
double tau = tGrid[i];
double sigsqtau =
Math.Pow(sigma[volatility - 1], 2) * tau;
double sqrt_sigsqtau = Math.Sqrt(sigsqtau);
double sigsq = sigma[volatility - 1] *
sigma[volatility - 1];
for (int j = 0; j < nv; j++)
{
// Values of the underlying price where
// evaluations are made:
double S = evaluationPoints[j];
double d1 = (Math.Log(S /
strikePrices[strikePrice - 1]) + (r - dividend)
* tau + 0.5 * sigsqtau) / sqrt_sigsqtau;
double n01pdf_d1 =
Math.Exp((-0.5) * d1 * d1) / sqrt2pi;
double nu = Math.Exp((-dividend) * tau) *
S * n01pdf_d1 * Math.Sqrt(tau);
if (iTDeriv == 0)
{
if (iSDeriv == 0)
{
double d2 = d1 - sqrt_sigsqtau;
if (iSigDeriv == 0)
{
if (iRDeriv == 0)
{
BSValuesOption[i][j] =
Math.Exp((-dividend) * tau) * S *
Cdf.Normal(d1) -
Math.Exp((-r) * tau) *
strikePrices[strikePrice - 1] *
Cdf.Normal(d2);
}
else
{
BSValuesOption[i][j] =
Math.Exp((-r) * tau) *
strikePrices[strikePrice - 1]
* tau * Cdf.Normal(d2);
}
}
else if (iSigDeriv == 1)
{
//greek = Vega
BSValuesOption[i][j] = nu;
}
else if (iSigDeriv == 2)
{
//greek = Volga
BSValuesOption[i][j] = nu * d1 * d2 /
sigma[volatility - 1];
}
}
else if (iSDeriv == 1)
{
//greek = delta
if (iSigDeriv == 0)
{
BSValuesOption[i][j] =
Math.Exp((-dividend) * tau) *
Cdf.Normal(d1);
}
else if (iSigDeriv == 1)
{
//greek = Vanna
BSValuesOption[i][j] = (nu / S) *
(1.0 - d1 / sqrt_sigsqtau);
}
}
else if (iSDeriv == 2)
{
if (ispeed == 0)
{
//greek = gamma
BSValuesOption[i][j] =
Math.Exp((-dividend) * tau) *
n01pdf_d1 /(S * sqrt_sigsqtau);
}
else
{
//greek = speed
BSValuesOption[i][j] =
(-Math.Exp((-dividend) * tau)) *
n01pdf_d1 * (1.0 + d1 / sqrt_sigsqtau)
/ (S * S * sqrt_sigsqtau);
}
}
else if (iSDeriv == 3 && ispeed == 2)
{
//greek = speed
BSValuesOption[i][j] =
(-Math.Exp((-dividend) * tau)) *
n01pdf_d1 * (1.0 + d1 / sqrt_sigsqtau) /
(S * S * sqrt_sigsqtau);
}
}
else if (iTDeriv == 1)
{
double d2 = d1 - sqrt_sigsqtau;
if (iSDeriv == 0)
{
//greek = theta
BSValuesOption[i][j] =
Math.Exp((-dividend) * tau) * S *
(dividend * Cdf.Normal(d1) - 0.5 * sigsq *
n01pdf_d1 / sqrt_sigsqtau) - r *
Math.Exp((-r) * tau) *
strikePrices[strikePrice - 1] *
Cdf.Normal(d2);
}
else if (iSDeriv == 1)
{
//greek = charm
BSValuesOption[i][j] =
Math.Exp((-dividend) * tau) * ((-dividend)
* Cdf.Normal(d1) + n01pdf_d1 *
((r - dividend) * tau - 0.5 * d2 *
sqrt_sigsqtau) / (tau * sqrt_sigsqtau));
}
else if (iSDeriv == 2)
{
//greek = color
BSValuesOption[i][j] =
(-Math.Exp((-dividend) * tau)) * n01pdf_d1 *
(2.0 * dividend * tau + 1.0 + d1 *
(2.0 * (r - dividend) * tau - d2 *
sqrt_sigsqtau) / sqrt_sigsqtau) /
(2.0 * S * tau * sqrt_sigsqtau);
}
}
}
}
double relerrmax = 0.0;
int gNi = 3 * iTDeriv + iSDeriv;
if (iSigDeriv == 1 && iSDeriv == 0)
gNi = 6; //vega
if (iSigDeriv == 2 && iSDeriv == 0)
gNi = 7; //volga
if (iSigDeriv == 1 && iSDeriv == 1)
gNi = 8; //vanna
if (iRDeriv == 1)
gNi = 9; //rho
if (ispeed >= 1)
gNi = 9 + ispeed; //speed
for (int i = 0; i < nv; i++)
{
for (int j = 0; j < nt; j++)
{
double sVo = splineValuesOption[j][i];
double BSVo = BSValuesOption[j][i];
//greeks(itd=1) ~ d/dtau = -d/dt for iSDeriv > 0:
if (iTDeriv == 1 && iSDeriv > 0)
sVo = -sVo;
double relerr = Math.Abs((sVo - BSVo) / BSVo);
if (relerr > relerrmax)
relerrmax = relerr;
reavg[gNi] += relerr;
irect[gNi] += 1;
}
}
if (relerrmax > rex[gNi])
rex[gNi] = relerrmax;
if ((volatility == 1) && (strikePrice == 1))
{
for (int i = 0; i < nv; i++)
{
Console.Out.Write(
" underlying price: {0,4:f1}; FK " +
greekName[gNi] + ": ", evaluationPoints[i]);
for (int j = 0; j < nt; j++)
{
double sVo = splineValuesOption[j][i];
//greeks(itd=1) ~ d/dtau = -d/dt for isd > 0:
if (iTDeriv == 1 && iSDeriv > 0)
sVo = -sVo;
Console.Out.Write("{0,13:f10} ", sVo);
}
Console.Out.WriteLine();
Console.Out.Write(" BS "
+ greekName[gNi] + ": ");
for (int j = 0; j < nt; j++)
{
Console.Out.Write("{0,13:f10} ",
BSValuesOption[j][i]);
}
Console.Out.WriteLine();
}
}
}
}
}
}
}
}
class MyCoefficients : FeynmanKac.IPdeCoefficients
{
const double zero = 0.0;
const double half = 0.5;
double sigma, strikePrice, interestRate;
double alpha, dividend;
public MyCoefficients(double[] myData)
{
this.strikePrice = myData[0];
this.sigma = myData[2];
this.alpha = myData[3];
this.interestRate = myData[4];
this.dividend = myData[5];
}
// The coefficient sigma(x)
public double Sigma(double x, double t)
{
return (sigma * Math.Pow(x, alpha * half));
}
// The coefficient derivative d(sigma) / dx
public double SigmaPrime(double x, double t)
{
return (half * alpha * sigma * Math.Pow(x, alpha * half - 1.0));
}
// 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 xMax, df, interestRate, strikePrice;
public MyBoundaries(double[] myData)
{
this.strikePrice = myData[0];
this.xMax = myData[1];
this.interestRate = myData[4];
}
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)
{
df = Math.Exp(interestRate * time);
bndCoeffs[0, 0] = 1.0;
bndCoeffs[0, 1] = 0.0;
bndCoeffs[0, 2] = 0.0;
bndCoeffs[0, 3] = xMax - df * strikePrice;
bndCoeffs[1, 0] = 0.0;
bndCoeffs[1, 1] = 1.0;
bndCoeffs[1, 2] = 0.0;
bndCoeffs[1, 3] = 1.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)
{
double zero = 0.0;
// The payoff function
double val = Math.Max(x - strikePrice, zero);
return val;
}
}
}
Constant Elasticity of Variance Model for Vanilla Call Option
Interest Rate: 0.050 Continuous Dividend: 0.000
Minimum and Maximum Prices of Underlying: 0.00 60.00
Number of equally spaced spline knots: 120
Number of unknowns: 363
Time in Years Prior to Expiration: 0.0833 0.3333 0.5833
Option valued at Underlying Prices: 19.00 20.00 21.00
3 point (parabolic) estimate of parameter derivatives;
epsfrac = 0.00100000
Strike=15.00, Sigma= 0.20, Alpha= 2.00:
years to expiration: 0.0833 0.3333 0.5833
underlying price: 19.0; FK Value: 4.0623732450 4.2575924184 4.4733805278
BS Value: 4.0623732509 4.2575929678 4.4733814062
underlying price: 20.0; FK Value: 5.0623700127 5.2505145764 5.4492418798
BS Value: 5.0623700120 5.2505143129 5.4492428547
underlying price: 21.0; FK Value: 6.0623699727 6.2485587059 6.4385718831
BS Value: 6.0623699726 6.2485585270 6.4385720688
underlying price: 19.0; FK Rho: 1.2447807022 4.8365676562 8.0884594648
BS Rho: 1.2447806658 4.8365650322 8.0884502627
underlying price: 20.0; FK Rho: 1.2448021850 4.8929216544 8.3041708173
BS Rho: 1.2448021908 4.8929245641 8.3041638392
underlying price: 21.0; FK Rho: 1.2448024992 4.9107294561 8.4114197621
BS Rho: 1.2448024996 4.9107310444 8.4114199038
underlying price: 19.0; FK Vega: 0.0003289870 0.3487168323 1.1153520921
BS Vega: 0.0003295819 0.3487535501 1.1153536190
underlying price: 20.0; FK Vega: 0.0000056652 0.1224632724 0.6032458218
BS Vega: 0.0000056246 0.1224675413 0.6033084039
underlying price: 21.0; FK Vega: 0.0000000623 0.0376974472 0.3028275296
BS Vega: 0.0000000563 0.0376857196 0.3028629419
underlying price: 19.0; FK Volga: 0.0286253687 8.3705172571 16.7944556484
BS Volga: 0.0286064650 8.3691191978 16.8219823169
underlying price: 20.0; FK Volga: 0.0007135181 4.2505026165 12.9315443687
BS Volga: 0.0007186004 4.2519372748 12.9612638820
underlying price: 21.0; FK Volga: 0.0000100364 1.7613083880 8.6626164020
BS Volga: 0.0000097963 1.7617504949 8.6676581034
underlying price: 19.0; FK Delta: 0.9999864098 0.9877532309 0.9652249945
BS Delta: 0.9999863811 0.9877520034 0.9652261127
underlying price: 20.0; FK Delta: 0.9999998142 0.9964646548 0.9842482622
BS Delta: 0.9999998151 0.9964644003 0.9842476147
underlying price: 21.0; FK Delta: 0.9999999983 0.9990831687 0.9932459040
BS Delta: 0.9999999985 0.9990834124 0.9932451927
underlying price: 19.0; FK Vanna: -0.0012418872 -0.3391850563 -0.6388552010
BS Vanna: -0.0012431594 -0.3391932673 -0.6387423326
underlying price: 20.0; FK Vanna: -0.0000244490 -0.1366771953 -0.3945466660
BS Vanna: -0.0000244825 -0.1367114682 -0.3945405194
underlying price: 21.0; FK Vanna: -0.0000002905 -0.0466333335 -0.2187406645
BS Vanna: -0.0000002726 -0.0466323413 -0.2187858632
underlying price: 19.0; FK Gamma: 0.0000543456 0.0144908955 0.0264849216
BS Gamma: 0.0000547782 0.0144911447 0.0264824761
underlying price: 20.0; FK Gamma: 0.0000008315 0.0045912854 0.0129288434
BS Gamma: 0.0000008437 0.0045925328 0.0129280372
underlying price: 21.0; FK Gamma: 0.0000000080 0.0012817012 0.0058860348
BS Gamma: 0.0000000077 0.0012818272 0.0058865489
underlying price: 19.0; FK Speed: -0.0002127758 -0.0157070513 -0.0181086989
BS Speed: -0.0002123854 -0.0156192867 -0.0179536520
underlying price: 20.0; FK Speed: -0.0000037305 -0.0056184183 -0.0098403706
BS Speed: -0.0000037568 -0.0055859333 -0.0097472434
underlying price: 21.0; FK Speed: -0.0000000385 -0.0017185470 -0.0048615664
BS Speed: -0.0000000378 -0.0017082128 -0.0048130214
underlying price: 19.0; FK Speed2: -0.0002310655 -0.0156276977 -0.0179516855
BS Speed2: -0.0002123854 -0.0156192867 -0.0179536520
underlying price: 20.0; FK Speed2: -0.0000043215 -0.0055923924 -0.0097502997
BS Speed2: -0.0000037568 -0.0055859333 -0.0097472434
underlying price: 21.0; FK Speed2: -0.0000000475 -0.0017117661 -0.0048153107
BS Speed2: -0.0000000378 -0.0017082128 -0.0048130214
underlying price: 19.0; FK Theta: -0.7472631891 -0.8301000450 -0.8845209253
BS Theta: -0.7472638978 -0.8301108199 -0.8844992143
underlying price: 20.0; FK Theta: -0.7468881086 -0.7706770630 -0.8152217385
BS Theta: -0.7468880640 -0.7706789470 -0.8152097697
underlying price: 21.0; FK Theta: -0.7468815742 -0.7479185416 -0.7728950748
BS Theta: -0.7468815673 -0.7479153725 -0.7728982104
underlying price: 19.0; FK Charm: -0.0014382828 -0.0879903285 -0.0843323992
BS Charm: -0.0014397520 -0.0879913927 -0.0843403333
underlying price: 20.0; FK Charm: -0.0000284881 -0.0364107814 -0.0547260337
BS Charm: -0.0000285354 -0.0364209077 -0.0547074804
underlying price: 21.0; FK Charm: -0.0000003396 -0.0126436426 -0.0313343015
BS Charm: -0.0000003190 -0.0126437838 -0.0313252716
underlying price: 19.0; FK Color: 0.0051622176 0.0685064195 0.0299871130
BS Color: 0.0051777484 0.0684737183 0.0300398444
underlying price: 20.0; FK Color: 0.0001188761 0.0355826975 0.0274292189
BS Color: 0.0001205713 0.0355891884 0.0274307898
underlying price: 21.0; FK Color: 0.0000015431 0.0143174419 0.0190897160
BS Color: 0.0000015141 0.0143247729 0.0190752019
Greek: Value; avg rel err: 0.000146170676; max rel err: 0.009030693731
Greek: Delta; avg rel err: 0.000035817236; max rel err: 0.001158483024
Greek: Gamma; avg rel err: 0.001088309906; max rel err: 0.044839095787
Greek: Theta; avg rel err: 0.000054196357; max rel err: 0.001412847204
Greek: Charm; avg rel err: 0.001213365580; max rel err: 0.064577964461
Greek: Color; avg rel err: 0.003323775096; max rel err: 0.136355625079
Greek: Vega; avg rel err: 0.001512805322; max rel err: 0.106097327138
Greek: Volga; avg rel err: 0.058535091480; max rel err: 1.639568432709
Greek: Vanna; avg rel err: 0.001061842143; max rel err: 0.065654941418
Greek: Rho; avg rel err: 0.000146868204; max rel err: 0.009438785575
Greek: Speed; avg rel err: 0.007252489626; max rel err: 0.123630427780
Greek: Speed2; avg rel err: 0.008430324710; max rel err: 0.255782408454
Link to C# source.