Example 3: European Option with two payoff strategies

This example evaluates the price of a European Option using two payoff strategies: Cash-or-Nothing and Vertical Spread. In the first case the payoff function is

 p(x)=\left\{ \begin{array}{lr} 0, & x \le K \\ B, & x \gt K \end{array} \right. .
The value B is regarded as the bet on the asset price, see Wilmott et al. (1995, p. 39-40). The second case has the payoff function
 p(x) = \max(x-K_1)-\max(x-K_2), \, K_2>K_1.
Both problems use the same boundary conditions. Each case requires a separate integration of the Black-Scholes differential equation, but only the payoff function evaluation differs in each case. The sets of parameters in the computation are:

  1. Strike and bet prices K_1 = {10.0},\; K_2={15.0}, and B={2.0}.
  2. Volatility  \sigma = {0.4}
  3. Times until expiration = {1/4, 1/2}
  4. Interest rate r = 0.1
  5. x_{\min}=0.0,\, x_{\max}=30.0
  6. nx = 61, \, n=3 \times nx = 183

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

public class FeynmanKacEx3 {

    public static void main(String args[]) throws Exception {
        int i;
        int nxGrid = 61;
        int ntGrid = 2;
        int nv = 12;
        // The strike price
        double strikePrice = 10.0;
        // The spread value
        double spreadValue = 15.0;
        // The Bet for the Cash-or-Nothing Call
        double bet = 2.0;
        // The sigma value
        double sigma = 0.4;
        // Time values for the options
        int nt = 2;
        double[] tGrid = {0.25, 0.5};
        // Values of the underlying where evaluations are made
        double[] evaluationPoints = new double[nv];
        // Value of the interest rate and continuous dividend
        double r = 0.1, dividend = 0.0;
        // Values of the min and max underlying values modeled
        double xMin = 0.0, xMax = 30.0;
        // Define parameters for the integration step.
        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;
        // Structure for the evaluation methods.
        int iData;
        double[] rData = new double[7];
        boolean[] timeDependence = new boolean[7];

        // Define an equally-spaced grid of points for the
        // underlying price
        dx = (xMax - xMin) / ((double) (nint));
        xGrid[0] = xMin;
        xGrid[nxGrid - 1] = xMax;
        for (i = 2; i <= nxGrid - 1; i++) {
            xGrid[i - 1] = xGrid[i - 2] + dx;
        }
        for (i = 1; i <= nv; i++) {
            evaluationPoints[i - 1] = 2.0 + (i - 1) * 2.0;
        }
        rData[0] = strikePrice;
        rData[1] = bet;
        rData[2] = spreadValue;
        rData[3] = xMax;
        rData[4] = sigma;
        rData[5] = r;
        rData[6] = dividend;
        // Flag the difference in payoff functions
        // 1 for the Bet, 2 for the Vertical Spread
        // In both cases, no time dependencies for
        // the coefficients and the left boundary
        // conditions
        timeDependence[5] = true;
        iData = 1;
        FeynmanKac betOption = new FeynmanKac(new MyCoefficients(rData));
        betOption.setTimeDependence(timeDependence);
        betOption.computeCoefficients(nlbcd, nrbcd,
                new MyBoundaries(rData, iData), xGrid, tGrid);
        double[][] betOptionCoefficients = betOption.getSplineCoefficients();
        double[][] splineValuesBetOption = new double[ntGrid][];

        iData = 2;
        FeynmanKac spreadOption = new FeynmanKac(new MyCoefficients(rData));

        spreadOption.setTimeDependence(timeDependence);
        spreadOption.computeCoefficients(nlbcd, nrbcd,
                new MyBoundaries(rData, iData), xGrid, tGrid);
        double[][] spreadOptionCoefficients
                = spreadOption.getSplineCoefficients();
        double[][] splineValuesSpreadOption = new double[ntGrid][];

        // Evaluate solutions at vector evaluationPoints, at each time value
        // prior to expiration.
        for (i = 0; i < ntGrid; i++) {
            splineValuesBetOption[i]
                    = betOption.getSplineValue(evaluationPoints,
                            betOptionCoefficients[i + 1], 0);
            splineValuesSpreadOption[i] = spreadOption.getSplineValue(
                    evaluationPoints, spreadOptionCoefficients[i + 1], 0);
        }

        System.out.printf("%2sEuropean Option Value for A Bet%n", " ");
        System.out.printf("%3sand a Vertical Spread, 3 and 6 Months "
                + "Prior to Expiry%n", " ");
        System.out.printf("%5sNumber of equally spaced spline knots:%4d%n",
                " ", nxGrid);
        System.out.printf("%5sNumber of unknowns:%4d\n", " ", n);
        System.out.printf(Locale.ENGLISH,
                "%5sStrike=%5.2f, Sigma=%5.2f, Interest Rate=" + "%5.2f%n", " ",
                strikePrice, sigma, r);
        System.out.printf(Locale.ENGLISH,
                "%5sBet=%5.2f, Spread Value=%5.2f%n%n",
                " ", bet, spreadValue);
        System.out.printf("%17s%18s%18s%n", "Underlying", "A Bet",
                "Vertical Spread");
        for (i = 0; i < nv; i++) {
            System.out.printf(Locale.ENGLISH,
                    "%8s%9.4f%9.4f%9.4f%9.4f%9.4f%n", " ",
                    evaluationPoints[i],
                    splineValuesBetOption[0][i],
                    splineValuesBetOption[1][i],
                    splineValuesSpreadOption[0][i],
                    splineValuesSpreadOption[1][i]);
        }
    }

    // These classes define the coefficients, payoff, boundary conditions
    // and forcing term for American and European Options.
    static class MyCoefficients implements FeynmanKac.PdeCoefficients {

        final double zero = 0.0;
        double sigma, strikePrice, interestRate;
        double spread, bet, dividend;
        int dataInt;
        double value = 0.0;

        public MyCoefficients(double[] rData) {
            this.strikePrice = rData[0];
            this.bet = rData[1];
            this.spread = rData[2];
            this.sigma = rData[4];
            this.interestRate = rData[5];
            this.dividend = rData[6];
        }

        // The coefficient sigma(x)
        public double sigma(double x, double t) {
            return (sigma * x);
        }

        // The coefficient derivative d(sigma) / dx
        public double sigmaPrime(double x, double t) {
            return sigma;
        }

        // 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 strikePrice, spread, bet, interestRate, df;
        private int dataInt;

        public MyBoundaries(double[] rData, int iData) {
            this.strikePrice = rData[0];
            this.bet = rData[1];
            this.spread = rData[2];
            this.interestRate = rData[5];
            this.dataInt = iData;
        }

        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) {
            // This is the discount factor using the risk-free
            // interest rate
            df = Math.exp(interestRate * time);
            // Use flag passed to decide on boundary condition
            switch (dataInt) {
                case 1:
                    bndCoeffs[0][0] = 1.0;
                    bndCoeffs[0][1] = 0.0;
                    bndCoeffs[0][2] = 0.0;
                    bndCoeffs[0][3] = bet * df;
                    break;
                case 2:
                    bndCoeffs[0][0] = 1.0;
                    bndCoeffs[0][1] = 0.0;
                    bndCoeffs[0][2] = 0.0;
                    bndCoeffs[0][3] = (spread - strikePrice) * df;
                    break;
            }
            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 double terminal(double x) {
            final double zero = 0.0;
            double value = 0.0;
            switch (dataInt) {
                // The payoff function - Use flag passed to decide which
                case 1:
                    // After reaching the strike price the payoff jumps
                    // from zero to the bet value.
                    value = zero;
                    if (x > strikePrice) {
                        value = bet;
                    }
                    break;
                case 2:
                    /* Function is zero up to strike price.
                     Then linear between strike price and spread.
                     Then has constant value Spread-Strike Price after
                     the value Spread. */
                    value = Math.max(x - strikePrice, zero)
                            - Math.max(x - spread, zero);
                    break;
            }
            return value;
        }
    }
}

Output

  European Option Value for A Bet
   and a Vertical Spread, 3 and 6 Months Prior to Expiry
     Number of equally spaced spline knots:  61
     Number of unknowns: 183
     Strike=10.00, Sigma= 0.40, Interest Rate= 0.10
     Bet= 2.00, Spread Value=15.00

       Underlying             A Bet   Vertical Spread
           2.0000  -0.0000   0.0000  -0.0000   0.0000
           4.0000   0.0000   0.0013   0.0000   0.0005
           6.0000   0.0112   0.0729   0.0038   0.0452
           8.0000   0.2686   0.4291   0.1486   0.3833
          10.0000   0.9948   0.9781   0.8898   1.1907
          12.0000   1.6103   1.4301   2.1904   2.2267
          14.0000   1.8650   1.6926   3.4267   3.1567
          16.0000   1.9335   1.8171   4.2274   3.8282
          18.0000   1.9477   1.8696   4.6261   4.2499
          20.0000   1.9501   1.8902   4.7903   4.4918
          22.0000   1.9505   1.8979   4.8493   4.6222
          24.0000   1.9506   1.9008   4.8685   4.6901
Link to Java source.