package com.imsl.test.example.math;

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

/**
 * <p>
 * Applies a diffusion model for options pricing.</p>
 * <p>
 * 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 <i>K</i>, 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:
 * <ol>
 * <li>power \(\alpha = {2.0, 1.0, 0.0}\)</li>
 * <li>strike price \( K = {15.0, 20.0, 25.0} \)</li>
 * <li>volatility \(\sigma = {0.2, 0.3, 0.4}\)</li>
 * <li>times until expiration = {1/12, 4 /12, 7 /12}</li>
 * <li>underlying prices = {19.0, 20.0, 21.0}</li>
 * <li>interest rate \(r = 0.05\)</li>
 * <li> \( x_{\min} = 0, x_{\max} = 60 \) </li>
 * <li>\( nx = 121, \, n = 3 \times nx = 363\)</li>
 * </ol>
 * </p>
 * <p>
 * With this model the Feynman-Kac differential equation $$ f_t +\mu(x,t) f_x
 * +\frac{\sigma^2(x,t)}{2}f_{xx}-\kappa(x,t)f=\phi(f,x,t), $$ is defined by
 * identifying:
 * <ul>
 * <li>\(x:\; S\)</li>
 * <li>\(\sigma(x,t):\; \sigma x^{\alpha/2};\; \frac{\partial \sigma}{\partial
 * x}= \frac{\alpha \sigma}{2}x^{\alpha/2-1} \)</li>
 * <li>\(\mu(x,t):\;rx\)</li>
 * <li>\(\kappa(x,t):\; r\)</li>
 * <li>\(\phi(f,x,t) \equiv 0\)</li>
 * </ul>
 * </p>
 * <p>
 * The payoff function is the "vanilla option", \(p(x) = \max(x-K, 0)\).
 * </p>
 *
 *
 * @see <a href="FeynmanKacEx2.java">Code</a>
 * @see <a href="doc-files/FeynmanKacEx2.out">Output</a>
 *
 */
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 methods.
                    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;
        }

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

        public double terminal(double x) {
            double zero = 0.0;
            // The payoff function
            double value = Math.max(x - strikePrice, zero);
            return value;
        }
    }
}