% Recall that net returns take values in [-1,infty)

muf=0.02;
mu=[0.08;0.03;0.05];
Omega = [0.3 0.02 0.01 ;
         0.02 0.15 0.03;
         0.01 0.03 0.18 ];

     

bone = ones(length(mu),1);
f=zeros(length(mu),1);
A=-eye(length(mu));
Aeq = [bone';mu'];
mup=0.04;
beq1 = [1; mup];
b=zeros(length(mu),1);
w = quadprog(Omega,f,A,b,Aeq,beq1);
w

ngrid=100;
mups=linspace(min(mu),max(mu),100);
varp=zeros(ngrid,1);

for i=1:ngrid
    beq = [1;mups(i)];
    w=quadprog(Omega,f,A,b,Aeq,beq);
    varp(i)=w'*Omega*w;
end

imin = find(varp==min(varp));

sharpratio = (mups'-ones(ngrid,1)*muf)./sqrt(varp);
sharpratio

Itangency = find(sharpratio==max(sharpratio));

sharpratio(Itangency)
mups(Itangency)

plot(sqrt(varp),mups)
