function [portStd, portRet, portWts, portIndx] = computebestportfolioPCT(expRet,expCov,portSize,targetRet)
% Oren Rosen
% 4/14/2009
% Copyright 2008 The MathWorks, Inc.
%
% This function implements a custom portfolio optimization routine:
%
% Find the weights that minimize portfolio standard deviation with
% constraints that require a given target portfolio return and a maximum
% number of non-zero weights (cardinality constraint).
%
% Inputs:
%
% expRet - Column vector of expected returns of equity universe.
% expCov - Square matrix of expected covariance of equity universe.
% portSize - Maximum number of non-zero weights. Less than the number of
% equities in the universe.
% targetRet - Target return for complete portfolio.
%
% Outputs:
%
% portStd - Standard deviation of optimal portfolio.
% portRet - Total return of optimal portfolio (should match targetRet).
% portWts - Vector of optimal weights. Length equal to universe size.
% PortIndx - Vector of integers dictating which equities are included in
% the portfolio.
%
% The implementation of this routine is done with nested optimization
% routines. The outer layer uses a custom genetic algorithm based on a bit
% string data type. The bit string has length equal to the number of
% equities in the universe, and a 1 (or 0) represents that an equity is
% included (not included) in the current portfolio. The portfolio
% cardinality constraint is built into the genetic algorithm evolution
% functions:
%
% generateinitpopNcK(...)
% crossoverNcK(...)
% mutationNcK(...)
%
% The genetic algorithm determines which equities are included in the best
% portfolio. Within the fitness function, the quadratic programming solver
% "quadprog" is used to find the optimal weights for the current choice of
% equities.
%
% NOTE: This version is meant to be used with a demonstration of the GA/PCT
% features. Note the addition of the last 2 options in the calls below to
% "gaoptimset". It is identical in all other ways to the version in the
% MLScript directory.
% Generate initial population for GA
univSize = length(expRet);
iPop = generateinitpopNcK(univSize,portSize);
% Start with default options
options = gaoptimset;
% Set population options
options = gaoptimset(options,'PopulationSize' ,univSize-portSize+1 );
options = gaoptimset(options,'InitialPopulation' , iPop);
% Set evolution options
options = gaoptimset(options,'CrossoverFraction' ,0.9 );
options = gaoptimset(options,'CrossoverFcn' ,@crossoverNcK);
options = gaoptimset(options,'MutationFcn' ,@mutationNcK);
options = gaoptimset(options,'Vectorize' ,'off');
% Set diaplay options
options = gaoptimset(options,'Display' ,'off');
% *** Set options for PCT Demo ***
options = gaoptimset(options,'UseParallel','always');
options = gaoptimset(options,'Generations',4); % Limit generations
% Maximum portfolio standard deviation is the largest individual stock
% standard deviation. This will be used in the GA fitness function below.
maxPortStd = sqrt(max(diag(expCov)));
% Run GA
[portBits,fval,gaExit] = ga(@rankport,univSize,[],[],[],[],[],[],[],options);
% Output of GA is logical vector of the best equities to
% include in the portfolio. Use this data to run quadprog
% one more time to insure that we get the correct weights.
% Extract subsets of expected return and covariance
portBits = logical(portBits);
subCov = expCov(portBits,portBits);
subRet = expRet(portBits);
% Calculate optimal weights, compute portfolio statistics
[subWts,subVar,mvExit] = minvar(subCov,subRet,targetRet);
portRet = subWts'*subRet;
portStd = sqrt(subVar);
% Build full weight vector and index vector for optimal portfolio
indices = 1:univSize;
portIndx = indices(portBits);
portWts = zeros(univSize,1);
portWts(portIndx) = subWts;
% *****************************************************************%
% Fitness function for GA
%
% The input "x" is a bit string that represents which stocks to include
% in the current portfolio. The fitness of this portfolio is measured
% as the minimum portfolio standard deviation possible given the target
% return constraint and the current choice of equities.
function risk = rankport(x)
% Use the current bit-string "x" to extract subsets of the equities
% expected returns and covariance.
x = logical(x);
subCov = expCov(x,x);
subRet = expRet(x);
% "minvar" is a user defined m-function. It is essentially a
% wrapper for "quadprog" and calculates the weights for the
% current portfolio that minimizes portfolio variance. If a
% different measure of risk is desired, replace this function.
[weights,fval,mvExit] = minvar(subCov,subRet,targetRet);
% Risk is measure in this routine as portfolio standard deviation.
% If the optimization routine in "minvar" failed, set risk to a
% high valuec so that this choice of portfolio is flagged as
% un-optimal. This could happen simply because the expected returns
% for the current choice of stocks are inadequete to achieve the
% target return.
if( mvExit == 1 )
risk = sqrt(fval);
else
risk = 10*maxPortStd;
end
end
% *****************************************************************%
end