clear all;
% cooperation: 0
% defect: 1
B=[1,0.8,0.7,0.65]% hyperbolic discount factors
Q=[0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1]% probability of forgiveness
N=200000;  % number of games
A=50; Z=12; D=25; AC=40;% payoff matrix
deltaE=0.5;

%in the model, the expected value of R(i+1) is (1-s)R(i)+AC*s
R(1)=AC; sen=0.1; DEV=8;
for i=2:N;
    r1(i)=rand;
    r2(i)=rand;
    R(i)=R(i-1)+sen*((AC+DEV)-R(i-1))*r1(i) - sen*(R(i-1)-(AC-DEV))*r2(i) ;
    
end

for bi=1:length(B);
    beta=B(bi);
    for qi=1:length(Q)
        q=Q(qi);
        beta
        q

    
    delta=deltaE/(beta + deltaE*(1-beta));




% threthold value:
WA = @(d,C)  (A*d*(1-d)*(1-q)^2+C*q*(1+d*(1-q))+Z*(1-d)*(1-q)) / ((1-d)*(1-d^2*(1-q)^2));
LE= @(d,C) C + beta*d*C/(1-d);
RE = @(d,C) A + beta*d*WA(d,C);
d=delta;
fun = @(C) LE(d,C)-RE(d,C);
x = fzero(fun,40);
%DR=sum(R < x)/length(R)
%De=(R<x);

r(1,1)=0;r(2,1)=0;
s(1,1)=0;s(2,1)=0;
payoff(1,1) = R(1); payoff(2,1) = R(1); u=payoff;
for i=2:N;
    
    
  % player 1's strategy
    if s(2,i-1)==0;
        if (1-s)*R(i)+40*s > x; s(1,i)=0;
        else s(1,i)=1; end
    else
        if rand<q; s(1,i)=0; %forgive
        else s(1,i)=1; end  % not forgive
    end
    
    % player 2's strategy
    if s(1,i-1)==0;
        s(2,i)=0;
    else
        if rand<q; s(2,i)=0; %forgive
        else s(2,i)=1; end  % not forgive
    end
   
  %  payoff(1,i) = (s(1,i)==0 & s(2,i)==0 )*R(i) + (s(1,i)==0 & s(2,i)==1 )*Z ...
   %     + (s(1,i)==1 & s(2,i)==0 )*A +(s(1,i)==1 & s(2,i)==1 )*D;
    %payoff(2,i) = (s(2,i)==0 & s(1,i)==0 )*AC + (s(2,i)==0 & s(1,i)==1 )*Z ...
     %   + (s(2,i)==1 & s(1,i)==0 )*A +(s(2,i)==1 & s(1,i)==1 )*D;
    
   % u(1,i)=(u(1,i-1)*(i-1)+payoff(1,i))/i;
   % u(2,i)=(u(2,i-1)*(i-1)+payoff(2,i))/i;
    
   % C1(i)=sum(s(1,:)==0)/length(s(1,:));
   % C2(i)=sum(s(2,:)==0)/length(s(2,:));
   %  P1(i)=mean(payoff(1,:));
    
end
%plot(C1);axis([0 N 0.9, 1])
%plot(P1);axis([0 N 30, 40])

%mean(payoff,2)
%sum(s(1,:)==0)/length(s(1,:))



MC(bi,qi)=sum(s(1,:)==0 & s(2,:)==0)/length(s(1,:))
    end
end

MC=MC.*100;

plot(Q,MC,'--s','LineWidth',1.5,...
                       'MarkerEdgeColor','k',...
                       'MarkerFaceColor','g',...
                       'MarkerSize',5);
                   xlabel('Forgiveness level (q)')
                   b=[cellstr(num2str(get(gca,'xtick')'*100))];   pct = char(ones(size(b,1),1)*'%'); 
                  new_xticks = [char(b),pct];set(gca,'xticklabel',new_xticks);
                   
                  ylabel('Cooperation rate')
                   
                  a=[cellstr(num2str(get(gca,'ytick')'))];   pct = char(ones(size(a,1),1)*'%'); 
                  new_yticks = [char(a),pct];set(gca,'yticklabel',new_yticks);
                  

txt = ['\beta = ',num2str(B(1))]; text(Q(length(Q)),MC(1,length(Q)),txt,'FontSize',16);
txt = ['\beta = ',num2str(B(2))]; text(Q(length(Q)),MC(2,length(Q)),txt,'FontSize',16);
txt = ['\beta = ',num2str(B(3))]; text(Q(length(Q)),MC(3,length(Q)),txt,'FontSize',16);
txt = ['\beta = ',num2str(B(4))]; text(Q(length(Q)),MC(4,length(Q)),txt,'FontSize',16)
       
                   

% Data
%MC=[98.3170   98.6495   98.7130   98.7295   98.5600   98.2350   97.8660   97.5440   96.5465   96.0455
%   90.2430   92.7655   92.3810   91.8635   91.0535   90.1285   88.9710   87.1785   85.4520   83.4115
%   74.9190   81.9195   83.8945   83.0730   81.6785   80.2900   77.7905   75.6035   72.2450   69.6645
%   66.3000   73.9740   75.9250   75.9300   75.2315   73.1125   69.8690   67.6670   64.0055   60.2590]
%