FIRST CODE: NO PRIORITIES
(NOTE: run each code separately)

clear all;
rand('state',sum(100*clock));
%%%%%%%%%%%%%%%%%%%%%%%SIMULATION BM-RTB CADA DA%%%%%%%%%%%%%%%%%%%
%%Parameters
nocount=4001;
noas=4;
mina=-.1;
noiter=50;
minnosc=3;
varscho=3;
mincapa=10;
varcapa=3;
%% Welfare and alert matrices
welfabm=zeros(noas,2,varscho,varcapa,noiter);
alertbm=zeros(noas,2,varscho,varcapa,noiter);
welfaca=zeros(noas,2,varscho,varcapa,noiter);
alertca=zeros(noas,2,varscho,varcapa,noiter);
welfada=zeros(noas,2,varscho,varcapa,noiter);
welnaibm=zeros(noas,2,varscho,varcapa,noiter);
welnaica=zeros(noas,2,varscho,varcapa,noiter);
welnaida=zeros(noas,2,varscho,varcapa,noiter);
%%%%%%%%%%%%%FOR EACH NUMBER OF SCHOOLS

for morescho=1:varscho
noscho=minnosc+morescho; %Regarded as J sometimes

%%%%%%%%%%%%%FOR EACH SCHOOL CAPACITY

for morecapa=1:varcapa
capac=mincapa+10*morecapa;

nostra=factorial(noscho);
nostud=noscho*capac;
nonaive=nostud/2; %%When naive mode is yes, half students naive

%%Creation permutation matrix (J! rows, J columns)
%%It is the strategy space
permu=perms(1:noscho);
[ordprob,ordperm]=sort(permu,2);
%ordprob just defines some matrix with proper dimensions I later use
%ordperm(i,j): in what position strategy i locates school j

%%%%%%%%%%%%%%%%%%%%%%%%%%FOR EACH ITERATION
for iter=1:noiter;

%%Generation valuations: vij=a·uj+(1-a)·uij
commval=rand(1,noscho);
commval2=commval(ones(nostud,1),:);
privval=rand(nostud,noscho);

%%%%%%%%%%%%%%%%%%%%%%%%%%%% FOR EACH VALUE OF a...
for aid=1:noas
a=mina+.2*aid;
vals=a*commval2+(1-a)*privval;

%%Valuations rescaled to the simplex
rescval=zeros(nostud,noscho);
for stu=1:nostud
rescval(stu,:)=vals(stu,:)/sum(vals(stu,:),2);
end

%%True ordinal preferences
[hale,ranking]=sort(rescval,2,'descend');
%Note: ranking(i,j), school student i puts in j-th postion

%% Bias for the truth
truth=zeros(nostra,nostud);
for estu=1:nostud
truth(:,estu)=sum((permu==ranking(estu*ones(nostra,1),:)),2);
end
biastru=.0000001*truth;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%DA%%%%%%%%%%%%%%%%%%%%%%
%%NOTE: DA does not depend on naive mode

%% Demands
demda=zeros(nostra,1);
for str=1:nostra
computo=ones(nostud,1);
for col=1:noscho
computo=computo.*(ranking(:,col)==permu(str,col));
end
demda(str)=sum(computo);
end

%% Cutoffs
cutoffs=zeros(noscho,1);
cutbuild=zeros(noscho,1);
cutdemda=zeros(nostra,1);
%Round 1 (most popular school)
for scda=1:noscho
deno=demda'*(permu(:,1)==scda);
cutbuild(scda)=min([1 capac/max([.0000001 deno])]);
end
[cutmin,posicut]=min(cutbuild);
cutoffs(posicut)=cutmin;
cutscs=[posicut];
%Rest of rounds
for roda=2:noscho
for scda=1:noscho
locscda=ordperm(:,scda);
for stda=1:nostra
cutdemda(stda)=max(cutoffs(permu(stda,1:locscda(stda))));
end
locrepe=locscda(:,ones(1,roda-1));
cutahead=(locrepe-ordperm(:,cutscs)>0);
criterio=(locscda-sum(cutahead,2)==1);
nume=capac+cutdemda'*(demda.*criterio);
deno=demda'*criterio;
cutbuild(scda)=min([1 nume/max([.0000001 deno])]);
end
cutbuild=cutbuild+50000*(cutoffs>0);
[cutmin,posicut]=min(cutbuild);
cutoffs(posicut)=cutmin;
cutscs=[cutscs posicut];
end
%%%% Calculation Welfare
% Probabilities for each student
probtool=zeros(nostud,1);
probalda=zeros(nostud,1);
for number=1:noscho
inccut=cutoffs(ranking(:,number))-probtool;
probalda(:,number)=(inccut>0).*inccut;
probtool=probtool+(inccut>0).*inccut;
end
%Rearrange probabilities (now prob(i,j): student i,
%school ranked j-th (in ranking), want j to identify school.
ordprda=zeros(nostud,noscho);
[whatever,ordrank]=sort(ranking,2);
payda=zeros(nostud,1);
for st3=1:nostud
prost=probalda(st3,:);
ordprda(st3,:)=prost(ordrank(st3,:));
payda(st3)=ordprda(st3,:)*rescval(st3,:)';
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NAIVE MODE NO/YES
for naimode=1:2
% 1=no naive; 2=naive
naive=naimode-1;

%%Welfare DA
welfada(aid,naimode,morescho,morecapa,iter)=mean(payda);
welnaida(aid,naimode,morescho,morecapa,iter)=mean(payda(1:nonaive));



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BOSTON MECHANISM RTB%%%%%%%%%%%%%%%%%%%
coincid=zeros(nocount-1,1);
reveal1=ranking;
count=1
while count<nocount
reveal0=reveal1;

%%Demands
demanda=zeros(nostra,1);
for str=1:nostra
computo=ones(nostud,1);
for col=1:noscho
computo=computo.*(reveal0(:,col)==permu(str,col));
end
demanda(str)=sum(computo);
end

%%Calculation probabilities per round
probrtb=zeros(noscho,noscho);
unass=demanda;
avail=capac*ones(noscho,1);
for ronda=1:noscho;
for esc=1:noscho;
pedido=sum(unass.*(permu(:,ronda)==esc));
probrtb(esc,ronda)=min([1 avail(esc)/max([0.000000000001 pedido])]);
avail(esc)=(avail(esc)-pedido)*(probrtb(esc,ronda)==1);
unass=unass.*(ones(nostra,1)-(permu(:,ronda)==esc)*probrtb(esc,ronda));
end
end

%%Assignment probabilities for each strategy
prob=zeros(nostra,noscho);
probunas=ones(nostra,1);
for colu=1:noscho
prob(:,colu)=probunas.*probrtb(permu(:,colu),colu);
probunas=probunas-prob(:,colu);
end

%Rearrange probabilities (now prob(i,j): strategy i,
%school ranked j-th (in permu), want j to identify school.
ordprob=zeros(nostra,noscho);
for st3=1:nostra
prost=prob(st3,:);
ordprob(st3,:)=prost(ordperm(st3,:));
end

% %%Demand bias
% dembias=.0000001*demanda(:,ones(1,nostud));

%%Best Responses and Welfare
payoff=ordprob*rescval'+biastru;
[maxpay,best]=max(payoff);
reveal1=permu(best',:);

%%%If naive mode...
if naive==1;
reveal1(1:nonaive,:)=reveal0(1:nonaive,:);
for naistu=1:nonaive
ranktool=ranking(naistu,:);
rankrepe=ranktool(ones(nostra,1),:);
location=find(sum((permu-rankrepe==0),2)==noscho);
maxpay(naistu)=ordprob(location,:)*rescval(naistu,:)';
end
end

welfabm(aid,naimode,morescho,morecapa,iter)=mean(maxpay);
welnaibm(aid,naimode,morescho,morecapa,iter)=mean(maxpay(1:nonaive));


%%Check if best response is equilibrium
coincid(count)=length(find(reveal1-reveal0==0));
if coincid(count)==noscho*nostud
break
end

% Auxilliary converge criterion
if (count>nocount/2) && (coincid(count)>=(.98+naive*.01)*noscho*nostud)
break
end

%% "Smoothing" to avoid cycles
lott=rand(nostud,1);
keep0=(lott<.975);
reveal1=reveal1+keep0(:,ones(1,noscho)).*(reveal0-reveal1);

%%Finally...
count=count+1
end
alertbm(aid,naimode,morescho,morecapa,iter)=(count==nocount)*coincid(nocount-1)/(noscho*nostud);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%CADA as DA-FCP%%%%%%%%%%%%%%%%%%%%%%%%%%

coincid=zeros(nocount-1,1);
reveal1=ranking;
count=1
while count<nocount
reveal0=reveal1;
%%Demands
demca=zeros(nostra,1);
for str=1:nostra
computo=ones(nostud,1);
for col=1:noscho
computo=computo.*(reveal0(:,col)==permu(str,col));
end
demca(str)=sum(computo);
end
%% FIRST: First choice assignment probabilities
remnext=capac*ones(noscho,1);
prfirst=zeros(noscho,1);
over=zeros(noscho,1);
for sc=1:noscho
deno=demca'*(permu(:,1)==sc);
prfirst(sc)=min([1 capac/max([.0000001 deno])]);
remnext(sc)=max([0 capac-deno]);
over(sc)=(capac<deno);
end
%% NEXT: cutoffs DA with remaining slots
% Overdemanded and underdemanded schools
overdem=find(over);
underdem=find(over==0);
nousc=length(underdem);
demnext=demca.*(ones(nostra,1)-prfirst(permu(:,1)));
% Cleaning rows with underdemanded targets
permurep=permu(:,1);
permurep=permurep(:,ones(1,nousc));
underrep=underdem(:,ones(1,nostra))';
demnext=demnext(find(min(abs(permurep-underrep),[],2)));
permnext=permu(find(min(abs(permurep-underrep),[],2)),:);
% Cleaning cells with overdemanded schools
for str=1:size(permnext,1)
for sc=overdem'
permnext(str,:)=[permnext(str,find(permnext(str,:)-sc)) 0];
end
end
permnext=permnext(:,1:nousc);
nostrane=size(permnext,1);
% Last before DA: changing school id's to 1:nousc
tool=zeros(nostrane,nousc);
for usc=underdem'
tool=tool+(permnext==usc)*find(underdem==usc);
end
permnext=tool;
[whatever,ordnext]=sort(permnext,2);
%%%%%%%%%%%%% At last, DA with remaining slots:
%% Cutoffs
if nousc==1
cutoffs=[1];
else
cutoffs=zeros(nousc,1);
cutbuild=zeros(nousc,1);
cutdemne=zeros(nostrane,1);
%Round 1 (most popular school)
for number=1:nousc
nume=remnext(underdem(number));
deno=demnext'*(permnext(:,1)==number);
cutbuild(number)=min([1 nume/max([.0000001 deno])]);
end
[cutmin,posicut]=min(cutbuild);
cutoffs(posicut)=cutmin;
cutscs=[posicut];
%Rest of rounds
for ronext=2:nousc
for number=1:nousc
locscne=ordnext(:,number);
for st=1:nostrane
cutdemne(st)=max(cutoffs(permnext(st,1:locscne(st))));
end
locrepe=locscne(:,ones(1,ronext-1));
cutahead=(locrepe-ordnext(:,cutscs)>0);
criterio=(locscne-sum(cutahead,2)==1);
nume=remnext(underdem(number))+cutdemne'*(demnext.*criterio);
deno=demnext'*criterio;
cutbuild(number)=min([1 nume/max([.0000001 deno])]);
end
cutbuild=cutbuild+50000*(cutoffs>0);
[cutmin,posicut]=min(cutbuild);
cutoffs(posicut)=cutmin;
cutscs=[cutscs posicut];
end
end
%% Building assignment probabilities
% Reordering cutoffs
cutnext=zeros(noscho,1);
for number=1:nousc
cutnext(underdem(number))=cutoffs(number);
end
% Probabilities for each strategy
probtar=prfirst(permu(:,1));
probalca=zeros(nostra,noscho);
probalca(:,1)=probtar;
prtarrep=probtar(:,ones(1,noscho-1));
probtool=zeros(nostra,1);
for number=2:noscho
inccut=cutnext(permu(:,number))-probtool;
probalca(:,number)=(inccut>0).*inccut;
probtool=probtool+(inccut>0).*inccut;
end
probalca(:,2:noscho)=(1-prtarrep).*probalca(:,2:noscho);
%Rearrange probabilities (now prob(i,j): strategy i,
%school ranked j-th (in permu), want j to identify school.
ordprca=zeros(nostra,noscho);
for st3=1:nostra
prost=probalca(st3,:);
ordprca(st3,:)=prost(ordperm(st3,:));
end

%%Best Responses and Welfare
payca=ordprca*rescval'+biastru;
[maxpayca,bestca]=max(payca);
reveal1=permu(bestca',:);

%%%If naive mode...
if naive==1;
reveal1(1:nonaive,:)=reveal0(1:nonaive,:);
for naistu=1:nonaive
ranktool=ranking(naistu,:);
rankrepe=ranktool(ones(nostra,1),:);
location=find(sum((permu-rankrepe==0),2)==noscho);
maxpayca(naistu)=ordprca(location,:)*rescval(naistu,:)';
end
end

welfaca(aid,naimode,morescho,morecapa,iter)=mean(maxpayca);
welnaica(aid,naimode,morescho,morecapa,iter)=mean(maxpayca(1:nonaive));

%%Check if best response is equilibrium
coincid(count)=length(find(reveal1-reveal0==0));
if coincid(count)==noscho*nostud
break
end

% Auxilliary converge criterion
if (count>nocount/2) && (coincid(count)>=(.98+naive*.01)*noscho*nostud)
break
end

%% "Smoothing" to avoid cycles
lott=rand(nostud,1);
keep0=(lott<.975);
reveal1=reveal1+keep0(:,ones(1,noscho)).*(reveal0-reveal1);

%%Finally...
count=count+1
end
alertca(aid,naimode,morescho,morecapa,iter)=(count==nocount)*coincid(nocount-1)/(noscho*nostud);



end
end
end
end
end

%%Average Welfares
avewelbm=mean(welfabm,5);
avewelca=mean(welfaca,5);
avewelda=mean(welfada,5);
avewnbm=mean(welnaibm,5);
avewnca=mean(welnaica,5);
avewnda=mean(welnaida,5);
save('results.mat','ave*','ale*') 


SECOND CODE: WEAK PRIORITIES
clear all;
rand('state',sum(100*clock));
%%%%%%%%%%%%%%%%%%%%%%%SIMULATION BM-RTB DA WEAK PRIORITIES%%%%%%%%%%%%%
%%Parameters
nocount=4001;
noas=2;
mina=-.2;
inca=.4;
corre=.2;
noiter=25;
minnosc=3;
varscho=3;
mincapa=10;
varcapa=3;
nodraw=500;
%% Welfare and alert matrices
welfabm=zeros(noas,2,2,varscho,varcapa,noiter);
alertbm=zeros(noas,2,2,varscho,varcapa,noiter);
welfaca=zeros(noas,2,2,varscho,varcapa,noiter);
alertca=zeros(noas,2,2,varscho,varcapa,noiter);
welfada=zeros(noas,2,2,varscho,varcapa,noiter);
welnaibm=zeros(noas,2,2,varscho,varcapa,noiter);
welnaida=zeros(noas,2,2,varscho,varcapa,noiter);
%%%%%%%%%%%%%FOR EACH NUMBER OF SCHOOLS

for morescho=1:varscho
noscho=minnosc+morescho; %Regarded as J sometimes

%%%%%%%%%%%%%FOR EACH SCHOOL CAPACITY

for morecapa=1:varcapa
capac=mincapa+10*morecapa;

nostra=factorial(noscho);
nostud=noscho*capac;
nonaive=nostud/2; %%When naive mode is yes, half students naive

%%Creation permutation matrix (J! rows, J columns)
%%It is the strategy space
permu=perms(1:noscho);
[ordprob,ordperm]=sort(permu,2);
%ordprob just defines some matrix with proper dimensions I later use
%ordperm(i,j): in what position strategy i locates school j

%%Generation random tiebreakers for DA
tiebreak=rand(nostud,nodraw);


%%%%%%%%%%%%%%%%%%%%%%%%%%FOR EACH ITERATION
for iter=1:noiter;

%%Generation valuations: vij=a·uj+(1-a)·uij
commval=rand(1,noscho);
commval2=commval(ones(nostud,1),:);
privval=rand(nostud,noscho);

%%Generation priorities
priolott=rand(nostud,noscho);
priority=(priolott<1/noscho); %the (independent) probability of having
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%priority at j is 1/noscho


%%%%%%%%%%%%%%%%%%%%%%%%%%FOR EACH CORRELATION-WITH-PRIORITY MODE
for corrmode=1:2
c=(corrmode-1)*corre;

%%%%%%%%%%%%%%%%%%%%%%%%%%%% FOR EACH VALUE OF a...
for aid=1:noas
a=mina+inca*aid;
vals=(a*commval2+(1-a)*privval)*(1-c)+c*priority;

%%Valuations rescaled to the simplex
rescval=zeros(nostud,noscho);
for stu=1:nostud
rescval(stu,:)=vals(stu,:)/sum(vals(stu,:),2);
end

%%True ordinal preferences
[hale,ranking]=sort(rescval,2,'descend');
%Note: ranking(i,j), school student i puts in j-th postion

%% Bias for the truth
truth=zeros(nostra,nostud);
for estu=1:nostud
truth(:,estu)=sum((permu==ranking(estu*ones(nostra,1),:)),2);
end
biastru=.0000001*truth;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%DA%%%%%%%%%%%%%%%%%%%%%%
%%NOTE: DA does not depend on naive mode

allocda=zeros(nostud,nodraw);
probalda=zeros(nostud,noscho);
for draw=1:nodraw
DRAW=draw
%ordering
orderno=zeros(nostud,noscho);
orden=zeros(nostud,noscho);
for school=1:noscho
orderno(:,school)=tiebreak(:,draw)+1.00001*(1-priority(:,school));
end
[whatever,orden]=sort(orderno);
%assignment
remda=capac*ones(1,noscho);
applic=ones(nostud,1);
asign=zeros(nostud,1);
while sum(remda)>0
expan=zeros(nostud,noscho);
remda=capac*ones(1,noscho);
for sc=1:noscho
for pos=1:nostud
who=orden(pos,sc);
applies=(ranking(who,applic(who))==sc);
expan(who,sc)=sc*applies*(remda(sc)>0);
remda(sc)=remda(sc)-(expan(who,sc)>0);
end
end
asign=sum(expan,2);
applic=applic+(asign==0);
end
%assignment given draw
allocda(:,draw)=asign;
end
%random assignment probalda(i,j): student i, school j
for studda=1:nostud
for escuda=1:noscho
probalda(studda,escuda)=length(find(allocda(studda,:)==escuda))/nodraw;
end
end

%Computar welfare
payda=diag(probalda*rescval');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%NAIVE MODE NO/YES
for naimode=1:2
% 1=no naive; 2=naive
naive=naimode-1;

%%Welfare DA
welfada(aid,corrmode,naimode,morescho,morecapa,iter)=mean(payda);
welnaida(aid,corrmode,naimode,morescho,morecapa,iter)=mean(payda(1:nonaive));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%BOSTON MECHANISM RTB%%%%%%%%%%%%%%%%%%%
coincid=zeros(nocount-1,1);
reveal1=ranking;
maxpay=zeros(nostud,1);
count=1
while count<nocount
reveal0=reveal1;

%%Calculation probabilities per round
probprio=zeros(noscho,noscho);
probnopr=zeros(noscho,noscho);
unass=ones(nostud,1);
assprob=zeros(nostud,1);
avail=capac*ones(noscho,1);
for ronda=1:noscho;
for esc=1:noscho;
pediprio=unass'*((reveal0(:,ronda)==esc).*priority(:,esc));
pedinopr=unass'*((reveal0(:,ronda)==esc).*(1-priority(:,esc)));
probprio(esc,ronda)=min([1 avail(esc)/max([0.000000000001 pediprio])]);
excess=max([0 avail(esc)-pediprio]);
probnopr(esc,ronda)=min([1 excess/max([0.000000000001 pedinopr])]);
avail(esc)=max([0 excess-pedinopr]);

assprob=(reveal0(:,ronda)==esc).*(priority(:,esc)*probprio(esc,ronda)+(1-priority(:,esc))*probnopr(esc,ronda));
unass=unass.*(1-assprob);
end
end

%%Assignment probabilities and welfare, for each sutdent
for stu=1:nostud
prob=zeros(nostra,noscho);
probunas=ones(nostra,1);
for colu=1:noscho
haveprio=priority(stu,permu(:,colu))';

prob(:,colu)=probunas.*(probprio(permu(:,colu),colu).*haveprio+(1-haveprio).*probnopr(permu(:,colu),colu));
probunas=probunas-prob(:,colu);
end
%Rearrange probabilities (now prob(i,j): strategy i,
%school ranked j-th (in permu), want j to identify school.
ordprob=zeros(nostra,noscho);
for st3=1:nostra
prost=prob(st3,:);
ordprob(st3,:)=prost(ordperm(st3,:));
end

payoffst=ordprob*rescval(stu,:)'+biastru(:,stu);

%%%If naive mode...
if (naive==1) && (stu<=nonaive);
reveal1(stu,:)=ranking(stu,:);
ranktool=ranking(stu,:);
rankrepe=ranktool(ones(nostra,1),:);
location=find(sum((permu-rankrepe==0),2)==noscho);
maxpay(stu)=payoffst(location);
else
%%If not, best response and welfare for this student
[maxpayst,best]=max(payoffst);
maxpay(stu)=maxpayst;
reveal1(stu,:)=permu(best,:);
end
end

welfabm(aid,corrmode,naimode,morescho,morecapa,iter)=mean(maxpay);
welnaibm(aid,corrmode,naimode,morescho,morecapa,iter)=mean(maxpay(1:nonaive));

%%Check if best response is equilibrium
coincid(count)=length(find(reveal1-reveal0==0));
if coincid(count)==noscho*nostud
break
end

% Auxilliary converge criterion
if (count>nocount/2) && (coincid(count)>=(.98+naive*.01)*noscho*nostud)
break
end

%% "Smoothing" to avoid cycles
lott=rand(nostud,1);
keep0=(lott<.975);
reveal1=reveal1+keep0(:,ones(1,noscho)).*(reveal0-reveal1);

%%Finally...
count=count+1
end
alertbm(aid,corrmode,naimode,morescho,morecapa,iter)=(count==nocount)*coincid(nocount-1)/(noscho*nostud);

end
end
end
end
end
end

%%Average Welfares
avewelbm=mean(welfabm,6);
avewelda=mean(welfada,6);
avewnbm=mean(welnaibm,6);
avewnda=mean(welnaida,6);
save('resprio.mat','ave*','ale*')