skip to main navigation skip to secondary navigation skip to content
Board of Governors of the Federal Reserve System
skip to content

Occasional Staff StudiesThe Anderson-Moore Algorithm (AMA)


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: The summation (from i equals negative tau through theta) of H sub i times x sub (t+1) equals psi times z sub t, where t equals 0, 1, ... to infinity

with initial conditions, if any, given by constraints of the form

Formula 2: x sub i equals x sub i superscript data , where i equals negative tau through -1

where both tau and theta are non-negative, and xt is an L dimensional vector of endogenous variables with

Formula 3: the limit as t goes to infinity of the norm of x sub t is less than infinity

and zt is a k dimensional vector of exogenous variables.

Solutions can be cast in the form

Formula 5: x sub t minus x star equals the summation from i equals negative tau through -1 of B sub i times (x sub t+1 minus x star)

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: formula 6: phi = the inverse of (H sub 0 + the product of H sub 1 and B sub 1; formula 7: F equals the product of phi and H sub 1 and B

We can write

Formula 8: x sub t minus x star equals B sub 1 times ((x sub t minus 1) minus x star) plus the summation where s equals 0 through infinity of (F to the s power) times phi times psi times z sub (t + s)

and when

Formulae 9-11: formula 9: z sub t plus 1 equals Upsilon times z sub t ; formula 10: vec(undercase upsilon) equals the inverse of (I minus the tensor product of Upsilon to the T and F) times vec(phi times psi)

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: Formula 12: V sub (t+1) equals (1+R) times V sub t, minus D sub (t+1); Formula 13: D equals (1-delta) times D sub (t minus 1)

INPUTS

Formula 14: The summation where i goes from negative tau to theta of H sub i times x sub (t+1) equals the product of phi and z sub t

Model Variable
Description
Dimensions
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: x sub t equals B times the column vector containing values x sub (t minus tau) through x sub (t minus 1) plus the row vector containing (0,..., 0, I) times the summation where s goes from zero to infinity of the product of F to the s and the column vector containing the elemtents 0 and phi times psi times z sub t + s

Model Variable
Description
Dimensions
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

H equals the 2 by 6 matrix containing the elements: 0, 0, -1.1, 0, 1, 1., 0, -.7, 0,1,0,0; psi equals the 2 by 2 matrix containing the elements 4., 1., 3., and -2; Upsilon equals the 2 by 2 matrix containing the elements .9, .1, .05, .2

produces output:

B equals the 2 by 2 matrix containing the elements 0., 1.225, 0., .7; F equals the 2 by 2 matrix containing elements .909091, .909091, 0., 0.; phi equals the 2 by 2 matrix containing the elements .909091, 1.75, 0., 1.; phi times psi equals the 2 by 2 matrix containing the elements 1.61364, -4.40909, 3., -2; undercase upsilon equals the 2 by 2 matrix containing elements 21.0857, -3.15714, 3., -2.

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: The product of Gamme sub zero and y sub t equals the product of Gamme sub 1 and y sub (t-1) plus C plus the product of psi and z sub t plus the product of pi and eta sub t

Model Variable
Description
Dimensions
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: y sub t equals the product of Theta sub 1 and y sub (t-1) plus theta sub c plus the product of theta sub 0 and z sub t plus (theta sub y) times the summation with s from 1 to infinity of (theta sub f superscript (s-1)) times theta sub z times E sub t times z sub (t + s)

Model Variable  
Description  
Dimensions
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: Gamma sub(0) equals the 4 by four matrix containing elements: 1, 0, 0, 0, 0, 1, 0, 0, -1.1, 0, 1, 1., 0, 1, 0, 0; Gamma sub(1) equals the 4 by 4 matrix containing elements: 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, .7, 0, 0; C equals the 4 by 1 matrix containing all zero elements; Psi equals the 4 by 2 element matrix containing elements 0, 0, 0, 0, 4., 1., 3., -2.; Pi equals the 4 by 2 matrix containing elements 1, 0, 0, 1, 0, 0, 0, 0

produces output: Theta sub(c) equals the 4 by 1 matrix containing the elements 0., 0., 0., 0.; Theta sub(0) equals the 4 by 2 matrix containing the elements 1.61364, -4.40909, 3., -2., 3.675, -2.45, 2.1, -1.4; Theta sub(f) equals the 2 by 2 matrix containing elements 0.909091, 1.23405, 0., 0.; Theta sub(y) equals the 4 by 2 matrix containing the elements 1.28794, 1.74832, -2.2204510 to the negative 16, -1.1102210 to the negative 16, 1.41673, 0.702491, -1.1102210 to the negative 16, 1.22066; Theta sub(z) equals the 2 by 2 matrix containing the elements -0.0796719, -2.29972, 2.4577, -1.63846; Theta sub(y) times Theta sub(z) equals the 4 by 2 matrix containing the elements 4.19421, -5.82645, -2.5516810 to the negative 16, 6.9254610 to the negative 16, 1.61364, -4.40909, 3., -2.

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.

the norm of (G1 sub gensys minus G1 sub AMA) is less than 1 times 10 to the -8; the norm of (CC sub gensys minus CC sub AMA) is less than 1 times 10 to the -8; the norm of (impact sub gensys minus impact sub AMA) is less than 1 times 10 to the -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

fmat sub gensys is similar to fmat sub AMA; the norm of (real(ywt sub gensys times fwt sub gensys) minus the product of ywt sub AM and fmat sub AMA and ywt sub AMA) is less than 1 times 10 to the -8


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 G0’s 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 z’y -------
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 G0’s 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 z’y -------
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 G0’s 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 z’y -------
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 .

F equals H sub plus times the inverse of (H sub 0 plus the product of H sub plus and B) equals phi times H sub plus; H equals the 2 by 6 matrix containing the elements g1, 0, g0, negative pi, 0, 0, 0, 0, 0, 0, 0, I; F equals the inverse of the the 2 by 2 matrix containing g0, pi, B sub pi*L, B sub pi*R and the 2 by 2 matrix containing the elements 0, 0, 0, I, equals the matrix containing 0, fU, and fL; F equals the product of the three following matrices: the matrix containing I and V superscript -1, the matrix containing zeros and f sub U and delta sub f, and the matrix containing I and V; F to the k equals the product of the following three matrices: the matrix containing I and V superscript -1, the matrix containing zero f sub U times delta sub f superscript (k-1) and delta sub f superscript k, and the matrix containing I and V; F to the k times phi times psi equals F to the k times psi times the column vector containing elements phi U and phi L; F to the k times phi times psi equals psi times the product of the following four matrices: the matrix containing the elements I and V superscript -1, the matrix containing the elements zero, f sub U times delta sub f superscript (k-1), and delta sub f superscript k, the matrix containing the elements I and V, and the matrix containing the elements phi sub U and phi sub L; F to the k times phi times psi equals F to the k times psi times the column vector containing elements phi sub U and phi sub L; F to the k times phi times psi equals psi times the matrix containing the two elements f sub U times delta sub negative times f superscript (k-1) times V times phi sub L and the inverse of V times delta sub f superscript k times V time phi sub L.

Last update: August 2, 2013