Stability of fixed point of label dynamics

1. Cooperation Indicator for Strategies i=1,...,16 for given label XY,’1’ stands for cooperation

"1-CCCC";
(Subscript[ind, CC])[1]=1;
(Subscript[ind, CD])[1]=1;
(Subscript[ind, DC])[1]=1;
(Subscript[ind, DD])[1]=1;
"2-CCCD";
(Subscript[ind, CC])[2]=1;
(Subscript[ind, CD])[2]=1;
(Subscript[ind, DC])[2]=1;
(Subscript[ind, DD])[2]=0;
"3-CCDC";
(Subscript[ind, CC])[3]=1;
(Subscript[ind, CD])[3]=1;
(Subscript[ind, DC])[3]=0;
(Subscript[ind, DD])[3]=1;
"4-CDCC";
(Subscript[ind, CC])[4]=1;
(Subscript[ind, CD])[4]=0;
(Subscript[ind, DC])[4]=1;
(Subscript[ind, DD])[4]=1;
"5-DCCC";
(Subscript[ind, CC])[5]=0;
(Subscript[ind, CD])[5]=1;
(Subscript[ind, DC])[5]=1;
(Subscript[ind, DD])[5]=1;
"6-CCDD";
(Subscript[ind, CC])[6]=1;
(Subscript[ind, CD])[6]=1;
(Subscript[ind, DC])[6]=0;
(Subscript[ind, DD])[6]=0;
"7-CDCD";
(Subscript[ind, CC])[7]=1;
(Subscript[ind, CD])[7]=0;
(Subscript[ind, DC])[7]=1;
(Subscript[ind, DD])[7]=0;
"8-DCCD";
(Subscript[ind, CC])[8]=0;
(Subscript[ind, CD])[8]=1;
(Subscript[ind, DC])[8]=1;
(Subscript[ind, DD])[8]=0;
"9-DCDC";
(Subscript[ind, CC])[9]=0;
(Subscript[ind, CD])[9]=1;
(Subscript[ind, DC])[9]=0;
(Subscript[ind, DD])[9]=1;
"10-DDCC";
(Subscript[ind, CC])[10]=0;
(Subscript[ind, CD])[10]=0;
(Subscript[ind, DC])[10]=1;
(Subscript[ind, DD])[10]=1;
"11-CDDC";
(Subscript[ind, CC])[11]=1;
(Subscript[ind, CD])[11]=0;
(Subscript[ind, DC])[11]=0;
(Subscript[ind, DD])[11]=1;
"12-CDDD";
(Subscript[ind, CC])[12]=1;
(Subscript[ind, CD])[12]=0;
(Subscript[ind, DC])[12]=0;
(Subscript[ind, DD])[12]=0;
"13-DCDD";
(Subscript[ind, CC])[13]=0;
(Subscript[ind, CD])[13]=1;
(Subscript[ind, DC])[13]=0;
(Subscript[ind, DD])[13]=0;
"14-DDCD";
(Subscript[ind, CC])[14]=0;
(Subscript[ind, CD])[14]=0;
(Subscript[ind, DC])[14]=1;
(Subscript[ind, DD])[14]=0;
"15-DDDC";
(Subscript[ind, CC])[15]=0;
(Subscript[ind, CD])[15]=0;
(Subscript[ind, DC])[15]=0;
(Subscript[ind, DD])[15]=1;
"16-DDDD";
(Subscript[ind, CC])[16]=0;
(Subscript[ind, CD])[16]=0;
(Subscript[ind, DC])[16]=0;
(Subscript[ind, DD])[16]=0;

2. Define (Subscript[p, sXY]) - Strategies s’s Probability to Cooperate Conditional on the Receiver’s Label

(Subscript[p, sXY])[s_,XY_]=((Subscript[ind, XY])[s](1-η)+Sum[(Subscript[ind, i])[s]η/3,{i,{CC,CD,DC,DD}}]-(Subscript[ind, XY])[s]η/3)(1-η);

3. Define (Subscript[p, XYs]) - Probability to carry label XY for someone who plays strategy s

(Subscript[p, CCs])[s_,t_]=(Subscript[p, sXY])[s,CC](Subscript[p, CC])[t-1]+(Subscript[p, sXY])[s,CD](Subscript[p, CD])[t-1]//Simplify;
(Subscript[p, CDs])[s_,t_]=(Subscript[p, sXY])[s,DC](Subscript[p, DC])[t-1]+(Subscript[p, sXY])[s,DD](Subscript[p, DD])[t-1]//Simplify;
(Subscript[p, DCs])[s_,t_]=(1-(Subscript[p, sXY])[s,CC])(Subscript[p, CC])[t-1]+(1-(Subscript[p, sXY])[s,CD])(Subscript[p, CD])[t-1]//Simplify;
(Subscript[p, DDs])[s_,t_]=(1-(Subscript[p, sXY])[s,DC])(Subscript[p, DC])[t-1]+(1-(Subscript[p, sXY])[s,DD])(Subscript[p, DD])[t-1]//Simplify;

4. Equations of Recurrence (degrees of freedom=3)

eq1=(Subscript[p, CC])[t]==  Collect[(Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}];
eq2=(Subscript[p, CD])[t]==Collect[ (Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}];
eq3=(Subscript[p, DC])[t]== Collect[((Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}];

5. Recurrence in Matrix Notation ((Subscript[p, CC])[t],(Subscript[p, CD])[t],(Subscript[p, DC])[t])^T=W*((Subscript[p, CC])[t-1],(Subscript[p, CD])[t-1],(Subscript[p, DC])[t-1])^T

(Subscript[w, 11])=Coefficient[Collect[(Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, CC])[-1+t]]//FullSimplify;
(Subscript[w, 12])=Coefficient[Collect[(Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, CD])[-1+t]]//FullSimplify;
(Subscript[w, 13])=Coefficient[Collect[(Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, DC])[-1+t]]//FullSimplify;
(Subscript[w, 21])=Coefficient[Collect[ (Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, CC])[-1+t]]//FullSimplify;
(Subscript[w, 22])=Coefficient[Collect[ (Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, CD])[-1+t]]//FullSimplify;
(Subscript[w, 23])=Coefficient[Collect[ (Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, DC])[-1+t]]//FullSimplify;
(Subscript[w, 31])=Coefficient[Collect[((Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, CC])[-1+t]]//FullSimplify;
(Subscript[w, 32])=Coefficient[Collect[((Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, CD])[-1+t]]//FullSimplify;
(Subscript[w, 33])=Coefficient[Collect[((Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}]/.{(Subscript[p, DD])[-1+t]->1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]})/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t]}],(Subscript[p, DC])[-1+t]]//FullSimplify;
W3={{(Subscript[w, 11]),(Subscript[w, 12]),(Subscript[w, 13])} ,{(Subscript[w, 21]),(Subscript[w, 22]),(Subscript[w, 23])},{(Subscript[w, 31]),(Subscript[w, 32]),(Subscript[w, 33])}};
MatrixForm[W3];
W3s=W3//FullSimplify
"rows"
(Subscript[w, 11])-(Subscript[w, 12])//Simplify
(Subscript[w, 21])-(Subscript[w, 22])//Simplify
(Subscript[w, 21])+(Subscript[w, 23])//Simplify
(Subscript[w, 31])-(Subscript[w, 32])//Simplify
"columns"
(Subscript[w, 21])+(Subscript[w, 11])//Simplify
(Subscript[w, 31])+(Subscript[w, 11])//Simplify
(Subscript[w, 21])-(Subscript[w, 31])//Simplify
(Subscript[w, 22])+(Subscript[w, 12])//Simplify
(Subscript[w, 32])+(Subscript[w, 12])//Simplify

6. Eigenvalues: the general case, RESULT: Trace(W)=(Subscript[λ, 1])+(Subscript[λ, 2])+(Subscript[λ, 3])=3(Subscript[λ, 1])

Tr[W3s]//Simplify
Det[W3s]//FullSimplify
e=Eigenvalues[W3s,Cubics->True];
"Eigenvalues"
e[[1,1]]//Simplify
e[[1,2]]
e[[1,3]]
"Numerator,Denominator"
Numerator[e[[1,2]]]//Simplify
Denominator[e[[1,2]]];
Numerator[e[[1,3]]];
Denominator[e[[1,3]]];
e[[1,2]]*e[[1,3]]//Simplify

7. Eigenvalues: for the case (Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 6])=1, RESULT: η>0 implies (Subscript[λ, max])<1

W3s1236=W3s/.{(Subscript[x, 4])->0,(Subscript[x, 5])->0,(Subscript[x, 6])->1-(Subscript[x, 1])-(Subscript[x, 2])-(Subscript[x, 3]),(Subscript[x, 7])->0,(Subscript[x, 8])->0,(Subscript[x, 9])->0,(Subscript[x, 10])->0,(Subscript[x, 11])->0,(Subscript[x, 12])->0,(Subscript[x, 13])->0,(Subscript[x, 14])->0,(Subscript[x, 15])->0}//Simplify;
MatrixForm[%]//FullSimplify
Tr[W3s1236]//Simplify
Det[W3s1236]
CP=CharacteristicPolynomial[W3s1236,λ]//FullSimplify
CoefficientList[CP,λ]//FullSimplify
e1236=Eigenvalues[W3s1236]//FullSimplify
e1236[[3]]/.{ (Subscript[x, 1])-> 0,  (Subscript[x, 2])->  1, (Subscript[x, 3])-> 0}//FullSimplify

8. Simulation: Example

eq1sim=(Subscript[p, CC])[t]==  Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}]//Simplify;
eq2sim=(Subscript[p, CD])[t]== Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}]//Simplify;
eq3sim=(Subscript[p, DC])[t]== Collect[(Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}]/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t]}];
eq4sim=(Subscript[p, DD])[t]== Collect[(Sum[(Subscript[x, s])*(Subscript[p, DDs])[s,t],{s,1,16}]/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
param={η->1/10,(Subscript[x, 1])->1/16,(Subscript[x, 2])->1/16,(Subscript[x, 3])->1/16,(Subscript[x, 4])->1/16,(Subscript[x, 5])->1/16,(Subscript[x, 6])->1/16,(Subscript[x, 7])->1/16,(Subscript[x, 8])->1/16,(Subscript[x, 9])->1/16,(Subscript[x, 10])->1/16,(Subscript[x, 11])->1/16,(Subscript[x, 12])->1/16,(Subscript[x, 13])->1/16,(Subscript[x, 14])->1/16,(Subscript[x, 15])->1/16};
param2={η->1/10,(Subscript[x, 1])->11/16,(Subscript[x, 2])->0,(Subscript[x, 3])->5/16,(Subscript[x, 4])->0,(Subscript[x, 5])->0,(Subscript[x, 6])->0,(Subscript[x, 7])->0,(Subscript[x, 8])->0,(Subscript[x, 9])->0,(Subscript[x, 10])->0,(Subscript[x, 11])->0,(Subscript[x, 12])->0,(Subscript[x, 13])->0,(Subscript[x, 14])->0,(Subscript[x, 15])->0};
T=10;
sol=RecurrenceTable[{eq1sim/.param,eq2sim/.param,eq3sim/.param,eq4sim/.param,(Subscript[p, CC])[0]==1/10,(Subscript[p, CD])[0]==2/10,(Subscript[p, DC])[0]==3/10,(Subscript[p, DD])[0]==4/10},{(Subscript[p, CC]),(Subscript[p, CD]),(Subscript[p, DC]),(Subscript[p, DD])},{t,0,T}];
ListPlot[{Table[sol[[n,1]],{n,1,T}],Table[sol[[n,2]],{n,1,T}],Table[sol[[n,3]],{n,1,T}],Table[sol[[n,4]],{n,1,T}]},PlotMarkers->Automatic,PlotLegends->{"pCC","pCD","pDC","pDD"}]
sol2=RecurrenceTable[{eq1sim/.param2,eq2sim/.param2,eq3sim/.param2,eq4sim/.param2,(Subscript[p, CC])[0]==1/10,(Subscript[p, CD])[0]==2/10,(Subscript[p, DC])[0]==3/10,(Subscript[p, DD])[0]==4/10},{(Subscript[p, CC]),(Subscript[p, CD]),(Subscript[p, DC]),(Subscript[p, DD])},{t,0,T}];
ListPlot[{Table[sol2[[n,1]],{n,1,T}],Table[sol2[[n,2]],{n,1,T}],Table[sol2[[n,3]],{n,1,T}],Table[sol2[[n,4]],{n,1,T}]},PlotMarkers->Automatic,PlotLegends->{"pCC","pCD","pDC","pDD"}]

9. Fixed Point

pcch=Collect[(Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
pcdh=Collect[ (Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
pdch=Collect[((Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
pddh=Collect[((Sum[(Subscript[x, s])*(Subscript[p, DDs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
w11test=Coefficient[pcch,(Subscript[p, CC])[-1+t]];
w12test=Coefficient[pcch,(Subscript[p, CD])[-1+t]];
w13test=Coefficient[pcch,(Subscript[p, DC])[-1+t]];
w14test=Coefficient[pcch,(Subscript[p, DD])[-1+t]];
w21test=Coefficient[pcdh,(Subscript[p, CC])[-1+t]];
w22test=Coefficient[pcdh,(Subscript[p, CD])[-1+t]];
w23test=Coefficient[pcdh,(Subscript[p, DC])[-1+t]];
w24test=Coefficient[pcdh,(Subscript[p, DD])[-1+t]];
w31test=Coefficient[pdch,(Subscript[p, CC])[-1+t]];
w32test=Coefficient[pdch,(Subscript[p, CD])[-1+t]];
w33test=Coefficient[pdch,(Subscript[p, DC])[-1+t]];
w34test=Coefficient[pdch,(Subscript[p, DD])[-1+t]];
w41test=Coefficient[pddh,(Subscript[p, CC])[-1+t]];
w42test=Coefficient[pddh,(Subscript[p, CD])[-1+t]];
w43test=Coefficient[pddh,(Subscript[p, DC])[-1+t]];
w44test=Coefficient[pddh,(Subscript[p, DD])[-1+t]];
W4test={{w11test, w12test, w13test, w14test},{w21test, w22test, w23test, w24test},{w31test, w32test, w33test, w34test},{w41test, w42test, w43test, w44test}};
"ZEROS";
zeros=w13test+w14test+w21test+w22test+w33test+w34test+w41test+w42test
"RELATIONS"
w11test-(1-w31test)//Simplify
w12test-(1-w32test)//Simplify
w23test-(1-w43test)//Simplify
w24test-(1-w44test)//Simplify
"FP";
(Subscript[FPp, CC])=(w12test/(1-w11test))(1/(2+w12test/(1-w11test)+(1-w23test)/w24test));
(Subscript[FPp, CD])=(1/(2+w12test/(1-w11test)+(1-w23test)/w24test));
(Subscript[FPp, DC])=(1/(2+w12test/(1-w11test)+(1-w23test)/w24test));
(Subscript[FPp, DD])=((1-w23test)/w24test)(1/(2+w12test/(1-w11test)+(1-w23test)/w24test));

10. Simulation: Loop (note: need to execute 1.-5. and 8.-9. first)

"Set T and number of simulations";
T=100;
M=10000;
"Element in the 15-Simplex";
Inix=FindInstance[0<(Subscript[x, 1])<1&&0<(Subscript[x, 2])<1&&0<(Subscript[x, 3])<1&&0<(Subscript[x, 4])<1&&0<(Subscript[x, 5])<1&&0<(Subscript[x, 6])<1&&0<(Subscript[x, 7])<1&&0<(Subscript[x, 8])<1&&0<(Subscript[x, 9])<1&&0<(Subscript[x, 10])<1&&0<(Subscript[x, 11])<1&&0<(Subscript[x, 12])<1&&0<(Subscript[x, 13])<1&&0<(Subscript[x, 14])<1&&0<(Subscript[x, 15])<1&&0<(Subscript[x, 16])<1&&(Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15])+(Subscript[x, 16])==1,{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]),(Subscript[x, 5]),(Subscript[x, 6]),(Subscript[x, 7]),(Subscript[x, 8]),(Subscript[x, 9]),(Subscript[x, 10]),(Subscript[x, 11]),(Subscript[x, 12]),(Subscript[x, 13]),(Subscript[x, 14]),(Subscript[x, 15]),(Subscript[x, 16])},100];
"Intitialisation for labels";
InLab=FindInstance[0<(Subscript[y, 1])<1&&0<(Subscript[y, 2])<1&&0<(Subscript[y, 3])<1&&0<(Subscript[y, 4])<1&&(Subscript[y, 1])+(Subscript[y, 2])+(Subscript[y, 3])+(Subscript[y, 4])==1,{(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]),(Subscript[y, 4])},100];
For[n=1,n≤ M,n++,
"POPULATION STATE: choose one of the 100 initials";
k[n]=RandomInteger[{1,100}];
"Initialisation=Reordering";
order[n]=RandomSample[Range[16],16];
"Parameter population state";
par[n]=N[{η-> RandomReal[{0,1/10}],(Subscript[x, order[n][[1]]])->((Subscript[x, 1])/.Inix[[k[n],1]]),(Subscript[x, order[n][[2]]])->((Subscript[x, 2])/.Inix[[k[n],2]]),(Subscript[x, order[n][[3]]])->((Subscript[x, 3])/.Inix[[k[n],3]]),(Subscript[x, order[n][[4]]])->((Subscript[x, 4])/.Inix[[k[n],4]]),(Subscript[x, order[n][[5]]])->((Subscript[x, 5])/.Inix[[k[n],5]]),(Subscript[x, order[n][[6]]])->((Subscript[x, 6])/.Inix[[k[n],6]]),(Subscript[x, order[n][[7]]])->((Subscript[x, 7])/.Inix[[k[n],7]]),(Subscript[x, order[n][[8]]])->((Subscript[x, 8])/.Inix[[k[n],8]]),(Subscript[x, order[n][[9]]])->((Subscript[x, 9])/.Inix[[k[n],9]]),(Subscript[x, order[n][[10]]])->((Subscript[x, 10])/.Inix[[k[n],10]]),(Subscript[x, order[n][[11]]])->((Subscript[x, 11])/.Inix[[k[n],11]]),(Subscript[x, order[n][[12]]])->((Subscript[x, 12])/.Inix[[k[n],12]]),(Subscript[x, order[n][[13]]])->((Subscript[x, 13])/.Inix[[k[n],13]]),(Subscript[x, order[n][[14]]])->((Subscript[x, 14])/.Inix[[k[n],14]]),(Subscript[x, order[n][[15]]])->((Subscript[x, 15])/.Inix[[k[n],15]]),(Subscript[x, order[n][[16]]])->((Subscript[x, 16])/.Inix[[k[n],16]])}];

"LABEL: choose one of the 100 initials";
l[n]=RandomInteger[{1,100}];
"Initialisation=Reordering";
order2[n]=RandomSample[Range[4],4];
"solve recurrence";
solu[n]=RecurrenceTable[{eq1sim/.par[n],eq2sim/.par[n],eq3sim/.par[n],eq4sim/.par[n],(Subscript[p, CC])[0]==((Subscript[y, order2[n][[1]]])/.InLab[[l[n],order2[n][[1]]]]),(Subscript[p, CD])[0]==((Subscript[y, order2[n][[2]]])/.InLab[[l[n],order2[n][[2]]]]),(Subscript[p, DC])[0]==((Subscript[y, order2[n][[3]]])/.InLab[[l[n],order2[n][[3]]]]),(Subscript[p, DD])[0]==((Subscript[y, order2[n][[4]]])/.InLab[[l[n],order2[n][[4]]]])},{(Subscript[p, CC]),(Subscript[p, CD]),(Subscript[p, DC]),(Subscript[p, DD])},{t,0,T}];
For[i=1,i≤ T,i++,dev[n,i]=Abs[((Subscript[FPp, CC])/.par[n])-solu[n][[i]][[1]]]+Abs[((Subscript[FPp, CD])/.par[n])-solu[n][[i]][[2]]]+Abs[((Subscript[FPp, DC])/.par[n])-solu[n][[i]][[3]]]+Abs[((Subscript[FPp, DD])/.par[n])-solu[n][[i]][[4]]]];
s[n]=Select[Table[dev[n,j],{j,1,T}],#< 0.10&,1];
p[n]=Position[Table[dev[n,j],{j,1,T}],s[n][[1]]]
]

"Substract 1 for Initialisation";
H=Histogram[Table[Flatten[p[k]][[1]]-1,{k,1,M-1}], {5},"Probability",AxesLabel->{"Number of Iterations","Frequency / Distribution"},ImageSize->Large,Ticks->{{0,10,20,30,40,50,60,70,80,90,100},{0,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90,1}},ChartStyle->Gray,AxesStyle->Directive[Black,14]];
For[k=1,k≤ T,k++,
S[0]=0;
"Number of trails which needed weakly less than k iterations for convergence";
S[k]=S[k-1]+Count[Table[Flatten[p[j]][[1]]-1,{j,1,T}],k];
]
L=ListPlot[ N[Table[S[j]/S[T],{j,1,T}]],AxesLabel->{"Number of Iterations","Frequency / Distribution"},ImageSize->Large, Ticks->{{0,10,20,30,40,50,60,70,80,90,100},{0,0.10,0.20,0.30,0.40,0.50,0.60,0.70,0.80,0.90,1}},Joined->{False, True},PlotMarkers->{"●",Small},PlotStyle->Black,AxesStyle->Directive[Black,14]];
Show[H,L]

Check Solution Presented in the Paper

11. Fixed Point of Recurrence (note: need to execute 1.-5. first)

pcch=Collect[(Sum[(Subscript[x, s])*(Subscript[p, CCs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
pcdh=Collect[ (Sum[(Subscript[x, s])*(Subscript[p, CDs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))}//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
pdch=Collect[((Sum[(Subscript[x, s])*(Subscript[p, DCs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
pddh=Collect[((Sum[(Subscript[x, s])*(Subscript[p, DDs])[s,t],{s,1,16}])/.{(Subscript[x, 16])->1-((Subscript[x, 1])+(Subscript[x, 2])+(Subscript[x, 3])+(Subscript[x, 4])+(Subscript[x, 5])+(Subscript[x, 6])+(Subscript[x, 7])+(Subscript[x, 8])+(Subscript[x, 9])+(Subscript[x, 10])+(Subscript[x, 11])+(Subscript[x, 12])+(Subscript[x, 13])+(Subscript[x, 14])+(Subscript[x, 15]))})//Simplify,{(Subscript[p, CC])[-1+t],(Subscript[p, CD])[-1+t],(Subscript[p, DC])[-1+t],(Subscript[p, DD])[-1+t]}];
eq1fph=(Subscript[p, CC])[t]==pcch;
eq2fph=(Subscript[p, CD])[t]==pcdh;
eq3fph=(Subscript[p, DC])[t]==pdch;
eq4fph=(Subscript[p, DC])[t]==pddh;
eq5fp=(Subscript[p, DD])[-1+t]==1-(Subscript[p, CC])[-1+t]-(Subscript[p, CD])[-1+t]-(Subscript[p, DC])[-1+t]/.{(Subscript[p, CC])[t]->(Subscript[p, CC]),(Subscript[p, CD])[t]->(Subscript[p, CD]),(Subscript[p, DC])[t]->(Subscript[p, DC]),(Subscript[p, DD])[t]->(Subscript[p, DD]),(Subscript[p, CC])[t-1]->(Subscript[p, CC]),(Subscript[p, CD])[t-1]->(Subscript[p, CD]),(Subscript[p, DC])[t-1]->(Subscript[p, DC]),(Subscript[p, DD])[t-1]->(Subscript[p, DD])};
eq1fp=eq1fph/.{(Subscript[p, CC])[t]->(Subscript[p, CC]),(Subscript[p, CD])[t]->(Subscript[p, CD]),(Subscript[p, DC])[t]->(Subscript[p, DC]),(Subscript[p, DD])[t]->(Subscript[p, DD]),(Subscript[p, CC])[t-1]->(Subscript[p, CC]),(Subscript[p, CD])[t-1]->(Subscript[p, CD]),(Subscript[p, DC])[t-1]->(Subscript[p, DC]),(Subscript[p, DD])[t-1]->(Subscript[p, DD])};
eq2fp=eq2fph/.{(Subscript[p, CC])[t]->(Subscript[p, CC]),(Subscript[p, CD])[t]->(Subscript[p, CD]),(Subscript[p, DC])[t]->(Subscript[p, DC]),(Subscript[p, DD])[t]->(Subscript[p, DD]),(Subscript[p, CC])[t-1]->(Subscript[p, CC]),(Subscript[p, CD])[t-1]->(Subscript[p, CD]),(Subscript[p, DC])[t-1]->(Subscript[p, DC]),(Subscript[p, DD])[t-1]->(Subscript[p, DD])};
eq3fp=eq3fph/.{(Subscript[p, CC])[t]->(Subscript[p, CC]),(Subscript[p, CD])[t]->(Subscript[p, CD]),(Subscript[p, DC])[t]->(Subscript[p, DC]),(Subscript[p, DD])[t]->(Subscript[p, DD]),(Subscript[p, CC])[t-1]->(Subscript[p, CC]),(Subscript[p, CD])[t-1]->(Subscript[p, CD]),(Subscript[p, DC])[t-1]->(Subscript[p, DC]),(Subscript[p, DD])[t-1]->(Subscript[p, DD])};
eq4fp=eq4fph/.{(Subscript[p, CC])[t]->(Subscript[p, CC]),(Subscript[p, CD])[t]->(Subscript[p, CD]),(Subscript[p, DC])[t]->(Subscript[p, DC]),(Subscript[p, DD])[t]->(Subscript[p, DD]),(Subscript[p, CC])[t-1]->(Subscript[p, CC]),(Subscript[p, CD])[t-1]->(Subscript[p, CD]),(Subscript[p, DC])[t-1]->(Subscript[p, DC]),(Subscript[p, DD])[t-1]->(Subscript[p, DD])};
fp=Solve[{eq1fp,eq2fp,eq3fp,eq4fp,eq5fp},{(Subscript[p, CC]),(Subscript[p, CD]),(Subscript[p, DC]),(Subscript[p, DD])}];

w11test=Coefficient[pcch,(Subscript[p, CC])[-1+t]];
w12test=Coefficient[pcch,(Subscript[p, CD])[-1+t]];
w13test=Coefficient[pcch,(Subscript[p, DC])[-1+t]];
w14test=Coefficient[pcch,(Subscript[p, DD])[-1+t]];
w21test=Coefficient[pcdh,(Subscript[p, CC])[-1+t]];
w22test=Coefficient[pcdh,(Subscript[p, CD])[-1+t]];
w23test=Coefficient[pcdh,(Subscript[p, DC])[-1+t]];
w24test=Coefficient[pcdh,(Subscript[p, DD])[-1+t]];
w31test=Coefficient[pdch,(Subscript[p, CC])[-1+t]];
w32test=Coefficient[pdch,(Subscript[p, CD])[-1+t]];
w33test=Coefficient[pdch,(Subscript[p, DC])[-1+t]];
w34test=Coefficient[pdch,(Subscript[p, DD])[-1+t]];
w41test=Coefficient[pddh,(Subscript[p, CC])[-1+t]];
w42test=Coefficient[pddh,(Subscript[p, CD])[-1+t]];
w43test=Coefficient[pddh,(Subscript[p, DC])[-1+t]];
w44test=Coefficient[pddh,(Subscript[p, DD])[-1+t]];
W4test={{w11test, w12test, w13test, w14test},{w21test, w22test, w23test, w24test},{w31test, w32test, w33test, w34test},{w41test, w42test, w43test, w44test}};
"ZEROS";
zeros=w13test+w14test+w21test+w22test+w33test+w34test+w41test+w42test
"RELATIONS"
w11test-(1-w31test)//Simplify
w12test-(1-w32test)//Simplify
w23test-(1-w43test)//Simplify
w24test-(1-w44test)//Simplify

12. Test of Symmetry and Solution presented in the paper (note: need to execute 1.-5. first)

eq1s=(Subscript[p, CC])==w11 (Subscript[p, CC])+w12 (Subscript[p, CD]);
eq2s=(Subscript[p, CD])==w23 (Subscript[p, DC])+w24 (Subscript[p, DD]) ;
eq3s=(Subscript[p, DC])==(1-w11) (Subscript[p, CC])+(1-w12) (Subscript[p, CD]) ;
eq4s=(Subscript[p, DD])==(1-w23) (Subscript[p, DC])+(1-w24) (Subscript[p, DD]) ;
eq5s=(Subscript[p, DD])==1- (Subscript[p, CC])- (Subscript[p, CD])- (Subscript[p, DC]) ;
fps1=Solve[eq1s&&eq2s &&eq3s &&eq4s,{(Subscript[p, CC]),(Subscript[p, CD]),(Subscript[p, DC]),(Subscript[p, DD])}]
fps2=Solve[eq1s&&eq2s &&eq3s &&eq4s&&eq5s,{(Subscript[p, CC]),(Subscript[p, CD]),(Subscript[p, DC]),(Subscript[p, DD])}]//FullSimplify
((Subscript[p, CD])/.fps2[[1,2]])-(1/(2+(w12/(1-w11))+((1-w23)/w24)))//Simplify
((Subscript[p, CD])/.fps2[[1,2]])-((Subscript[p, DC])/.fps2[[1,3]])//Simplify
((Subscript[p, CC])/.fps2[[1,1]])-w12/(1-w11)((Subscript[p, CD])/.fps2[[1,2]])//Simplify