A diffusion model for Call Options

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:

power $\alpha = {2.0, 1.0, 0.0}$
strike price
volatility $\sigma = {0.2, 0.3, 0.4}$
times until expiration = {1/12, 4 /12, 7 /12}
underlying prices = {19.0, 20.0, 21.0}
interest rate
$x_{\min} = 0, x_{\max} = 60$
$nxGrid = 121, \, n = 3 \times nxGrid = 363$

With this model the Feynman-Kac differential equation is defined by identifying:

$x:\; S$
$\sigma(x,t):\; \sigma x^{\alpha/2};\; \frac{\partial \sigma}{\partial x}= \frac{\alpha \sigma}{2}x^{\alpha/2-1}$
$\mu(x,t):\;rx$
$\kappa(x,t):\; r$
$\phi(f,x,t) \equiv 0$

The payoff function is the "vanilla option", $p(x) = \max(x-K, 0)$ .

using System;
using Imsl.Math;

public class FeynmanKacEx2
{
   public static void Main(String[] args)
   {
      // 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
      bool[] timeDependence = new bool[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;
      }

      Console.Out.WriteLine("  Constant Elasticity of Variance Model" +
                            " for Vanilla Call");
      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}", nxGrid - 1);
      Console.Out.WriteLine("       Number of unknowns:{0,4:d}", n);
      Console.Out.WriteLine("       Time in Years Prior to Expiration: " +
                            "{0,7:f4}{1,7:f4}{2,7:f4}",
                            tGrid[0], tGrid[1], tGrid[2]);
      Console.Out.WriteLine("       Option valued at Underlying Prices:" +
                            "{0,7:f2}{1,7:f2}{2,7:f2}\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[] optionTimeCoeffs = new double[n];
               double[][] splineValuesOption = new double[ntGrid][];

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

               Console.Out.WriteLine("  Strike={0,5:f2}, Sigma={1,5:f2}, "
                                     + "Alpha={2,5:f2}",
                                       strikePrices[i3 - 1],
                                       sigma[i2 - 1],
                                       alpha[i1 - 1]);

               for (int i = 0; i < nv; i++)
               {
                  Console.Out.Write("                       Call " +
                                    "Option Values  ");
                  for (int j = 0; j < nt; j++)
                  {
                     Console.Out.Write("{0,7:f4} ", splineValuesOption[j][i]);
                  }
                  Console.Out.WriteLine();
               }
               Console.Out.WriteLine();
            }
         }
      }
   }

   internal class MyCoefficients : FeynmanKac.IPdeCoefficients
   {
      const double zero = 0.0;
      const double 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;
      }
   }

   internal class MyBoundaries : FeynmanKac.IBoundaries
   {

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

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
       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 C# source.