1function [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
2    SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
3%  [b,rts,ia,nexact,nnumeric,lgroots,aimcode] = ...
4%                       SPAmalg(h,neq,nlag,nlead,condn,uprbnd)
5%
6%  Solve a linear perfect foresight model using the matlab eig
7%  function to find the invariant subspace associated with the big
8%  roots.  This procedure will fail if the companion matrix is
9%  defective and does not have a linearly independent set of
10%  eigenvectors associated with the big roots.
11%
12%  Input arguments:
13%
14%    h         Structural coefficient matrix (neq,neq*(nlag+1+nlead)).
15%    neq       Number of equations.
16%    nlag      Number of lags.
17%    nlead     Number of leads.
18%    condn     Zero tolerance used as a condition number test
19%              by numeric_shift and reduced_form.
20%    uprbnd    Inclusive upper bound for the modulus of roots
21%              allowed in the reduced form.
22%
23%  Output arguments:
24%
25%    b         Reduced form coefficient matrix (neq,neq*nlag).
26%    rts       Roots returned by eig.
27%    ia        Dimension of companion matrix (number of non-trivial
28%              elements in rts).
29%    nexact    Number of exact shiftrights.
30%    nnumeric  Number of numeric shiftrights.
31%    lgroots   Number of roots greater in modulus than uprbnd.
32%    aimcode     Return code: see function aimerr.
33
34% Original author: Gary Anderson
35% Original file downloaded from:
36% http://www.federalreserve.gov/Pubs/oss/oss4/code.html
37% Adapted for Dynare by Dynare Team, in order to deal
38% with infinite or nan values.
39%
40% This code is in the public domain and may be used freely.
41% However the authors would appreciate acknowledgement of the source by
42% citation of any of the following papers:
43%
44% Anderson, G. and Moore, G.
45% "A Linear Algebraic Procedure for Solving Linear Perfect Foresight
46% Models."
47% Economics Letters, 17, 1985.
48%
49% Anderson, G.
50% "Solving Linear Rational Expectations Models: A Horse Race"
51% Computational Economics, 2008, vol. 31, issue 2, pages 95-113
52%
53% Anderson, G.
54% "A Reliable and Computationally Efficient Algorithm for Imposing the
55% Saddle Point Property in Dynamic Models"
56% Journal of Economic Dynamics and Control, 2010, vol. 34, issue 3,
57% pages 472-489
58
59b=[];rts=[];ia=[];nexact=[];nnumeric=[];lgroots=[];aimcode=[];
60if(nlag<1 || nlead<1)
61    error('Aim_eig: model must have at least one lag and one lead');
62end
63% Initialization.
64nexact=0;nnumeric=0;lgroots=0;iq=0;aimcode=0;b=0;qrows=neq*nlead;qcols=neq*(nlag+nlead);
65bcols=neq*nlag;q=zeros(qrows,qcols);rts=zeros(qcols,1);
66[h,q,iq,nexact]=SPExact_shift(h,q,iq,qrows,qcols,neq);
67if (iq>qrows)
68    aimcode = 61;
69    return
70end
71[h,q,iq,nnumeric]=SPNumeric_shift(h,q,iq,qrows,qcols,neq,condn);
72if (iq>qrows)
73    aimcode = 62;
74    return
75end
76[a,ia,js] = SPBuild_a(h,qcols,neq);
77if (ia ~= 0)
78    if any(any(isnan(a))) || any(any(isinf(a)))
79        disp('A is NAN or INF')
80        aimcode=63;
81        return
82    end
83    [w,rts,lgroots,flag_trouble]=SPEigensystem(a,uprbnd,min(length(js),qrows-iq+1));
84    if flag_trouble==1
85        disp('Problem in SPEIG');
86        aimcode=64;
87        return
88    end
89    q = SPCopy_w(q,w,js,iq,qrows);
90end
91test=nexact+nnumeric+lgroots;
92if (test > qrows)
93    aimcode = 3;
94elseif (test < qrows)
95    aimcode = 4;
96end
97if(aimcode==0)
98    [nonsing,b]=SPReduced_form(q,qrows,qcols,bcols,neq,condn);
99    if ( nonsing && aimcode==0)
100        aimcode =  1;
101    elseif (~nonsing && aimcode==0)
102        aimcode =  5;
103    elseif (~nonsing && aimcode==3)
104        aimcode = 35;
105    elseif (~nonsing && aimcode==4)
106        aimcode = 45;
107    end
108end