package com.imsl.test.example.math;

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

/**
 * <p>
 * Evaluates the price of a European option with two payoff strategies.</p>
 * <p>
 * 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:
 * <ol>
 * <li> Strike and bet prices \(K_1 = {10.0},\; K_2={15.0}, and B={2.0}.\)</li>
 * <li> Volatility \( \sigma = {0.4}\)</li>
 * <li> Times until expiration = {1/4, 1/2}</li>
 * <li> Interest rate \(r = 0.1\)</li>
 * <li> \( x_{\min}=0.0,\, x_{\max}=30.0\)</li>
 * <li>\( nx = 61, \, n=3 \times nx = 183\)</li>
 * </ol>
 * </p>
 * <p>
 * </P>
 *
 * @see <a href="FeynmanKacEx3.java">Code</a>
 * @see <a href="doc-files/FeynmanKacEx3.out">Output</a>
 *
 */
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;
        }
    }
}