function [PRisk, PRoR, PWts] = frontera(ERet, ECov, NPts, RTarget, ACLUBounds, ACGroups, ACGLUBounds)
% FRONTCON  Mean-variance efficient frontier with portfolio constraints.
%   This function returns the mean-variance efficient frontier with user specified
%   asset constraints, covariance and returns. Given a collection of NASSETS risky 
%   assets, computes a portfolio of asset investment weights which minimize the
%   risk for given values of the expected return.  The portfolio risk is minimized
%   subject to constraints on the asset weights or on groups of asset weights.
%
%   [PortRisk, PortReturn, PortWts] = frontcon(ExpReturn, ExpCovariance, ...
%      NumPorts, PortReturn, AssetBounds, Groups, GroupBounds)
%
%   Inputs: 
%     ExpReturn is a 1xNASSETS vector specifying the expected (mean) 
%     return of each asset.
%    
%     ExpCovariance is an NASSETSxNASSETS matrix specifying the covariance of
%     the asset returns.
%    
%     NumPorts is the number of efficient portfolios, NPORTS.  The default value
%     is 10 if NumPorts is entered as the empty matrix [] or not specified. 
%     NumPorts should be entered as [] if PortReturn is to be entered.
%
%     PortReturn is a vector of length NPORTS containing the target return
%     values on the frontier.
%     You can specify the portfolio return values in two ways.  If PortReturn
%     is not entered or is empty, NumPorts equally spaced returns between the
%     minimum and maximum possible values will be used.  
%
%     AssetBounds is a 2xNASSETS matrix containing the lower and upper
%     bounds on the weight to be allocated in each asset in the portfolio.
%     The default lower bound is all zeros (no short-selling), the default 
%     upper bound is all ones (any asset may comprise the entire portfolio). 
%    
%     Groups is an NGROUPSxNASSETS matrix specifying NGROUPS asset  
%     groups or classes.  Each row specifies a group:  Groups(i,j) = 1 if 
%     the j'th asset belongs to the i'th group.  Groups(i,j) = 0 if the j'th
%     asset is not in the i'th group.
%    
%     GroupBounds is a NGROUPSx2 matrix specifying, for each group, the lower
%     and upper bounds of the total weights of all assets in that group.  The
%     default lower bound is all zeros, the default upper bound is all ones.   
%
%
%   Outputs: 
%     PortRisk is an NPORTSx1 vector of the standard deviation of return for
%     each portfolio. 
%    
%     PortReturn is an NPORTSx1 vector of the expected return of each portfolio.
%    
%     PortWts is an NPORTSxNASSETS matrix of weights allocated to each
%     asset. Each row represents a different portfolio. The total of all 
%     weights in a portfolio is 1.
%            
%
%   Notes: 
%     A plot of the efficient frontier is returned if the function is
%     invoked without output arguments.
%
%     The asset returns are assumed to be jointly normal, with expected mean
%     returns of ExpReturn and return covariance ExpCovariance.  The variance
%     of a portfolio with 1 by NASSET weights PortWts is given by
%     PortVar = PortWts*ExpCovariance*PortWts'.  The portfolio expected
%     return is PortReturn = dot(ExpReturn, PortWts).
%
if (nargin < 2)
     error('finance:frontcon:missingInputs','You must enter ExpReturn and ExpCovariance.');
end

% Make sure that the number of returns entered matches the number of rows/columns in the  
% covariance matrix (which represents the number of assets).
ERet = ERet(:);
NASSETS = length(ERet);
   
if any( size(ECov)~=NASSETS )
    error('finance:frontcon:nonSquareCovMatrix','The covariance matrix must be NxN, where N = number of assets');
end
   
EC = eig(ECov);   
if(min(min(EC)) < (-1E-14 * max(max(abs(EC)))))
    warning('Covariance matrix must be positive semi-definite.');
end
clear EC;
   
% Determine which optional arguments were entered as non-empty.
if (nargin < 3 | isempty(NPts)) 
   NPtsEntered = 0; 
   NPts = [];
else
   NPtsEntered = 1; 
   NPts = double(NPts);
end

if (nargin < 4 | isempty(RTarget))    
   RTargetEntered = 0; 
   RTarget = [];
else
   RTargetEntered = 1; 
end

%When entering the target rate of return, enter NPts as an empty matrix.
if (NPtsEntered ==1) & (RTargetEntered == 1)
  if (length(RTarget(:)) ~= NPts )
    error('finance:frontcon:needEmptyNumPorts','When entering PortReturn, enter NumPorts as an empty matrix')
  end
end

if (nargin < 5 | isempty(ACLUBounds))  
   ACLUBoundsEntered = 0; 
else
   ACLUBoundsEntered = 1; 
end

if (nargin < 6 | isempty(ACGroups))  
   ACGroupsEntered  = 0; 
else
   ACGroupsEntered  = 1; 
end

if (nargin < 7 | isempty(ACGLUBounds))  
   ACGLUBoundsEntered  = 0; 
else
   ACGLUBoundsEntered  = 1; 
end

if(ACGroupsEntered & ~ACGLUBoundsEntered)
   % Set default values for group bounds
   NGroups = size(ACGroups, 1);
	ACGLUBounds = [zeros(NGroups,1), ones(NGroups,1)];   
end


%----------------------------------------------------------------
% Generate Constraint array needed to call PORTOPT
%----------------------------------------------------------------

% ACLUBounds, ACGroups, ACGLUBounds

if(ACLUBoundsEntered)
   % Allow short-selling, ineqparse applies the most restrictive constraint!
   ConSet = pcpval(1, NASSETS); 
   ConSet = [ConSet; portcons('AssetLims', ACLUBounds(1,:), ACLUBounds(2,:))];
else
   ConSet = portcons('Default', NASSETS);
end

if(ACGroupsEntered)
   ConSet = [ConSet; portcons('GroupLims', ACGroups, ACGLUBounds(:,1), ACGLUBounds(:,2))];
end

%----------------------------------------------------------------
% Generate output
%----------------------------------------------------------------

if nargout == 0
   portopt(ERet, ECov, NPts, RTarget, ConSet);
else
   [PRisk, PRoR, PWts] = portopt(ERet, ECov, NPts, RTarget, ConSet);
end
function ConSet = portcons(varargin)
% PORTCONS Generates a matrix of constraints for a portfolio
%   of asset investments by using linear inequalities. The 
%   inequalities are of the type A*Wts' <= b, where Wts is the 
%   matrix of weights. The matrix ConSet is defined as 
%   ConSet = [A b].
%
%   ConSet = portcons('ConstType', Data1, ..., DataN) creates a matrix 
%   ConSet, based on the constraint type ConstType, and the constraint 
%   parameters Data1, ..., DataN. See below for a list of constraint 
%   types and their corresponding constraint parameters.
%   
%   ConSet = portcons('ConstType1', Data11, ..., Data1N, 'ConstType2', Data21, ..., Data2N, ...) 
%   creates a matrix ConSet, based on the on the constraint types ConstTypeN, 
%   and the corresponding constraint parameters DataN1, ..., DataNN. See 
%   below for a list of constraint types and their corresponding constraint 
%   parameters.
%
%   Constraint types and corresponding constraint data:
%
%     Constraint type: 'Default'
%     All allocations are >= 0; no short selling is allowed. The combined
%     value of the portfolio allocations is normalized to 1.
%
%     Constraint data: NumAssets (required)
%     NumAssets is a scalar representing the number of assets in the 
%     portfolio.
%
%
%     Constraint type: 'PortValue'
%     Require the total value of the portfolio to be fixed to PVal.
%
%     Constraint data: PVal (required), NumAssets(required)
%     PVal is a scalar representing the total value of the portfolio.
%     NumAssets is a scalar representing the number of assets in the portfolio.
%     See the function PCPVAL for more information.
%     
%
%     Constraint type: 'AssetLims'
%     Specify the minimum and maximum allocation per asset.
%
%     Constraint data: AssetMin (required), AssetMax (required), NumAssets (optional)
%     AssetMin is a scalar or vector of length NASSETS, specifying the  
%     minimum allocation per asset.
%     AssetMax is a scalar or vector of length NASSETS, specifying the  
%     maximum allocation per asset.
%     NumAssets is a scalar representing the number of assets in the portfolio.
%     See the function PCALIMS for more information.     
% 
%     
%     Constraint type: 'GroupLims'
%     Specify minimum and maximum allocations to groups of assets.
%
%     Constraint data: Groups (required), GroupMin (required), GroupMax (required)
%     Groups  is an NGROUPS by NASSETS matrix, specifying which assets 
%     belong to each group.
%     GroupMin is a scalar or a vector of length NGROUPS, specifying  
%     the minimum combined allocations in each group.
%     GroupMax is a scalar or a vector of length NGROUPS, specifying  
%     the maximum combined allocations in each group.
%     See the function PCGLIMS for more information.     
%
%
%     Constraint type: 'GroupComparison' 
%     Specify the asset group-to-group comparison constraints. 
%
%     Constraint data: GroupA (required), AtoBmin (required), AtoBmax (required), GroupB (required)
%     GroupA and GroupB are NGROUPS by NASSETS matrices, specifying the 
%     groups to compare.
%     AtoBmin is a scalar or a vector of length NGROUPS, specifying the  
%     minimum ratios of allocations in group A to allocations in group B.
%     AtoBmax is a scalar or a vector of length NGROUPS, specifying the  
%     maximum ratios of allocations in group A to allocations in group B.
%     See the function PCGCOMP for more information.
%
%
%     Constraint type: 'Custom'
%     Custom linear inequality constraints A*Wts' <= b.  
%
%     Constraint data: A (required), b (required)
%     A is an NCONSTRAINTS by NASSETS matrix, specifying the weights for 
%     each asset in each inequality equation.
%     b is a vector of length NCONSTRAINTS, specifying the right hand sides  
%     of the inequalities.
%
%
%   See also PORTOPT, PCPVAL, PCGLIMS, PCGCOMP, PCALIMS
%

%   Copyright 1995-2006 The MathWorks, Inc.
%   $Revision: 1.9.2.2 $   $ Date: 1998/01/30 13:45:34 $

% Find the types in the argument list
TypeIndex  = [];
for i=1:nargin
   if(ischar(varargin{i}))
      TypeIndex  = [TypeIndex i];  % Create a vector with the positions (indexes) of 
   end                             %the types within varargin
end

NumTypes = length(TypeIndex);
if(NumTypes == 0)
   error('finance:portcons:invalidInputList','No types found in the input list');
end


% TypeIndex contains the index (in varargin) to each type. Now, for each type, 
% check its value and call the appropriate function to generate its corresponding 
% ConSet. Then, concatenate as we move through the cell array.
ConSet=[];
for i=1:NumTypes
   
   % Find the number of data arguments passed for this type
   if(i == NumTypes)
      NumArgs = nargin - TypeIndex(i); % We are looking at the last type in varargin
   else
      NumArgs = TypeIndex(i+1) - TypeIndex(i) - 1;
   end
   
   switch(lower(varargin{TypeIndex(i)}))
   case 'default'
      % Make sure that one data argument is being passed in.
      if(NumArgs < 1)
         error('finance:portcons:missingDefaultTypeArg','Type "Default" requires at least one data argument.')
      end      
      
      PcpvalConSet = pcpval(1, varargin{TypeIndex(i)+1});
   	PcalimsConSet = pcalims(0, NaN, varargin{TypeIndex(i)+1});
   
      % Make sure that the number of columns is adequate for concatenation.
		if(isempty(ConSet))   
         ConSet = [PcpvalConSet; PcalimsConSet];               
      else
         if(size(ConSet,2) == size(PcpvalConSet,2) & size(ConSet,2) == size(PcalimsConSet,2))
            ConSet = [ConSet; PcpvalConSet; PcalimsConSet];
         else
            error('finance:portcons:invalidNumAssets','Number of assets implied in type "Default" is not consistent.');
         end   
      end         
      
      
   case 'assetlims'
      % Make sure that at least one data argument is passed in.
      if(NumArgs < 1)
         error('finance:portcons:missingAssetLimsArg','Type "AssetLims" requires at least one data argument.')
      end
      
      if(NumArgs == 1)
         PcalimsConSet = pcalims(varargin{TypeIndex(i)+1});
      elseif(NumArgs == 2)
         PcalimsConSet = pcalims(varargin{TypeIndex(i)+1:TypeIndex(i)+2});
      elseif (NumArgs == 3)
         PcalimsConSet = pcalims(varargin{TypeIndex(i)+1:TypeIndex(i)+3});
      else
         warning('More than three data arguments passed to type "AssetLims".');
      end
      
 		% Make sure that the number of columns is adequate for concatenation. 
      if(isempty(ConSet))   
         ConSet = PcalimsConSet;               
      else
         if(size(ConSet,2) == size(PcalimsConSet,2))
            ConSet = [ConSet; PcalimsConSet];
         else
            error('finance:portcons:invalidNumAssets','Number of assets implied in type "AssetLims" is not consistent.');
         end   
      end  
      
      
   case 'portvalue'
      % Make sure that at least two data arguments are passed in.
      if(NumArgs < 2)
         error('finance:portcons:missingPortValueArg','Type "PortValue" requires two data arguments');
      end
      
      PcpvalConSet = pcpval(varargin{TypeIndex(i)+1:TypeIndex(i)+2});
      
 		% Make sure that the number of columns is adequate for concatenation. 
      if(isempty(ConSet))   
         ConSet = PcpvalConSet;               
      else
         if(size(ConSet,2) == size(PcpvalConSet,2))
            ConSet = [ConSet; PcpvalConSet];
         else
            error('finance:portcons:invalidNumAssets','Number of assets implied in type "PortValue" is not consistent.');
         end   
      end  
      
      
	case 'grouplims'
      % Make sure that at least three data arguments are passed in
      if(NumArgs < 3)
         error('finance:portcons:missingGroupLimsArg','Type "GroupLims" requires three data arguments.');
      end  
      
      PcglimsConSet = pcglims(varargin{TypeIndex(i)+1:TypeIndex(i)+3});
      
     	% Make sure that the number of columns is adequate for concatenation. 
      if(isempty(ConSet))   
         ConSet = PcglimsConSet;               
      else
         if(size(ConSet,2) == size(PcglimsConSet,2))
            ConSet = [ConSet; PcglimsConSet];
         else
            error('finance:portcons:invalidAssetLims','Number of assets implied in type "GroupLims" is not consistent.');
         end   
      end  
      
      
   case 'groupcomparison'
      % Make sure that at least four data arguments are passed in
      if(NumArgs < 4)
         error('finance:portcons:missingGroupComparisonArg','Type "GroupComparison" requires four data arguments.');
      end   
      
      PcgcompConSet = pcgcomp(varargin{TypeIndex(i)+1:TypeIndex(i)+4});
      
      % Make sure that the number of columns is adequate for concatenation. 
      if(isempty(ConSet))   
         ConSet = PcgcompConSet;               
      else
         if(size(ConSet,2) == size(PcgcompConSet,2))
            ConSet = [ConSet; PcgcompConSet];
         else
            error('finance:portcons:invalidNumAssets','Number of assets implied in type "GroupComparison" is not consistent.');
         end   
      end  
      
      
   case 'custom'
      % Make sure that at least two data arguments are passed in
      if(NumArgs < 2)
         error('finance:portcons:missingCustomArg','Type "Custom" requires two data arguments');
      end
      
      A = varargin{TypeIndex(i)+1};
      b = varargin{TypeIndex(i)+2};
      b = b(:);
      
      if(length(b) ~= size(A,1))
         error('finance:portcons:mismatchElementsBRowsA','Number of elements in b must correspond to number of rows in A.');
      end
      
      NewConSet = [A b];
      
      % Make sure that the number of columns is adequate for concatenation. 
      if(isempty(ConSet))   
         ConSet = NewConSet;               
      else
         if(size(ConSet,2) == size(NewConSet,2))
            ConSet = [ConSet; NewConSet];
         else
            error('finance:portcons:invalidNumAssets','Number of assets implied in type "Custom" is not consistent.');
         end   
      end  

                  
   otherwise
      warning(['Type "' varargin{TypeIndex(i)} '" not recognized. Ignoring it.']);
   end
end

function [A,b] = pcpval(PortValue, NumAssets)
%PCPVAL Linear inequalities for fixing total portfolio value.
%   This function scales the total value of a portfolio of NumAssets assets 
%   to PortValue.  All portfolio weights, bounds, return, and risk values 
%   except ExpReturn and ExpCovariance (see PORTOPT) are in terms of PortValue.  
%  
%   [A,b] = pcpval(PortValue, NumAssets)
%
%   Inputs:
%     PortValue : Total value of asset portfolio (scalar).  PortValue is the
%     sum of the allocations in all the assets.  Use PortValue = 1 to specify
%     weights as fractions of the portfolio and return and risk numbers as
%     rates instead of value.
%
%     NumAssets : Number of available asset investments.
%
%   Outputs:
%     Matrix A and vector b such that A*Pwts' <= b sets the value,
%     where Pwts is a 1 by NumAssets vector of asset allocations.
%
%     Alternate Usage:  If called with fewer than 2 output arguments, 
%     returns A and b concatenated together: Cons = [A, b];
%     Cons = pcpval(PortValue, NumAssets)
% 
%   See also PORTOPT, PCALIMS, PCGLIMS, PCGCOMP, PORTCONS.

%   Copyright 1995-2006 The MathWorks, Inc.
%   $Revision: 1.8.2.2 $   $ Date: 1998/01/30 13:45:34 $

%----------------------------------------------------------------------------
% Process arguments and defaults
% PortValue [scalar]
% NumAssets [scalar]
%----------------------------------------------------------------------------

if nargin<2,
  error('finance:pcpval:missingInputs','Must specify PortValue and NumAssets');
end

%----------------------------------------------------------------------------
% build inequalities
%----------------------------------------------------------------------------

% Constrain the portfolio value from both sides
A = [ones(1,NumAssets); -ones(1,NumAssets)];
b = [PortValue; -PortValue];

%----------------------------------------------------------------------------
% Concatenation feature 
%----------------------------------------------------------------------------
if nargout<2,
  A = [A b];
end

%----------------------------------------------------------------------------
% end of function PCPVAL
%----------------------------------------------------------------------------
function [A,b] = pcalims(AssetMin, AssetMax, NumAssets)
%PCALIMS Linear inequalities for individual asset min and max allocation.
%   This function specifies the lower and upper bounds of portfolio
%   allocations in each of NASSETS available asset investments.
%
%   [A,b] = pcalims(AssetMin, AssetMax, NumAssets)  
% 
%   Inputs:
%     AssetMin, AssetMax : scalar or NASSETS long vector of minimum and
%     maximum allocations in each asset.  The entry NaN indicates no
%     constraint for that asset in that direction.  Scalar bounds are
%     applied to all the assets. 
%
%     NumAssets : (optional) number of assets NASSETS.  If not specified,
%     NumAssets is the length of AssetMin or AssetMax.
%
%   Outputs:
%     Matrix A and vector b such that A*Pwts' <= b enforces the constraints,
%     where Pwts is a 1 by NASSETS vector of asset allocations.
%
%     Alternate Usage:  If called with fewer than 2 output arguments, 
%     returns A and b concatenated together: Cons = [A, b];
%     Cons = pcalims(AssetMin, AssetMax, NumAssets)
% 
%   See also PORTOPT, PCPVAL, PCGLIMS, PCGCOMP, PORTCONS.
%

%   Copyright 1995-2006 The MathWorks, Inc.
%   $Revision: 1.8.2.2 $   $ Date: 1998/01/30 13:45:34 $

%----------------------------------------------------------------------------
% Process arguments and defaults
% AssetMin  [NumAssets by 1]
% AssetMax  [NumAssets by 1]
% NumAssets [scalar]
%----------------------------------------------------------------------------

% change min and max values to column vectors
if nargin<1
  error('finance:pcalims:missingInputs','Enter at least a vector of minimum asset values')
else
  AssetMin = AssetMin(:);
end

if nargin<2,
  AssetMax = NaN;
else
  AssetMax = AssetMax(:);
end

% Get the number of assets and expand scalar arguments
if nargin<3,
  NumAssets = max( length(AssetMax), length(AssetMin) );
end

if length(AssetMin)==1,
  AssetMin = AssetMin(ones(NumAssets,1));
elseif length(AssetMin)~= NumAssets,
  error('finance:pcalims:mismatchAssetMinAssetMax','Incompatible number of assets specified: AssetMin, AssetMax');
end

if length(AssetMax)==1,
  AssetMax = AssetMax(ones(NumAssets,1));
elseif length(AssetMax)~= NumAssets,
  error('finance:pcalims:mismatchAssetMinAssetMax','Incompatible number of assets specified: AssetMin, AssetMax');
end


%----------------------------------------------------------------------------
% build inequalities
%----------------------------------------------------------------------------

% Build all the max/min equations
A = [eye(NumAssets); -eye(NumAssets)];
b = [AssetMax; -AssetMin];

% Remove equations which are not specified
Mask = isnan(b);
A(Mask,:) = [];
b(Mask)   = [];

%----------------------------------------------------------------------------
% Concatenation feature 
%----------------------------------------------------------------------------
if nargout<2,
  A = [A b];
end

%----------------------------------------------------------------------------
% end of function PCALIMS
%----------------------------------------------------------------------------

function [PortRisk, PortReturn, PortWts] = portopt(ExpReturn, ExpCovariance, ... 
   NumPorts, PortReturn, ConSet)
%PORTOPT  Portfolios on constrained efficient frontier.
%  Returns portfolios on the mean-variance efficient frontier given asset
%  properties and a user specified set of portfolio constraints (ConSet).
%  Among a collection of NASSETS risky assets, computes a portfolio of asset
%  investment weights which minimize the risk for given values of the  
%  expected return.  The portfolio weights satisfy constraints specified in
%  linear inequality equations.  See PORTCONS for the generation of
%  portfolio constraints.
%
%  [PortRisk, PortReturn, PortWts] = portopt(ExpReturn, ExpCovariance, ...
%                                        NumPorts, PortReturn, ConSet) 
%
%  Inputs: 
%     ExpReturn is a 1xNASSETS vector specifying the expected (mean) 
%     return of each asset.
%
%     ExpCovariance is an NASSETSxNASSETS matrix specifying the 
%     covariance of the asset returns.  ExpCovariance must be symmetric with
%     no negative eigenvalues (positive semi-definite).
%
%     Either NumPorts or PortReturn specifies the set of efficient
%     portfolios computed.
%
%     NumPorts is the number of portfolios generated along the efficient
%     frontier when specific portfolio return values are not requested.
%     The default is 10 portfolios equally spaced between the minimum risk
%     point and the maximum possible return.  Enter NumPorts as an empty
%     matrix [], when specifying PortReturn.
%
%     PortReturn is a vector of length NPORTS containing the target return
%     values along the frontier.  If PortReturn is not entered or is empty,
%     NumPorts equally spaced returns between the minimum and maximum
%     possible values will be used.   
%
%     ConSet is a matrix of constraints for a portfolio of asset
%     investments.  An eligible 1 by NASSETS vector of asset allocation
%     weights, PortWts, satisfies the inequalities A*PortWts' <= b, 
%     where A = ConSet(:,1:end-1) and b = ConSet(:,end).  ConSet should
%     include at least an equation bounding the total value of the
%     portfolio below.
% 
%     See PORTCONS for a list of portfolio constraint types and their
%     corresponding financial parameters.
%
%     If the variable ConSet is not specified, a default constraint set 
%     will be used.  The default constraints scale the total value of the
%     portfolio to 1, and place the minimum weight of every asset at 0 to
%     prevent short-selling.
%
%  Outputs: 
%     PortRisk is an NPORTSx1 vector of the standard deviation of return 
%     for each portfolio. 
%            
%     PortReturn is an NPORTSx1 vector of the expected return of each 
%     portfolio.
%
%     PortWts is an NPORTSxNASSETS matrix of weights allocated to each
%     asset. Each row represents a different portfolio. 
%            
%
%  Notes: 
%     A plot of the efficient frontier is returned if the function is
%     invoked without output arguments.
%             
%
%  See also PORTCONS, PORTSTATS, EWSTATS, FRONTCON

%
%  Copyright 1995-2006 The MathWorks, Inc.
%  $Revision: 1.23.2.5 $   $ Date: 1997/08/15 14:07:00 $

%----------------------------------------------------------------
% Input argument validation
%----------------------------------------------------------------

% Check for input errors
if (nargin < 2)
     error('finance:portopt:missingInputs','You must enter ExpReturn and ExpCovariance.');
end

% Make sure that the number of returns entered matches the number of rows/columns in the  
% covariance matrix (which represents the number of assets).
ExpReturn = ExpReturn(:);
NASSETS = length(ExpReturn);
   
[covRows, covCols] = size(ExpCovariance);
if(covRows ~= covCols)
    error('finance:portopt:invalidCovMatrix','The covariance matrix must be NxN, where N = number of assets');
end
 
if norm(ExpCovariance - ExpCovariance',inf) > eps
	warning('finance:portopt:asymmetricCovariance', ...
		'Non-symmetric covariance matrix will be made symmetric.');
	ExpCovariance = 0.5*(ExpCovariance + ExpCovariance');
end
[L, D] = ldl(ExpCovariance);
if any(diag(D) < 0)
	warning('finance:portopt:nonposdefCovMatrix', ...
		'Non-positive-semidefinite covariance input.');
end
clear L
clear D

if size(ExpCovariance, 1) ~= NASSETS
  error('finance:portopt:invalidExpReturns','The number of expected returns does not equal the number of assets.');
end

% Determine which optional arguments were entered as non-empty.
if (nargin < 3 | isempty(NumPorts))  
   NumPortsEntered = 0; 
else
   NumPortsEntered = 1;
   NumPorts = double(NumPorts);
end

if (nargin < 4 | isempty(PortReturn))  
   PortReturnEntered = 0; 
else
   PortReturnEntered = 1; 
end

% When entering the target rate of return (PortReturn), 
% enter NumPorts as an empty matrix.
if (NumPortsEntered ==1) & (PortReturnEntered == 1)
   error('finance:portopt:emptyNumPorts','When entering the target rate of return, PortReturn, enter NumPorts as an empty matrix.');
end


% if 'default' problem, set a flag to trim negative weights to 0 (rmt4)
if (nargin < 5 | isempty(ConSet))
	% A default constraint matrix ConSet will be created if this is not entered.
	ConSet = portcons('default', NASSETS);
	trimweights = true;
else
	trimweights = false;
end


%----------------------------------------------------------------
% Constraint array construction
%----------------------------------------------------------------
% A = ConSet(:,1:end-1);
% b = ConSet(:,end);
% 
% % Call the function EQPARSE in order to find equality equations implied 
% % by the linear inequalities A* Wts <= b.
% % The first Neq rows of Aset and Bset will be interpreted as equalitiy
% % equations 
% [Aset, Bset, Neq]= eqparse(A,b);
% 
% Aeq = []; Beq = [];
% if(Neq > 0)
%    Aeq  = Aset(1:Neq, :);        %equalities constraints
%    Beq  = Bset(1:Neq, :);        %equalities constraints
% end
% 
% Aineq = []; Bineq = [];
% if(Neq < size(Aset,1))
%    Aineq = Aset(Neq+1:end, :);   %inequalities constraints
%    Bineq = Bset(Neq+1:end, :);   %inequalities constraints
% end

% Call INEQPARSE to find inequalites, equalites, and bounds implied by
% the linear inequalities in ConSet.
[Aineq, Bineq, Aeq, Beq, LB, UB] = ineqparse(ConSet);

% If PortReturn has been entered, check them for legality relative to
% constraints. If it has not been entered, construct the default 
% return range.
W0 = ones(NASSETS, 1)/NASSETS;


%----------------------------------------------------------------
% Maximum return calculation
%----------------------------------------------------------------

% Find the maximum expected return achievable, given the individual asset
% expected returns and all the other constraints. 

% Set the options: ('LargeScale' mode os turned off because it can cause 
%                   some warnings to be thrown.)

% set a larger value for MaxIter (rmt4)
MaxIter = 100000;
options = optimset('LargeScale','off','Display','off','MaxIter',MaxIter);

[MaxReturnWeights, Fval, ErrorFlag] = linprog(-ExpReturn, Aineq, Bineq, Aeq, Beq, LB, UB, W0, options);

if ErrorFlag ~= 1
   error('finance:portopt:noPortfolios','No portfolios satisfy all the input constraints');
end

MaxReturn = MaxReturnWeights'*ExpReturn;

% Find the minimum variance return.
F = zeros(NASSETS, 1);

[MinVarWeights,  Fval, ErrorFlag] = quadprog(ExpCovariance, F, Aineq, Bineq, Aeq, Beq, LB, UB, W0, options);

if ErrorFlag ~= 1
   error('finance:portopt:noMinVarPortfolio','A solution was not feasible for the minimum variance portfolio.');
end

MinVarReturn = MinVarWeights'*ExpReturn;


%----------------------------------------------------------------
% Calculate return corresponding to minimum risk (variance)
%----------------------------------------------------------------
if PortReturnEntered
  % check the requested returns against points on the frontier
  % use a small numerical fudge factor and assume returns are postitive
  if ( min(PortReturn) < MinVarReturn*(1 - 1000*eps) )
    error('finance:portopt:missingMinPortfolios',sprintf( ...
        ['One or more requested returns are less than the return %f ',...
         'of the least risky portfolio'], MinVarReturn));
  end
  if ( max(PortReturn) > MaxReturn*(1 + 1000*eps) )
    error('finance:portopt:missingMaxPortfolios',sprintf( ...
        ['One or more requested returns are greater than the maximum ',...
         'achievable return of %f'], MaxReturn));
  end
 
  NumFrontPoints = length(PortReturn);
  PortfOptResults = zeros(NumFrontPoints, 2 + NASSETS);
  StartPoint = 1;
  EndPoint = NumFrontPoints;

else
   % If NumPorts is not entered, set a default value.
   if NumPortsEntered
      NumFrontPoints = NumPorts;
   else
   	NumFrontPoints = 10; 
	end

	MinVarStd = sqrt(MinVarWeights' * ExpCovariance * MinVarWeights);
      
   if MaxReturn > MinVarReturn % This is just the algorithm from MATLAB function, LINSPACE
      PortReturn = [MinVarReturn+(0:NumFrontPoints-2)*(MaxReturn-MinVarReturn)/(NumFrontPoints-1) MaxReturn];
   else
      PortReturn = MaxReturn;
      NumFrontPoints = 1;
   end
		
	PortfOptResults = zeros(NumFrontPoints, 2 + NASSETS);
	PortfOptResults(1, :) = [MinVarReturn MinVarStd MinVarWeights(:)'];
   StartPoint = 2;
   EndPoint = NumFrontPoints-1;
end

FrontPointConstraint = -ExpReturn';
Aeq = [FrontPointConstraint; Aeq ]; % Add a new equality constraint
Beq = [0; Beq];                    % Add a new equality constraint

W0 = MaxReturnWeights;
% Set the options with larger MaxIter (rmt4)
options = optimset('LargeScale','off','Display','off','MaxIter',MaxIter);
   
for Point = StartPoint:EndPoint
   Beq(1) = -PortReturn(Point);
   
   [Weights, Fval, ErrorFlag] = quadprog(ExpCovariance, F, Aineq, Bineq, Aeq, Beq, LB, UB, W0, options);  
   if ErrorFlag ~=1
		PortfOptResults(Point, :) = [Beq(2) nan*ones(1, NASSETS+1)];
	else
		Return = dot(Weights, ExpReturn);
		Std = sqrt(Weights'*ExpCovariance*Weights);
		PortfOptResults(Point, :) = [Return Std Weights(:)'];
	end

end

if ~PortReturnEntered
   % The last row corresponds to the point of maximum return, which has already been calculated
   Std = sqrt(MaxReturnWeights'*ExpCovariance*MaxReturnWeights);
   PortfOptResults(end, :) = [MaxReturn, Std, MaxReturnWeights(:)'];
end


%----------------------------------------------------------------
% Validate results and generate output
%----------------------------------------------------------------

ErrorIndex = find(isnan(PortfOptResults(:, 2)));
if ~isempty(ErrorIndex)
	NumErrors = num2str(length(ErrorIndex));
	warning(['A solution was not feasible for ', NumErrors, ' expected return(s).']);
end

PortReturn = PortfOptResults(:, 1);
PortRisk = PortfOptResults(:, 2);
PortWts = PortfOptResults(:, 3:size(PortfOptResults, 2));

% if lower-bound weights is zero, force all negative weights to be zero (rmt4)
if trimweights
	PortWts(PortWts < 0) = 0;
end

%A plot of the efficient frontier is returned if the function is
%invoked without output arguments.
if nargout == 0
 		FrontWin = figure;

      set(FrontWin,'NumberTitle','off');
		set(FrontWin,'Name','Efficient Frontier');
		set(FrontWin,'Resize','on');
	
		set(FrontWin,'Tag','FrontWin');
	
      plot(PortRisk, PortReturn);
      title('Mean-Variance-Efficient Frontier', 'Color', 'k');
	  xlabel('Risk (Standard Deviation)');
      ylabel('Expected Return');
      grid on;
      
      % prevent any output
      clear PortRisk;
end

%----------------------------------------------------------------
% Auxiliary function(s)
%----------------------------------------------------------------

function [A,b,Aeq,beq,LB,UB] = ineqparse(Ain, bin)
%INEQPARSE Find inequalities, equalities, and bounds implied by Ain*x <= bin.
%  Identifies equalities specified as Arow*x <= bval, -Arow*x <= -bval.
%  Parses out duplicate entries and all zero equations.
%  Finds bound constraints among inequalities.
%
%  [A,b,Aeq,beq,LB,UB] = ineqparse(Ain, bin)
%  [A,b,Aeq,beq,LB,UB] = ineqparse([Ain, bin])
%
%  The function does not catch linear combinations of inequalities which
%  together imply an equality constraint.
%
%  See also PORTCONS.
%
%----------------------------------------------------------------------
%
% % Test the logic to parse out inequalites with single-entry equations
% % 1 : redundant equaltiy
% % 2 : equality
% % 3 : redundant upper bound
% % 4 : lower bound
% % 5 : upper bound
% Astart = [1 2 -1 1 3 -1 -2 -4 5 3]'
% Ain = full(sparse(1:length(Astart),abs(Astart),Astart))
% [A,b,Aeq,beq,LB,UB] = ineqparse(Ain, zeros(length(Astart),1))
%
% % Catch rows which are multiples
% m = 1:length(Astart)
% Ain = diag(m)*Ain
% [A,b,Aeq,beq,LB,UB] = ineqparse(Ain, zeros(length(Astart),1))
%
% % Degenerate case with equality, lower, upper bounds
% C = portcons('default',3,'AssetLims',0,[0.5 0.6 0.7],3)
% [A,b,Aeq,beq,LB,UB] = ineqparse(C)
%
% % Case with a general inequality constraint
% C = portcons('default',3,'AssetLims',0,[0.5 0.6 0.7],3, ...
%              'Custom',[0.1 0.2 0.3],0.40)
% [A,b,Aeq,beq,LB,UB] = ineqparse(C)
%

% J. Akao 8/22/99

% find usage ineqparse(ConSet)
if nargin==1
  bin = Ain(:,end);
  Ain = Ain(:,1:end-1);
end

[NumEquations, NumVars] = size(Ain);
if any( size(bin)~=[NumEquations, 1] )
  error('finance:portopt:mismatchAandB','dimensions of A and b are inconsistent')
end

% Pull out degenerate rows
I = all(Ain==0,2);
if(any(I))
   warning('Degenerate rows found in constraint matrix. Eliminating these contraints');
   Ain(I,:) = [];
   bin(I) = [];
end

% Constraint rows
ConRows = [Ain, bin];

% Form numerator and denominator dot products.
%
% row I and row J are the same direction when:
%   rowI*rowJ' == sqrt(rowI*rowI')*sqrt(rowJ*rowJ') 
%        numIJ == denIJ
%
% row I and row J are the opposite direction when:
%   rowI*rowJ' == - sqrt(rowI*rowI')*sqrt(rowJ*rowJ') 
%        numIJ == - denIJ

% square (rowI*rowJ') but keep the sign
numIJsqrt = ConRows*ConRows';
numIJ = sign(numIJsqrt).*(numIJsqrt.*numIJsqrt);

% form (rowI*rowI') times (rowJ*rowJ')
rowKdot = dot(ConRows, ConRows, 2);
[rowIdot, rowJdot] = meshgrid(rowKdot, rowKdot);
denIJ = rowIdot .* rowJdot;

% record which equations are negations or duplicates
% denIJ is always positive
% take the upper triangular part only
%
% isdupIJ [NumEqs x NumEqs] row I is a positive multiple of row J
% isnegIJ [NumEqs x NumEqs] row I is a negative multiple of row J
reltol = 1000*eps;
isdupIJ = ( denIJ*(1-reltol) <  numIJ ) & (  numIJ < denIJ*(1+reltol) );
isnegIJ = ( denIJ*(1-reltol) < -numIJ ) & ( -numIJ < denIJ*(1+reltol) );

isdupIJ = triu(isdupIJ, 1);
isnegIJ = triu(isnegIJ, 1);

% search through the equations and clean out equalities and duplicates.
% store the equalities separately.
%
% ConEqs  [NumEqs   x NumVars+1] : [Aeq, beq]
% ConRows [NumInEqs x NumVars+1] : [A, b]
ConEqs = zeros(0, NumVars+1);

i=1;
while (i < size(ConRows,1) )
  % find negations and duplicates of this row
  RowIsNeg = isnegIJ(i,:);
  RowIsDup = isdupIJ(i,:);
  
  % negations and duplicates should be removed from the inequality list
  IndRemove = RowIsNeg | RowIsDup;
  
  if any(RowIsNeg)
    % add the row to the equality list
    ConEqs = [ConEqs; ConRows(i,:)];
    
    % remove the row from the inequality list along with negs and dups
    IndRemove(i) = 1;
  else
    % equation i has been left in
    i = i + 1;
  end
  
  % remove equations from the inequality list
  ConRows(IndRemove,:) = [];
  isnegIJ = isnegIJ(~IndRemove, ~IndRemove);
  isdupIJ = isdupIJ(~IndRemove, ~IndRemove);
end

% Break up into left and right hand sides
Aeq = ConEqs(:,1:NumVars);
beq = ConEqs(:,NumVars+1);
A = ConRows(:,1:NumVars);
b = ConRows(:,NumVars+1);

% search through the inequalities and find bounds
% SingleValue * x(Ind) <= b(Ind)
%
% IndSingle   [NumInEqs x 1] true if only 1 non-zero value in row of A
% SingleValue [NumInEqs x 1] only valid for IndSingle == 1
%
% VarNum       [NumInEqs x NumVars] column of each entry of A
% SingleVarNum [NumInEqs x 1] column of first non-zero entry in A
%
IndSingle   = sum(A~=0 , 2) == 1;
SingleValue = sum(A    , 2);

IndLower = IndSingle & ( SingleValue < 0 );
IndUpper = IndSingle & ( SingleValue > 0 );

VarNum = (1:NumVars);
VarNum = VarNum(ones(size(A,1),1),:);
VarNum(A==0) = Inf;
SingleVarNum = min(VarNum,[],2);

if any(IndLower)
  LB = -Inf*ones(NumVars,1);

  % find the variable and the bound value
  VarNum = SingleVarNum(IndLower);
  BVal = b(IndLower)./SingleValue(IndLower);

  % apply the most restrictive bound to each variable
  UniqVar = unique(VarNum);
  if length(UniqVar)==length(VarNum)
    % no variables have multiple parallel bounds
    LB(SingleVarNum(IndLower)) = BVal;
  else
    for Var=UniqVar(:)'
      LB(Var) = max( BVal( VarNum==Var ) );
    end
  end
    
else
  LB = [];
end

if any(IndUpper)
  UB = Inf*ones(NumVars,1);
  
  % find the variable and the bound value
  VarNum = SingleVarNum(IndUpper);
  BVal = b(IndUpper)./SingleValue(IndUpper);

  % apply the most restrictive bound to each variable
  UniqVar = unique(VarNum);
  if length(UniqVar)==length(VarNum)
    % no variables have multiple parallel bounds
    UB(SingleVarNum(IndUpper)) = BVal;
  else
    for Var=UniqVar(:)'
      UB(Var) = min( BVal( VarNum==Var ) );
    end
  end
    
else
  UB = [];
end

% remove lower or upper bound inequalities
A(IndSingle,:) = [];
b(IndSingle) = [];


%----------------------------------------------------------------
% Auxiliary function(s)
%----------------------------------------------------------------

function [Aset, Bset, Neq]= eqparse(A, b)
% EQPARSE find equality equations implied by linear inequalities A*x <= b
% Orders equality equations to the top of the list and returns the number of
% equations to be treated as equalities.
% 
% [Aset, Bset, Neq] = eqparse(A,b)


[NumEquations, NumVars] = size(A);
if any( size(b)~=[NumEquations, 1] )
  error('finance:portopt:mismatchAandB','dimensions of A and b are inconsistent')
end

% Pull out degenerate rows
I = all(A==0,2);
if(any(I))
   warning('Degenerate rows found in constraint matrix. Eliminating these contraints');
   A(I,:) = [];
end


ConLE = [A, b]; % equations which are less than or equal

% Normalize the rows so that dot(row,row) = 1
RowNorm = sqrt(dot( ConLE, ConLE, 2));
ConLE = ConLE./RowNorm(:,ones(1,NumVars+1));

% start a variable for any equality equations found
ConEq = zeros(0,NumVars+1);

% pull out equalities from the inequalities
i = 1; % number of inequalities found + 1
while ( i < size(ConLE,1) )
  % candidate row
  ConRow = ConLE(i,:);
  
  % find which remaining rows are negatives of ConRow
  % dot(ConRow, ConRow) == -1
  % where dot(Row, Row) = 1 because of normalization
  DotWithRows = ConLE(i+1:end,:)*ConRow';
  RowIsNeg = (-1-1000*eps) < DotWithRows & DotWithRows < (-1+1000*eps);
  
  % find which rows are the same as ConRow
  RowIsSame = (1-1000*eps) < DotWithRows & DotWithRows < (1+1000*eps);
  
  if any(RowIsNeg) 
    % add Row i to the equality list
    ConEq = [ConEq; ConRow];

    % build list to delete from the inequalities
    IndRemove = [i; i + find(RowIsNeg | RowIsSame)];
  else
    % build the delete list from the duplicate rows
    IndRemove = i + find(RowIsSame);
    
    % equation i has been left in
    i = i + 1;
  end
  ConLE(IndRemove,:) = [];

end

Aset = [ConEq(:,1:NumVars); ConLE(:,1:NumVars)];
Bset = [ConEq(:,NumVars+1); ConLE(:,NumVars+1)];
Neq = size(ConEq,1);



