Example 5: Calculating the Greeks using the Feynman-Kac Class

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. \alpha = 2. The second method calculates the Greeks using the closed-form BS evaluations which can be found at http://en.wikipedia.org/wiki/The_Greeks.

Example 5 calculates FK and BS solutions V(S, t) to the BS problem and the following Greeks:

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 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 Example 5), 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 \{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 Example 5, \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 Example 5 to calculate the Greek Speed . (For comparison purposes, Speed is also calculated directly in Example 5 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 Example 5 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.

import com.imsl.math.*;
import com.imsl.stat.*;
import java.util.*;
import com.imsl.IMSLException;

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 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][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;

         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 value = Math.max(x - strikePrice, zero);
         return value;
      }
   }
}

Output

  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.8929216545  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.3028275297 
                          BS   Vega:  0.0000000563  0.0376857196  0.3028629419 
 underlying price: 19.0;  FK  Volga:  0.0286253243  8.3705172571 16.7944557372 
                          BS  Volga:  0.0286064650  8.3691191978 16.8219823169 
 underlying price: 20.0;  FK  Volga:  0.0007137846  4.2505027498 12.9315444575 
                          BS  Volga:  0.0007186004  4.2519372748 12.9612638820 
 underlying price: 21.0;  FK  Volga:  0.0000100364  1.7613084768  8.6626165796 
                          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.0048153106 
                          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.009030693732
 Greek:  Delta;  avg rel err:  0.000035817236;  max rel err:  0.001158483024
 Greek:  Gamma;  avg rel err:  0.001088461868;  max rel err:  0.044851604910
 Greek:  Theta;  avg rel err:  0.000054196357;  max rel err:  0.001412847201
 Greek:  Charm;  avg rel err:  0.001213382637;  max rel err:  0.064579338827
 Greek:  Color;  avg rel err:  0.003323753774;  max rel err:  0.136355681844
 Greek:   Vega;  avg rel err:  0.001513788862;  max rel err:  0.106176227011
 Greek:  Volga;  avg rel err:  0.058530439717;  max rel err:  1.639567620281
 Greek:  Vanna;  avg rel err:  0.001062052211;  max rel err:  0.065672253096
 Greek:    Rho;  avg rel err:  0.000146868204;  max rel err:  0.009438786906
 Greek:  Speed;  avg rel err:  0.007251538925;  max rel err:  0.123630425554
 Greek: Speed2;  avg rel err:  0.008429577530;  max rel err:  0.255722275986
Link to Java source.