package com.imsl.test.example.math;

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

/**
 * <p>
 * Evaluates the price of a convertible bond.</p>
 *
 * <p>
 * This example evaluates the price of a convertible bond. Here, convertibility
 * means that the bond may, at any time of the bond holder's choosing, be
 * converted to a multiple of the specified asset. Thus a convertible bond with
 * price \(x\) returns an amount \(K\) at time \(T\) \(unless\) the owner has
 * converted the bond to \(\nu x,\, \nu \ge 1\) units of the asset at some time
 * prior to \(T\). This definition, the differential equation and boundary
 * conditions are given in Chapter 18 of Wilmott et al. (1996). Using a constant
 * interest rate and volatility factor, the parameters and boundary conditions
 * are:
 * <ol>
 * <li>Bond face value \(K = {1}\) , conversion factor \( \nu = 1.125\)</li>
 * <li>Volatility \( \sigma = {0.25}\)</li>
 * <li>Times until expiration = {1/2, 1}</li>
 * <li>Interest rate \(r\) = 0.1, dividend \(D\) = 0.02</li>
 * <li>\(x_{\min} = 0,\, x_{\max} = 4\)</li>
 * <li>\(nx = 61,\, n = 3 \times nx = 183\)</li>
 * <li>Boundary conditions \(f(0,t) = K \exp(r-(T-t)),\, f(x_{\max},t)=\nu
 * x_{\max}\)</li>
 * <li>Terminal data \(f(x,T)=\max(K,\nu x)\)</li>
 * <li>Constraint for bond holder \(f(x,t) \ge \nu x \)</li>
 * </ol>
 * </p>
 * <p>
 * Note that the error tolerance is set to a pure absolute error of value
 * \(10^{-3}\). The free boundary constraint \(f(x,t) \ge \nu x\) is achieved by
 * use of a non-linear forcing term in interface <code>ForcingTerm</code>. The
 * coefficient values of the Hermite quintic spline representing the approximate
 * solution of the differential equation at the initial time point are provided
 * with the interface <code>InitialData</code>.
 * </p>
 *
 *
 * @see <a href="FeynmanKacEx4.java">Code</a>
 * @see <a href="doc-files/FeynmanKacEx4.out">Output</a>
 *
 */
public class FeynmanKacEx4 {

    public static void main(String args[]) throws Exception {
        int i;
        int nxGrid = 61;
        int ntGrid = 2;
        int nv = 13;

        // Compute value of a Convertible Bond
        // The face value
        double KS = 1.0;
        // The sigma or volatility value
        double sigma = 0.25;
        // Time values for the options
        double[] tGrid = {0.5, 1.0};
        // Values of the underlying where evaluation are made
        double[] evaluationPoints = new double[nv];
        // Value of the interest rate, continuous dividend and factor
        double r = 0.1, dividend = 0.02, factor = 1.125;
        // Values of the min and max underlying values modeled
        double xMin = 0.0, xMax = 4.0;
        // Define parameters for the integration step.
        int nint = nxGrid - 1, n = 3 * nxGrid;
        double[] xGrid = new double[nxGrid];
        // Array for user-defined data
        double[] rData = new double[8];
        double dx, atol;
        // Number of left/right boundary conditions.
        int nlbcd = 3, nrbcd = 3;
        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] = (i - 1) * 0.25;
        }
        // Pass the data for evaluation
        rData[0] = KS;
        rData[1] = xMax;
        rData[2] = sigma;
        rData[3] = r;
        rData[4] = dividend;
        rData[5] = factor;
        // Use a pure absolute error tolerance for the integration
        atol = 1.0e-3;
        rData[6] = atol;
        // Compute value of convertible bond
        FeynmanKac convertibleBond = new FeynmanKac(new MyCoefficients(rData));

        MyForcingTerm forceTerm = new MyForcingTerm(rData);
        MyInitialData initialData = new MyInitialData(rData);

        convertibleBond.setForcingTerm(forceTerm);
        convertibleBond.setInitialData(initialData);
        //convertibleBond.setErrorTolerances(1.0e-3,0.0);
        convertibleBond.setAbsoluteErrorTolerances(1.0e-3);
        convertibleBond.setRelativeErrorTolerances(0.0);

        // Only the left boundary conditions are time dependent
        timeDependence[4] = true;
        convertibleBond.setTimeDependence(timeDependence);

        convertibleBond.computeCoefficients(nlbcd, nrbcd,
                new MyBoundaries(rData), xGrid, tGrid);
        double[][] bondCoefficients = convertibleBond.getSplineCoefficients();

        double[][] bondSplineValues = new double[ntGrid + 1][];

        /*
         * Evaluate and display solutions at vector of points XS(:),
         * at each time value prior to expiration.
         */
        for (i = 0; i <= ntGrid; i++) {
            bondSplineValues[i]
                    = convertibleBond.getSplineValue(evaluationPoints,
                            bondCoefficients[i], 0);
        }

        System.out.printf("%2sConvertible Bond Value, 0+, 6 and 12 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%n",
                " ", KS, sigma);
        System.out.printf(Locale.ENGLISH, "%5sInterest Rate=%5.2f, "
                + "Dividend=%5.2f, Factor=%6.3f%n%n",
                " ", r, dividend, factor);
        System.out.printf("%15s%18s%n", "Underlying", "Bond Value");
        for (i = 0; i < nv; i++) {
            System.out.printf(Locale.ENGLISH, "%7s%8.4f%8.4f%8.4f%8.4f%n",
                    " ", evaluationPoints[i],
                    bondSplineValues[0][i],
                    bondSplineValues[1][i],
                    bondSplineValues[2][i]);
        }
    }

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

        private double sigma, strikePrice, interestRate;
        private double dividend, factor, zero = 0.0;
        private double value = 0.0;

        public MyCoefficients(double[] rData) {
            this.strikePrice = rData[0];
            this.sigma = rData[2];
            this.interestRate = rData[3];
            this.dividend = rData[4];
            this.factor = rData[5];
        }

        // 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 interestRate, strikePrice, dp, factor, xMax;

        public MyBoundaries(double[] rData) {
            this.strikePrice = rData[0];
            this.xMax = rData[1];
            this.interestRate = rData[3];
            this.factor = rData[5];
        }

        public void leftBoundaries(double time, double[][] bndCoeffs) {
            dp = strikePrice * Math.exp(time * interestRate);
            bndCoeffs[0][0] = 1.0;
            bndCoeffs[0][1] = 0.0;
            bndCoeffs[0][2] = 0.0;
            bndCoeffs[0][3] = dp;
            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) {
            bndCoeffs[0][0] = 1.0;
            bndCoeffs[0][1] = 0.0;
            bndCoeffs[0][2] = 0.0;
            bndCoeffs[0][3] = factor * xMax;
            bndCoeffs[1][0] = 0.0;
            bndCoeffs[1][1] = 1.0;
            bndCoeffs[1][2] = 0.0;
            bndCoeffs[1][3] = factor;
            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) {
            // The payoff function
            double value = Math.max(factor * x, strikePrice);
            return value;
        }
    }

    static class MyForcingTerm implements FeynmanKac.ForcingTerm {

        private final double zero = 0.0, one = 1.0;
        private double value, strikePrice, interestRate;
        private double rt, mu, factor;

        public MyForcingTerm(double[] rData) {
            this.value = rData[6];
            this.strikePrice = rData[0];
            this.interestRate = rData[3];
            this.factor = rData[5];
        }

        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;
            int ndeg = xlocal.length;
            double[] yl = new double[6];
            double[] bf = new double[6];
            for (int i = 0; i < local; i++) {
                yl[i] = y[3 * interval - 3 + i];
                phi[i] = zero;
            }
            for (int i = 0; i < local; i++) {
                for (int j = 0; j < local; j++) {
                    dphi[i][j] = zero;
                }
            }

            mu = 2.0;
            /*
             * This is the local definition of the forcing term -
             * It "forces" the constraint f >= factor*x.
             */
            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 - factor * 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 * factor * strikePrice;
            }
            /*
             * This is the local derivative matrix for the forcing term
             */
            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 - factor * xlocal[l - 1]);
                        dphi[j - 1][i - 1] += qw[l - 1] * bf[i - 1]
                                * bf[j - 1] * Math.pow(value * rt, mu) * rt;
                    }
                }
            }
            for (int i = 0; i < local; i++) {
                for (int j = 0; j < local; j++) {
                    dphi[i][j] = -mu * dphi[i][j] * width
                            * factor * strikePrice;
                }
            }
        }
    }

    static class MyInitialData implements FeynmanKac.InitialData {

        private double[] data;

        public MyInitialData(double[] rData) {
            data = new double[rData.length];
            for (int i = 0; i < rData.length; i++) {
                data[i] = rData[i];
            }
        }

        public void init(double[] xGrid, double[] tGrid, double tp,
                double[] yprime, double[] y, double[] atol, double[] rtol) {
            int nxGrid = xGrid.length;
            if (tp == 0.0) { // Set initial data precisely.
                for (int i = 1; i <= nxGrid; i++) {
                    if (xGrid[i - 1] * data[5] < data[0]) {
                        y[3 * i - 3] = data[0];
                        y[3 * i - 2] = 0.0;
                        y[3 * i - 1] = 0.0;
                    } else {
                        y[3 * i - 3] = xGrid[i - 1] * data[5];
                        y[3 * i - 2] = data[5];
                        y[3 * i - 1] = 0.0;
                    }
                }
            }
        }
    }
}