(* Replication files for Gibbs (2016) - This notebook calculated the Tmaps (unrestricted and restricted) for Lucas-type monetary model explored in Section 3 of the paper. It provides numerical equilibrium solutions and E-stability results. The final portion of the notebook provides code to simulate the models under real-time learning. Figure 6 of the paper is demonstrated.*)

(*Unrestriced Regression Weights with two misspecified models*)
T[M1_List,M2_List,w_List]:=Module[{A,B,CC,R,G,Mhold,Ma,Mb,TM1,TM2,TW},
A=μ+α (w[[1]]M1[[1]]+w[[2]]M2[[1]]);
B=α w[[1]]M1[[2]]+ζ1;
CC=α w[[2]] M2[[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}};
R={{A},{M1[[1]]},{M2[[1]]}};
G={{B,CC},{M1[[2]],0},{0,M2[[2]]}};
Mhold=R.R\[Transpose]+G.{{z1,s},{s,z2}}.G\[Transpose];

TW={{(Mhold[[1]][[2]]-w[[2]] Mhold[[2]][[3]])/Mhold[[2]][[2]]},{(Mhold[[1]][[3]]-w[[1]]Mhold[[2]][[3]])/Mhold[[3]][[3]]}};
{TM1[[1]],TM1[[2]],TM2[[1]],TM2[[2]],TW[[1]],TW[[2]]}
];
T[{a1,b1},{a2,b2},{γ1,γ2}]//Simplify

(* Equilibrium and Stability *)

Stability[αh_,ζ1h_,ζ2h_,sh_,z1h_,z2h_,μh_,eq_]:=Module[{rule,temp,forecast,realforecast,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}];
Print["Equilibria - Including non-real valued"];
Print[rule //MatrixForm];
temp=Table[γ1(a1+b1)+γ2(a2+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},{{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[.7,-1,-.9,0,1,1,0,6]//MatrixForm

(* Restricted Optimal Weights*)

TR[M1_List,M2_List,w_List]:=Module[{A,B,CC,R,G,Mhold,Ma,Mb,TM1,TM2,TW},
A=μ+α (w[[1]]M1[[1]]+(1-w[[1]])M2[[1]]);
B=α w[[1]]M1[[2]]+ζ1;
CC=α (1-w[[1]]) M2[[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}};
R={{A-M2[[1]]},{M1[[1]]-M2[[1]]}};
G={{B,CC-M2[[2]]},{M1[[2]],-M2[[2]]}};
Mhold= R.R\[Transpose]+G.{{z1,s},{s,z2}}.G\[Transpose];
TW=Mhold[[1]][[2]]/Mhold[[2]][[2]];
Flatten@{TM1[[1]],TM1[[2]],TM2[[1]],TM2[[2]],TW}
];

(* Constant gain learning simulations *)
s=0;
μ=0;
(* This simulates optimal weights forecast combination. It takes in parameter values α, ζ1, ζ2; a value for the gain g;  standard deviations for shocks; a seed for the random number generator; a length of the simulation; and the equilibrium from which to start beliefs. *)
Simul1[αh_,ζ1h_,ζ2h_,g_,σ1_,σ2_,σ3_,seed_,tt_,eq_]:=Module[{rule,rule2,rule3,rule4,y,theta,theta1,y1,y2,zz1,zz2,zz3,R,xx,yy,jj,ll,MA},
SeedRandom[seed];
rule=NSolve[{T[{a1,b1},{a2,b2},{γ1,γ2}]=={a1,b1,a2,b2,γ1,γ2}}/.{α->αh,ζ1->ζ1h,ζ2->ζ2h,z1->σ1,z2->σ2},{a1,b1,a2,b2,γ1,γ2}];
theta = Table[{a1,b1,a2,b2,γ1,γ2}/.rule[[k]],{k,1,7}];
(* Check if real*)
rule2=DeleteCases[theta,_Complex,Infinity];
rule3 =Table[If[Length[rule2[[k]]]==6,rule2[[k]],"Null"],{k,1,7}];
rule4=DeleteCases[rule3,"Null"];
y={rule4[[1]][[5]]*rule4[[1]][[1]]+rule4[[1]][[6]]*rule4[[1]][[3]]};
(* Simulate Economy *)
ll =Length[rule4];
Print[ll];
theta1={rule4[[eq]]};
R=IdentityMatrix[6];
For[i=1,i<tt,i++,
zz1=RandomVariate[NormalDistribution[0,σ1]];
zz2=RandomVariate[NormalDistribution[0,σ2]] ;
zz3=RandomVariate[NormalDistribution[0,σ3]];
y1 = theta1[[i,1]]+theta1[[i,2]]*zz1;
y2 = theta1[[i,3]] + theta1[[i,4]]*zz2;
y=Insert[y,αh*(theta1[[i,5]]*y1 + theta1[[i,6]]*y2)+zz1*ζ1h +zz2*ζ2h+zz3,i+1];
yy = {{y[[i+1]]},{y[[i+1]]},{y[[i+1]]}};
xx = {{1, zz1,0,0,0,0},{0,0,1,zz2,0,0},{0,0,0,0,y1,y2}};
jj= theta1[[i]];
jj = jj + g*Inverse[R].Transpose[xx].(yy-xx.jj);
R = R + g *(Transpose[xx].xx-R);
theta1=Insert[theta1,Flatten@Transpose[jj],i+1];
];
Print[rule4];
MA = ExponentialMovingAverage[y,1/10];

{Table[y[[k]],{k,1,tt,1}],Table[MA[[k]],{k,1,tt,1}],Table[theta1[[k]][[1]],{k,1,tt,1}],Table[theta1[[k]][[4]],{k,1,tt,1}],Table[theta1[[k]][[5]],{k,1,tt,1}],Table[theta1[[k]][[6]],{k,1,tt,1}]}

];
(*Simulation for Restricted Optimal Weights Case*)
Simul2[αh_,ζ1h_,ζ2h_,g_,σ1_,σ2_,σ3_,seed_,tt_,eq_]:=Module[{rule,rule2,rule3,rule4,y,theta,theta1,y1,y2,zz1,zz2,zz3,R,xx,yy,jj,ll,MA,MV},
SeedRandom[seed];
rule=NSolve[{TR[{a1,b1},{a2,b2},{γ1}]=={a1,b1,a2,b2,γ1}}/.{α->αh,ζ1->ζ1h,ζ2->ζ2h,z1->σ1,z2->σ2},{a1,b1,a2,b2,γ1}];
theta = Table[{a1,b1,a2,b2,γ1}/.rule[[k]],{k,1,3}];
(* Check if real*)
rule2=DeleteCases[theta,_Complex,Infinity];
rule3 =Table[If[Length[rule2[[k]]]==5,rule2[[k]],"Null"],{k,1,3}];
rule4=DeleteCases[rule3,"Null"];
Print[rule4];
y={rule4[[1]][[5]]*rule4[[1]][[1]]+(1-rule4[[1]][[5]])*rule4[[1]][[3]]};
(* Simulate Economy *)
ll =Length[rule4];
Print[ll];
theta1={rule4[[eq]]};
R=IdentityMatrix[5];
For[i=1,i<tt,i++,
zz1=RandomVariate[NormalDistribution[0,σ1]];
zz2=RandomVariate[NormalDistribution[0,σ2]] ;
zz3=RandomVariate[NormalDistribution[0,σ3]];
y1 = theta1[[i,1]]+theta1[[i,2]]*zz1;
y2 = theta1[[i,3]] + theta1[[i,4]]*zz2;
y=Insert[y,αh*(theta1[[i,5]]*y1 + (1-theta1[[i,5]])*y2)+zz1*ζ1h +zz2*ζ2h+zz3,i+1];
yy = {{y[[i+1]]},{y[[i+1]]},{y[[i+1]]-y2}};
xx = {{1, zz1,0,0,0},{0,0,1,zz2,0},{0,0,0,0,y1-y2}};
jj= theta1[[i]];
jj = jj + g*Inverse[R].Transpose[xx].(yy-xx.jj);
R = R + g *(Transpose[xx].xx-R);
theta1=Insert[theta1,Flatten@Transpose[jj],i+1];
];
MA = ExponentialMovingAverage[y,1/10];
{Table[y[[k]],{k,1,tt,1}],Table[MA[[k]],{k,1,tt,1}],Table[theta1[[k]][[1]],{k,1,tt,1}],Table[theta1[[k]][[4]],{k,1,tt,1}],Table[theta1[[k]][[5]],{k,1,tt,1}]}
];

(*Dynamic Predictor Selection Case*)
Simul3[αh_,ζ1h_,ζ2h_,β_,g_,λ_,σ1_,σ2_,σ3_,seed_,tt_]:=Module[{rule,rule2,rule3,rule4,y,theta,theta1,y1,y2,zz1,zz2,zz3,R,xx,yy,jj,ll,MA,MV,w1,MSE1,MSE2,mm1,mm2},
SeedRandom[seed];
theta1={{0,1,0,1}};
w1={.5};
y={0};
MSE1={0};
MSE2 ={0};
R=IdentityMatrix[4];
For[i=1,i<tt,i++,
zz1=RandomVariate[NormalDistribution[0,σ1]];
zz2=RandomVariate[NormalDistribution[0,σ2]] ;
zz3=RandomVariate[NormalDistribution[0,σ3]];
y1 = theta1[[i,1]]+theta1[[i,2]]*zz1;
y2 = theta1[[i,3]] + theta1[[i,4]]*zz2;
w1=Insert[w1,Exp[-β*MSE1[[i]]]/(Exp[-β*MSE1[[i]]]+Exp[-β*MSE2[[i]]]),i+1];
(*Print[Exp[-β*MSE1[[i]]]/(Exp[-β*MSE1[[i]]]+Exp[-β*MSE2[[i]]])];*)
y=Insert[y,αh*(w1[[i]]*y1 + (1-w1[[i]])*y2)+zz1*ζ1h +zz2*ζ2h+zz3,i+1];
mm1=MSE1[[i]];
mm2=MSE2[[i]];
mm1=mm1 +λ ((y[[i+1]]-y1)^2-mm1);
mm2=mm2 +λ ((y[[i+1]]-y2)^2-mm2);
MSE1=Insert[MSE1,mm1,i+1];
MSE2=Insert[MSE2,mm2,i+1];
yy = {{y[[i+1]]},{y[[i+1]]}};
xx = {{1, zz1,0,0},{0,0,1,zz2}};
jj= theta1[[i]];
jj = jj + g*Inverse[R].Transpose[xx].(yy-xx.jj);
R = R + g *(Transpose[xx].xx-R);
theta1=Insert[theta1,Flatten@Transpose[jj],i+1];
];
MA = ExponentialMovingAverage[y,1/10];
(*ListLinePlot[{Table[y[[k]],{k,500,tt,1}],Table[MA[[k]],{k,500,tt,1}]},PlotStyle->{{Opacity[0.4]},{Opacity[1]}},PlotRange->{{0,tt-500},{-3,5}}]
ListLinePlot[{Table[w1[[k]],{k,1,tt}]},PlotRange-> {{0,tt-500},{-1,2}}]*)
{Table[y[[k]],{k,1,tt,1}],Table[MA[[k]],{k,1,tt,1}],Table[theta1[[k]][[1]],{k,1,tt,1}],Table[theta1[[k]][[4]],{k,1,tt,1}],w1}
];

(* Generates Figure 6 from Gibbs (2016) *)
AllG[seed_,g_,start_,range_]:=Module[{L1,L2,L3},
L1=Line[{{850,-4},{950,-4}}];
L2 =Line[{{850,-4.75},{950,-4.75}}];
L3=Line[{{850,-5.5},{950,-5.5}}];
Show[ListLinePlot[{Simul1[.9,.5 ,.75,g,.1,.1,.25,seed,range,1][[3]],Simul2[.9,.5,.75,g,.1,.1,.25,seed,range,1][[3]],Simul3[.9,.5,.75,50,g,.35,.1,.1,.25,seed,range][[3]]},PlotStyle->{{Opacity[1],RGBColor[0,.5,.1],Dashed},{Opacity[1],Blue},{Opacity[1],Purple,Thick,Dotted}},PlotRange->{{start,range},{-7.5,4}},Frame->True,FrameStyle->Directive[{15,15}]],Graphics[{RGBColor[0,.5,.1],Dashed,L1}],Graphics[{Blue,L2}],Graphics[{Purple,Thick,Dotted,L3}]]];
AllG[31,.08,800,1800]