(* This notebook calculated the Tmap for the Fisherian Model example explored in Section 3.3. It gives numerical values of the equilibria and calculates E-stability. It also demonstrates how to calculate a version of the bifurcation diagrams shown in the paper.*)
(* Optimal Weights with two misspecified models in the Fisherian Model *)
T[M1_List,M2_List,w_List]:=Module[{A,B,CC,Ma,Mb,TM1,TM2,TW1,TW2,theta, thetaj,Ez,u1,u2,rho},
A=μ+α (w[[1]]M1[[1]]+w[[2]]M2[[1]]);
B=α w[[1]]M1[[2]]ρ1+ζ1;
CC=α w[[2]] M2[[2]]ρ2+ζ2;
Ma={{1,0},{0,z1}}//Inverse;
Mb={{1,0},{0,z2}}//Inverse;

TM1=Ma.{{A},{B z1+CC s}};
TM2=Mb.{{A},{B s+CC z2}};
theta ={{M1[[1]]},{M1[[2]]}};
thetaj ={{M2[[1]]},{M2[[2]]}};
Ez ={{1,0,0},{0,z1,s},{0,s,z2}};
u1 ={{1,0,0},{0,1,0}};
u2 ={{1,0,0},{0,0,1}};
rho ={{1,0,0},{0,ρ1,0},{0,0,ρ2}};
TW1 = Inverse[-Transpose[theta].u1.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u1].theta]*(Transpose[theta].u1.rho.Ez.rho.{{μ},{ζ1},{ζ2}} + w[[2]]*Transpose[theta].u1.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u2].thetaj);
TW2 = Inverse[-Transpose[thetaj].u2.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u2].thetaj]*(Transpose[thetaj].u2.rho.Ez.rho.{{μ},{ζ1},{ζ2}} + w[[1]]*Transpose[thetaj].u2.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u1].theta);
{TM1[[1]],TM1[[2]],TM2[[1]],TM2[[2]],Flatten@TW1, Flatten@TW2}



];
T[{a1,b1},{a2,b2},{γ1,γ2}]/.{μ->0,s->0}//Simplify

Stability[αh_,ζ1h_,ζ2h_,sh_,z1h_,z2h_,r1h_,r2h_,μh_,eq_]:=Module[{rule,temp,forecast,rule2},
(*rule=NSolve[{T[{a1,b1},{a2,b2},{γ1,γ2}]=={a1,b1,a2,b2,γ1,γ2}}/.{α->αh,ζ1->ζ1h,ζ2->ζ2h,s->sh,z1->z1h,z2->z2h,μ->μh},{a1,b1,a2,b2,γ1,γ2}];*)
rule=NSolve[{T[{a1,b1},{a2,b2},{γ1,γ2}]=={a1,b1,a2,b2,γ1,γ2}}/.{α->αh,ζ1->ζ1h,ζ2->ζ2h,s->sh,z1->z1h,z2->z2h,μ->μh,ρ1-> r1h,ρ2->r2h},{a1,b1,a2,b2,γ1,γ2}];
Print["Equilibria - Including non-real valued"];
Print[rule //MatrixForm];
temp=Table[γ1(a1+r1h b1)+γ2(a2+r2h b2)/.rule[[k]],{k,1,7}];
forecast=DeleteCases[temp,_Complex];
Print["Eigenvalues governing E-stability for all solutions"];
rule2=Table[D[Flatten@T[{a1,b1},{a2,b2},{γ1,γ2}]/. {α->αh, ζ1 -> ζ1h, ζ2 -> ζ2h, s->sh, z1 ->z1h, z2->z2h, μ->μh,ρ1-> r1h,ρ2->r2h},{{a1,b1,a2,b2,γ1,γ2}}]/.{rule[[k]][[1]],  rule[[k]][[2]], rule[[k]][[3]], rule[[k]][[4]], rule[[k]][[5]],rule[[k]][[6]]} //Eigenvalues,{k,1,7}];
Print[rule2//MatrixForm];
Print["Equilibrium Forecasts for real solutions if X_{t-1}=(1 1)'"];
forecast 
];
Stability[.95,-.1,.9,0,1,1,.7,.9,0,6]//MatrixForm

(* Pseudo-Bifurcation Diagram *)

EqForecast2[αh_,ζ1h_,ζ2h_,sh_,z1h_,z2h_,r1h_,r2h_,μh_,x1_,x2_]:=Module[{rule,rule1,rule2,ll,y1,y2,y3,y4,y5,y6,y7,ya1,ya2,ya3,ya4,ya5,ya6,ya7,REE},
rule=NSolve[{T[{a1,b1},{a2,b2},{γ1,γ2}]=={a1,b1,a2,b2,γ1,γ2}}/.{α->αh,ζ1->ζ1h,ζ2->ζ2h,s->sh,z1->z1h,z2->z2h,ρ1->r1h,ρ2->r2h,μ->μh},{a1,b1,a2,b2,γ1,γ2}];
y1=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[1,1]],  rule[[1,2]], rule[[1,3]], rule[[1,4]], rule[[1,5]],rule[[1,6]]};
y2=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[2,1]],  rule[[2,2]], rule[[2,3]], rule[[2,4]], rule[[2,5]],rule[[2,6]]};
y3=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[3,1]],  rule[[3,2]], rule[[3,3]], rule[[3,4]], rule[[3,5]],rule[[3,6]]};
y4=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[4,1]],  rule[[4,2]], rule[[4,3]], rule[[4,4]], rule[[4,5]],rule[[4,6]]};
y5=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[5,1]],  rule[[5,2]], rule[[5,3]], rule[[5,4]], rule[[5,5]],rule[[5,6]]};
y6=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[6,1]],  rule[[6,2]], rule[[6,3]], rule[[6,4]], rule[[6,5]],rule[[6,6]]};
y7=γ1(a1+r1h*b1)*x1+γ2(a2+r2h*b2)*x2/.{rule[[7,1]],  rule[[7,2]], rule[[7,3]], rule[[7,4]], rule[[7,5]],rule[[7,6]]};

rule1={y1,y2,y3,y4,y5,y6,y7};
rule2=DeleteCases[rule1,_Complex];
ll=Length[rule2];
Table[{αh,rule2[[k]]},{k,1,ll}]
(*If[size==1,{{αh,y1},{αh,ya1}},If[size==3,{{αh,y1},{αh,y2},{αh,y3},{αh,ya1},{αh,ya2},{αh,ya3}},{{αh,y1},{αh,y2},{αh,y3},{αh,y4},{αh,y5},{αh,y6},{αh,y7},{αh,ya1},{αh,ya2},{αh,ya3},{αh,ya4},{αh,ya5},{αh,ya6},{αh,ya7}}]]*)
];
EqForecast2[.9,.9,-.9,0,1,1,.9,.9,0,1,1]

output=Flatten[Table[EqForecast2[αh,.9,-.9,0,1,1,.9,.9,0,1,1],{αh,0.0001,1.5,0.0002}],1];

ListPlot[output]