function fstar=MRtreat(bounds)

%evaluates propositon 1 in J. Stoye (2005):
%Minimax Regret Treatment Choice with Incomplete Data and Many Treatments.
%syntax: fstar=MRtreat(bounds)
%bounds should be (2xT); for further documentation see
%http://homepages.nyu.edu/~js3909/

T=length(bounds(1,:))

deltas=bounds(1,:)-bounds(2,:);

%find f, lambda, and t=max{t:ft>0}
success=0;
t=1;
A=eye(T+1);
A(T+1,T+1)=0;
A(1,1)=bounds(1,1)-bounds(2,1);
A(1,T+1)=1;
A(T+1,:)=[ones(1,T),0];
B=zeros(T,1);
B=[B;1];
B(1)=bounds(1,1);
bounds(1,T+1)=-inf;
while success==0
    t=t+1;
    A(t,t)=deltas(t);
    A(t,T+1)=1;
    B(t)=bounds(1,t);
    guess=A\B;
    lambda=guess(T+1);
    success=(lambda>=bounds(1,t+1));
end
bounds=bounds(:,1:T);

%case distinction
nu=max(bounds(2,:))
nu-bounds(2,:);
r=(nu-bounds(2,:))./deltas;
for u=1:T
    uu(u)=sum(r(1:u));
end
tstar=sum(uu<1)+1
if tstar>t
    plan=1
elseif uu(tstar)~=1
    plan=2
else
    plan=3
end

%case (i)
if plan==1
    fstar=guess(1:T)'
end

if plan>1

%common between (ii) and (iii)
lambda=bounds(1,tstar);
f=(bounds(1,:)-lambda)./deltas;
f=zeros(1,T)+[f(1,1:tstar),zeros(1,T-tstar)];
maxmins=find(bounds(2,:)==nu);
for count=1:length(maxmins)
    fstar(count,:)=f;
    fstar(count,maxmins(count))=fstar(count,maxmins(count))+1-sum(f);
end

%case (ii)
if plan==2
    fstar
end

%case (iii)
if plan==3
    if tstar>=T
        lambda=guess(T+1);
    else
        lambda=max(guess(T+1),bounds(1,tstar+1));
    end
    f=(bounds(1,:)-lambda)./deltas;
    f=zeros(1,T)+[f(1,1:tstar),zeros(1,T-tstar)];
    for count=1:length(maxmins)
        fstar(length(maxmins)+count,:)=f;
        fstar(length(maxmins)+count,maxmins(count))=fstar(length(maxmins)+count,maxmins(count))+1-sum(f);
    end
    fstar
end
        
end

%construction of sigma^star
for count=1:length(fstar(:,1))
    for n=1:T
        regret(count,n)=bounds(1,n)-sum(fstar(count,:).*bounds(2,:))-fstar(count,n).*deltas(n);
    end
end
regret

if plan==1
    A=zeros(t+1,t+1);
    for count=1:t
        A(count,count)=deltas(count);
    end
    A(t+1,1:t)=ones(1,t);
    A(1:t,t+1)=-ones(t,1);
    B=[-bounds(2,1:t)';1];
    sigmastar=zeros(1,T);
    solution=A\B;
    sigmastar(1:t)=solution(1:t)'
    mu=bounds(1,:).*sigmastar+bounds(2,:).*(1-sigmastar)
end

if plan>1
    sigmastar=zeros(1,T);
    count=1;
    while sum(sigmastar)<1
        sigmastar(count)=min((nu-bounds(2,count))/deltas(count),1-sum(sigmastar));
        count=count+1;
    end
    sigmastar
end