/*
* -------------------------------------------------------------------------
* $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);
}
}