%December 6, 2003
%Author: Francisco M. Torralba
%Description: This matlab program simulates a 10,000-period history of x(t), as described by "Ecological rationality, Bayesian learning and rational expectations" 
clc;
clear all;
P=[.6 .2 .2;              %transition probability matrix
    .4 .3 .3
    .3 .4 .3];
P1=zeros(3,3);
P2=zeros(3,3);
q=0;                      %proportion of agents with good estimates of P 

T=10000;                  %length of the simulation
J=100;                    %number of workers
k=100;                    %length of observed history
alpha=0.01;               %fraction of the population which updates their decision every period (for the rational expectations case)
s0=3;                     %initial state 
m2=0.8;                   %popularity weight   
V=[-1 1 3];               %values of the stochastic payoff in sector A
y=zeros(T,1);             %stores the history of the stochastic payoff
h=round(rand(T,1)*100+1); %generates the sequence of individuals who are allowed to switch
predictedy=zeros(T+1,1);
expectedy=zeros(T+1,1);
x=zeros(T+1,1);
z=zeros(T+1,1);
u=zeros(T+1,1);
y1=zeros(100,1);
y2=zeros(10,1);
gamma=1/3;
b=2/3;

y=markov(P,T+1,s0,V);      %generates the sequence of stochastic payoffs (see markov.m)
y1=markov(P,101,s0,V)';    
                            %stores the initial history of stochastic payoffs observed by each agent the first time they are allowed to switch
y2=markov(P,11,s0,V)';
tran1=zeros(3,3);
tran2=zeros(3,3);

%Generating initial histories for Bayesian learners

%for type 1
for t=1:99;
  %If the state is 1
  if y1(t)==V(1);
     if y1(t+1)==V(1);
         tran1(1,1)=tran1(1,1)+1;
     elseif y1(t+1)==V(2);
             tran1(1,2)=tran1(1,2)+1;
     else tran1(1,3)=tran1(1,3)+1;
     end        
  else g=1;
  end
  %If the state is 2
  if y1(t)==V(2);
     if y1(t+1)==V(1);
         tran1(2,1)=tran1(2,1)+1;
     elseif y1(t+1)==V(2);
             tran1(2,2)=tran1(2,2)+1;
     else tran1(2,3)=tran1(2,3)+1;
     end
  else g=1;
  end
  %If the state is 3
  if y1(t)==V(3);
     if y1(t+1)==V(1);
         tran1(3,1)=tran1(3,1)+1;
     elseif y1(t+1)==V(2);
             tran1(3,2)=tran1(3,2)+1;
     else tran1(3,3)=tran1(3,3)+1;
     end
  else g=1;
  end
end
for r=1:3;
    for c=1:3;
    P1(r,c)=tran1(r,c)/sum(tran1(r,:));
    end
end
tran2=zeros(3,3);
%for type 2 
for t=1:9;
  %If the state is 1
  if y2(t)==V(1);
     if y2(t+1)==V(1);
         tran2(1,1)=tran2(1,1)+1;
     elseif y2(t+1)==V(2);
         tran2(1,2)=tran2(1,2)+1;
     else tran2(1,3)=tran2(1,3)+1;
     end
  else c=1;
  end
  %If the state is 2
  if y2(t)==V(2);
     if y2(t+1)==V(1);
         tran2(2,1)=tran2(2,1)+1;
     elseif y2(t+1)==V(2);
         tran2(2,2)=tran2(2,2)+1;
     else tran2(2,3)=tran2(2,3)+1;
     end
  else c=1;
  end
  %If the state is 3
  if y2(t)==V(3);
     if y2(t+1)==V(1);
         tran2(3,1)=tran2(3,1)+1;
     elseif y2(t+1)==V(2);
         tran2(3,2)=tran2(3,2)+1;
     else tran2(3,3)=tran2(3,3)+1;
     end
  else c=1;
  end
end  
if tran2(1,1)+tran2(1,2)+tran2(1,3)==0;
    P2(1,1)=0;
    P2(1,2)=0;
    P2(1,3)=0;
else
P2(1,1)=tran2(1,1)/(tran2(1,1)+tran2(1,2)+tran2(1,3));
P2(1,2)=tran2(1,2)/(tran2(1,1)+tran2(1,2)+tran2(1,3));
P2(1,3)=tran2(1,3)/(tran2(1,1)+tran2(1,2)+tran2(1,3));
end
if tran2(2,1)+tran2(2,2)+tran2(2,3)==0;
    P2(2,1)=0;
    P2(2,2)=0;
    P2(2,3)=0;
else
P2(2,1)=tran2(2,1)/(tran2(2,1)+tran2(2,2)+tran2(2,3));
P2(2,2)=tran2(2,2)/(tran2(2,1)+tran2(2,2)+tran2(2,3));
P2(2,3)=tran2(2,3)/(tran2(2,1)+tran2(2,2)+tran2(2,3));
end
if tran2(3,1)+tran2(3,2)+tran2(3,3)==0;
    P2(3,1)=0;
    P2(3,2)=0;
    P2(3,3)=0;
else
P2(3,1)=tran2(3,1)/(tran2(3,1)+tran2(3,2)+tran2(3,3));
P2(3,2)=tran2(3,2)/(tran2(3,1)+tran2(3,2)+tran2(3,3));
P2(3,3)=tran2(3,3)/(tran2(3,1)+tran2(3,2)+tran2(3,3));
end

Q=P1;
Tran=tran1;
Pest=zeros(3,T*3);
for j=2:101;
    if j<=q*J;
        Qadd=P1;
        tranadd=tran1;
    else Qadd=P2;
         tranadd=tran2; 
    end
    Q=[Q Qadd];
    Tran=[Tran tranadd];
end

x=zeros(T,1);              %Stores the evolution of x under rational expectations
x(1)=0.3;
expectedy=zeros(T+1,1);    %Stores the history of expected y under rational expectations
expectedy(1)=0;
z=zeros(T,1);
w1=zeros(T+1,1);
z(1)=x(1);
u=zeros(T,1);             %Store the evolution of x under the heuristic rule  
u(1)=x(1);
predictedy=zeros(T,1);    %Stores the history of expected y under Bayesian learning
v=zeros(T,1);
v(1)=.2532;
for t=1:(T-1);
    v(t+1)=.2532;            %Expectation of x(t) under rational expectations
    %Rational expectations
        if y(t)==V(1);
        expectedy(t+1)=P(1,:)*V';
            if expectedy(t+1)>=exp(gamma*(x(t)-b));
               x(t+1)=(1-alpha)*x(t)+alpha;
            else
               x(t+1)=(1-alpha)*x(t);
            end
    elseif y(t)==V(2);
        expectedy(t+1)=P(2,:)*V';
            if expectedy(t+1)>=exp(gamma*(x(t)-b));
               x(t+1)=(1-alpha)*x(t)+alpha;
            else
               x(t+1)=(1-alpha)*x(t);
            end
    else
        expectedy(t+1)=P(3,:)*V';
            if expectedy(t+1)>=exp(gamma*(x(t)-b));
               x(t+1)=(1-alpha)*x(t)+alpha;
            else
               x(t+1)=(1-alpha)*x(t);
            end
    end

    %Rule-of-thumb learning
      
    predictedy2(t+1)=y(t);
    
    if m2*(2*u(t)-1)+(1-m2)*(predictedy2(t+1)-exp(gamma*(u(t)-b)))>=0;
        u(t+1)=(1-alpha)*u(t)+alpha;
    else
        u(t+1)=(1-alpha)*u(t);
    end
      
    %Learners 
    Pest(:,t*3-2:t*3)=Q(:,((h(t)*3-2):(h(t)*3)));
       
    if y(t)==V(1);
        predictedy(t+1)=Pest(1,t*3-2:t*3)*V';
            if (predictedy(t+1)-exp(gamma*(z(t)-b)))>=0;
               z(t+1)=(1-alpha)*z(t)+alpha;
            else
               z(t+1)=(1-alpha)*z(t);
            end
            %Update the estimate of the transition probability matrix for individual h(t)
            if y(t+1)==V(1);
            Tran(1,((h(t)*3-2)))=Tran(1,((h(t)*3-2)))+1;
            Q(1,((h(t)*3-2)))=Tran(1,((h(t)*3-2)))/sum(Tran(1,((h(t)*3-2)):((h(t)*3))));
            elseif y(t+1)==V(2);
            Tran(1,((h(t)*3-1)))=Tran(1,((h(t)*3-1)))+1;
            Q(1,((h(t)*3-1)))=Tran(1,((h(t)*3-1)))/sum(Tran(1,((h(t)*3-2)):((h(t)*3))));
            else Tran(1,((h(t)*3)))=Tran(1,((h(t)*3)))+1;
            Q(1,((h(t)*3)))=Tran(1,((h(t)*3)))/sum(Tran(1,((h(t)*3-2)):((h(t)*3))));    
            end   
          
                
    elseif y(t)==V(2);
        predictedy(t+1)=Pest(2,t*3-2:t*3)*V';
            if (predictedy(t+1)-exp(gamma*(z(t)-b)))>=0;
               z(t+1)=(1-alpha)*z(t)+alpha;
            else
               z(t+1)=(1-alpha)*z(t);
            end
            %Update the estimate of the transition probability matrix for individual h(t)
            if y(t+1)==V(1);
            Tran(2,((h(t)*3-2)))=Tran(2,((h(t)*3-2)))+1;
            Q(2,((h(t)*3-2)))=Tran(2,((h(t)*3-2)))/sum(Tran(2,((h(t)*3-2)):((h(t)*3))));
            elseif y(t+1)==V(2);
            Tran(2,((h(t)*3-1)))=Tran(2,((h(t)*3-1)))+1;
            Q(2,((h(t)*3-1)))=Tran(2,((h(t)*3-1)))/sum(Tran(2,((h(t)*3-2)):((h(t)*3))));
            else Tran(2,((h(t)*3)))=Tran(2,((h(t)*3)))+1;
            Q(2,((h(t)*3)))=Tran(2,((h(t)*3)))/sum(Tran(2,((h(t)*3-2)):((h(t)*3))));
            end
                
            
    else
        predictedy(t+1)=Pest(3,t*3-2:t*3)*V';
            if (predictedy(t+1)-exp(gamma*(z(t)-b)))>=0;
               z(t+1)=(1-alpha)*z(t)+alpha;
            else
               z(t+1)=(1-alpha)*z(t);
            end
            %Update the estimate of the transition probability matrix for individual h(t)
            if y(t+1)==V(1);
            Tran(3,((h(t)*3-2)))=Tran(3,((h(t)*3-2)))+1;
            Q(3,((h(t)*3-2)))=Tran(3,((h(t)*3-2)))/sum(Tran(3,((h(t)*3-2)):((h(t)*3))));
            elseif y(t+1)==V(2);
            Tran(3,((h(t)*3-1)))=Tran(3,((h(t)*3-1)))+1;
            Q(3,((h(t)*3-1)))=Tran(3,((h(t)*3-1)))/sum(Tran(3,((h(t)*3-2)):((h(t)*3))));
            else Tran(3,((h(t)*3)))=Tran(3,((h(t)*3)))+1;
            Q(3,((h(t)*3)))=Tran(3,((h(t)*3)))/sum(Tran(3,((h(t)*3-2)):((h(t)*3))));    
            end   
               
  end
  end 
 
end   

plot(v,'--')
hold on;
plot(x,'blue')
hold on;
plot(z,'red')
ylim([0 1])
plot(u,'green')
set(gca,'YTick',0:0.1:1)
set(gca,'Xtick',0:1000:T)
set(gca,'YTickLabel',{'0','0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9','1'})

hold off;