function [x,mixe,mixsd] = qqp(rt,e,C,lbd,ubd,x0);
%  [x,mixe,mixsd]  = qqp(rt,e,C,lbd,ubd,x0)
%        computes an optimal asset mix
%     Copyright William F. Sharpe
%     This version: April 13, 1995
%        inputs:
%           rt: risk tolerance
%            e: expected return vector
%            C: covariance matrix
%          lbd: lower bound vector
%          ubd: upper bound vector
%           x0: initial feasible mix vector
%
%       outputs:
%            x: optimal mix vector
%         mixe: x'*e
%        mixsd: sqrt(x'*C*x);
%
%       maximizes:    rt*(x'*e) - x'*C*x
%         subject to: sum(x) = sum(x0)
%                     lbd <= x <= ubd  
% algorithm based on:
%    William F. Sharpe, "An Algorithm for Portfolio Improvement,"
%    in Advances in Mathematical Programming and Financial Planning
%    K.D. Lawrence, J.B. Guerard, Jr., and Gary D. Reeves, Editors
%    JAI Press, Inc., 1987, pp. 155-170.

% maximum number of iterations
    maxit = 500;

% convert any row vectors to column vectors
   if size(e,2)  >1; e=e';     end;
   if size(lbd,2)>1; lbd=lbd'; end;   
   if size(ubd,2)>1; ubd=ubd'; end;
   if size(x0,2) >1; x0=x0';   end;

% return if x0 is not feasible
   if (x0>=lbd) & (x0<=ubd)
       % ok
     else
       disp('******** ERROR: x0 is not feasible *****');
       disp('         press any key to return');
       return;
   end;

% set minimum MU change to continue
    minMUchg = 0.0001;

% initialize
    x = x0;  
    n = size(C,1);
  
% continue to improve portfolio until further improvement impossible
%   when done, return

  iterations = 0;

  while 1==1;
    % compute marginal utilities
        mu =  rt*e - 2*C*x;
   % find best variable to add
        [MUadd,Aadd] = max(mu - 1E200*(x>=ubd));
    % find best variable to subtract
        [MUsub,Asub] =  min(mu + 1E200*(x<=lbd));
    % terminate and return if change in mu is less than minimum
        if (MUadd - MUsub) <= minMUchg
           % compute mix e and sd
              mixe  = x'*e;
              mixsd = sqrt(x'*C*x);
           % terminate
              return
        end
    % set up delta vector
        d = zeros(n,1);
        d(Aadd) = 1;
        d(Asub) = -1;
    % compute step size
       k = zeros(3,1);
       % optimal unconstrained step size
          k(1) = ((rt*d'*e)-2*(x'*C*d))/(2*(d'*C*d));
       % maximum step size based on upper bounds
          Slack = ubd - x;
          AUp = find(d>0);
          k(2) = min(Slack(AUp)./d(AUp));
       % maximum step size based on lower bounds
          Slack = x-lbd;
          ADown = find(d<0);
          k(3) = min(Slack(ADown)./ (-d(ADown)));
       % minimum step size
          kmin = min(k);
       % terminate and return if minumum step size is zero
         if (kmin == 0)
            % compute mix e and sd
               mixe  = x'*e;
               mixsd = sqrt(x'*C*x);
            % terminate 
               return;  
         end
      % count and terminate if maximum iterations exceeded 
         iterations = iterations+1;
         if iterations > maxit
            % compute mix e and sd
               mixe  = x'*e;
               mixsd = sqrt(x'*C*x);
            % terminate 
               return;  
         end
    % change mix
       x = x + ( kmin*d) ;
  end;