clear all;
global sigma;
global r;
global target_return;

%Plots efficient frontier for a portfolio of stocks
%based on Markowitz model

%------------- Enter parameters ----------------------
sigma=    [ .08 -.05 -.05 -.05 
		-.05  .16 -.02 -.02 
		-.05 -.02 .35 .06 
		-.05 -.02 .06 .35] ; %covariance matrix
r= [.05 ; .10 ; .15 ; .30];  %Expected return for each stock
tgvec=0.01:0.01:.2; %vector of target returns for plotting
%--------------------------------------------

%Adds 2 parameters w_n-1 and w_n to the input vector using the constraints
% 1. w1+w2+......wn=1
% 2. r1w1+r2w2+......rnwn=target_return
% where w1,..wn are stock weigths in portfolio
% r1,..rn are expected returns for each stock
function ret=GetParamvec(x)
	global r;
	global target_return;
	n=size(r,1);
	rvec=r(1:n-2);
	rhsmat=zeros(2,1);
	rhsmat(1)=1-sum(x,1);
	rhsmat(2)=target_return-r(1:n-2)'*x;
	lhsmat=ones(2,2);
	lhsmat(2,1)=r(n-1);
	lhsmat(2,2)=r(n);
	unk=inv(lhsmat)*rhsmat;
	test=inv(lhsmat);
	ret=[x;unk];
endfunction

%returns variance of portfolio
function ret=GetVariance(x)
	global sigma;
	w=GetParamvec(x);
	ret=w'*sigma*w;
endfunction

n=size(r,1);
%remove last 2 parameters as they will be determined by constraints
initialparam=ones(n-2,1)*(1/n); %initial parameter vector of weights

tgvec=tgvec';
psigvec=tgvec*0;
pretvec=tgvec*0;
resmat=zeros(size(tgvec,1),n+2);

for i=1:size(tgvec,1),
	target_return=tgvec(i);
	res=fmins('GetVariance',initialparam);
	w=GetParamvec(res);
	psigvec(i)=w'*sigma*w;
	pretvec(i)=w'*r;
	resmat(i,:)=[psigvec(i) pretvec(i) w'];
end

%results matrix 
%1st col->variance 
%2ndcol->returns
%3..n cols ->weights
resmat  

plot(psigvec,pretvec);
title('Efficient frontier for portfolio');
xlabel('Risk');
ylabel('Return');