/*
 * -------------------------------------------------------------------------
 *      $Id: EfficientFrontier.java,v 1.4 2005/03/11 16:48:20 estewart Exp $
 * -------------------------------------------------------------------------
 *      Copyright (c) 2004 Visual Numerics Inc. All Rights Reserved.
 *
 *      This software is confidential information which is proprietary to
 *      and a trade secret of Visual Numerics, Inc.  Use, duplication or
 *      disclosure is subject to the terms of an appropriate license
 *      agreement.
 *
 *      VISUAL NUMERICS MAKES NO REPRESENTATIONS OR WARRANTIES ABOUT THE
 *      SUITABILITY OF THE SOFTWARE, EITHER EXPRESS OR IMPLIED, INCLUDING
 *      BUT NOT LIMITED TO THE IMPLIED WARRANTIES OF MERCHANTABILITY,
 *      FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. VISUAL
 *      NUMERICS SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED BY LICENSEE
 *      AS A RESULT OF USING, MODIFYING OR DISTRIBUTING THIS SOFTWARE OR
 *      ITS DERIVATIVES.
 *--------------------------------------------------------------------------
 */

package com.imsl.demo.portfolio;
import com.imsl.stat.Covariances;
import com.imsl.math.LinearProgramming;
import com.imsl.math.QuadraticProgramming;

/** Compute an Efficient Frontier. User may input a set of investment time series
 * and various constraints on the portfolio content. Expected returns and a
 * variance-covariance matrix may be input instead of the historical time series.
 *
 * @author  ed
 * @created February 27, 2004
 */
public class EfficientFrontier {

    private double[][] totalReturns, fractionalReturns, varcov, cw, bw, solution;
    private double[] expectedReturns;
    private int m;
    private final int numPortfolios = 100;
    private boolean haveTotalReturns;
    private boolean haveExpectedReturns = false, haveVarCov = false;

    /** Creates a new instance of EfficientFrontier and sets the default constraints.
     *
     * @param m_ The number of elements in the portfolio
     */
    public EfficientFrontier(int m_) {
        this.m = m_;

        // set up the default linear constraints on the weights
        // for 3 input time series, this matrix looks like:
        // 1  0  0
        // 0  1  0
        // 0  0  1
        // 1  1  1
        cw = new double[m+1][m];
        for (int i=0; i<m; i++) {
            cw[i][i] = 1;
            cw[m][i] = 1;
        }

        // set up the default bounds on the weights
        // for 3 input time series, this matrix looks like:
        // 0  1
        // 0  1
        // 0  1
        // 1  1
        // indicating the weights must be between 0 and 1 for each
        // and that the total of the weights must be 1
        bw = new double[m+1][2];
        for (int i=0; i<m+1; i++) {
            bw[i][1] = 1;
        }
        bw[m][0] = 1;
    }

    /** Compute the Efficient Frontier
     * @throws IMSLException If there is a problem with LinearProgramming or
     * QuadraticProgramming, this will be thrown.
     */
    public void compute() throws com.imsl.IMSLException {
        if (haveTotalReturns) {
            computeFractionalReturns();
        }

        if (!haveExpectedReturns) {
            computeExpectedReturns();
        }

        if (!haveVarCov) {
            Covariances cov = new Covariances(com.imsl.math.Matrix.transpose(fractionalReturns));
            varcov = cov.compute(Covariances.VARIANCE_COVARIANCE_MATRIX);
        }

        // prepare inputs into LinearProgramming
        int[] cType = new int[m+1];
        double[] b = new double[m+1];
        double[] upperLimit = new double[m+1];
        double[] negativeExpectedReturns = new double[m];
        for (int i=0; i<m+1; i++) {
            cType[i] = 3;
            b[i] = bw[i][0];
            upperLimit[i] = bw[i][1];
            if (i<m) {
                negativeExpectedReturns[i] = -1*expectedReturns[i];
            }
        }

        //f = linprog( c(*,0:p-1), c(*,p), x(*), bu=c(*,p+1), irt=0*c(*,p)+3, obj=y )
        LinearProgramming lp = new LinearProgramming(cw, b, expectedReturns);
        lp.setUpperLimit(upperLimit);
        lp.setConstraintType(cType);
        lp.solve();
        double y = lp.getOptimalValue();

        //f = linprog( c(*,0:p-1), c(*,p), -x(*), bu=c(*,p+1), irt=0*c(*,p)+3, obj=z )
        lp = new LinearProgramming(cw, b, negativeExpectedReturns);
        lp.setUpperLimit(upperLimit);
        lp.setConstraintType(cType);
        lp.solve();
        double z = lp.getOptimalValue();

        //f = rebin( interpol([y+1e-4,-z-1e-4],n), n, p+2 )
        double y0 = y+1e-4;
        double z0 = -1*z-1e-4;
        double step = (z0-y0)/(numPortfolios-1.0);
        double[] frontier = new double[numPortfolios];
        for (int i=0; i<numPortfolios; i++) {
            frontier[i] = i*step+y0;
        }

        // now build the inputs into QuadraticProgramming
        //b = [ x, -x, c(*,0:p-1), -c(*,0:p-1) ]
        double[][] qa = new double[2*(m+1)+2][m];
        for (int i=0; i<m; i++) {
            qa[0][i] = expectedReturns[i];
            qa[1][i] = -1*expectedReturns[i];
            for (int j=0; j<m+1; j++) {
                qa[j+2][i] = cw[j][i];
                qa[j+m+3][i] = -1*cw[j][i];
            }
        }

        //d = [ 0, 0, c(*,p), -c(*,p+1) ]
        double[] qb = new double[2*(m+1)+2];
        // the first two elements are dynamic and change inside the quadprog loop
        for (int i=0; i<m+1; i++) {
            qb[i+2] = bw[i][0];
            qb[i+m+3] = -1*bw[i][1];
        }

        // now loop through and call QuadraticProgramming for each point
        //  for i = 0, n-1 do begin
        //    d(0:1) = [ f(i), -f(i) ]
        //    f(i,2:*) = quadprog( b, d, g, v, meq=0, obj=h )
        //    f(i) = sqrt( 2*h )
        //  endfor
        QuadraticProgramming qp = null;
        double[] qg = new double[m];
        solution = new double[numPortfolios][m+2];
        for (int i=0; i<numPortfolios; i++) {
            qb[0] = frontier[i];
            qb[1] = -1.0*frontier[i];

            try {
                qp = new QuadraticProgramming(varcov, qg, null, null, qa, qb);

                double[] distWeights = qp.getSolution();
                // optimal object function found
                // CNL has this as an optional output, but JMSL doesn't
                //      tmp = 0.0;
                //      for (i = 1; i <= n; i++) {
                //          k = (i - 1) * n;
                //          tmp += (*lv_x)[i - 1] * imsl_sdot (n, &h[k], 1, *lv_x, 1);
                //        }
                //      *obj = imsl_sdot (n, g, 1, *lv_x, 1) + .5 * tmp;

                double tmp=0.0;
                for (int j=0; j<m; j++) {
                    tmp += distWeights[j] * dot(m, varcov[j], distWeights);
                }
                double obj = 0.5*tmp + dot(m, qg, distWeights);

                solution[i][0] = Math.sqrt(2.0*obj);
                solution[i][1] = frontier[i];
                for (int j=0; j<m; j++) {
                    solution[i][j+2] = distWeights[j];
                }
            } catch (QuadraticProgramming.InconsistentSystemException ise){
                for (int j=0; j<m+2; j++) {
                    solution[i][j] = Double.NaN;
                }
            };
        }
    }

    // compute a dot product of two vectors
    private double dot(int n, double x[], double y[]) {
        double    sum = 0.0;
        for (int k = 0;  k < n;  k++) {
            sum += x[k]*y[k];
        }
        return sum;
    }

    // Convert total returns into fractional returns
    private void computeFractionalReturns() {
        fractionalReturns = new double[m][];
        for (int i=0; i<m; i++) {
            int p = totalReturns[i].length;
            fractionalReturns[i] = new double[p-1];
            for (int j=0; j<p-1; j++) {
                //u = (r(1:*,*)-r(0:m-2,*))/r(0:m-2,*)
                fractionalReturns[i][j] = (totalReturns[i][j+1]-totalReturns[i][j])/totalReturns[i][j];
            }
        }
    }

    //compute expected return from fractional returns, the geometric mean
    private void computeExpectedReturns() {
        //x = exp(total(alog(1+u),dimens=0)/m)-1
        expectedReturns = new double[m];
        double[] totalLog = new double[m];
        for (int i=0; i<m; i++) {
            for (int j=0; j<fractionalReturns[i].length; j++) {
                totalLog[i] += Math.log(1+fractionalReturns[i][j]);
            }
        }
        for (int i=0; i<m; i++) {
            expectedReturns[i] = Math.exp(totalLog[i]/fractionalReturns[i].length)-1;
        }
    }

    /** Set the total returns of the m stocks to consider in the portfolios.
     * These are converted to fractional returns internally.
     * @param r The total returns (i.e., stock prices), an [m][p] array.
     */
    public void setTotalReturns(double[][] r) {
        this.totalReturns = r;
        this.haveTotalReturns = true;
    }

    /** Set the fractional returns of the stocks to consider in the portfolios
     * @param r The fractional returns, an [m][p] array.
     */
    public void setFractionalReturns(double[][] r) {
        this.fractionalReturns = r;
        this.haveTotalReturns = false;
    }

    /** Set the constraints on the weights
     *  This isn't as general as it coupld be, and the last row is
     *  added on automatically.
     * @param b A [m][2] array of lower and upper bounds for the weights of
     * each of the m stocks in the portfolio.
     */
    public void setBounds(double[][] b) {
        bw = new double[m+1][2];
        for (int i=0; i<m; i++) {
            bw[i][0] = b[i][0];
            bw[i][1] = b[i][1];
        }
        bw[m][0] = 1;
        bw[m][1] = 1;
    }

    /** Set the coefficient matrix for the constraints
     *  This isn't as general as it coupld be, and the last row is
     *  added on automatically.
     * @param c A [m][m] array; defaults to the identity matrix.
     */
    public void setCoefficientMatrix(double[][] c) {
        cw = new double[m+1][m];
        for (int i=0; i<m; i++) {
            for (int j=0; j<m; j++) {
                cw[i][j] = c[i][j];
            }
            cw[m][i] = 1;
        }
    }

    /** Set the variance-covariance matrix. If a user also sets the
     *  expected returns, there is no need to input time series.
     * @param cov A [m][m] array of covariances for the m stocks. If time series
     * of returns are input, this is computed as:
     * Covariances cov = new Covariances(com.imsl.math.Matrix.transpose(fractionalReturns));
* varcov = cov.compute(Covariances.VARIANCE_COVARIANCE_MATRIX); */
public void setCovarianceMatrix(double[][] cov) { this.haveVarCov = true; this.varcov = cov; } /** Get the variance-covariance matrix computed from input time series. * @return c A [m][m] array of covariances for the m stocks. */ public double[][] getCovarianceMatrix() { return this.varcov; } /** Set the expected return of each component of the portfolio * @param e A vector of the expected returns for each of the m investments. */ public void setExpectedReturns(double[] e) { this.haveExpectedReturns = true; this.expectedReturns = e; } /** Get the expected return of each component of the portfolio, computed as the geometric * mean of the input time series. * @return e A vector of the expected returns for each of the m investments. */ public double[] getExpectedReturns() { return this.expectedReturns; } /** Return the solution, which contains the risk, return and portfolio ratios * @return The entire solution. Solution[0][*] are the values of risk; * Solution[1][*] are the values of return; Solution[2 to 2+m][*] are the * weights of each of the m investments that make up the portfolio. */ public double[][] getSolution() { return com.imsl.math.Matrix.transpose(this.solution); } }