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;