Code for Bear, Kagan, & Rand (2017)

AA_BearKaganRand_master_file.m

%   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %
%                 CODE FOR BEAR, KAGAN & RAND (2017) PROCB
%                                                                       
% This MATLAB script runs the low mutation limit (steady state) calculations     
% shown in Figures S1-S12
%                                                                       
% This script calls the following functions:  
%
% calculate_payoff_YDYI(strategy1,strategy2,b,c,d,p,YD,YI)
%       - returns the payoff of strategy1 against strategy2 given   
%         the specified values of b, c, d, p, YD, and YI
%
% steady_state_calc(payoffTable,N,w,a)
%       - returns the steady state distribution of the population in
%         the limit of low mutation, where payoffTable is an MxM matrix
%         with payoffTable(i,j) being the payoff of strategy i against 
%         strategy j; N is the population size; w is the selection 
%         strength; and a is the amount of assortment
%
% vary_YD(N,b,c,d,w,a)
%       - returns and plots the average results from steady_state_calc 
%         when varying YD and p (with YI fixed at .5)
%
% vary_YI(N,b,c,d,w,a)
%       - returns and plots the average results from steady_state_calc 
%         when varying YI and p (with YD fixed at 1)
%   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %   %


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% STEADY STATE CALCULATION %% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% parameters used in Figures S1-S12
N = 50; 
w = 3;
a = 0;

bs = [4,2,8,4,4,4,4];
cs = [1,1,1,2,.5,1,1];
ds = [1,1,1,1,1,.75,2];


for i=1:length(bs) 
    vary_YD(N,b,c,d,w,a);
    vary_YI(N,b,c,d,w,a);
end
% Note: These functions may take a long time to run.

calculate_payoff_YDYI.m

function payoff = calculate_payoff_YDYI(s1,s2,b,c,d,p,YD,YI)
%% SETUP %%

% s1 and s2 are agent 1's and agent 2's respective strategy profiles
%   where a strategy profile consists of [pR,p1,pIR,pI1,T]

% STRATEGY PARAMETERS %
% pR is probability of cooperating when agent deliberates and realizes the
% game is repeated/coordination game
% p1 is probability of cooperating when agent deliberates and realizes the
% game is 1-shot/social dilemma
% pIR is probability of cooperating when agent uses intuition in a repeated
% game
% pI1 is probability of cooperating when agent uses intuition in a 1 shot
% game
% T is (threshold) max cost in [0,d] agent is willing to pay to deliberate

% ENVIRONMENT/GAME PARAMETERS % 
% b is benefit of cooperation
% c is cost of cooperation
% d is max possible cost of deliberation (such that the cost of
%   deliberation varies uniformly in the interval [0,d])
% p is probability of repeated game
% YD is the sensitivity of deliberation 
% YI is the sensitivity of intuition

% OUTPUTS %
% payoff is s1's payoff


%%%%%%%%%%%%%%%%%%%%%%%
%% CALCULATE PAYOFFS %%
%%%%%%%%%%%%%%%%%%%%%%%

% Rename each agent's threshold to T1 and T2, respectively
T1 = s1(5);
T2 = s2(5);

%%%%
% Calculate repeated game payoffs

% Repeated Game Payoff Matrix 
repPayoffMat=[0	0
              0 (b-c)];
          
%   Rows correspond to agent 1's strategy 
%       where row 1 = agent 1 defects, row 2 = agent 1 cooperates
%   Columns correspond to agent 2's strategy
%       where col 1 = agent 2 defects, col 2 = agent 2 cooperates
%   Values correspond to agent 1's payoffs
%       e.g., if agent 1 cooperates and agent 2 defects, 
%       agent 1 gets a payoff of 0 because repPayoffMat(2,1) = 0 
 
% We consider repeated-game payoffs in 4 situations:
%   rPayDD: both agents deliberate 
%   rPayDI: only agent 1 deliberates
%   rPayID: only agent 2 deliberates
%   rPayII: neither agent deliberates

%Given the possibility of  innaccurate deliberation and sensitivity of intuition
%   rPayDD:
%     rPayDtDt: Both agents deliberate, both deliberations track game type
%     rPayDtDf: Both agents deliberate, agent 1's deliberation tracks, agent
%                 2's fails
%     rPayDfDt: Both agents deliberate, agent 1's deliberation fails,
%     agent 2's tracks
%     rPayDfDf: Both agents deliberate, both deliberations fail
%   rPayDI:
%     rPayDtIt: 1 deliberates (tracks) 1 intuits (tracks)
%     rPayDtIf: 1 deliberates (tracks) 1 intuits (fails)
%     rPayDfIt: 1 deliberates (fails) 1 intuits (tracks)
%     rPayDfIf: 1 deliberates (fails) 1 intuits (fails)
%   rPayID:
%     rPayItDt: agent 2 deliberates (tracks), agent 1 intuits (tracks)
%     rPayIfDt: agent 2 deliberates (tracks), agent 1 intuits (fails)
%     rPayItDf: agent 2 deliberates (fails), agent 1 intuits (tracks)
%     rPayIfDf: agent 2 deliberates (fails), agent 1 intuits (fails)
%   rPayII:
%     rPayItIt: both intuit (both track)
%     rPayItIf: both intuit (one tracks, 1 fails)
%     rPayIfIt: both intuit (one tracks, 1 fails)
%     rPayIfIf: both intuit (both fail)

% E.g., to calculate rPayDfIf, 
%   We know agent 1 deliberates but fails to reailze they're in a repeated game => 
%       their prob. of cooperating is p1, i.e., s1(2)
%   We know agent 2 uses their intuitive response but fails to realize they're in a repeated game =>
%       their prob. of cooperating is pI1, i.e., s2(4)
%   Thus, rPayDfIf is a weighted average of the payoff agent 1 gets when both
%       cooperate (repPayoffMat(2,2)), when only player 1 cooperates
%       (repPayoffMat(2,1)), and so on
%   where the prob. both cooperate is given by s1(2)*s2(4), prob. only
%   player 1 cooperates is given by s1(2)*(1-s2(4)), and so on

rPayDtDt=s1(1)*s2(1)*repPayoffMat(2,2) + s1(1)*(1-s2(1))*repPayoffMat(2,1) + ...
    (1-s1(1))*s2(1)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(1))*repPayoffMat(1,1);
rPayDtDf=s1(1)*s2(2)*repPayoffMat(2,2) + s1(1)*(1-s2(2))*repPayoffMat(2,1) + ...
    (1-s1(1))*s2(2)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(2))*repPayoffMat(1,1);
rPayDfDt=s1(2)*s2(1)*repPayoffMat(2,2) + s1(2)*(1-s2(1))*repPayoffMat(2,1) + ...
    (1-s1(2))*s2(1)*repPayoffMat(1,2) + (1-s1(2))*(1-s2(1))*repPayoffMat(1,1);
rPayDfDf=s1(2)*s2(2)*repPayoffMat(2,2) + s1(2)*(1-s2(2))*repPayoffMat(2,1) + ...
    (1-s1(2))*s2(2)*repPayoffMat(1,2) + (1-s1(2))*(1-s2(2))*repPayoffMat(1,1);

rPayDtIt=s1(1)*s2(3)*repPayoffMat(2,2) + s1(1)*(1-s2(3))*repPayoffMat(2,1) + ...
    (1-s1(1))*s2(3)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(3))*repPayoffMat(1,1);
rPayDtIf=s1(1)*s2(4)*repPayoffMat(2,2) + s1(1)*(1-s2(4))*repPayoffMat(2,1) + ...
    (1-s1(1))*s2(4)*repPayoffMat(1,2) + (1-s1(1))*(1-s2(4))*repPayoffMat(1,1);
rPayDfIt=s1(2)*s2(3)*repPayoffMat(2,2) + s1(2)*(1-s2(3))*repPayoffMat(2,1) + ...
    (1-s1(2))*s2(3)*repPayoffMat(1,2) + (1-s1(2))*(1-s2(3))*repPayoffMat(1,1);
rPayDfIf=s1(2)*s2(4)*repPayoffMat(2,2) + s1(2)*(1-s2(4))*repPayoffMat(2,1) + ...
    (1-s1(2))*s2(4)*repPayoffMat(1,2) + (1-s1(2))*(1-s2(4))*repPayoffMat(1,1);

rPayItDt=s1(3)*s2(1)*repPayoffMat(2,2) + s1(3)*(1-s2(1))*repPayoffMat(2,1) + ...
    (1-s1(3))*s2(1)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(1))*repPayoffMat(1,1);
rPayIfDt=s1(4)*s2(1)*repPayoffMat(2,2) + s1(4)*(1-s2(1))*repPayoffMat(2,1) + ...
    (1-s1(4))*s2(1)*repPayoffMat(1,2) + (1-s1(4))*(1-s2(1))*repPayoffMat(1,1);
rPayItDf=s1(3)*s2(2)*repPayoffMat(2,2) + s1(3)*(1-s2(2))*repPayoffMat(2,1) + ...
    (1-s1(3))*s2(2)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(2))*repPayoffMat(1,1);
rPayIfDf=s1(4)*s2(2)*repPayoffMat(2,2) + s1(4)*(1-s2(2))*repPayoffMat(2,1) + ...
    (1-s1(4))*s2(2)*repPayoffMat(1,2) + (1-s1(4))*(1-s2(2))*repPayoffMat(1,1);

rPayItIt=s1(3)*s2(3)*repPayoffMat(2,2) + s1(3)*(1-s2(3))*repPayoffMat(2,1) + ...
    (1-s1(3))*s2(3)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(3))*repPayoffMat(1,1);
rPayItIf=s1(3)*s2(4)*repPayoffMat(2,2) + s1(3)*(1-s2(4))*repPayoffMat(2,1) + ...
    (1-s1(3))*s2(4)*repPayoffMat(1,2) + (1-s1(3))*(1-s2(4))*repPayoffMat(1,1);
rPayIfIt=s1(4)*s2(3)*repPayoffMat(2,2) + s1(4)*(1-s2(3))*repPayoffMat(2,1) + ...
    (1-s1(4))*s2(3)*repPayoffMat(1,2) + (1-s1(4))*(1-s2(3))*repPayoffMat(1,1);
rPayIfIf=s1(4)*s2(4)*repPayoffMat(2,2) + s1(4)*(1-s2(4))*repPayoffMat(2,1) + ...
    (1-s1(4))*s2(4)*repPayoffMat(1,2) + (1-s1(4))*(1-s2(4))*repPayoffMat(1,1);

%Thus: 
rPayDD=YD*YD*rPayDtDt + YD*(1-YD)*rPayDtDf + (1-YD)*YD*rPayDfDt + (1-YD)*(1-YD)*rPayDfDf;
rPayDI=YD*YI*rPayDtIt + (1-YD)*YI*rPayDfIt + YD*(1-YI)*rPayDtIf +(1-YD)*(1-YI)*rPayDfIf;
rPayID=YI*YD*rPayItDt + (1-YI)*(1-YD)*rPayIfDf + (1-YI)*YD*rPayIfDt +YI*(1-YD)*rPayItDf;
rPayII=YI*YI*rPayItIt + (1-YI)*YI*rPayIfIt + YI*(1-YI)*rPayItIf + (1-YI)*(1-YI)*rPayIfIf;

%%%
% Calculate 1-shot game payoffs
%   (these work the same way as the repeated game payoffs)

% 1-shot Game Payoff Matrix
onesPayoffMat=[0 b
              -c (b-c)]; 

oPayDtDt=s1(2)*s2(2)*onesPayoffMat(2,2) + s1(2)*(1-s2(2))*onesPayoffMat(2,1) + ...
    (1-s1(2))*s2(2)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(2))*onesPayoffMat(1,1);
oPayDtDf=s1(2)*s2(1)*onesPayoffMat(2,2) + s1(2)*(1-s2(1))*onesPayoffMat(2,1) + ...
    (1-s1(2))*s2(1)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(1))*onesPayoffMat(1,1);
oPayDfDt=s1(1)*s2(2)*onesPayoffMat(2,2) + s1(1)*(1-s2(2))*onesPayoffMat(2,1) + ...
    (1-s1(1))*s2(2)*onesPayoffMat(1,2) + (1-s1(1))*(1-s2(2))*onesPayoffMat(1,1);
oPayDfDf=s1(1)*s2(1)*onesPayoffMat(2,2) + s1(1)*(1-s2(1))*onesPayoffMat(2,1) + ...
    (1-s1(1))*s2(1)*onesPayoffMat(1,2) + (1-s1(1))*(1-s2(1))*onesPayoffMat(1,1);
oPayDtIt=s1(2)*s2(4)*onesPayoffMat(2,2) + s1(2)*(1-s2(4))*onesPayoffMat(2,1) + ...
    (1-s1(2))*s2(4)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(4))*onesPayoffMat(1,1);
oPayDtIf=s1(2)*s2(3)*onesPayoffMat(2,2) + s1(2)*(1-s2(3))*onesPayoffMat(2,1) + ...
    (1-s1(2))*s2(3)*onesPayoffMat(1,2) + (1-s1(2))*(1-s2(3))*onesPayoffMat(1,1);
oPayDfIt=s1(1)*s2(4)*onesPayoffMat(2,2) + s1(1)*(1-s2(4))*onesPayoffMat(2,1) + ...
    (1-s1(1))*s2(4)*onesPayoffMat(1,2) + (1-s1(1))*(1-s2(4))*onesPayoffMat(1,1);
oPayDfIf=s1(1)*s2(3)*onesPayoffMat(2,2) + s1(1)*(1-s2(3))*onesPayoffMat(2,1) + ...
    (1-s1(1))*s2(3)*onesPayoffMat(1,2) + (1-s1(1))*(1-s2(3))*onesPayoffMat(1,1);
oPayItDt=s1(4)*s2(2)*onesPayoffMat(2,2) + s1(4)*(1-s2(2))*onesPayoffMat(2,1) + ...
    (1-s1(4))*s2(2)*onesPayoffMat(1,2) + (1-s1(4))*(1-s2(2))*onesPayoffMat(1,1);
oPayIfDt=s1(3)*s2(2)*onesPayoffMat(2,2) + s1(3)*(1-s2(2))*onesPayoffMat(2,1) + ...
    (1-s1(3))*s2(2)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(2))*onesPayoffMat(1,1);
oPayItDf=s1(4)*s2(1)*onesPayoffMat(2,2) + s1(4)*(1-s2(1))*onesPayoffMat(2,1) + ...
    (1-s1(4))*s2(1)*onesPayoffMat(1,2) + (1-s1(4))*(1-s2(1))*onesPayoffMat(1,1);
oPayIfDf=s1(3)*s2(1)*onesPayoffMat(2,2) + s1(3)*(1-s2(1))*onesPayoffMat(2,1) + ...
    (1-s1(3))*s2(1)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(1))*onesPayoffMat(1,1);

oPayItIt=s1(4)*s2(4)*onesPayoffMat(2,2) + s1(4)*(1-s2(4))*onesPayoffMat(2,1) + ...
    (1-s1(4))*s2(4)*onesPayoffMat(1,2) + (1-s1(4))*(1-s2(4))*onesPayoffMat(1,1);
oPayItIf=s1(4)*s2(3)*onesPayoffMat(2,2) + s1(4)*(1-s2(3))*onesPayoffMat(2,1) + ...
    (1-s1(4))*s2(3)*onesPayoffMat(1,2) + (1-s1(4))*(1-s2(3))*onesPayoffMat(1,1);
oPayIfIt=s1(3)*s2(4)*onesPayoffMat(2,2) + s1(3)*(1-s2(4))*onesPayoffMat(2,1) + ...
    (1-s1(3))*s2(4)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(4))*onesPayoffMat(1,1);
oPayIfIf=s1(3)*s2(3)*onesPayoffMat(2,2) + s1(3)*(1-s2(3))*onesPayoffMat(2,1) + ...
    (1-s1(3))*s2(3)*onesPayoffMat(1,2) + (1-s1(3))*(1-s2(3))*onesPayoffMat(1,1);

oPayDD=YD*YD*oPayDtDt + YD*(1-YD)*oPayDtDf + (1-YD)*YD*oPayDfDt + (1-YD)*(1-YD)*oPayDfDf;
oPayDI=YD*YI*oPayDtIt + (1-YD)*YI*oPayDfIt + YD*(1-YI)*oPayDtIf +(1-YD)*(1-YI)*oPayDfIf;
oPayID=YI*YD*oPayItDt + (1-YI)*(1-YD)*oPayIfDf + (1-YI)*YD*oPayIfDt +YI*(1-YD)*oPayItDf;
oPayII=YI*YI*oPayItIt + YI*(1-YI)*oPayItIf + (1-YI)*YI*oPayIfIt + (1-YI)*(1-YI)*oPayIfIf;

%%%
% Calculate agent 1's final payoff based on
%   Prob. of repeated game (p) vs. 1-shot game (1-p)
%   Prob. agents 1 and 2 deliberate (T1/d and T2/d) vs. use intuition (1-T1/d and 1-T2/d)
%   Average costs agent 1 and 2 pay when deliberating (T1/2 and T2/2)
%       (because the cost of deliberation varies uniformly, an agent who's  
%       willing to pay max cost of deliberation T pays an average cost of
%       T/2 to deliberate)

% We take a weighted avg. of agent 1's payoff in 4 situations:
%   1. payoff when both agents deliberate = (p*rPayDD + (1-p)*oPayDD - T1/2)
%       this happens with probability (T1/d)*(T2/d)
%   2. payoff when only agent 1 deliberates = (p*rPayDI + (1-p)*oPayDI - T1/2)
%       this happens with probability (T1/d)*(1 - T2/d)
%   ... and so on

payoff = (T1/d)*(T2/d)*(p*rPayDD + (1-p)*oPayDD - T1/2) + ... 
    (T1/d)*(1 - T2/d)*(p*rPayDI + (1-p)*oPayDI - T1/2) + ...  
    (1 - T1/d)*(T2/d)*(p*rPayID + (1-p)*oPayID) + ...         
    (1 - T1/d)*(1 - T2/d)*(p*rPayII + (1-p)*oPayII);          


end

steady_state_calc.m

function [ssDist,transitionPs]=steady_state_calc(payoffTable,N,w,a)
% Calculate the steady state distribution of the Moran process with 
% exponential fitness, assuming the limit of low mutation
%
% INPUTS
% payoffTable is an MxM matrix of payoffs, where M is the number of strategies
% N is the population size
% w is the selection strength (0=neutral drift, ->Inf=deterministic)
% a is assortment (probability of interacting with a player having the same
%                 strategy as you, rather than a randomly selected player)
%
% OUTPUTS
% ssDist is an Mx1 vector with the steady state frequency of each strategy
% transitionPs is an MxM vector with transition probilities between
%              homogeneous states of the population

numStrats=length(payoffTable);
transitionPs=zeros(numStrats);

% To do the steady state calc, we need to assemble a transition matrix
% betewen homogeneous states of the population
for s1=1:numStrats
    for s2=1:numStrats
        
        if s1~=s2
            
        % for each combination of strategies s1 and s2, the transition prob
        % going from s1 to s2 is the product of the probability a s1 mutant
        % arises, and the probability that mutant goes to fixation

                            %prob s1 mutant arises           %prob that s1 goes to fixation
                            %(assuming uniform mutation)
            transitionPs(s1,s2)=(1/(numStrats-1))            *calc_invasionprob(payoffTable,N,w,s1,s2,a);

        end
    end
end

% The diagonal entries of the transition matrix are the probability you
% stay in the given state - which is just 1 minus all the transition
% probabilities out of that state
for i=1:numStrats
    transitionPs(i,i)=1-sum(transitionPs(:,i));
end

% Now that you've got the transition matrix, find the eigenvector with
% eigenvalue 1, and that's your steady state distribution!
[vec,val]=eig(transitionPs);
[~, idxofEigWithVal1]=min(abs(diag(val)-1));  % For some reason, sometimes Matlab gets the eigenvalue wrong s/t its not literally 1, just very close

% normalize the eigenvector corresponding to eigenvalue so that it sums to 1, and you're done! 
ssDist=vec(:,idxofEigWithVal1)./sum(vec(:,idxofEigWithVal1));

end


function rhoBA=calc_invasionprob(payoffTable,N,w,sA,sB,a)
% here we calculate rhoAB, the probability that a sA mutant goes to fixation in
% a sB population

% this probability is the inverse of (1 + a sum of products) so I
% start by intitalizing rhoInv to 1, and then add on the sum of series of
% products, and then take the inverse
rhoInv=1;

% loop through each possible number of steps the sA mutant would need to go through in order to take over
% (this ranges from 1, if there are already N-1 sA players, 
% to N-1, if there is only the 1 initial sA mutant)
for k=1:N-1
    
    % same thing - initialize the eventual product of various terms as 1,
    % then multiply it by stuff
    thisProd=1;
    
    % loop through each possible make-up of the population (indexed by k=the
    % number of sA players) that would need to be traversed to get from 
    % there being only 1 sA player to there being k sA players
    for i=1:k        
        
        % for each population make-up, multiply by the ratio of the fitness 
        % of an sA player in the current pop to the fitness of an sB player
        % in the current pop (where there are i sA players)
        
        % fitness = exp( w * expectedpayoff)
        % expectedpayoff = a * payoff against your own strategy
        %                  + (1-a) *( (payoff against sB * fraction of other players who are sB) 
        %                            +(payoff against sA * fraction of other players who are sA)
        %                           )

        fi=exp(w* (a*payoffTable(sA,sA)+(1-a)*(payoffTable(sA,sB)*(N-i)/(N-1)   + payoffTable(sA,sA)*(i-1)/(N-1))));
        gi=exp(w* (a*payoffTable(sB,sB)+(1-a)*(payoffTable(sB,sB)*(N-i-1)/(N-1) + payoffTable(sB,sA)*(i)/(N-1)  )));

        % multiple the current state of the product by the ratio for this value of i
        thisProd=thisProd*(gi/fi);
    end
    % now that we've finished with the product for k, add it onto the sum
    rhoInv=rhoInv+thisProd;
end

% now that we've finished with the sum, take the inverse
rhoBA=1/rhoInv;
        
end 

vary_YD.m

function [ssDists,avgStrats] = vary_YD(N,b,c,d,w,a)
% Calculates and plots the average results from steady_state_calc 
% when varying YD and p (with YI fixed at .5)

% INPUTS
% N is the population size
% b is benefit of cooperation
% c is cost of cooperation
% d is max possible cost of deliberation (such that the cost of
%   deliberation varies uniformly in the interval [0,d])
% w is the selection strength (0=neutral drift, ->Inf=deterministic)
% a is assortment (probability of interacting with a player having the same
%                 strategy as you, rather than a randomly selected player)
%
% OUTPUTS
% ssDists stores the ssDist steady state frequencies of each strategy for
%   each p x YD combination.
% avgStrats stores the average strategy of of the steady state populations
%   for each p x YD combination.

ps_steadyState = 0:.01:1;
YD_steadyState= 0.5:0.01:1;

%% Create strat space
% generate all possible combinations of the first three strat parameters (sR,s1,sI)
% set to be 0 or 1
[x1,x2,x3] = ind2sub([2,2,2], find(zeros(2,2,2)==0));
strats = [x1 x2 x3] - 1;

% we now combine the above with the 11 possible deliberation thresholds,
% running from 0 to d in increments of d/10
Tvals = 0:d/10:d;
strats = repmat(strats,length(Tvals),1);
Tvals = sort(repmat(Tvals,1,length(strats)/length(Tvals)));
strats = [strats Tvals'];

% because we're using the YIYD calc payoff function, have to add in another
% column (4th column) for second intuitive response
strats = [strats(:,1:3) strats(:,3) strats(:,4)];

% ... to create a total of 2*2*2*11 (88) possible strategies
M = size(strats,1); % define M to be number of strats


%% Calculate steady state distribution

% initialize matrix to store steady state distributions for each p value
% and each YD value
ssDists = zeros(M,length(ps_steadyState),length(YD_steadyState));
%ssDists = zeros(M ,length(ps_steadyState));

% initialize matrix to store average strat values for each p value and each
% YD value
avgStrats = zeros(length(ps_steadyState),5,length(YD_steadyState));
%avgStrats = zeros(length(ps_steadyState),4); 

for YD=YD_steadyState
    YD % display YD value that's running
    
   for p=ps_steadyState     
        p % display p value that's running

        % initialize payoff table 
        payoffTable=zeros(M);

        % create payoff table
        for i=1:M
            for j=1:M
                payoffTable(i,j)=calculate_payoff_YDYI(strats(i,:),strats(j,:),b,c,d,p,YD,.5);
            end
        end

        % run steady state using this payoff table and store result in ssDists
        ssDists(:,p==ps_steadyState,YD==YD_steadyState)=steady_state_calc(payoffTable,N,w,a);

        % use ssDists to get average strategy values by taking a weighted average of
        % the strat values based on the distribution in ssDists
        avgStrats(p==ps_steadyState,:,YD==YD_steadyState) = sum(repmat(ssDists(:,p==ps_steadyState,YD==YD_steadyState),1,5)...
            .*strats);      
     end
end


%% Create Surface Plots

% create a surface showing avg. amount of deliberation in p x YD space
T = squeeze(avgStrats(:,5,:)); % stores this average value
figure
surf(YD_steadyState,ps_steadyState,T)
xlabel('Sensitivity of Deliberation y_d')
ylabel('Probability of Coordination p')
colormap(parula)
view(2); view(-90,90) % change the orientation of the plot...
set(gca,'ydir','reverse') % ... and reverse axis direction for clearer presentation

% create a surface showing closest intuitive profile in p x YD space
intuitiveProfile = round(squeeze(avgStrats(:,3,:))).*round(squeeze(avgStrats(:,4,:))); % stores the intuitive responses as 1,2
figure
surf(YD_steadyState,ps_steadyState,intuitiveProfile)
xlabel('Sensitivity of Deliberation y_d')
ylabel('Probability of Coordination p')
colormap([.8 0 0; 1 1 0]) % we create a custom colormap where...
% red always intuitively defects and yellow always intuitively cooperates
view(2); view(-90,90) % change the orientation of the plot...
set(gca,'ydir','reverse') % ... and reverse axis direction for clearer presentation

vary_YI.m

function [ssDists,avgStrats] = vary_YI(N,b,c,d,w,a)
% Calculates and plots the average results from steady_state_calc 
% when varying YI and p (with YD fixed at 1)

% INPUTS
% N is the population size
% b is benefit of cooperation
% c is cost of cooperation
% d is max possible cost of deliberation (such that the cost of
%   deliberation varies uniformly in the interval [0,d])
% w is the selection strength (0=neutral drift, ->Inf=deterministic)
% a is assortment (probability of interacting with a player having the same
%                 strategy as you, rather than a randomly selected player)
%
% OUTPUTS
% ssDists stores the ssDist steady state frequencies of each strategy for
%   each p x YI combination.
% avgStrats stores the average strategy of of the steady state populations
%   for each p x YI combination.


ps_steadyState = 0:.01:1;
YI_steadyState= 0.5:0.01:1;

%% Generate Strategy Space 

% For the steady state calculation, we must generate a discrete strategy
% set. To do so:
%   pR, p1, pIR and pI1 are all restricted to being either 0 or 1;  
%   we discretize the deliberation threshold (T) in intervals of d/10

% generate all possible combinations of the first four strat parameters (pR,p1,pIR SI1)
% set to be 0 or 1
[x1,x2,x3,x4] = ind2sub([2,2,2,2], find(zeros(2,2,2,2)==0));
strats = [x1 x2 x3 x4] - 1;

% we now combine the above with the 11 possible deliberation thresholds,
% running from 0 to d in increments of d/10
Tvals = 0:d/10:d;
strats = repmat(strats,length(Tvals),1);
Tvals = sort(repmat(Tvals,1,length(strats)/length(Tvals)));
strats = [strats Tvals'];

% ... to create a total of 2*2*2*2*11 (176) possible strategies
M = size(strats,1); % define M to be number of strats


%% Main Steady State Loop 

% initialize matrix to store steady state distributions for each p value
% each 1 value and each YI value
ssDists = zeros(M,length(ps_steadyState),length(YI_steadyState));

% initialize matrix to store average strat values for each p value, YI val and 1 val
avgStrats = zeros(length(ps_steadyState),5,length(YI_steadyState));

for YI=YI_steadyState
    YI % display YI value that's running

   for p=ps_steadyState     
        p % display p value that's running

        % initialize payoff table 
        payoffTable=zeros(M);

        % create payoff table
        for i=1:M
            for j=1:M
                payoffTable(i,j)=calculate_payoff_YDYI(strats(i,:),strats(j,:),b,c,d,p,1,YI);
            end
        end

        % run steady state using this payoff table and store result in ssDists
        ssDists(:,p==ps_steadyState,YI==YI_steadyState)=steady_state_calc(payoffTable,N,w,a);

        % use ssDists to get average strategy values by taking a weighted average of
        % the strat values based on the distribution in ssDists
        avgStrats(p==ps_steadyState,:,YI==YI_steadyState) = sum(repmat(ssDists(:,p==ps_steadyState,YI==YI_steadyState),1,5).*strats);

   end
end

%% Create Surface Plots

% create a surface showing avg. amount of deliberation in p x YI space
T = squeeze(avgStrats(:,5,:)); % stores this average value
figure
surf(YI_steadyState,ps_steadyState,T)
xlabel('Sensitivity of Intuition y_i')
ylabel('Probability of Coordination p')
colormap(parula)
view(2); view(-90,90) % change the orientation of the plot...
set(gca,'ydir','reverse') % ... and reverse axis direction for clearer presentation

% create a surface showing closest intuitive profile in p x YI space
intuitiveProfile = round(squeeze(avgStrats(:,3,:))).*round(squeeze(avgStrats(:,4,:))) ...
    + 2*(round(squeeze(avgStrats(:,3,:))) - round(squeeze(avgStrats(:,4,:)))); % stores the intuitive responses as 1,2,3
figure
surf(YI_steadyState,ps_steadyState,intuitiveProfile)
xlabel('Sensitivity of Intuition y_i')
ylabel('Probability of Coordination p')
colormap([.8 0 0; 1 1 0; 0 .8 0]) % we create a custom colormap where...
% red always intuitively defects, yellow always intuitively cooperates, and
% green discriminates its response (defect in SD; cooperate in CG)
view(2); view(-90,90) % change the orientation of the plot...
set(gca,'ydir','reverse') % ... and reverse axis direction for clearer presentation