(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)

BeginPackage["ratexpmodsolver`",{"SingularValueDecompositionUnstandard`","LinearAlgebra`MatrixManipulation`",
    "NumericalMath`","ControlSystems`","NumericalControlMethods`"}]

(*Copyright \[Copyright] Frank Hespeler 2008*)



ratexpmodsolver::usage=
  "Calculates the coefficients for the simulation formula from C. A. Sims, Solving linear rational expectations models, Computational Economics, vol. 20, no. 1-2, pp. 1\[Dash]20, October 2002, p.12 using the results from the generalized schur decomposition and the original system matrices. The algorithm offer also the alterntive solution according to the method presented in F. Hespeler, Solution algprithm to a class of monetary rational equilibrium macromodels with optimal monetary policy, Computational Economics, vol. 31, no. 3, pp. 207\[Dash]223, April 2008. Inputs: The matrices obtained fom the system's matrix pencil \!\(\*
StyleBox[\"ss\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"tt\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"qq\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"zz\",\nFontSlant->\"Italic\"]\), the system's coefficienct matrix for the expectational error \!\(\*
StyleBox[\"gga\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)the system's coefficienct matrix for the shock term \!\(\*
StyleBox[\"schockcoff\",\nFontSlant->\"Italic\"]\), the number of predetermined variables \!\(\*
StyleBox[\"predet\",\nFontSlant->\"Italic\"]\), the numbers of state and nonstate variables \!\(\*
StyleBox[\"n1\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)and\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"n2\",\nFontSlant->\"Italic\"]\), the value for a switch from syable to unstable eigenvalues \!\(\*
StyleBox[\"cutoff\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\), the value for the rate of persistence in the shock \!\(\*
StyleBox[\"inert\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)and two ptional arguments:\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)the measure \!\(\*
StyleBox[\"precise\",\nFontSlant->\"Italic\"]\) for the internally minimally required precision  and \!\(\*
StyleBox[\"tol\",\nFontSlant->\"Italic\"]\) for the tolerance level in the computation of pseudoinverses. Output: the coefficient linking the stable variables of t+1 to the stable variables in t \!\(\*
StyleBox[\"coffx1x1\",\nFontSlant->\"Italic\"]\), the coefficient linking the stable variables of t+1 to the unstable variables in t \!\(\*
StyleBox[\"coffx1x2\",\nFontSlant->\"Italic\"]\), the coefficient linking the unstable variables of t+1 to the stable variables in t \!\(\*
StyleBox[\"coffx2x1\",\nFontSlant->\"Italic\"]\), the coefficient linking the unstable variables of t+1 to the un stable variables in t \!\(\*
StyleBox[\"coffx2x2\",\nFontSlant->\"Italic\"]\), the coefficient linking the stable variables of t+1 to the system's shock in t \!\(\*
StyleBox[\"coffx1z\",\nFontSlant->\"Italic\"]\), the coefficient linking the stable variables of t+1 to the system's correction of expectations in t \!\(\*
StyleBox[\"coffx1zfut\",\nFontSlant->\"Italic\"]\), the coefficient linking the unstable variables of t+1 to the system's shock in t \!\(\*
StyleBox[\"coffx2z\",\nFontSlant->\"Italic\"]\), the coefficient linking the unstable variables of t+1 to the system's correction of expectations in t \!\(\*
StyleBox[\"coffx2zfut\",\nFontSlant->\"Italic\"]\), seven messages giving some information on the chosen type of the solution and the number of stable eigenvalues \!\(\*
StyleBox[\"stabcon\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\".\",\nFontSlant->\"Italic\"]\)"

Begin["`Private`"]

ratexpmodsolver[ss_,tt_,qq_,zz_,gga_,schockcoff_,predet_,n1_,n2_,cutoff_,
      inert_,precise_:10^(-15),
      tol_:Automatic]:=(test=Table[tt[[i,i]]/ss[[i,i]],{i,1,Length[tt]}];
      test=Replace[test,ComplexInfinity\[Rule]10^10,5];
      stabcon=
        Count[Table[TrueQ[Abs[test[[i]]]<cutoff],{i,1,Length[test]}],True];
      s11=SubMatrix[ss,{1,1},{stabcon,stabcon}];Off[Inverse::"luc"];
             s11in=Inverse[SubMatrix[ss,{1,1},{stabcon,stabcon}]];
      s12=SubMatrix[ss,{1,stabcon+1},{stabcon,Length[tt]-stabcon}];
             s21=SubMatrix[ss,{stabcon+1,1},{Length[tt]-stabcon,stabcon}];
      s22=
        SubMatrix[
          ss,{stabcon+1,stabcon+1},{Length[tt]-stabcon,Length[tt]-stabcon}];
      t11=SubMatrix[tt,{1,1},{stabcon,stabcon}];
      t12=SubMatrix[tt,{1,stabcon+1},{stabcon,Length[tt]-stabcon}];
      t21=SubMatrix[tt,{stabcon+1,1},{Length[tt]-stabcon,stabcon}];
      t22=
        SubMatrix[
          tt,{stabcon+1,stabcon+1},{Length[tt]-stabcon,Length[tt]-stabcon}];
      t22in=Inverse[
          SubMatrix[
            tt,{stabcon+1,stabcon+1},{Length[tt]-stabcon,
              Length[tt]-stabcon}]];On[Inverse::"luc"];
      z11=SubMatrix[zz,{1,1},{predet,stabcon}];
      z12=SubMatrix[zz,{1,stabcon+1},{predet,Length[tt]-stabcon}];
      z21=SubMatrix[zz,{predet+1,1},{Length[tt]-predet,stabcon}];
      z22=
        SubMatrix[
          zz,{predet+1,stabcon+1},{Length[tt]-predet,Length[tt]-stabcon}];
      q11con=SubMatrix[ConjugateTranspose[qq],{1,1},{stabcon,stabcon}];
      q12con=SubMatrix[
          ConjugateTranspose[qq],{1,stabcon+1},{stabcon,Length[tt]-stabcon}];
      q21con=SubMatrix[
          ConjugateTranspose[qq],{stabcon+1,1},{Length[tt]-stabcon,stabcon}];
      q22con=
        SubMatrix[
          ConjugateTranspose[qq],{stabcon+1,stabcon+1},{Length[tt]-stabcon,
            Length[tt]-stabcon}];
      z11con=SubMatrix[ConjugateTranspose[zz],{1,1},{stabcon,predet}];
      z12con=SubMatrix[
          ConjugateTranspose[zz],{1,predet+1},{stabcon,Length[tt]-predet}];
      z21con=SubMatrix[
          ConjugateTranspose[zz],{stabcon+1,1},{Length[tt]-stabcon,predet}];
      z22con=
        SubMatrix[
          ConjugateTranspose[zz],{stabcon+1,predet+1},{Length[tt]-stabcon,
            Length[tt]-predet}];
      z1=SubMatrix[zz,{1,1},{stabcon,Length[tt]}];
      q1con=SubMatrix[ConjugateTranspose[qq],{1,1},{stabcon,Length[tt]}];
      q2con=SubMatrix[
          ConjugateTranspose[qq],{stabcon+1,1},{Length[tt]-stabcon,
            Length[tt]}];
      {test1,test2,test3}=SingularValueDecompositionUnstandard[q2con.gga];
      result1=-s22.t22in.q2con.schockcoff;
      Do[result1=
          AppendRows[
            result1,-inert^(i-1)*
              t22.MatrixPower[t22in.s22,i].t22in.q2con.schockcoff],{i,2,
          Length[tt]-stabcon}];
      (*result1=-s22.Inverse[
                IdentityMatrix[Length[t22in]]-
                  inert*t22in.s22].t22in.q2con.schockcoff*);
      {tes1,tes2,tes3}=SingularValueDecompositionUnstandard[result1];
      If[Norm[(IdentityMatrix[Length[test1]]-
                  test1.ConjugateTranspose[test1]).tes1]\[LessEqual]precise*
            Length[gga],mess1=ToString["solution existent"];
        helpa1=LinearSolve[
            q2con.gga,(-s22.Inverse[
                    Chop[IdentityMatrix[Length[t22in]]-
                        inert*t22in.s22]].t22in.q2con.schockcoff)];
        helpa=PseudoInverse[-q2con.gga,
              Tolerance\[Rule]tol].(-inert*
                s22.Inverse[
                    Chop[IdentityMatrix[Length[t22in]]-
                        inert*t22in.s22]].t22in.q2con.schockcoff);ifcon5=1,
        mess1={ToString[
              "no solution existent, absolue value of biggest deviation:"],
            Max[Abs[Chop[(IdentityMatrix[Length[test1]]-
                        test1.ConjugateTranspose[test1]).tes1,precise]]]};
        ifcon5=0];
      If[ifcon5\[NotEqual]1,{coffx1x1,coffx1x2,coffx2x1,coffx2x2,coffx1z,
            coffx1zfut,coffx2z,coffx2zfut,mess1,mess2,mess3,mess4,mess5,mess6,
            mess7,stabcon}={{},{},{},{},{},{},{},{},mess1,mess2,mess3,mess4,
            mess5,mess6,mess7,stabcon},
        {test11b,test22b,test33b}=
          SingularValueDecompositionUnstandard[q2con.gga];{tes11b,tes22b,
            tes33b}=SingularValueDecompositionUnstandard[q1con.gga];
        If[
          MatrixRank[q2con.gga]\[Equal]MatrixRank[
              AppendColumns[q2con.gga,q1con.gga]],
          {test11,test22,test33}=
            SingularValueDecompositionUnstandard[
              ConjugateTranspose[q2con.gga]];{tes11a,tes22a,tes33a}=
            SingularValueDecompositionUnstandard[
              ConjugateTranspose[q1con.gga]];
          If[Norm[(IdentityMatrix[Length[test11]]-
                      test11.ConjugateTranspose[test11]).tes11a]<=
              precise*Length[gga],
            help=q1con.gga.ConjugateTranspose[
                  ConjugateTranspose[test33].Inverse[
                      test22].ConjugateTranspose[test11]];
            mess2=ToString["unique solution"];ifcon1=1,
            mess2=ToString[
                "no unique solution, but solution to sims procedure existent"]\
;ifcon1=1;help=(q1con.gga).PseudoInverse[q2con.gga]],
          mess2=ToString[
              "no unique solution, and no solution using sims procedure"];
          ifcon1=0];
        If[(q1con.gga).PseudoInverse[q2con.gga]\[Equal]ZeroMatrix[
              Dimensions[q1con][[1]],Dimensions[gga][[2]]],
          mess3=ToString["Cf.mess 2"] ,
          mess3=ToString[
              "no influence of expectation error on stable part"]];
        If[MatrixRank[z11]\[Equal]Length[z11],
          mess4=ToString[
              "x2 can be expressed as linear function of x1, Cf. Paul Klein"],
          mess4=ToString[
              "x2 can not be expressed as linear function of x1, Cf. Paul Klein"]\
];
        If[
          q11con.SubMatrix[gga,{1,1},{stabcon,predet}]==
            q12con.SubMatrix[gga,{stabcon+1,1},{Length[tt]-stabcon,predet}],
          
          mess5=ToString[{"influence of endogenous expectations error on stable part is zero",\
(stabcon-ifcon3)*(Length[tt]-stabcon) ,"degrees of freedom"}],
          mess5=ToString[
              "influence of endogenous expectations error on stable does not cancel out"]\
];
        mess6=
          ToString[{"Matrix rank q2con times gga:" ,
              MatrixRank[q2con.gga,Tolerance\[Rule]tol]}];
        mess7=
          ToString[{"matrix rank q1con times gga:", 
              MatrixRank[q1con.gga,Tolerance\[Rule]tol]}];
        If[ifcon1==1,coffx1x1=z11.s11in.(t11.z11con+(t12-help.t22).z21con);
          coffx1x2=z11.s11in.(t11.z12con+(t12-help.t22).z22con);
          coffx2x1=z21.s11in.(t11.z11con+(t12-help.t22).z21con);
          coffx2x2=z21.s11in.(t11.z12con+(t12-help.t22).z22con);
          coffx1z=z11.s11in.(q1con-help.q2con).schockcoff;
          coffx2z=z21.s11in.(q1con-help.q2con).schockcoff;
          
          coffx1zfut={(z11.s11in.(s12-help.s22)-z12),t22in.s22,
              t22in.q2con.schockcoff};
          
          coffx2zfut={(z21.s11in.(s12-help.s22)-z22),t22in.s22,
              t22in.q2con.schockcoff},Null];
        If[ifcon5\[Equal]1&&ifcon1\[Equal]0,
          coffx1x1=z11.s11in.(t11.z11con+(t12-s12.t22).z21con);
          coffx1x2=z11.s11in.(t11.z12con+(t12-s12.t22).z22con);
          coffx2x1=z21.s11in.(t11.z11con+(t12-s12.t22).z21con);
          coffx2x2=z21.s11in.(t11.z12con+(t12-s12.t22).z22con);
          coffx1z=z11.s11in.q1con.schockcoff+z11.s11in.q1con.gga.helpa;
          coffx2z=z21.s11in.q1con.schockcoff+z21.s11in.q1con.gga.helpa;
          
          coffx1zfut={(z11.s11in.(s12-s12.s22)-z12),t22in.s22,
              t22in.q2con.schockcoff};
          
          coffx2zfut={(z21.s11in.(s12-s12.s22)-z22),t22in.s22,
              t22in.q2con.schockcoff};
          mess1=ToString["solution existent: general solution"],
          Null]];{coffx1x1,coffx1x2,coffx2x1,coffx2x2,coffx1z,coffx1zfut,
        coffx2z,coffx2zfut,mess1,mess2,mess3,mess4,mess5,mess6,mess7,
        stabcon});
Clear["ratexpmodsolver`*`*"]
End[]
EndPackage[]