package com.imsl.test.example.math;

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

/**
 * <p>
 * Compares American vs European options on a vanilla put.
 * </p>
 * <p>
 * The value of the American Option on a Vanilla Put can be no smaller than its
 * European counterpart. That is due to the American Option providing the
 * opportunity to exercise at any time prior to expiration. This example
 * compares this difference - or premium value of the American Option - at two
 * time values using the Black-Scholes model. The example is based on Wilmott et
 * al. (1996, p. 176), and uses the non-linear forcing or weighting term
 * described in Hanson, R. (2008), for evaluating the price of the American
 * Option. The coefficients, payoff, boundary conditions and forcing term for
 * American and European options are defined through interfaces
 * <code>PdeCoefficients</code>, <code>Boundaries</code> and
 * <code>ForcingTerm</code>, respectively. One breakpoint is set exactly at the
 * strike price. The sets of parameters in the computation are:
 * <ol>
 * <li> Strike price \(K = {10.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>
 * The payoff function is the "vanilla option", \(p(x) = \max(K-x,0)\) .
 * </p>
 *
 *
 * @see <a href="FeynmanKacEx1.java">Code</a>
 * @see <a href="doc-files/FeynmanKacEx1.out">Output</a>
 *
 */
public class FeynmanKacEx1 {

    public static void main(String args[]) throws Exception {
        // Compute American Option Premium for Vanilla Put
        // The strike price
        double KS = 10.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 = {0.0, 2.0, 4.0, 6.0, 8.0, 10.0,
            12.0, 14.0, 16.0};
        // Value of the interest rate
        double r = 0.1;
        // Values of the min and max underlying values modeled
        double xMin = 0.0, xMax = 30.0;
        // Define parameters for the integration step.
        int nxGrid = 61;
        int nv = 9;
        int nint = nxGrid - 1, n = 3 * nxGrid;
        double[] xGrid = new double[nxGrid];
        double dx;
        int nlbcd = 2, nrbcd = 3;
        // Define an 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;
        }

        FeynmanKac.Boundaries boundaries = new FeynmanKac.Boundaries() {

            public void leftBoundaries(double time, double[][] bndCoeffs) {
                bndCoeffs[0][0] = 0.0;
                bndCoeffs[0][1] = 1.0;
                bndCoeffs[0][2] = 0.0;
                bndCoeffs[0][3] = -1.0;
                bndCoeffs[1][0] = 0.0;
                bndCoeffs[1][1] = 0.0;
                bndCoeffs[1][2] = 1.0;
                bndCoeffs[1][3] = 0.0;
            }

            @Override
            public void rightBoundaries(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 double terminal(double x) {
                double zero = 0.0;
                // Strike price
                double strikePrice = 10.0;
                // The payoff function
                double value = Math.max(strikePrice - x, zero);

                return value;
            }
        };

        FeynmanKac.PdeCoefficients pdeCoefficients
                = new FeynmanKac.PdeCoefficients() {
            // The coefficient sigma(x)

            public double sigma(double x, double t) {
                double sigma = 0.4;
                return (sigma * x);
            }

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

            // The coefficient mu(x)
            public double mu(double x, double t) {
                double interestRate = 0.1;
                double dividend = 0.0;
                return ((interestRate - dividend) * x);
            }

            // The coefficient kappa(x)
            public double kappa(double x, double t) {
                double interestRate = 0.1;
                return interestRate;
            }
        };

        FeynmanKac.ForcingTerm forceTerm = new FeynmanKac.ForcingTerm() {

            public void force(int interval, double[] y, double time,
                    double width, double[] xlocal, double[] qw,
                    double[][] u, double[] phi, double[][] dphi) {
                final int local = 6;
                final double zero = 0.0, one = 1.0;
                double[] yl = new double[local];
                double[] bf = new double[local];
                double value, strikePrice, interestRate;
                double rt, mu;
                int ndeg = xlocal.length;

                for (int i = 0; i < local; i++) {
                    yl[i] = y[3 * interval - 3 + i];
                    phi[i] = zero;
                }
                strikePrice = 10.0;
                interestRate = 0.1;
                value = 1.0e-5;
                mu = 2.0;
                // This is the local definition of the forcing term
                for (int j = 1; j <= local; j++) {
                    for (int l = 1; l <= ndeg; l++) {
                        bf[0] = u[0][l - 1];
                        bf[1] = u[1][l - 1];
                        bf[2] = u[2][l - 1];
                        bf[3] = u[6][l - 1];
                        bf[4] = u[7][l - 1];
                        bf[5] = u[8][l - 1];
                        rt = 0.0;
                        for (int k = 0; k < local; k++) {
                            rt += yl[k] * bf[k];
                        }
                        rt = value / (rt + value
                                - (strikePrice - xlocal[l - 1]));
                        phi[j - 1] += qw[l - 1] * bf[j - 1] * Math.pow(rt, mu);
                    }
                }
                for (int i = 0; i < local; i++) {
                    phi[i] = -phi[i] * width * interestRate * strikePrice;
                }
                // This is the local derivative matrix for the forcing term
                for (int i = 0; i < local; i++) {
                    for (int j = 0; j < local; j++) {
                        dphi[i][j] = zero;
                    }
                }
                for (int j = 1; j <= local; j++) {
                    for (int i = 1; i <= local; i++) {
                        for (int l = 1; l <= ndeg; l++) {
                            bf[0] = u[0][l - 1];
                            bf[1] = u[1][l - 1];
                            bf[2] = u[2][l - 1];
                            bf[3] = u[6][l - 1];
                            bf[4] = u[7][l - 1];
                            bf[5] = u[8][l - 1];
                            rt = 0.0;
                            for (int k = 0; k < local; k++) {
                                rt += yl[k] * bf[k];
                            }
                            rt = one / (rt + value
                                    - (strikePrice - xlocal[l - 1]));
                            dphi[j - 1][i - 1] += qw[l - 1] * bf[i - 1]
                                    * bf[j - 1] * Math.pow(rt, mu + 1.0);
                        }
                    }
                }
                for (int i = 0; i < local; i++) {
                    for (int j = 0; j < local; j++) {
                        dphi[i][j] = mu * dphi[i][j] * width
                                * Math.pow(value, mu) * interestRate
                                * strikePrice;
                    }
                }
            }
        };

        FeynmanKac european = new FeynmanKac(pdeCoefficients);
        FeynmanKac american = new FeynmanKac(pdeCoefficients);

        american.setForcingTerm(forceTerm);

        european.computeCoefficients(nlbcd, nrbcd, boundaries, xGrid, tGrid);
        american.computeCoefficients(nlbcd, nrbcd, boundaries, xGrid, tGrid);

        // Evaluate solutions at vector of points evaluationPoints, at each
        // time value prior to expiration.
        double[][] europeanCoefficients = european.getSplineCoefficients();
        double[][] americanCoefficients = american.getSplineCoefficients();

        double[][] splineValuesEuropean = new double[nt][];
        double[][] splineValuesAmerican = new double[nt][];

        for (int i = 0; i < nt; i++) {
            splineValuesEuropean[i]
                    = european.getSplineValue(evaluationPoints,
                            europeanCoefficients[i + 1], 0);
            splineValuesAmerican[i]
                    = american.getSplineValue(evaluationPoints,
                            americanCoefficients[i + 1], 0);
        }

        System.out.printf("%nAmerican Option Premium for Vanilla Put, "
                + "3 and 6 Months Prior to Expiry%n");
        System.out.printf("%7sNumber of equally spaced spline knots:"
                + "%4d%n", " ", nxGrid);
        System.out.printf("%7sNumber of unknowns:%4d%n", " ", n);
        System.out.printf(Locale.ENGLISH, "%7sStrike=%6.2f, sigma=%5.2f,"
                + " Interest Rate=%5.2f%n%n", " ", KS, sigma, r);
        System.out.printf("%7s%10s%20s%20s%n", " ", "Underlying",
                "European", "American");
        for (int i = 0; i < nv; i++) {
            System.out.printf(Locale.ENGLISH, "%7s%10.4f%10.4f%10.4f"
                    + "%10.4f%10.4f%n", " ",
                    evaluationPoints[i],
                    splineValuesEuropean[0][i],
                    splineValuesEuropean[1][i],
                    splineValuesAmerican[0][i],
                    splineValuesAmerican[1][i]);
        }
    }
}