Example 2: A diffusion model for Call Options

In Beckers (1980) there is a model for a Stochastic Differential Equation of option pricing. The idea is a "constant elasticity of variance diffusion (or CEV) class"

 dS = \mu S dt + \sigma S^{\alpha/2}dW,\; 0 \le \alpha \lt 2.
The Black-Scholes model is the limiting case \alpha \to 2. A numerical solution of this diffusion model yields the price of a call option. Various values of the strike price K , time values, \sigma and power coefficient \alpha are used to evaluate the option price at values of the underlying price. The sets of parameters in the computation are:

  1. power \alpha = {2.0, 1.0, 0.0}
  2. strike price  K = {15.0, 20.0, 25.0}
  3. volatility \sigma = {0.2, 0.3, 0.4}
  4. times until expiration = {1/12, 4 /12, 7 /12}
  5. underlying prices = {19.0, 20.0, 21.0}
  6. interest rate r = 0.05
  7. x_{\min} = 0, x_{\max} = 60
  8. nx = 121, \, n = 3 \times nx = 363
With this model the Feynman-Kac differential equation is defined by identifying: The payoff function is the "vanilla option", p(x) = \max(x-K, 0).

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

public class FeynmanKacEx2 {

   public static void main(String args[]) throws Exception {
      // Compute Constant Elasticity of Variance Model for Vanilla Call
      // The set of strike prices
      double[] strikePrices = {15.0, 20.0, 25.0};
      // The set of sigma values
      double[] sigma = {0.2, 0.3, 0.4};
      // The set of model diffusion powers
      double[] alpha = {2.0, 1.0, 0.0};
      // Time values for the options
      int nt = 3;
      double[] tGrid = {1.0 / 12.0, 4.0 / 12.0, 7.0 / 12.0};
      // Values of the underlying where evaluations are made
      double[] evaluationPoints = {19.0, 20.0, 21.0};
      // Value of the interest rate and continuous dividend
      double r = 0.05, dividend = 0.0;
      // Values of the min and max underlying values modeled
      double xMin = 0.0, xMax = 60.0;
      // Define parameters for the integration step. */
      int nxGrid = 121;
      int ntGrid = 3;
      int nv = 3;
      int nint = nxGrid - 1, n = 3 * nxGrid;
      double[] xGrid = new double[nxGrid];
      double dx;
      // Number of left/right boundary conditions
      int nlbcd = 3, nrbcd = 3;
      // Time dependency
      boolean[] timeDependence = new boolean[7];
      double[] userData = new double[6];
      int j;
      // Define equally-spaced grid of points for the underlying price
      dx = (xMax - xMin) / ((double) nint);
      xGrid[0] = xMin;
      xGrid[nxGrid - 1] = xMax;
      for (int i = 2; i <= nxGrid - 1; i++) {
         xGrid[i - 1] = xGrid[i - 2] + dx;
      }

      System.out.printf("%2sConstant Elasticity of Variance Model" +
            " for Vanilla Call%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",
            " ", nxGrid - 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]);
      System.out.printf(Locale.ENGLISH,
            "%7sOption valued at Underlying Prices:%7.2f%7.2f%7.2f%n%n",
            " ", evaluationPoints[0], evaluationPoints[1],
            evaluationPoints[2]);

      for (int i1 = 1; i1 <= 3; i1++) /* Loop over power */ {
         for (int i2 = 1; i2 <= 3; i2++) /* Loop over volatility */ {
            for (int i3 = 1; i3 <= 3; i3++) /* Loop over strike price */ {
               // Pass data through into evaluation routines.
               userData[0] = strikePrices[i3 - 1];
               userData[1] = xMax;
               userData[2] = sigma[i2 - 1];
               userData[3] = alpha[i1 - 1];
               userData[4] = r;
               userData[5] = dividend;

               FeynmanKac callOption = new FeynmanKac(
                     new MyCoefficients(userData));

               // Right boundary condition is time-dependent
               timeDependence[5] = true;
               callOption.setTimeDependence(timeDependence);
               callOption.computeCoefficients(nlbcd, nrbcd,
                     new MyBoundaries(userData), xGrid, tGrid);
               double[][] optionCoefficients=callOption.getSplineCoefficients();
               double[][] splineValuesOption = new double[ntGrid][];

               // Evaluate solution at vector evaluationPoints, at each time
               // value prior to expiration.
               for (int i = 0; i < ntGrid; i++) {
                  splineValuesOption[i] = callOption.getSplineValue(
                        evaluationPoints, optionCoefficients[i + 1], 0);
               }

               System.out.printf(Locale.ENGLISH, "%2sStrike=%5.2f, Sigma=" +
                     "%5.2f, Alpha=%5.2f%n", " ", strikePrices[i3 - 1],
                     sigma[i2 - 1], alpha[i1 - 1]);
               for (int i = 0; i < nv; i++) {
                  System.out.printf("%23sCall Option Values%2s", " ", " ");
                  for (j = 0; j < nt; j++) {
                     System.out.printf(Locale.ENGLISH, "%7.4f ",
                           splineValuesOption[j][i]);
                  }
                  System.out.printf("%n");
               }
               System.out.printf("%n");
            }
         }
      }
   }

   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
       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

  Strike=15.00, Sigma= 0.20, Alpha= 2.00
                       Call Option Values   4.0624  4.2576  4.4734 
                       Call Option Values   5.0624  5.2505  5.4492 
                       Call Option Values   6.0624  6.2486  6.4386 

  Strike=20.00, Sigma= 0.20, Alpha= 2.00
                       Call Option Values   0.1312  0.5951  0.9693 
                       Call Option Values   0.5024  1.0880  1.5093 
                       Call Option Values   1.1980  1.7478  2.1745 

  Strike=25.00, Sigma= 0.20, Alpha= 2.00
                       Call Option Values   0.0000  0.0111  0.0751 
                       Call Option Values   0.0000  0.0376  0.1630 
                       Call Option Values   0.0006  0.1036  0.3150 

  Strike=15.00, Sigma= 0.30, Alpha= 2.00
                       Call Option Values   4.0636  4.3405  4.6627 
                       Call Option Values   5.0625  5.2951  5.5794 
                       Call Option Values   6.0624  6.2712  6.5248 

  Strike=20.00, Sigma= 0.30, Alpha= 2.00
                       Call Option Values   0.3107  1.0261  1.5479 
                       Call Option Values   0.7317  1.5404  2.0999 
                       Call Option Values   1.3762  2.1674  2.7362 

  Strike=25.00, Sigma= 0.30, Alpha= 2.00
                       Call Option Values   0.0005  0.1124  0.3564 
                       Call Option Values   0.0035  0.2184  0.5565 
                       Call Option Values   0.0184  0.3869  0.8230 

  Strike=15.00, Sigma= 0.40, Alpha= 2.00
                       Call Option Values   4.0755  4.5143  4.9673 
                       Call Option Values   5.0660  5.4210  5.8328 
                       Call Option Values   6.0633  6.3588  6.7306 

  Strike=20.00, Sigma= 0.40, Alpha= 2.00
                       Call Option Values   0.5109  1.4625  2.1260 
                       Call Option Values   0.9611  1.9934  2.6915 
                       Call Option Values   1.5807  2.6088  3.3202 

  Strike=25.00, Sigma= 0.40, Alpha= 2.00
                       Call Option Values   0.0081  0.3302  0.7795 
                       Call Option Values   0.0287  0.5178  1.0656 
                       Call Option Values   0.0820  0.7690  1.4097 

  Strike=15.00, Sigma= 0.20, Alpha= 1.00
                       Call Option Values   4.0624  4.2479  4.4312 
                       Call Option Values   5.0624  5.2479  5.4312 
                       Call Option Values   6.0624  6.2479  6.4312 

  Strike=20.00, Sigma= 0.20, Alpha= 1.00
                       Call Option Values   0.0000  0.0219  0.1051 
                       Call Option Values   0.1497  0.4107  0.6484 
                       Call Option Values   1.0832  1.3314  1.5773 

  Strike=25.00, Sigma= 0.20, Alpha= 1.00
                       Call Option Values  -0.0000 -0.0000  0.0000 
                       Call Option Values  -0.0000 -0.0000  0.0000 
                       Call Option Values  -0.0000 -0.0000  0.0000 

  Strike=15.00, Sigma= 0.30, Alpha= 1.00
                       Call Option Values   4.0624  4.2479  4.4312 
                       Call Option Values   5.0624  5.2479  5.4312 
                       Call Option Values   6.0624  6.2479  6.4312 

  Strike=20.00, Sigma= 0.30, Alpha= 1.00
                       Call Option Values   0.0010  0.0786  0.2208 
                       Call Option Values   0.1993  0.4997  0.7539 
                       Call Option Values   1.0835  1.3444  1.6022 

  Strike=25.00, Sigma= 0.30, Alpha= 1.00
                       Call Option Values  -0.0000  0.0000  0.0000 
                       Call Option Values  -0.0000  0.0000  0.0000 
                       Call Option Values  -0.0000  0.0000  0.0004 

  Strike=15.00, Sigma= 0.40, Alpha= 1.00
                       Call Option Values   4.0624  4.2479  4.4312 
                       Call Option Values   5.0624  5.2479  5.4312 
                       Call Option Values   6.0624  6.2479  6.4312 

  Strike=20.00, Sigma= 0.40, Alpha= 1.00
                       Call Option Values   0.0072  0.1540  0.3446 
                       Call Option Values   0.2498  0.5950  0.8728 
                       Call Option Values   1.0868  1.3795  1.6586 

  Strike=25.00, Sigma= 0.40, Alpha= 1.00
                       Call Option Values   0.0000  0.0000  0.0000 
                       Call Option Values   0.0000  0.0000  0.0005 
                       Call Option Values   0.0000  0.0002  0.0057 

  Strike=15.00, Sigma= 0.20, Alpha= 0.00
                       Call Option Values   4.0624  4.2479  4.4311 
                       Call Option Values   5.0624  5.2479  5.4312 
                       Call Option Values   6.0624  6.2479  6.4312 

  Strike=20.00, Sigma= 0.20, Alpha= 0.00
                       Call Option Values   0.0001  0.0001  0.0002 
                       Call Option Values   0.0816  0.3316  0.5748 
                       Call Option Values   1.0817  1.3308  1.5748 

  Strike=25.00, Sigma= 0.20, Alpha= 0.00
                       Call Option Values   0.0000 -0.0000 -0.0000 
                       Call Option Values   0.0000 -0.0000 -0.0000 
                       Call Option Values  -0.0000  0.0000 -0.0000 

  Strike=15.00, Sigma= 0.30, Alpha= 0.00
                       Call Option Values   4.0624  4.2479  4.4312 
                       Call Option Values   5.0624  5.2479  5.4312 
                       Call Option Values   6.0624  6.2479  6.4312 

  Strike=20.00, Sigma= 0.30, Alpha= 0.00
                       Call Option Values   0.0000 -0.0000  0.0026 
                       Call Option Values   0.0894  0.3326  0.5753 
                       Call Option Values   1.0826  1.3306  1.5749 

  Strike=25.00, Sigma= 0.30, Alpha= 0.00
                       Call Option Values   0.0000 -0.0000  0.0000 
                       Call Option Values   0.0000 -0.0000  0.0000 
                       Call Option Values   0.0000 -0.0000 -0.0000 

  Strike=15.00, Sigma= 0.40, Alpha= 0.00
                       Call Option Values   4.0624  4.2479  4.4312 
                       Call Option Values   5.0624  5.2479  5.4312 
                       Call Option Values   6.0624  6.2479  6.4312 

  Strike=20.00, Sigma= 0.40, Alpha= 0.00
                       Call Option Values  -0.0000  0.0001  0.0108 
                       Call Option Values   0.0985  0.3383  0.5781 
                       Call Option Values   1.0830  1.3306  1.5749 

  Strike=25.00, Sigma= 0.40, Alpha= 0.00
                       Call Option Values   0.0000  0.0000  0.0000 
                       Call Option Values   0.0000 -0.0000  0.0000 
                       Call Option Values   0.0000 -0.0000  0.0000 

Link to Java source.