package com.imsl.test.example.math;
import com.imsl.math.*;
import com.imsl.stat.*;
import java.util.*;
/**
*
* Solves for the "Greeks" of mathematical finance.
*
*
* In this example, the FeynmanKac class is used to calculate the
* so-called "Greeks" of mathematical finance. The "Greeks" are defined in terms
* of various derivatives of Feynman-Kac solutions and have special
* interpretations in 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 Feynman-Kac (FK)
* solution to the diffusion model for call options as described in example 2
* {@link FeynmanKacEx2}, for the Black-Scholes case (\( \alpha \rightarrow
* 2\)). The second method calculates the Greeks using the closed-form
* Black-Scholes (BS) evaluations, which can be found here:
* Wikipedia:
* The Greeks.
*
*
*
* Example 5 calculates FK and BS solutions \(V(S, t)\) to the BS problem and
* the following Greeks:
*
* - \(Delta\;=\;\frac{\partial V}{\partial S}\) is the first derivative of
* the Value, \(V(S,t)\), of a portfolio of derivative security derived
* from an underlying instrument with respect to the price, S of the
* underlying instrument.
*
* - \(Gamma\;=\;\frac{\partial^2 V}{\partial S^2}\)
* - \(Theta\;=\;-\frac{\partial V}{\partial t}\) is the negative first
* derivative of V with respect to time t
* - \(Charm\;=\;\frac{\partial^2 V}{\partial S \partial t}\)
* - \(Color\;=\;\frac{\partial^3 V}{\partial S^2 \partial t}\)
* - \(Rho\;=\;-\frac{\partial V}{\partial r}\) is the first derivative of
* V with respect to risk free rate r
* - \(Vega\;=\;\frac{\partial V}{\partial \sigma}\) measures sensitivity to
* volatility parameter \(\sigma\) of the underlying S
* - \(Volga\;=\;\frac{\partial^2 V}{\partial \sigma^2}\)
* - \(Vanna\;=\;\frac{\partial^2 V}{\partial S \partial \sigma}\)
* - \(Speed\;=\;\frac{\partial^3 V}{\partial S^3}\)
*
*
*
* 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 \(V(S, t)\) 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 \(\sigma\). 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 \(f(\sigma + \epsilon)\) truncated to
* m+1 terms (to allow an m-degree polynomial fit of m+1
* data points):
*
*
* $$ f(\sigma + \epsilon)\;=\;\sum_{n = 0}^m {\frac{f^{(n)}(\sigma)}{n!}
* \epsilon^n} $$
*
*
* we are able to derive the following parabolic (3 point) estimation of first
* and second derivatives \(f^{(1)}(\sigma)\) and \(f^{(2)}(\sigma)\) in terms
* of the three values: \(f(\sigma - \epsilon)\), \(f(\sigma)\), and \(f(\sigma
* + \epsilon)\), where \(\epsilon = \epsilon_{frac} \sigma\) and \(0 \lt
* \epsilon_{frac} \lt\lt 1\):
*
*
* $$ 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} $$
*
*
* $$ 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} $$
*
*
* Similarly, the quartic (5 point) estimation of \(f^{(1)}(\sigma)\) and
* \(f^{(2)}(\sigma)\) in terms of \(f(\sigma - 2\epsilon)\), \(f(\sigma -
* \epsilon)\), \(f(\sigma)\), \(f(\sigma + \epsilon)\), and \(f(\sigma +
* 2\epsilon)\) is:
*
*
* $$ f^{(1)}(\sigma)\;\approx\;\frac{4}{3} f^{[1]}(\sigma, \epsilon) \;-\;
* \frac{1}{3}f^{[1]}(\sigma, 2\epsilon) $$
*
*
* $$ f^{(2)}(\sigma)\;\approx\;\frac{4}{3} f^{[2]}(\sigma, \epsilon) \;-\;
* \frac{1}{3}f^{[2]}(\sigma, 2\epsilon) $$
*
*
* 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 Hanson, R. (2008), 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 the 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 \(\{S,\; f^{(2)}(S)\}\) and the second and
* third solutions for solution grid and evaluation point sets \(\{S+\epsilon,\;
* f^{(2)}(S+\epsilon)\}\) and \(\{S-\epsilon,\; f^{(2)}(S-\epsilon)\}\), where
* the solution grid and evaluation point sets are shifted up and down by
* \(\epsilon\). In this example, \(\epsilon\) is set to \(\epsilon_{frac}
* \bar{S}\), where \(\bar{S}\) is the average value of S over the range
* of grid values and \(0 \lt \epsilon_{frac} \lt\lt 1\). The third derivative
* solution can then be obtained using the parabolic estimate:
*
*
* $$ f^{(3)}(S)\;=\;\frac{\partial f^{(2)}(S)}{\partial
* S}\;\approx\;\frac{f^{(2)}(S + \epsilon)-f^{(2)}(S - \epsilon)}{2 \epsilon}
* $$
*
*
* This procedure is implemented in this example to calculate the Greek
* Speed. (For comparison purposes, Speed
* is also calculated directly in this example 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 method
* setInitialStepsize), by choosing more closely spaced S
* and t grid points (by adjusting method
* computeCoefficients input parameter arrays xGrid
* and tGrid), and by adjusting \(\epsilon_{frac}\) so that the
* central differences used to calculate the derivatives are not too small to
* compromise accuracy.
*
*
* @see Code
* @see Output
*
*/
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, xMax = 60.0;
// Value of the interest rate and continuous dividend
private static double r = 0.05, dividend = 0.0;
// Define parameters for the integration step.
private static int nXGgrid = 121, nTGrid = 3, nv = 3;
private static int nint = nXGgrid - 1, n = 3 * nXGgrid;
private static double[] xGrid = new double[nXGgrid];
private static double dx;
// Time dependency
private static boolean[] timeDependence = new boolean[7];
// Number of left/right boundary conditions
private static int nlbcd = 3, 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. * Math.PI);
private static String greekName[] = {
" Value", " Delta", " Gamma", " Theta", " Charm",
" Color", " Vega", " Volga", " Vanna", " Rho", " Speed", "Speed2"};
// Time values for the options
private static int nt = 3;
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[]) throws Exception {
// 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;
}
System.out.printf(" Constant Elasticity of Variance Model"
+ " for Vanilla Call Option%n");
System.out.printf(Locale.ENGLISH,
"%7sInterest Rate:%7.3f Continuous Dividend:%7.3f%n",
" ", r, dividend);
System.out.printf(Locale.ENGLISH,
"%7sMinimum and Maximum Prices of Underlying:%7.2f%7.2f%n",
" ", xMin, xMax);
System.out.printf("%7sNumber of equally spaced spline knots:%4d%n",
" ", nXGgrid - 1);
System.out.printf("%7sNumber of unknowns:%4d%n%n", " ", n);
System.out.printf(Locale.ENGLISH,
"%7sTime in Years Prior to Expiration: %7.4f%7.4f%7.4f%n",
" ", tGrid[0], tGrid[1], tGrid[2]); // tGrid[] = tau == T - t
System.out.printf(Locale.ENGLISH,
"%7sOption valued at Underlying Prices:%7.2f%7.2f%7.2f%n%n",
" ", evaluationPoints[0], evaluationPoints[1],
evaluationPoints[2]);
/*
* 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) {
System.out.printf(Locale.ENGLISH,
" 3 point (parabolic) estimate of "
+ "parameter derivatives;%n");
} else {
System.out.printf(Locale.ENGLISH,
" 5 point (quartic) estimate of "
+ "parameter derivatives%n");
}
System.out.printf(Locale.ENGLISH, " epsfrac = %11.8f%n", 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);
}
}
System.out.println();
for (int ig = 0; ig < 12; ig++) {
reavg[ig] /= irect[ig];
System.out.printf(Locale.ENGLISH, "%n Greek: " + greekName[ig]
+ "; avg rel err: %15.12f" + "; max rel err: %15.12f",
reavg[ig], rex[ig]);
}
System.out.println();
} // end main
private static void calculateGreeks(int volatility, int strikePrice,
int iquart) throws Exception {
int i1 = 1;
int nt = 3;
if ((volatility == 1) && (strikePrice == 1)) {
System.out.printf(Locale.ENGLISH, "%n"
+ " Strike=%5.2f, Sigma="
+ "%5.2f, Alpha=%5.2f:",
strikePrices[strikePrice - 1], sigma[volatility - 1],
alpha[i1 - 1]);
System.out.printf(Locale.ENGLISH, "%n"
+ " years to expiration: "
+ " %7.4f %7.4f %7.4f"
+ "%n",
tGrid[0], tGrid[1], tGrid[2]);
}
/* 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 methods.
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][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[][] optionCoefficients
= new double[tGrid.length + 1][3
* xGrid.length];
if (iTDeriv == 0) {
optionCoefficients
= callOption.
getSplineCoefficients();
} else {
optionCoefficients
= callOption.
getSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++) {
// FK option values for tau = tGrid[i]:
splineValuesOption[i]
= callOption.getSplineValue(
evaluationPoints,
optionCoefficients[i + 1],
iSDeriv);
}
}
if (iSigDeriv > 0 || iRDeriv > 0 || ispeed == 1) {
System.arraycopy(xGrid, 0, xGridP, 0, nXGgrid);
System.arraycopy(xGrid, 0, xGridM, 0, nXGgrid);
System.arraycopy(evaluationPoints, 0,
evaluationPointsP, 0, nv);
System.arraycopy(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. + epsfrac);
} else if (iRDeriv > 0) {
userData[4] = r * (1. + 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
= new double[tGrid.length + 1][3
* xGrid.length];
if (iTDeriv == 0) {
optionCoefficientsP
= callOptionP.
getSplineCoefficients();
} else {
optionCoefficientsP
= callOptionP.
getSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++) {
// FK option values for tau = tGrid[i]:
splineValuesOptionP[i]
= callOptionP.getSplineValue(
evaluationPointsP,
optionCoefficientsP[i + 1],
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. - epsfrac);
} else {
userData[4] = r * (1. - 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
= new double[tGrid.length + 1][3
* xGrid.length];
if (iTDeriv == 0) {
optionCoefficientsM
= callOptionM.
getSplineCoefficients();
} else {
optionCoefficientsM
= callOptionM.
getSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++) {
// FK option values for tau = tGrid[i]:
splineValuesOptionM[i]
= callOptionM.
getSplineValue(
evaluationPointsM,
optionCoefficientsM[i + 1],
iSDeriv);
}
if (iquart == 1) {
System.arraycopy(xGrid, 0, xGridPP, 0,
nXGgrid);
System.arraycopy(xGrid, 0, xGridMM, 0,
nXGgrid);
System.arraycopy(evaluationPoints, 0,
evaluationPointsPP, 0, nv);
System.arraycopy(evaluationPoints, 0,
evaluationPointsMM, 0, nv);
if (ispeed == 1) {
for (int i = 0; i < nXGgrid; i++) {
xGridPP[i] = xGrid[i] + 2. * dx2;
xGridMM[i] = xGrid[i] - 2. * dx2;
}
for (int i = 0; i < nv; i++) {
evaluationPointsPP[i]
= evaluationPoints[i]
+ 2. * dx2;
evaluationPointsMM[i]
= evaluationPoints[i]
- 2. * dx2;
}
xMaxPP = xMax + 2. * dx2;
xMaxMM = xMax - 2. * dx2;
}
userData[1] = xMaxPP;
if (iSigDeriv > 0) {
// calculate spline values for
// sigmaPP = sigma[i2-1]
// *(1. + 2.*epsfrac):
userData[2]
= sigma[volatility - 1]
* (1. + 2. * epsfrac);
} else if (iRDeriv > 0) {
// calculate spline values for
// rPP = r*(1. + 2.*epsfrac):
userData[4] = r * (1. + 2. * 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[][] optionCoefficientsPP
= new double[tGrid.length + 1][3
* xGrid.length];
if (iTDeriv == 0) {
optionCoefficientsPP
= callOptionPP.
getSplineCoefficients();
} else {
optionCoefficientsPP
= callOptionPP.
getSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++) {
// FK option values for tau = tGrid[i]:
splineValuesOptionPP[i]
= callOptionPP.getSplineValue(
evaluationPointsPP,
optionCoefficientsPP[i + 1],
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. - 2. * epsfrac);
} else if (iRDeriv > 0) {
userData[4] = r * (1. - 2. * 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[][] optionCoefficientsMM
= new double[tGrid.length + 1][3
* xGrid.length];
if (iTDeriv == 0) {
optionCoefficientsMM
= callOptionMM.
getSplineCoefficients();
} else {
optionCoefficientsMM
= callOptionMM.
getSplineCoefficientsPrime();
}
for (int i = 0; i < nTGrid; i++) {
// FK option values for tau = tGrid[i]:
splineValuesOptionMM[i]
= callOptionMM.getSplineValue(
evaluationPointsMM,
optionCoefficientsMM[i + 1],
iSDeriv);
}
}
}
if (iSigDeriv == 1 || iRDeriv == 1 || ispeed == 1) {
double eps = 0., f11 = 0., f12 = 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. * eps);
if (iquart == 0) {
splineValuesOption[i][j] = f11;
} else {
f12 = (splineValuesOptionPP[i][j]
- splineValuesOptionMM[i][j])
/ (4. * eps);
splineValuesOption[i][j] = (4.
* f11 - f12) / 3.;
}
}
}
}
if (iSigDeriv == 2) {
double eps = sigma[volatility - 1] * epsfrac;
double f21 = 0.;
double f22 = 0.;
for (int i = 0; i < nTGrid; i++) {
for (int j = 0; j < nv; j++) {
f21 = (splineValuesOptionP[i][j]
+ splineValuesOptionM[i][j] - 2.
* splineValuesOption[i][j])
/ (eps * eps);
if (iquart == 0) {
splineValuesOption[i][j] = f21;
} else {
f22 = (splineValuesOptionPP[i][j]
+ splineValuesOptionMM[i][j]
- 2.
* splineValuesOption[i][j])
/ (4. * eps * eps);
splineValuesOption[i][j]
= (4. * f21 - f22) / 3.;
}
}
}
}
// Evaluate BS solution at vector evaluationPoints,
// at each time value prior to expiration.
double[][] BSValuesOption = new double[nTGrid][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. - 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. + 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. + 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. * dividend * tau + 1.
+ d1 * (2. * (r - dividend)
* tau - d2 * sqrt_sigsqtau)
/ sqrt_sigsqtau)
/ (2. * S * tau
* sqrt_sigsqtau);
}
}
}
}
double relerrmax = 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], 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++) {
System.out.printf(" underlying price:"
+ " %4.1f; 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;
}
System.out.printf(Locale.ENGLISH,
"%13.10f ", sVo);
}
System.out.println();
System.out.printf(" "
+ " BS "
+ greekName[gNi] + ": ");
for (int j = 0; j < nt; j++) {
System.out.printf(Locale.ENGLISH,
"%13.10f ",
BSValuesOption[j][i]);
}
System.out.println();
}
}
}
}
}
}
}
}
static class MyCoefficients implements FeynmanKac.PdeCoefficients {
final double zero = 0.0, half = 0.5;
private double sigma, strikePrice, interestRate;
private 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;
}
}
static class MyBoundaries implements FeynmanKac.Boundaries {
private 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;
}
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;
}
public double terminal(double x) {
double zero = 0.0;
// The payoff function
double value = Math.max(x - strikePrice, zero);
return value;
}
}
}