SPAmalg

  [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
                       SPAmalg(h,neq,nlag,nlead,condn,uprbnd)

  Solve a linear perfect foresight model using the matlab eig
  function to find the invariant subspace associated with the big
  roots.  This procedure will fail if the companion matrix is
  defective and does not have a linearly independent set of
  eigenvectors associated with the big roots.
 
  Input arguments:
 
    h         Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
    neq       Number of equations.
    nlag      Number of lags.
    nlead     Number of leads.
    condn     Zero tolerance used as a condition number test
              by numeric_shift and reduced_form.
    uprbnd    Inclusive upper bound for the modulus of roots
              allowed in the reduced form.
 
  Output arguments:
 
    b         Reduced form coefficient matrix (neq,neq*nlag).
    rts       Roots returned by eig.
    ia        Dimension of companion matrix (number of non-trivial
              elements in rts).
    nexact    Number of exact shiftrights.
    nnumeric  Number of numeric shiftrights.
    lgroots   Number of roots greater in modulus than uprbnd.
    aimcode     Return code: see function aimerr.

0001 %  [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
0002 %                       SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
0003 %
0004 %  Solve a linear perfect foresight model using the matlab eig
0005 %  function to find the invariant subspace associated with the big
0006 %  roots.  This procedure will fail if the companion matrix is
0007 %  defective and does not have a linearly independent set of
0008 %  eigenvectors associated with the big roots.
0009 %
0010 %  Input arguments:
0011 %
0012 %    h         Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
0013 %    neq       Number of equations.
0014 %    nlag      Number of lags.
0015 %    nlead     Number of leads.
0016 %    condn     Zero tolerance used as a condition number test
0017 %              by numeric_shift and reduced_form.
0018 %    uprbnd    Inclusive upper bound for the modulus of roots
0019 %              allowed in the reduced form.
0020 %
0021 %  Output arguments:
0022 %
0023 %    b         Reduced form coefficient matrix (neq,neq*nlag).
0024 %    rts       Roots returned by eig.
0025 %    ia        Dimension of companion matrix (number of non-trivial
0026 %              elements in rts).
0027 %    nexact    Number of exact shiftrights.
0028 %    nnumeric  Number of numeric shiftrights.
0029 %    lgroots   Number of roots greater in modulus than uprbnd.
0030 %    aimcode     Return code: see function aimerr.
0031 
0032 % Original author: Gary Anderson
0033 % Original file downloaded from:
0034 % http://www.federalreserve.gov/Pubs/oss/oss4/code.html
0035 % Adapted for Dynare by Dynare Team, in order to deal
0036 % with infinite or nan values.
0037 %
0038 % This code is in the public domain and may be used freely.
0039 % However the authors would appreciate acknowledgement of the source by
0040 % citation of any of the following papers:
0041 %
0042 % Anderson, G. and Moore, G.
0043 % "A Linear Algebraic Procedure for Solving Linear Perfect Foresight
0044 % Models."
0045 % Economics Letters, 17, 1985.
0046 %
0047 % Anderson, G.
0048 % "Solving Linear Rational Expectations Models: A Horse Race"
0049 % Computational Economics, 2008, vol. 31, issue 2, pages 95-113
0050 %
0051 % Anderson, G.
0052 % "A Reliable and Computationally Efficient Algorithm for Imposing the
0053 % Saddle Point Property in Dynamic Models"
0054 % Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
0055 % pages 472-489
0056 
0057 function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
0058                         SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
0059 b=[];rts=[];ia=[];nexact=[];nnumeric=[];lgroots=[];aimcode=[];
0060 if(nlag<1 || nlead<1) 
0061     error('Aim_eig: model must have at least one lag and one lead');
0062 end
0063 % Initialization.
0064 nexact=0;nnumeric=0;lgroots=0;iq=0;aimcode=0;b=0;qrows=neq*nlead;qcols=neq*(nlag+nlead);
0065 bcols=neq*nlag;q=zeros(qrows,qcols);rts=zeros(qcols,1);
0066 [h,q,iq,nexact]=SPExact_shift(h,q,iq,qrows,qcols,neq);
0067 if (iq>qrows)
0068     aimcode = 61;
0069     return;
0070 end
0071 [h,q,iq,nnumeric]=SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn);
0072 if (iq>qrows)
0073     aimcode = 62;
0074     return;
0075 end
0076 [a,ia,js] = SPBuild_a(h,qcols,neq);
0077 if (ia ~= 0)
0078     if any(any(isnan(a))) || any(any(isinf(a))) 
0079         disp('A is NAN or INF')
0080         aimcode=63; 
0081         return 
0082     end 
0083     [w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1));
0084     if flag_trouble==1; 
0085         disp('Problem in SPEIG'); 
0086         aimcode=64;
0087         return
0088     end 
0089     q = SPCopy_w(q,w,js,iq,qrows);
0090 end
0091 test=nexact+nnumeric+lgroots;
0092 if (test > qrows)
0093     aimcode = 3;
0094 elseif (test < qrows)
0095     aimcode = 4;
0096 end
0097 if(aimcode==0)
0098     [nonsing,b]=SPReduced_form(q,qrows,qcols,bcols,neq,condn);
0099     if ( nonsing && aimcode==0)
0100         aimcode =  1;
0101     elseif (~nonsing && aimcode==0)
0102         aimcode =  5;
0103     elseif (~nonsing && aimcode==3)
0104         aimcode = 35;
0105     elseif (~nonsing && aimcode==4)
0106         aimcode = 45;
0107     end
0108 end