gensysToAMA
A Matlab Implementation of the Anderson-Moore Algorithm
Using gensys Input and Output Matrices
Gary S. Anderson
April 24, 2007
Abstract
This note describes a Matlab program for solving linear rational expectation problems. The gensysToAMA program provides a version of the Anderson-Moore algorithm (AMA) that has a matrix interface exactly the same as the gensys program. The code allows the user to invoke the AMA solution code, a copy of their own version of gensys, or a copy of gensys that was available in early 2007. The code can also verify that the solutions obtained using gensys and AMA are equivalent. Timing tests reveal that, except for problems of small dimension, gensysToAMA computes solutions much more quickly than gensys.
Contents
1 Introduction and Summary
2 Usage
2.1 Installation
2.2 Examples using gensysToAMA
2.2.1 Run both algorithms on a small model
2.2.2 Large models much faster with AMA
3 Problem Statement and Notation
4 Algorithmic Solution Concepts
4.1 Anderson-Moore
4.2 Sims
5 Comparing Output and Computation Time
A Appendices
A.1 An early 2007 gensys implementation
A.1.1 gensys2007
A.1.2 qzdiv2007
A.1.3 qzswitch2007
A.2 convertFromGensysIn
A.3 convertToGensysOut
A.4 Linear Algebra for Comparisons
1 Introduction and Summary
This note describes a Matlab program for solving linear rational expectation problems. The gensysToAMA program provides a version of the Anderson-Moore algorithm (AMA) that has a matrix interface exactly the same as the gensys program. The code allows the user to invoke the AMA solution code, a copy of their own version of gensys, or a copy of gensys that was available in early 2007. The code can also verify that the solutions obtained using gensys and AMA are equivalent. Timing tests reveal that, except for problems of small dimension, gensysToAMA computes solutions much more quickly than gensys.
2 Usage
2.1 Installation
1. unzip the files into a directory (someDir) accessible by Matlab. This will create a directory, gen- sysToAMADist, containing the gensysToAMA programs and some example .mat input matrix files.
2. start matlab
3. place the gensysToAMA directory on the Matlab path using
addpath someDir/gensysToAMADist
4. during a matlab session, you can run a quick test of the installation:
type
>>isGensysToAMAOK
to verify the gensysToAMA program functions correctly. After a few seconds, you should get a
SUCCESS message.
2.2 Examples using gensysToAMA
The installation directory provides a number of .mat files that contain input matrices for the gensys (or gensysToAMA) program. These models vary in size and required computation time. For example, the two canada .mat files characterize the largest models and require the most computation time. Typing help gensysToAMA provides information on the gensysToAMA inputs and outputs.
>> help gensysToAMA
function [G1,CC,impact,fmat,fwt,ywt,gev,eu]=gensysToAMA(g0,g1,cc,psi,pi,div,
optionalArg)
gensys interface to both gensys and the Anderson-Moore algorithm.
Just as with gensys, system given as
g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
with z an exogenous variable process and eta being endogenously determined
one-step-ahead expectational errors. Returned system is
y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
If z(t) is i.i.d., the last term drops out.
If div is omitted from argument list, a div>1 is calculated.
eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
By Christopher A. Sims
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
when called with no optional args,
program first tries to find gensys on the matlab path
if that fails, the program runs gensys2007 in the gensysToAMA directory
an optional string argument may follow the original gensys arguments
gensys run gensys
program first tries to find gensys on the matlab path
if that fails, the program runs gensys2007 in the gensysToAMA directory
gensys2007 run gensys
the program runs gensys2007 in the gensysToAMA directory
ama run the anderson-moore algorithm with gensys inputs and outputs
both run the anderson-moore algorithm and the gensys program
verify that solutions are equivalent print out execution times
2.2.1 Run both algorithms on a small model
gensysToAMA Example
>> load inp
>> who
Your variables are:
cc div g0 g1 hmat nlags pi psi
>> [ggg,ccc,iii,fff,ffw,yyw,gev,eue]=gensysToAMA(g0,g1,cc,psi,pi,1.0,both);
gensysToAMA:running both ama and gensys for comparison
problem dimensions: g0:10 x 10, psi:10 x 2, pi:10 x 5
gensysToAMA:running ama
gensysToAMA:converting ama output to gensys format
gensysToAMA:running gensys
gensysToAMA: trying gensys on your matlab path
gensysToAMA: that failed, using gensys2007 in gensysToAMA dir
gensysToAMA:runs complete
no difference in sims and AMA results
AMATime=1.562500e-002 AMAFTime=0 convertToTime=0 convertFromTime=0 genSysTime=0
>>
2.2.2 Large models much faster with AMA
gensysToAMA Example
>> load canada2lagas1.mat
>> [ggg,ccc,iii,fff,ffw,yyw,gev,eue]=gensysToAMA(g0,g1,cc,psi,pi,1.0,both);
gensysToAMA:running both ama and gensys for comparison
problem dimensions: g0:222 x 222, psi:222 x 1, pi:222 x 111
gensysToAMA:running ama
gensysToAMA:converting ama output to gensys format
gensysToAMA:running gensys
gensysToAMA: trying gensys on your matlab path
gensysToAMA: that failed, using gensys2007 in gensysToAMA dir
gensysToAMA:runs complete
no difference in sims and AMA results
AMATime=8.593750e-001 AMAFTime=4.687500e-002
convertToTime=3.125000e-002 convertFromTime=9.375000e-002
genSysTime=4.273438e+001
>>
3 Problem Statement and Notation
These algorithms compute solutions for models of the form
Formula 1: 
with initial conditions, if any, given by constraints of the form
Formula 2: 
where both tau and theta are non-negative, and xt is an L dimensional vector of endogenous variables with
Formula 3: 
and zt is a k dimensional vector of exogenous variables.
Solutions can be cast in the form
Formula 5: 
Given any algorithm that computes the Bi, one can easily compute other quantities useful for character- izing the impact of exogenous variables. For models with tau = theta = 1 the formulae are especially simple. Let
Formulae 6-7: 
We can write
Formula 8: 
and when
Formulae 9-11: 
Consult Anderson [1997] for other useful formulae concerning rational expectations model solutions.
4 Algorithmic Solution Concepts
The following sections present the inputs and outputs for each of the algorithms for the following simple example:
Formulae 12-13: 
INPUTS
Formula 14: 
|
|
|
|
| x sub(t-tau), ..., x sub(t), ... x sub(t+theta) | Model Variables | L(tau + theta) x 1 |
| z sub t | Exogenous Variables | M x 1 |
| theta | Longest Lead | 1 x 1 |
| tau | Longest Lag | 1 x 1 |
| H sub i | Structural Coefficients Matrix | (L x L)(tau + theta + 1) |
| Psi | Exogenous Shock Coefficients Matrix | L x M |
| Upsilon | Optional Exogenous VAR Coefficients Matrix (z sub(t+1) = Upsilon*z sub(t) ) | M x M |
Formula 16: 
|
|
|
|
| B | reduced form coefficients matrix | L x L(tau + theta) |
| Phi | exogenous shock scaling matrkx | L x L |
| F | exogenous shock transfer matrix | L*theta x L*theta |
| undercase upsilon | autoregressive shock transfer matrix when z sub(t+1) equals uppercase Upsilon times z sub(t) the infinite sum simplifies to give x sub(t) equals B times the column vector containing elements x sub(t-tau) to x sub(t-1), plus undercase upsilon times z sub(t) | L x M |
Anderson-Moore input:
AIM Modeling Language Input
Parameter File Input
produces output:
Usage Notes for Anderson-Moore Algorithm
1. "Align" model variables so that the data history (without applying model equations), completely determines all of x sub (t-1), but none of xt.
2. Develop a "model file" containing the model equations written in the "AIM modeling language"
3. Apply the model pre-processor to create MATLAB programs for initializing the algorithm's input matrix,(H). Create Psi and, optionally, Upsilon matrices.
4. Execute the MATLAB programs to generate B, phi ,F and optionally (lower case) upsilon.
Users can obtain code for the algorithm and the preprocessor from the author (at http://www.bog.frb.fed.us/pubs/oss/oss4/aimindex.html, July, 1999.)
4.2 Sims
Formula 21: 
|
|
|
|
| y sub(t) | State Variables | L x 1 |
| z sub(t) | Exogenous Variables | M sub(1) x 1 |
| Eta sub(t) | Expectational Error | M sub(2) x 1 |
| Gamma sub(0) | Structural Coefficients Matrix | L x L |
| Gamma sub(1) | Structural Coefficients Matrix | L x L |
| C | Constants | L x 1 |
| Psi | Structural Exogenous Variables Coefficients Matrix | L x M sub(1) |
| Pi | Structural Exogenous Variables Coefficients Matrix | L x M sub(2) |
Outputs
Formula 23: 
|
|
|
|
| Theta sub(1) | L x L | |
| Theta sub(c) | L x 1 | |
| Theta sub(0) | L x M sub(1) | |
| Theta sub(y) | L x M sub(2) | |
| Theta sub(f) | M sub(2) x M sub(2) | |
| Theta sub(z) | M sub(2) x M sub(1) |
Example Input/Output:
Sims input: 
produces output: 
5 Comparing Output and Computation Time
Anderson and Moore describe their algorithm in Anderson and Moore [1983]. Sims describes his algorithm in Sims [1996], Although they attack a the same class of problems their notation is different. Their computer codes reflect these notational differences. Section A.2 presents the code for converting gensys style input into the form expected by the Anderson-Moore algorithm. style input Section A.3 presents the code for converting Anderson-Moore style output into the form generated by the gensys program. When running the code using the 'both' option, the code computes the two-norm of the difference between the gensys style output of the two programs. (The code does not compare the outputs when the rational expectations model does not have a unique saddle point solution.) In particular, the code verifies that the following two-norms are smaller than 1E-8.
Comparing fmat, fwt, ywt must account for that fact that applying any similarity transformation to fmat and correspondingly adjusting ywt and fwt produces an equivalent solution. Consequently, the code verifies that
References
Gary Anderson. A reliable and computationally efficient algorithm for imposing the saddle point property in dynamic models. URL http://www.federalreserve.gov/pubs/oss/oss4/aimindex.html. Unpublished Manuscript, Board of Governors of the Federal Reserve System., 1997.
Gary Anderson and George Moore. An efficient procedure for solving linear perfect foresight models. Un- published Manuscript, Board of Governors of the Federal Reserve System., 1983.
Christopher A. Sims. Solving linear rational expectations models. Seminar paper, 1996.
A Appendices
A.1 An early 2007 gensys implementation
A.1.1 gensys2007
1 function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys2007(g0,g1,c,psi,pi,div)
2 % function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys2007(g0,g1,c,psi,pi,div)
3 % frozen copy of gensys for gensysToAMA
4 %for use in case gensys not found on user path
5 % System given as
6 % g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
7 % with z an exogenous variable process and eta being endogenously determined
8 % one-step-ahead expectational errors. Returned system is
9 % y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
10 % If z(t) is i.i.d., the last term drops out.
11 % If div is omitted from argument list, a div>1 is calculated.
12 % eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
13 % existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
14 % By Christopher A. Sims
15 % Corrected 10/28/96 by CAS
16 eu=[0;0];
17 realsmall=1e-6;
18 fixdiv=(nargin==6);
19 n=size(g0,1);
20 [a b q z v]=qz(g0,g1);
21 if ~fixdiv, div=1.01; end
22 nunstab=0;
23 zxz=0;
24 for i=1:n
25 % ------------------div calc------------
26 if ~fixdiv
27 if abs(a(i,i)) > 0
28 divhat=abs(b(i,i))/abs(a(i,i));
29 % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
30 % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
31 % and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
32 if 1+realsmall<="div" 33 div=.5*(1+divhat);
34 end
35 end
36 end
37 % ----------------------------------------
38 nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
39 if abs(a(i,i))
41 end
42 end
43 div ;
44 nunstab;
45 if ~zxz
46 %alejandro indicates ordqz faster [a b q z]=qzdiv2007(div,a,b,q,z);
47 [a b q z]=qzdiv2007(div,a,b,q,z);
48 % [a b q z]=ordqz(div,a,b,q,z);
49 end
50 gev=[diag(a) diag(b)];
51 if zxz
52 disp(Coincident zeros. Indeterminacy and/or nonexistence.)
53 eu=[-2;-2];
54 % correction added 7/29/2003. Otherwise the failure to set output
55 % arguments leads to an error message and no output (including eu).
56 G1=[];C=[];impact=[];fmat=[];fwt=[];ywt=[];gev=[];
57 return
58 end
59 q1=q(1:n-nunstab,:);
60 q2=q(n-nunstab+1:n,:);
61 z1=z(:,1:n-nunstab);
62 z2=z(:,n-nunstab+1:n);
63 a2=a(n-nunstab+1:n,n-nunstab+1:n);
64 b2=b(n-nunstab+1:n,n-nunstab+1:n);
65 etawt=q2*pi;
66 % zwt=q2*psi;
67 [ueta,deta,veta]=svd(etawt);
68 md=min(size(deta));
69 bigev=find(diag(deta(1:md,1:md))>realsmall);
70 ueta=ueta(:,bigev);
71 veta=veta(:,bigev);
72 deta=deta(bigev,bigev);
73 % ------ corrected code, 3/10/04
74 eu(1) = length(bigev)>=nunstab;
75 % ------ Code below allowed "existence" in cases where the initial lagged state was free to take on v
76 % ------ inconsistent with existence, so long as the state could w.p.1 remain consistent with a stabl
77 % ------ if its initial lagged value was consistent with a stable solution. This is a mistake, thoug
78 % ------ are situations where we would like to know that this "existence for restricted initial state
79 %% [uz,dz,vz]=svd(zwt);
80 %% md=min(size(dz));
81 %% bigev=find(diag(dz(1:md,1:md))>realsmall);
82 %% uz=uz(:,bigev);
83 %% vz=vz(:,bigev);
84 %% dz=dz(bigev,bigev);
85 %% if isempty(bigev)
86 %% exist=1;
87 %% else
88 %% exist=norm(uz-ueta*ueta*uz) < realsmall*n;
89 %% end
90 %% if ~isempty(bigev)
91 %% zwtx0=b2\zwt;
92 %% zwtx=zwtx0;
93 %% M=b2\a2;
94 %% for i=2:nunstab
95 %% zwtx=[M*zwtx zwtx0];
96 %% end
97 %% zwtx=b2*zwtx;
98 %% [ux,dx,vx]=svd(zwtx);
99 %% md=min(size(dx));
100 %% bigev=find(diag(dx(1:md,1:md))>realsmall);
101 %% ux=ux(:,bigev);
102 %% vx=vx(:,bigev);
103 %% dx=dx(bigev,bigev);
104 %% existx=norm(ux-ueta*ueta*ux) < realsmall*n;
105 %% else
106 %% existx=1;
107 %% end
108 % ----------------------------------------------------
109 % Note that existence and uniqueness are not just matters of comparing
110 % numbers of roots and numbers of endogenous errors. These counts are
111 % reported below because usually they point to the source of the problem.
112 % ------------------------------------------------------
113 [ueta1,deta1,veta1]=svd(q1*pi);
114 md=min(size(deta1));
115 bigev=find(diag(deta1(1:md,1:md))>realsmall);
116 ueta1=ueta1(:,bigev);
117 veta1=veta1(:,bigev);
118 deta1=deta1(bigev,bigev);
119 %% if existx | nunstab==0
120 %% %disp(solution exists);
121 %% eu(1)=1;
122 %% else
123 %% if exist
124 %% %disp(solution exists for unforecastable z only);
125 %% eu(1)=-1;
126 %% %else
127 %% %fprintf(1,No solution. %d unstable roots. %d endog errors.\n,nunstab,size(ueta1,2));
128 %% end
129 %% %disp(Generalized eigenvalues)
130 %% %disp(gev);
131 %% %md=abs(diag(a))>realsmall;
132 %% %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
133 %% %disp(ev)
134 %% % return;
135 %% end
136 if isempty(veta1)
137 unique=1;
138 else
139 unique=norm(veta1-veta*veta*veta1)140 end
141 if unique
142 %disp(solution unique);
143 eu(2)=1;
144 else
145 fprintf(1,Indeterminacy. %d loose endog errors.\n,size(veta1,2)-size(veta,2));
146 %disp(Generalized eigenvalues)
147 %disp(gev);
148 %md=abs(diag(a))>realsmall;
149 %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
150 %disp(ev)
151 % return;
152 end
153 tmat = [eye(n-nunstab) -(ueta*(deta\veta)*veta1*deta1*ueta1)];
154 G0= [tmat*a; zeros(nunstab,n-nunstab) eye(nunstab)];
155 G1= [tmat*b; zeros(nunstab,n)];
156 % ----------------------
157 % G0 is always non-singular because by construction there are no zeros on
158 % the diagonal of a(1:n-nunstab,1:n-nunstab), which forms G0s ul corner.
159 % -----------------------
160 G0I=inv(G0);
161 G1=G0I*G1;
162 usix=n-nunstab+1:n;
163 C=G0I*[tmat*q*c;(a(usix,usix)-b(usix,usix))\q2*c];
164 impact=G0I*[tmat*q*psi;zeros(nunstab,size(psi,2))];
165 fmat=b(usix,usix)\a(usix,usix);
166 fwt=-b(usix,usix)\q2*psi;
167 ywt=G0I(:,usix);
168 % -------------------- above are output for system in terms of zy -------
169 G1=real(z*G1*z);
170 C=real(z*C);
171 impact=real(z*impact);
172 % Correction 10/28/96: formerly line below had real(z*ywt) on rhs, an error.
173 ywt=z*ywt;
A.1.2 qzdiv2007
1 function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys2007(g0,g1,c,psi,pi,div)
2 % function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys2007(g0,g1,c,psi,pi,div)
3 % frozen copy of gensys for gensysToAMA
4 %for use in case gensys not found on user path
5 % System given as
6 % g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
7 % with z an exogenous variable process and eta being endogenously determined
8 % one-step-ahead expectational errors. Returned system is
9 % y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
10 % If z(t) is i.i.d., the last term drops out.
11 % If div is omitted from argument list, a div>1 is calculated.
12 % eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
13 % existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
14 % By Christopher A. Sims
15 % Corrected 10/28/96 by CAS
16 eu=[0;0];
17 realsmall=1e-6;
18 fixdiv=(nargin==6);
19 n=size(g0,1);
20 [a b q z v]=qz(g0,g1);
21 if ~fixdiv, div=1.01; end
22 nunstab=0;
23 zxz=0;
24 for i=1:n
25 % ------------------div calc------------
26 if ~fixdiv
27 if abs(a(i,i)) > 0
28 divhat=abs(b(i,i))/abs(a(i,i));
29 % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
30 % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
31 % and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
32 if 1+realsmall<="div" 33 div=.5*(1+divhat);
34 end
35 end
36 end
37 % ----------------------------------------
38 nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
39 if abs(a(i,i))
41 end
42 end
43 div ;
44 nunstab;
45 if ~zxz
46 %alejandro indicates ordqz faster [a b q z]=qzdiv2007(div,a,b,q,z);
47 [a b q z]=qzdiv2007(div,a,b,q,z);
48 % [a b q z]=ordqz(div,a,b,q,z);
49 end
50 gev=[diag(a) diag(b)];
51 if zxz
52 disp(Coincident zeros. Indeterminacy and/or nonexistence.)
53 eu=[-2;-2];
54 % correction added 7/29/2003. Otherwise the failure to set output
55 % arguments leads to an error message and no output (including eu).
56 G1=[];C=[];impact=[];fmat=[];fwt=[];ywt=[];gev=[];
57 return
58 end
59 q1=q(1:n-nunstab,:);
60 q2=q(n-nunstab+1:n,:);
61 z1=z(:,1:n-nunstab);
62 z2=z(:,n-nunstab+1:n);
63 a2=a(n-nunstab+1:n,n-nunstab+1:n);
64 b2=b(n-nunstab+1:n,n-nunstab+1:n);
65 etawt=q2*pi;
66 % zwt=q2*psi;
67 [ueta,deta,veta]=svd(etawt);
68 md=min(size(deta));
69 bigev=find(diag(deta(1:md,1:md))>realsmall);
70 ueta=ueta(:,bigev);
71 veta=veta(:,bigev);
72 deta=deta(bigev,bigev);
73 % ------ corrected code, 3/10/04
74 eu(1) = length(bigev)>=nunstab;
75 % ------ Code below allowed "existence" in cases where the initial lagged state was free to take on v
76 % ------ inconsistent with existence, so long as the state could w.p.1 remain consistent with a stabl
77 % ------ if its initial lagged value was consistent with a stable solution. This is a mistake, thoug
78 % ------ are situations where we would like to know that this "existence for restricted initial state
79 %% [uz,dz,vz]=svd(zwt);
80 %% md=min(size(dz));
81 %% bigev=find(diag(dz(1:md,1:md))>realsmall);
82 %% uz=uz(:,bigev);
83 %% vz=vz(:,bigev);
84 %% dz=dz(bigev,bigev);
85 %% if isempty(bigev)
86 %% exist=1;
87 %% else
88 %% exist=norm(uz-ueta*ueta*uz) < realsmall*n;
89 %% end
90 %% if ~isempty(bigev)
91 %% zwtx0=b2\zwt;
92 %% zwtx=zwtx0;
93 %% M=b2\a2;
94 %% for i=2:nunstab
95 %% zwtx=[M*zwtx zwtx0];
96 %% end
97 %% zwtx=b2*zwtx;
98 %% [ux,dx,vx]=svd(zwtx);
99 %% md=min(size(dx));
100 %% bigev=find(diag(dx(1:md,1:md))>realsmall);
101 %% ux=ux(:,bigev);
102 %% vx=vx(:,bigev);
103 %% dx=dx(bigev,bigev);
104 %% existx=norm(ux-ueta*ueta*ux) < realsmall*n;
105 %% else
106 %% existx=1;
107 %% end
108 % ----------------------------------------------------
109 % Note that existence and uniqueness are not just matters of comparing
110 % numbers of roots and numbers of endogenous errors. These counts are
111 % reported below because usually they point to the source of the problem.
112 % ------------------------------------------------------
113 [ueta1,deta1,veta1]=svd(q1*pi);
114 md=min(size(deta1));
115 bigev=find(diag(deta1(1:md,1:md))>realsmall);
116 ueta1=ueta1(:,bigev);
117 veta1=veta1(:,bigev);
118 deta1=deta1(bigev,bigev);
119 %% if existx | nunstab==0
120 %% %disp(solution exists);
121 %% eu(1)=1;
122 %% else
123 %% if exist
124 %% %disp(solution exists for unforecastable z only);
125 %% eu(1)=-1;
126 %% %else
127 %% %fprintf(1,No solution. %d unstable roots. %d endog errors.\n,nunstab,size(ueta1,2));
128 %% end
129 %% %disp(Generalized eigenvalues)
130 %% %disp(gev);
131 %% %md=abs(diag(a))>realsmall;
132 %% %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
133 %% %disp(ev)
134 %% % return;
135 %% end
136 if isempty(veta1)
137 unique=1;
138 else
139 unique=norm(veta1-veta*veta*veta1)140 end
141 if unique
142 %disp(solution unique);
143 eu(2)=1;
144 else
145 fprintf(1,Indeterminacy. %d loose endog errors.\n,size(veta1,2)-size(veta,2));
146 %disp(Generalized eigenvalues)
147 %disp(gev);
148 %md=abs(diag(a))>realsmall;
149 %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
150 %disp(ev)
151 % return;
152 end
153 tmat = [eye(n-nunstab) -(ueta*(deta\veta)*veta1*deta1*ueta1)];
154 G0= [tmat*a; zeros(nunstab,n-nunstab) eye(nunstab)];
155 G1= [tmat*b; zeros(nunstab,n)];
156 % ----------------------
157 % G0 is always non-singular because by construction there are no zeros on
158 % the diagonal of a(1:n-nunstab,1:n-nunstab), which forms G0s ul corner.
159 % -----------------------
160 G0I=inv(G0);
161 G1=G0I*G1;
162 usix=n-nunstab+1:n;
163 C=G0I*[tmat*q*c;(a(usix,usix)-b(usix,usix))\q2*c];
164 impact=G0I*[tmat*q*psi;zeros(nunstab,size(psi,2))];
165 fmat=b(usix,usix)\a(usix,usix);
166 fwt=-b(usix,usix)\q2*psi;
167 ywt=G0I(:,usix);
168 % -------------------- above are output for system in terms of zy -------
169 G1=real(z*G1*z);
170 C=real(z*C);
171 impact=real(z*impact);
172 % Correction 10/28/96: formerly line below had real(z*ywt) on rhs, an error.
173 ywt=z*ywt;
A.1.3 qzswitch2007
1 function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys2007(g0,g1,c,psi,pi,div)
2 % function [G1,C,impact,fmat,fwt,ywt,gev,eu]=gensys2007(g0,g1,c,psi,pi,div)
3 % frozen copy of gensys for gensysToAMA
4 %for use in case gensys not found on user path
5 % System given as
6 % g0*y(t)=g1*y(t-1)+c+psi*z(t)+pi*eta(t),
7 % with z an exogenous variable process and eta being endogenously determined
8 % one-step-ahead expectational errors. Returned system is
9 % y(t)=G1*y(t-1)+C+impact*z(t)+ywt*inv(I-fmat*inv(L))*fwt*z(t+1) .
10 % If z(t) is i.i.d., the last term drops out.
11 % If div is omitted from argument list, a div>1 is calculated.
12 % eu(1)=1 for existence, eu(2)=1 for uniqueness. eu(1)=-1 for
13 % existence only with not-s.c. z; eu=[-2,-2] for coincident zeros.
14 % By Christopher A. Sims
15 % Corrected 10/28/96 by CAS
16 eu=[0;0];
17 realsmall=1e-6;
18 fixdiv=(nargin==6);
19 n=size(g0,1);
20 [a b q z v]=qz(g0,g1);
21 if ~fixdiv, div=1.01; end
22 nunstab=0;
23 zxz=0;
24 for i=1:n
25 % ------------------div calc------------
26 if ~fixdiv
27 if abs(a(i,i)) > 0
28 divhat=abs(b(i,i))/abs(a(i,i));
29 % bug detected by Vasco Curdia and Daria Finocchiaro, 2/25/2004 A root of
30 % exactly 1.01 and no root between 1 and 1.02, led to div being stuck at 1.01
31 % and the 1.01 root being misclassified as stable. Changing < to <= below fixes this.
32 if 1+realsmall<="div" 33 div=.5*(1+divhat);
34 end
35 end
36 end
37 % ----------------------------------------
38 nunstab=nunstab+(abs(b(i,i))>div*abs(a(i,i)));
39 if abs(a(i,i))
41 end
42 end
43 div ;
44 nunstab;
45 if ~zxz
46 %alejandro indicates ordqz faster [a b q z]=qzdiv2007(div,a,b,q,z);
47 [a b q z]=qzdiv2007(div,a,b,q,z);
48 % [a b q z]=ordqz(div,a,b,q,z);
49 end
50 gev=[diag(a) diag(b)];
51 if zxz
52 disp(Coincident zeros. Indeterminacy and/or nonexistence.)
53 eu=[-2;-2];
54 % correction added 7/29/2003. Otherwise the failure to set output
55 % arguments leads to an error message and no output (including eu).
56 G1=[];C=[];impact=[];fmat=[];fwt=[];ywt=[];gev=[];
57 return
58 end
59 q1=q(1:n-nunstab,:);
60 q2=q(n-nunstab+1:n,:);
61 z1=z(:,1:n-nunstab);
62 z2=z(:,n-nunstab+1:n);
63 a2=a(n-nunstab+1:n,n-nunstab+1:n);
64 b2=b(n-nunstab+1:n,n-nunstab+1:n);
65 etawt=q2*pi;
66 % zwt=q2*psi;
67 [ueta,deta,veta]=svd(etawt);
68 md=min(size(deta));
69 bigev=find(diag(deta(1:md,1:md))>realsmall);
70 ueta=ueta(:,bigev);
71 veta=veta(:,bigev);
72 deta=deta(bigev,bigev);
73 % ------ corrected code, 3/10/04
74 eu(1) = length(bigev)>=nunstab;
75 % ------ Code below allowed "existence" in cases where the initial lagged state was free to take on
76 % ------ inconsistent with existence, so long as the state could w.p.1 remain consistent with a stabl
77 % ------ if its initial lagged value was consistent with a stable solution. This is a mistake, thoug
78 % ------ are situations where we would like to know that this "existence for restricted initial state
79 %% [uz,dz,vz]=svd(zwt);
80 %% md=min(size(dz));
81 %% bigev=find(diag(dz(1:md,1:md))>realsmall);
82 %% uz=uz(:,bigev);
83 %% vz=vz(:,bigev);
84 %% dz=dz(bigev,bigev);
85 %% if isempty(bigev)
86 %% exist=1;
87 %% else
88 %% exist=norm(uz-ueta*ueta*uz) < realsmall*n;
89 %% end
90 %% if ~isempty(bigev)
91 %% zwtx0=b2\zwt;
92 %% zwtx=zwtx0;
93 %% M=b2\a2;
94 %% for i=2:nunstab
95 %% zwtx=[M*zwtx zwtx0];
96 %% end
97 %% zwtx=b2*zwtx;
98 %% [ux,dx,vx]=svd(zwtx);
99 %% md=min(size(dx));
100 %% bigev=find(diag(dx(1:md,1:md))>realsmall);
101 %% ux=ux(:,bigev);
102 %% vx=vx(:,bigev);
103 %% dx=dx(bigev,bigev);
104 %% existx=norm(ux-ueta*ueta*ux) < realsmall*n;
105 %% else
106 %% existx=1;
107 %% end
108 % ----------------------------------------------------
109 % Note that existence and uniqueness are not just matters of comparing
110 % numbers of roots and numbers of endogenous errors. These counts are
111 % reported below because usually they point to the source of the problem.
112 % ------------------------------------------------------
113 [ueta1,deta1,veta1]=svd(q1*pi);
114 md=min(size(deta1));
115 bigev=find(diag(deta1(1:md,1:md))>realsmall);
116 ueta1=ueta1(:,bigev);
117 veta1=veta1(:,bigev);
118 deta1=deta1(bigev,bigev);
119 %% if existx | nunstab==0
120 %% %disp(solution exists);
121 %% eu(1)=1;
122 %% else
123 %% if exist
124 %% %disp(solution exists for unforecastable z only);
125 %% eu(1)=-1;
126 %% %else
127 %% %fprintf(1,No solution. %d unstable roots. %d endog errors.\n,nunstab,size(ueta1,2));
128 %% end
129 %% %disp(Generalized eigenvalues)
130 %% %disp(gev);
131 %% %md=abs(diag(a))>realsmall;
132 %% %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
133 %% %disp(ev)
134 %% % return;
135 %% end
136 if isempty(veta1)
137 unique=1;
138 else
139 unique=norm(veta1-veta*veta*veta1)140 end
141 if unique
142 %disp(solution unique);
143 eu(2)=1;
144 else
145 fprintf(1,Indeterminacy. %d loose endog errors.\n,size(veta1,2)-size(veta,2));
146 %disp(Generalized eigenvalues)
147 %disp(gev);
148 %md=abs(diag(a))>realsmall;
149 %ev=diag(md.*diag(a)+(1-md).*diag(b))\ev;
150 %disp(ev)
151 % return;
152 end
153 tmat = [eye(n-nunstab) -(ueta*(deta\veta)*veta1*deta1*ueta1)];
154 G0= [tmat*a; zeros(nunstab,n-nunstab) eye(nunstab)];
155 G1= [tmat*b; zeros(nunstab,n)];
156 % ----------------------
157 % G0 is always non-singular because by construction there are no zeros on
158 % the diagonal of a(1:n-nunstab,1:n-nunstab), which forms G0s ul corner.
159 % -----------------------
160 G0I=inv(G0);
161 G1=G0I*G1;
162 usix=n-nunstab+1:n;
163 C=G0I*[tmat*q*c;(a(usix,usix)-b(usix,usix))\q2*c];
164 impact=G0I*[tmat*q*psi;zeros(nunstab,size(psi,2))];
165 fmat=b(usix,usix)\a(usix,usix);
166 fwt=-b(usix,usix)\q2*psi;
167 ywt=G0I(:,usix);
168 % -------------------- above are output for system in terms of zy -------
169 G1=real(z*G1*z);
170 C=real(z*C);
171 impact=real(z*impact);
172 % Correction 10/28/96: formerly line below had real(z*ywt) on rhs, an error.
173 ywt=z*ywt;
A.2 convertFromGensysIn
1 function [theHM,theH0,theHP]=convertFromGensysIn(g0,g1,pi)
2 %function [theHM,theH0,theHP]=convertFromGensysIn(g0,g1,pi)
3 gDim=size(g0,1);
4 piCol=size(pi,2);
5 theHM=sparse([...
6 -g1,zeros(gDim,piCol)
7 zeros(piCol,gDim+piCol)]);
8 theH0=sparse([...
9 g0,-pi;...
10 zeros(piCol,gDim+piCol)]);
11 theHP=sparse([...
12 zeros(gDim,gDim+piCol);...
13 zeros(piCol,gDim),eye(piCol)]);
A.3 convertToGensysOut
1 function [CC,G1,impact,ywt,fmat,fwt]=...
2 convertToGensysOut(bb,phi,theF,cc,g0,g1,psi,ncpi)
3 %function [CC,G1,impact,ywt,fmat,fwt]=...
4 %convertToGensysOut(bb,phi,theF,cc,g0,g1,psi,ncpi)
5
6 [nr,nc]=size(g1);
7 [nrpsi,ncpsi]=size(psi);
8 stateDim=size(bb,2)-ncpi;
9 G1=bb(1:nr,1:nc);
10
11 ststate=(g0-g1)\cc;
12 CC=(eye(nr)-G1)*ststate;
13
14 thePsi=[psi;zeros(ncpi,ncpsi)];
15 aa=phi*thePsi;
16 impact=aa(1:nr,:);
17 %no unique way to represent these components
18 [ywt,fmat,fwt]=smallF(theF,aa,stateDim);
19
20 function [onLeft,inMiddle,onRight]=smallF(anF,bigPhi,nn)
21 [fRows,fCols]=size(anF);
22 lilFL=anF(nn+1:fRows,nn+1:fRows);
23 uu=null(full(lilFL));
24 theNull=size(uu,2);
25 eOpts.disp=0;
26 lilFU=anF(1:nn,nn+1:fRows);
27 onLeft=lilFU;
28 onRight=bigPhi(nn+1:fRows,:);
29 inMiddle=lilFL;
A.4 Linear Algebra for Comparisons Under construction. Check code in gensysToAMA .
