(* This Mathematica file carries out some symbolic computations to check or provide some results in the Online Appendix of my paper "Stabilization Policy and Lags". The notations are essentially the same as in the paper and the Online Appendix. The run time is about 39 minutes on my computer. O. Loisel (CREST), November 2022. *)

Delete All Output;
Clear["Global`*"]
Model={"Basic New Keynesian model:","Model of Svensson (1997) and Ball (1999):","Model of Smets and Wouters (2007):","Model of Schmitt-Grohé and Uribe (1997):"};
For[m=1,m≤Length[Model],m++,
Clear@@DeleteCases[Names@"`*","Model"|"m"];
Print[" "];
Print[" "];
Print[Style[Model[[m]],Bold]];
If[m==1,
A={{σ(1-X),1},{ X,(β-X)/κ}};
B={{-X},{0}}];
If[m==2,
A={{1-(Subscript[λ, 2])X,-(Subscript[μ, 2])X},{-(Subscript[χ, 2])X,1-X}};
B={{(Subscript[μ, 2])X},{0}}];
If[m==3,
η=α (Subscript[ϕ, p]) + (Subscript[z, 1]) (Subscript[z, y]) ;
(Subscript[i, 1])=(Subscript[w, 1]) ;
(Subscript[w, 2])=(Subscript[w, 1]) + (1-(Subscript[w, 1]))(Subscript[w, 3])/(Subscript[w, 1]) ;
λ=γ (Subscript[c, 1])/(1-(Subscript[c, 1])) ;
(Subscript[π, 2])=((1-(Subscript[w, 1]))/(Subscript[w, 1]))(1-((1-(Subscript[w, 1]))/(Subscript[w, 1]))(Subscript[π, 1])) ;
(Subscript[q, 1])=((1-(Subscript[w, 1]))/(Subscript[w, 1]))(Subscript[k, 1]) ;
(Subscript[z, y])=((Subscript[i, y])/(1-(Subscript[k, 1])))(((Subscript[w, 1])/(1-(Subscript[w, 1])))- (Subscript[k, 1])) ;
(Subscript[A, 1])={{α (Subscript[ϕ, p])(Subscript[c, y]),α (Subscript[ϕ, p])(Subscript[i, y]),α (Subscript[ϕ, p])(Subscript[z, 1])(Subscript[z, y]),(1-(1+(Subscript[ϕ, p])-η)α)(Subscript[ϕ, p])},{(1+(Subscript[z, 1]))(Subscript[k, 1])(Subscript[c, y]),(1+(Subscript[z, 1]))(Subscript[k, 1])(Subscript[i, y])+(1-(Subscript[k, 1]))η,((Subscript[z, y])-α (Subscript[ϕ, p]))(Subscript[k, 1])(Subscript[z, 1]),(Subscript[k, 1])η-(1+(Subscript[z, 1]))(Subscript[k, 1])(Subscript[ϕ, p])},
{(Subscript[c, y]),(Subscript[i, y]),(Subscript[z, 1])(Subscript[z, y]),η-(Subscript[ϕ, p])},{-(Subscript[z, 1])(Subscript[c, y]),-(Subscript[z, 1])(Subscript[i, y]),α (Subscript[z, 1])(Subscript[ϕ, p]),(Subscript[z, 1])(Subscript[ϕ, p])},{-(Subscript[c, y]),-(Subscript[i, y]),α (Subscript[ϕ, p]),(Subscript[ϕ, p])},{α (Subscript[c, y]),α (Subscript[i, y]),-(1-α (Subscript[z, 1])(Subscript[z, y])),-(1+(Subscript[ϕ, p])-η)α},{-(1-λX/γ)/(1-λ/γ),0,1,-(Subscript[σ, l])}}/η;
(Subscript[A, 2])={{(1-(Subscript[c, 1]))-X+ (Subscript[c, 1]) (X^2),0},{0,(1-(Subscript[i, 1]))-X+ (Subscript[i, 1]) (X^2)},{-(1-(Subscript[q, 1]))(Subscript[c, y])/η,-(1-(Subscript[q, 1]))(Subscript[i, y])/η},{(1+(Subscript[z, 1]))(Subscript[k, 1]) (Subscript[c, y]) (1- (Subscript[k, 1]) X),(1+(Subscript[z, 1]))(Subscript[k, 1]) (Subscript[i, y]) -((1+(Subscript[z, 1]))(Subscript[k, 1]) (Subscript[i, y]) + (1-(Subscript[k, 1]))η)(Subscript[k, 1]) X},{-α (Subscript[π, 3]) (Subscript[c, y]) X,-α (Subscript[π, 3]) (Subscript[i, y]) X},{((Subscript[w, 4])/(1-λ/γ))(X-(λ/γ)(X^2)),0}};
(Subscript[A, 3])={{0,-(Subscript[c, 2]) + (Subscript[c, 2]) X},{0,0},{(1-(Subscript[q, 1]))α (Subscript[ϕ, p])/η,(1-(Subscript[q, 1]))(Subscript[ϕ, p])/η},{((Subscript[z, y])-α (Subscript[ϕ, p]))(Subscript[k, 1]) (Subscript[z, 1]) (1 - (Subscript[k, 1]) X),-((1+(Subscript[z, 1]))(Subscript[ϕ, p]) - η) (Subscript[k, 1]) (1 - (Subscript[k, 1]) X)},{(Subscript[π, 3]) (η - α (Subscript[z, 1]) (Subscript[z, y]))X,α (Subscript[π, 3]) (Subscript[ϕ, p]) X},{(1-(Subscript[w, 1]))-(1+(Subscript[w, 4]))X  + (Subscript[w, 1]) (X^2),(Subscript[w, 4]) (Subscript[σ, l]) X}};
(Subscript[A, 4])={{0,(Subscript[c, 3])},{(Subscript[i, 2]) X,0},{(Subscript[q, 1])-X,1},{0,0},{0,η((Subscript[π, 2]) - X + (Subscript[π, 1]) (X^2))},{0,(1+(Subscript[w, 1]))- (Subscript[w, 2]) X + (Subscript[w, 3]) (X^2)}};
A=ArrayFlatten[{{ArrayFlatten[{{IdentityMatrix[7],-(Subscript[A, 1]),ConstantArray[0,{7,2}]}}]},{ArrayFlatten[{{ConstantArray[0,{6,7}],(Subscript[A, 2]),(Subscript[A, 3]),(Subscript[A, 4])}}]}}];
B={{0},{0},{0},{0},{0},{0},{0},{-(Subscript[c, 3])},{0},{-1},{0},{0},{0}}];
If[m==4,
η=(1+γ)(Subscript[s, c]) +(1-(Subscript[s, k]))(σ-(Subscript[s, c]));
φ=(1-(Subscript[s, k]))σ - (Subscript[s, c]) ;
χ= 1 - β(1-δ) ;
(Subscript[A, 1])={{(1-(Subscript[s, k]))(Subscript[σs, i])+((1+γ)(Subscript[δs, c])(Subscript[s, k])-(1-δ)(1-(Subscript[s, k]))(Subscript[σs, i]))X},{(Subscript[σs, i])-((σ-(Subscript[s, c]))(Subscript[δs, k])+(1-δ)(Subscript[σs, i]))X},{((Subscript[φσs, i])/(Subscript[s, c]))+((1+γ)δσ (Subscript[s, k])-(1-δ)((Subscript[φσs, i])/(Subscript[s, c])))X},
{(1-(Subscript[s, k]))(Subscript[σs, i])-(((Subscript[γs, c])+σ)δ+(1-δ)(Subscript[σs, i]))(1-(Subscript[s, k]))X},{((Subscript[γs, c])+φ)((Subscript[σs, i])/(Subscript[s, c]))+(((Subscript[γs, c])+σ)(Subscript[δs, k])-((Subscript[γs, c])+φ)(1-δ)((Subscript[σs, i])/(Subscript[s, c])))X}}/(δ η);
(Subscript[A, 2])=(φ σ-(1-(Subscript[s, k]))χ (Subscript[s, c]))(Subscript[s, i])/(δ η (Subscript[s, c])) + ((1+γ)δ (σ^2) (Subscript[s, c]) (Subscript[s, k]) - (2-δ)φ (σ^2) (Subscript[s, i]) + ((γ (Subscript[s, c]) +σ)δ + (1-δ)σ (Subscript[s, i]))(1- (Subscript[s, k]))χ (Subscript[s, c]))X/(δ η σ (Subscript[s, c])) + ((1-δ)φ (Subscript[s, i]) -(1+γ)δ (Subscript[s, c]) (Subscript[s, k]))σ (X^2) /(δ η (Subscript[s, c])) ;
(Subscript[A, 3])= ((1-(Subscript[s, k]))σ (Subscript[s, i]) + ((1+γ)δ (Subscript[s, c]) (Subscript[s, k])  - (1-δ)(1-(Subscript[s, k]))σ (Subscript[s, i]) ) X)(1 - (Subscript[s, k]) + ω (Subscript[s, k]))τ/(δ η) ;
(Subscript[A, 4])= 1 - (1/β)X ;
(Subscript[B, 1])={{(1-(Subscript[s, k]))(Subscript[s, c])},{(Subscript[s, c])},{(1-(Subscript[s, k]))σ},{(1-(Subscript[s, k]))(Subscript[s, c])+ηω},{(Subscript[γs, c])+(1-(Subscript[s, k]))σ}}(- τ )/ ((1-τ)η );
(Subscript[B, 2])= (((1-(Subscript[s, k]))(Subscript[s, c]) + η ω)χ - (1-(Subscript[s, k]))(σ^2)  + (1-(Subscript[s, k]))(σ^2) X)τ / ((1-τ)η σ) ;
(Subscript[B, 3])= ((1-τ)η - (1-(Subscript[s, k]))τ (Subscript[s, c]))(1 - (Subscript[s, k]) + ω (Subscript[s, k]))τ / ((1-τ)η) ;
A=ArrayFlatten[{{ArrayFlatten[{{IdentityMatrix[5],-(Subscript[A, 1]),ConstantArray[0,{5,1}]}}]},{ArrayFlatten[{{0,0,0,0,0,(Subscript[A, 2]),0}}]},{ArrayFlatten[{{0,0,0,0,0,(Subscript[A, 3]),(Subscript[A, 4])}}]}}];
B=ArrayFlatten[{{-(Subscript[B, 1])},{(Subscript[B, 2])},{(Subscript[B, 3])}}]];
AX[X_]:=Evaluate[A];
If[m==4,
DetA0=Simplify[Det[AX[0]],(ω==0||ω==1)],
DetA0=Simplify[Det[AX[0]]]];
If[PossibleZeroQ[DetA0],
Print[•,"  ",Style["The assumption of non-zero det[A(0)]",Bold]," is not satisfied:"],
Print[•,"  ",Style["The assumption of non-zero det[A(0)]",Bold]," is satisfied:"]];
Print["   ","det[A(0)] = ",DetA0,"."];
n=Length[B];
If[m==4,
ψ=Table[Simplify[Det[Drop[ArrayFlatten[{{A,B}}],0,{index}]],(ω==0||ω==1)],{index,n+1}],
ψ=Table[Det[Drop[ArrayFlatten[{{A,B}}],0,{index}]],{index,n+1}]];
(* ψ denotes the (n+1)-element vector composed of ±(Subscript[W, 1])(X), ..., ±(Subscript[W, n])(X), and Q(X). *)
W=Drop[ψ,{n+1}];
(* W denotes the n-element vector composed of ±(Subscript[W, 1])(X), ..., and ±(Subscript[W, n])(X). *)
If[PossibleZeroQ[Times@@W],
Print[•,"  ",Style["The assumption of non-zero (Subscript[W, 1])(X), ..., (Subscript[W, n])(X)",Bold]," is not satisfied: "],
Print[•,"  ",Style["The assumption of non-zero (Subscript[W, 1])(X), ..., (Subscript[W, n])(X)",Bold]," is satisfied: "]];
Print["   ","W(X) = ",MatrixForm[W],"."];
If[m==4,
DJ=Table[Simplify[PolynomialGCD[ψ[[index]],ψ[[n+1]]],(ω==0||ω==1)],{index,1,n}],
DJ=Table[Simplify[PolynomialGCD[ψ[[index]],ψ[[n+1]]]],{index,1,n}]];
If[AllTrue[PossibleZeroQ[Exponent[DJ,X]],TrueQ],
Print[•,"  ",Style["Condition (22)",Bold]," is satisfied for any non-empty set J:"];
Print["   ","(Subscript[W, j])(X) and Q(X) have no common roots for any j in {1,...,",n,"}."],Print[•,"  ",Style["Condition (22)",Bold]," might not be satisfied for some non-empty set J:"];
Print["   ","(Subscript[W, j])(X) and Q(X) have at least one common root for some j in {1,...,",n,"}, further analysis is required."]];
If[Length[Position[Map[#==1&,Map[Length,Exponent[W,X,List],{1}],{-1}],True]]==0,
Print[•,"  ",Style["Condition (24)",Bold]," is violated for any singleton set J:"];
Print["   ","(Subscript[W, j])(X) has at least one non-zero root for any j in {1,...,",n,"}."],
Print[•,"  ",Style["Condition (24)",Bold]," is violated for any singleton set J except the set(s) in ",Position[Map[#==1&,Map[Length,Exponent[W,X,List],{1}],{-1}],True],"."];
];
GJ=Table[0,{index1,1,n},{index2,1,n}];
If[m==4,
Do[GJ[[index1,n]]=Simplify[PolynomialGCD[ψ[[index1]],ψ[[n]]],(ω==0||ω==1)],{index1,1,n-1}],
Do[GJ[[index1,index2]]=Simplify[PolynomialGCD[ψ[[index1]],ψ[[index2]]]],{index1,1,n-1},{index2,index1+1,n}]];
If[Length[Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True]]==0,
If[m==4,
Print[•,"  ",Style["Condition (24)",Bold]," is satisfied for any two-element set J including 7:"];
Print["   ","(Subscript[W, j])(X) and (Subscript[W, 7])(X) have no common non-zero roots for any j in {1,...,",n-1,"}."],
Print[•,"  ",Style["Condition (24)",Bold]," is satisfied for any two-element set J:"];
Print["   ","(Subscript[W, i])(X) and (Subscript[W, j])(X) have no common non-zero roots for any subset {i,j} of {1,...,",n,"}."]],
If[m==4,
Print[•,"  ",Style["Condition (24)",Bold]," is satisfied for any two-element set J including 7 except the set(s) in ",Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True],":"];
Print["   ","(Subscript[W, j])(X) and (Subscript[W, 7])(X) have no common non-zero roots for any set {j,7} with j in {1,...,",n-1,"} except the set(s) in ",Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True],", "];
Print["   ","and the greatest common divisor(s) of (Subscript[W, j])(X) and (Subscript[W, 7])(X) for the set(s) {j,7} in ",Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True]," is (are) "];
Print[Flatten[ReplacePart[GJ,Position[Map[#≤1&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True]->Nothing]],"."],
Print[•,"  ",Style["Condition (24)",Bold]," is satisfied for any two-element set J except the set(s) in ",Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True],":"];
Print["   ","(Subscript[W, i])(X) and (Subscript[W, j])(X) have no common non-zero roots for any subset {i,j} of {1,...,",n,"} except the set(s) in ",Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True],", "];
Print["   ","and the greatest common divisor(s) of (Subscript[W, i])(X) and (Subscript[W, j])(X) for the set(s) {i,j} in ",Position[Map[#≥2&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True]," is (are) "];
Print[Flatten[ReplacePart[GJ,Position[Map[#≤1&,Map[Length,Exponent[GJ,X,List],{2}],{-1}],True]->Nothing]],"."]]];
];