A Reliable and Computationally Efficient Algorithm for Imposing the Saddle Point Property in Dynamic Models

Finance and Economics Discussion Series Divisions of Research & Statistics and Monetary Affairs Federal Reserve Board, Washington, D.C. A Reliable and Computationally Efficient Algorithm for Imposing the Saddle Point Property in Dynamic Models Gary S. Anderson 2010-13 NOTE: Staff working papers in the Finance and Economics Discussion Series (FEDS) are preliminary materials circulated to stimulate discussion and critical comment. The analysis and conclusions set forth are those of the authors and do not indicate concurrence by other members of the research staff or the Board of Governors. References in publications to the Finance and Economics Discussion Series (other than acknowledgement) should be cleared with the author(s) to protect the tentative character of these papers.

A Reliable and Computationally Efficient Algorithm for Imposing the Saddle Point Property in Dynamic Models Gary S. Anderson Abstract This paper describes a set of algorithms for quickly and reliably solving linear rational expectations models. The utility, reliability and speed of these algorithms are a consequence of 1) the algorithm for computing the minimal dimension state space transition matrix for models with arbitrary numbers of lags or leads, 2) the availability of a simple modeling language for characterizing a linear model and 3) the use of the QR Decomposition and Arnoldi type eigenspace calculations. The paper also presents new formulae for computing and manipulating solutions for arbitrary exogenous processes. 1 Introduction and Summary Economists at the Board have an operational need for tools that are useful for building, estimating and simulating moderate to large scale rational expectations models. This context dictates a need for careful attention to computational efficiency and numerical stability of the algorithms. These algorithms have proved very durable and useful for staff at the central bank.ManyeconomistsattheFederalReserveBoardhaveusedthealgorithms in their daily work and their research.2 With the exception of researchers at other European central banks(Zagaglia, 2002), few economists outside the US 1 First, I would like to thank George Moore, my now deceased mentor, friend and coauthor of(Anderson & Moore, 1985). I also wish to thank Brian Madigan, Robert Tetlow, Andrew Levin, Jeff Fuhrer and Hoyt Bleakley for helpful comments. I am responsible for any remaining errors. The views expressed herein are mine and do notnecessarilyrepresenttheviewsoftheBoardofGovernorsoftheFederalReserve System. 2 See,forexample,(Bomfim,1996;Fuhrer&Moore,1995a;Fuhrer&Moore,1995b; Fuhrer & Moore, 1995; Fuhrer & Madigan, 1997; Fuhrer, 1997b; Fuhrer, 1997a; Preprint submitted to Elsevier 28 August 2009

central bank seem to know about the method. This paper attempts to make the method and approach more widely available by describing the underlying theory along with a number of improvements on the original algorithm.3 The most distinctive features of this approach are: • its algorithm for computing the minimal dimension state space transition matrix • its use of bi-orthogonality to characterize the asymptotic constraints that guarantee stability (See Section 3.1.2). • It’s reliance on QR Decomposition and the real Schur Decomposition for speed and accuracy. This unique combination of features makes the algorithm especially effective for large models. See (Anderson, 2006) for a systematic comparison of this algorithm with the alternatives procedures. The remainder of the paper is organized as follows. Section 2 states the saddle point problem. Section 3 describes the algorithms for solving the homogeneous and inhomogeneous versions of the problem and describes several implementations. Section 4 shows how to compute matrices often found useful for manipulating rational expectations model solutions: the observable structure(Section 4.1) and stochastic transition matrices(Section 4.2). The Appendices contain proofs for the linear algebra underlying the algorithm and the solution of a simple example model. 2 Saddle Point Problem Statement Consider linear models of the form: θ X H x = Ψz , t = 0,...,∞ (1) i t+i t i=−τ Fuhrer, 1996; Fuhrer et al., 1995; Fuhrer & Hooker, 1993; Fuhrer, 1994; Fuhrer, 1997c; Orphanides et al., 1997; Levin et al., 1998; Orphanides, 1998; Orphanides & Wieland, 1998; Edge et al., 2003; Orphanides & Williams, 2002). 3 At the Board, economists commonly refer to this family of algorithms as the AIM algorithm.Ametaphorrelatingourapproachtothe“shootingmethod”inspiredthe name. 2

with initial conditions, if any, given by constraints of the form x = xdata,i = −τ,...,−1 (2) i i where both τ and θ are non-negative, and x is an L dimensional vector of t endogenous variables with lim kx k < ∞ (3) t t→∞ and z is a k dimensional vector of exogenous variables. t Section 3 describes computationally efficient algorithms for determining the existence and uniqueness of solutions to this problem. 3 The Algorithms The uniqueness of solutions to system 1 requires that any transition matrix characterizing the dynamics of the linear system have an appropriate number of explosive and stable eigenvalues(Blanchard & Kahn, 1980), and that a certain set of asymptotic linear constraints are linearly independent of explicit and certain other auxiliary initial conditions(Anderson & Moore, 1985). The solution methodology entails (1) Manipulating the left hand side of equation 1 to obtain a state space transition matrix, A, along with a set of auxiliary initial conditions, Z for the homogeneous solution.       x x x −τ −τ+1 −τ        .   .   .  Z . .  = 0 and  . .  = A . .  (4)                   x x x θ θ θ−1 See Section 3.1.1. (2) Computing the eigenvalues and vectors spanning the left invariant space associated with large eigenvalues. See Section 3.1.2. VA = MV (5) with the eigenvalues of M all greater than one in absolute value. (3) Assembling asymptotic constraints, Q, (See Section 3.1.2.) by combining the: (a) auxiliary initial conditions identified in the computation of the transition matrix and 3

(b) the invariant space vectors   Z Q =   (6)   V (4) Investigating the rank of the linear space spanned by these asymptotic constraints and, when a unique solution exists, (a) Computing the auto-regressive representation, B. See Section 3.1.3. (b) Computing matrices, φ,F,ϑ for characterizing the impact of the inhomogeneous right hand side term. See Section 3.2.1. Figure 1 presents a flowchart of the algorithm. 3.1 Homogeneous Solution Suppose, for now, that Ψ = 0:4 θ X H x = 0,t ≥ 0 (8) i t+i i=−τ lim kx k < ∞ (9) t t→∞ The homogeneous specification 8 is not restrictive. Since the procedure can handleinhomogeneousversionsofequation1byrecastingtheprobleminterms of deviations from a steady state value. However, the next section provides a moreintuitive,flexibleandcomputationallyefficientalternativeforcomputing inhomogeneous solutions.5 4 Note that there is no unique steady state requirement. Steady state solutions, x∗ satisfying θ X (H )x∗ = 0 (7) i i=−τ lie in a linear subspace of RL. We will develop conditions that guarantee solutions that evolve from a given set of initial conditions to a single point in this subspace. Asaresult,onecanapplytheseroutinestomodelswithunitroots,seasonalfactors, cointegrating vectors and error correction terms. 5 The original algorithmic description and software implementation of these algorithmsdevelopedhomogeneoussolutions.Researchersobtainedsolutionsformodels withinhomogeneoussystemsbyaddinganequationoftheformi = i withinitial t t−1 condition i = 1 to the system. t−1 4

’ $ Begin & % ? Compute Unconstrained Autoregressive Representation(H],∗) and Auxiliary Initial Conditions(Z],∗), (See Section 3.1.1) H −−−→ {H],∗,Z],∗} ? Compute Convergence Constraint(V) (See Section 3.1.2) {H],∗,Z],∗} −−−→ V ? Compute Asymptotic Constraints (Q), (See Section 3.1.2) {Z],∗,V} −−−→ Q ? Investigate Rank of Q R (See Section A.3) {Q} −−−→ {B,φ,F,ϑ} ’ ? $ End & % Fig. 1. Algorithm Components 3.1.1 State Space Transition Matrix and Auxiliary Initial Conditions: A,Z Thissectiondescribeshowtodetermineafirst-orderstatespacerepresentation oftheequationsystem8.Themethodisanextensionoftheshufflingalgorithm developed in(Luenberger, 1978; Luenberger, 1977). If H is non-singular, we θ 5

can immediately obtain x in terms of x ...x t+θ t−τ t+θ−1 (cid:20) (cid:21) −H−1 H ... H (10) θ −τ θ−1 However, the natural specification of many economic models has singular H . θ This,first,algorithmappliesfullranklineartransformationstoequationsfrom the original linear system in order to express x in terms of x ...x . It t+θ t−τ t+θ−1 produces an unconstrained, typically explosive, autoregressive representation for the evolution of the components of the state space vectors and a set of vectors that provide important auxiliary initial conditions. Section A.1 presents a proof that repeating this process of annihi- Fig. 2. Initial Tableau lating and regrouping rows ulti- τ +θ H♯,0 mately produces an H],k = H],∗ with H],∗ non-singular. The proof θ identifies a benign rank condition that guarantees that the algorithm will successfully compute the unconstrained autoregression and the auxiliary initial conditions. Fig. 3. Forward Row Annihilation Figures 2-4 provide a graphical H♯,1 characterization of the linear al- Z♯,1 =F♯,1 gebraic transformations characterizing the algorithm. Figure 2 presents a graphical characterization of the relevant set of linear constraints for t = 0...(τ + θ). The figure represents the regions where the coefficients are poten- Fig. 4. Full Rank Leading Block tially non-zero as shaded gray. If H],0 is singular, one can find a lin- θ H♯,k with H♯,k non-singular ear combination of its rows which θ Z♯,k preserves the rank of H],0 but θ which annihilates one or more of F♯,k its rows. Consequently,onecanpre-multiply the matrices presented in Figure 2 by a unitary matrix to get the distribution of zeros displayed in Figure 3. Since the matrices repeat over time, one need only investigate the rank of the square matrix H],0. U H],0 = θ 1 6

(cid:20) (cid:21) F],1U H],0 . With some of the rows of U H],0 all zeros. One can regroup 1 θ 1 θ the rows in the new tableau to get H],1. By construction, the rank of H],1 can θ be no less than the rank of H],0. θ Note that changing the order of the equations in H],0 will not affect the rank θ of H],0 or the space spanned by the nonzero rows that will enter H],1. Since θ θ H],0 is not full rank, a subset of the rows of H],1 will span the same space as θ θ H],0. The other rows in H],1 will come from linear combinations of the original θ θ system of equations in effect “shifted forward” in time. One can think of the “Full Rank Leading Block” matrix as the result of premultiplications of the “Initial Tableau” by a sequence of unitary matrices. Implementations of the algorithm can take advantage of the fact that the rows of the matrices repeat. The regrouping can be done by “shifting equations forward” in time in an L×L(τ +1+θ) version of the tableau. Section A.1 presents a proof that repeating this process of annihilating and regrouping rows ultimately produces an H],k = H],∗ with H],∗ non-singular. θ Algorithm 1 presents pseudo code for an algorithm for computing the components of the state space transition matrix and the auxiliary initial conditions. Algorithm 1 Given H, 1 compute the unconstrained autoregression. 2 funct F (H) ≡ 3 1 k := 0 4 Z0 := ∅ 5 H0 := H 6 Γ = ∅ 7 while Hk is singular ∩rows(Zk) < L(τ +θ) 8 θ do 9 10 Determine a Non-singular matrix that annihilates L−r(Hk) Rows of Hk θ θ   Uk Uk =  Z := rowAnnihilator(Hk) 11   θ Uk N   0 UkHk ... UkHk Hk+1 :=  Z τ Z θ−1 12   UkHk ... UkHk N τ N θ   Zk Zk+1 :=   13   UkHk ... UkHk Z τ Z θ−1 k := k +1 14 od 15 7

(cid:20) (cid:21) 16 Γ = −H−1 H ... H θ −τ θ−1   0 I A =   17   Γ (cid:20) (cid:21) 18 return{ Hkτ ... Hk ,A,Zk} − θ . 19 The algorithm terminates with:   x −τ   (cid:20) (cid:21) .  H]∗ ... H]∗  . .  = 0 (11) −τ θ       x θ with H]∗ non singular. Let θ (cid:20) (cid:21) Γ] = −(H]∗)−1 H]∗ ... H]∗ (12) θ −τ θ−1 Then   x −τ    .  x = Γ] . .  (13) θ       x θ−1 This unconstrained auto-regression in x provides exactly what one needs to t construct the state space transition matrix.   0 I A] =   (14)   Γ] so that     x x −τ+1 −τ      .   .   . .  = A . .  (15)             x x θ θ−1 8

3.1.2 Asymptotic Linear Constraint Matrix: Q In order to compute solutions to equation 8 that converge, one must rule out explosive trajectories. Blanchard and Kahn(Blanchard & Kahn, 1980) used eigenvalue and eigenvector calculations to characterize the space in which the solutions must lie. In contrast, our approach uses an orthogonality constraint to characterize regions which the solutions must avoid. Each left eigenvector associated with a given eigenvalue is orthogonal to each right eigenvector associated with roots associated with different eigenvalues. Since vectors in the left invariant space associated with roots outside the unit circle are orthogonal to right eigenvectors associated with roots inside the unit circle, a given state vector that is part of a convergent trajectory must be orthogonal to each of these left invariant space vectors. See theorem 4 on page 20. Thus, the algorithm can exploit bi-orthogonality and a less burdensome computation of vectors spanning the left invariant space in order to rule out explosive trajectories. If the vectors in V span the invariant space associated with explosive roots, trajectories satisfying equation 8 are non-explosive if and only if VA = MV (16) with the eigenvalues of M similar to the Jordan blocks of A associated with all eigenvalues greater than one in absolute value.   x t−τ    .  V  . .  = 0 (17)       x t+θ−1 for some t.6 Combining V and Z] completely characterizes the space of stable solutions satisfying the linear system 8.   Z] Q =   (18)   V 6 If A has roots with magnitude 1 then trajectories can converge to either a limit cycle or a non-zero fixed point. Otherwise, non-explosive trajectories will converge to the origin. 9

The first set of equations come from the equations in equation system 8 which do not appear in the transformed system of Equation 11 but must nevertheless be satisfied. The second set of equations come from the constraint that the solutions are non-explosive. Algorithm 2 provides pseudo code for computing Q. Algorithm 2 Given A,Z],∗, 1 funct F (A,Z],∗) 2 2 Compute V, the vectors spanning the left 3 invariant space of A associated with eigenvalues 4 greater than one in magnitude 5   Z],∗ Q :=   6   V returnQ 7 . 8 3.1.3 Convergent Autoregressive Representation: B The first two algorithms together produce a matrix Q characterizing constraints guaranteeing that trajectories are not explosive. See theorem 5 and corollary 2 for a proof. (Hallet & McAdam, 1999) describes how to use the matrix Q from Section 3.1.2 to impose saddle point stability in non linear perfect foresight models. However, for linear models with unique saddle point solutions it is often useful to employ an autoregressive representation of the solution. Theorem 6 in Section A.3 provides a fully general characterization of the existence and uniqueness of a saddle point solution. A summary for typical applications of the algorithm follows. Partition Q = (cid:20) (cid:21) Q Q where Q has Lτ columns. When η = Lθ, Q is square. If Q is L R L R R 10

non-singular, the system has a unique solution7   B         x t−τ x t−τ B        .   .     2 . . .    = Q− R 1Q L and solutions are of the form x t = B   . .    , x t+k = B k    . .            x x   t−1 t−1 B θ (19) Algorithm 3 provides pseudo code for computing B. Algorithm 3 Given Q, 1 funct F (Q) 2 3 cnt = noRows(Q) 3   { { Q Q, ,∞ 0} } c c n n t t < > L L θ θ return 4  { { Q B , = ∞ − } Q− R 1Q L ,1} ( o Q th R e s r i w n i g s u e lar) . 5 3.2 Inhomogeneous Solution Now, suppose θ X H x = Ψz ,t ≥ 0 (20) i t+i t i=−τ lim kx k < ∞ (21) t t→∞ 3.2.1 Inhomogeneous Factor Matrices: φ,F Theorem 1 Given structural model matrices, H ,i = −τ,...,θ and i Ψ, convergent autoregression matrix B there exist inhomogeneous 7 If Q is singular, the system has an infinity of solutions. When η < Lθ, The R system has an infinity of solutions. When Q has more than Lθ rows, The system has a unique nontrivial solution only for specific values of x . data 11

factor matrices, φ and F such that with (cid:20) (cid:21) B ... B = B (22) −τ −1   −1 B     (cid:20) (cid:21) .  φ = H + H ...H  . .  (23)  0 1 θ         B θ   −1 (cid:20) (cid:21) ∞ 0 x = X B x + 0 ...0 I X (Fs ) (24) t i t+i   i=−τ s=0 φΨz t+s will satisfy the linear inhomogeneous system equation 20. See Section A.4 for derivations and formulae. Algorithm 4 provides pseudo code for computing φ and F. Algorithm 4 Given H,Q 1 funct F (H,Q) 2 4 Where 3   B B L R    . .  4 B =   . . . .   = Q− R 1Q L     Bθ Bθ L R   B R    .  6 φ = (H 0 +H +   . .   )−1     B Rθ   0 I     . . . ...          0 I           7 F =   0 0 I           . .  . .     .  .   B    −φH +     −φH +     ... −φH +   . R           .   0  I   .                I B B R θ−1 return(φ,F) 8 . 9 12

3.2.2 Exogenous VAR Impact Matrix, ϑ Modelers can augment the homogeneous linear perfect foresight solutions with particular solutions characterizing the impact of exogenous vector autoregressive variables. Theorem 2 When z = Υz (25) t+1 t one can show that   x t−τ   (cid:20) (cid:21) .  x = B B  . . +ϑz (26) t L R   t     x t−1 where   0    . .  (cid:20) (cid:21)  .  vec(ϑ) = 0 ...0 I (I −ΥT ⊗F)−1vec     (27)    0      φΨ See Section A.5 for derivations and formulae. 3.3 Implementations This set of algorithms has been implemented in a wide variety of languages. Three implementations, a Matlab, a “C”, and a Mathematica implementation, are available from the author.Each implementation avoids using the large tableauofFigure2.Theyeachshiftelementsintherowsofasinglecopyofthe matrix H. Each implementation eliminates inessential lags from the autoregressive representation before constructing the state space transition matrix for invariant space calculations. The most widely used version is written in MATLAB. The MATLAB version has a convenient modeling language front end for specifying the model equations and generating the H matrices. i The “C” version, designed especially for solving large scale models is by far the fastest implementation and most frugal with memory. It uses sparse lin- 13

ear algebra routines from SPARSKIT(Saad, 1994) and HARWELL(Numerical Analysis Group, 1995) to economize on memory. It avoids costly eigenvector calculations by computing vectors spanning the left invariant space using ARPACK(Lehoucq et al., 1996). For small models, one can employ a symbolic algebra version of the algorithms written in Mathematica. On present day computers this, code can easily construct symbolic state space transition matrices and compute symbolic expressions for eigenvalues for models with 5 to 10 equations. The code can often obtainsymbolicexpressionsfortheinvariantspacevectorswhenthetransition matrix is of dimension 10 or less. 4 Other Useful Rational Expectations Solution Calculations Economists use linear rational expectations models in a wide array of application. The following sections describe calculations which are useful for optimal policy design, model simulation and estimation exercises. 4.1 Observable Structure: S To compute the error term for estimation of the coefficients of these models, one must commit to a particular information set. Two typical choices are t and t-1 period expectations. Given structural model matrices, H ,i = −τ,...,θ and convergent autorei gression matrices B ,i = −τ,−1 there exists an observable structure matrix, i S   x t−τ    .  (cid:15) = S . .  (28) t       x t See Section A.6 for a derivation and formula for S. Algorithm 5 provides pseudo code for computing S for a given lag, k∗ in the availability of information. Algorithm 5 14

Given B,k∗ 1 funct F (B,k∗) 2 5   0 I     B ˜ =   3       B (cid:20) (cid:21) 4 S = 0 L×Lmax(0,k∗−1) H −τ ... H 0 +     B     (cid:20) (cid:21) .   5   H 1 ...H θ   . .   B ˜k∗ 0 L×Lmax(0,k∗−1)           B θ return(S) 6 . 7 4.2 Stochastic Transition Matrices: A,B To compute covariances, practitioners will find it useful to construct the stochastic transition matrices A and B. Given structural model matrices, H ,i = −τ,...,θ and convergent autoregresi sion matrices B ,i = −τ,−1 there exist stochastic transition matrices B, A i such that     x x t−τ+1 t−τ      .   .  (cid:20) (cid:21)  . .  = A . . +B (cid:15) +Ψ(E[z |I ]−E[z |I ]) (29)     t t t t t−1         x x t t−1 SeeSectionA.7forderivationsandformulaeforAandB.Algorithm6provides pseudo code for computing A and B Algorithm 6 Given H,Ψ,S 1 funct F (H,S) 2 6 15

  0 I     . . . ...     3 A =      0 I      S−1S ... ... S−1S 0 t−τ−max(k∗−1,0)+1 0 −1   0    . .   .    4 B =      0      S−1 0 return(A,B) 5 . 6 5 Conclusions This paper describes a set of algorithms that have proved very durable and useful for staff at the central bank. The most distinctive features of this approach are: • its algorithm for computing the minimal dimension state space transition matrix • its use of bi-orthogonality to characterize the asymptotic constraints that guarantee stability. • It’s reliance on QR Decomposition and the real Schur Decomposition for speed and accuracy. The unique combination of features makes the algorithm more efficient than all the alternatives—especially for large models. Staff at the Federal Reserves have developed a large scale linear rational expectations model consisting of 421 equations with one lead and one lag. This model provides an extreme example of the speed advantage of the Anderson-Moore Algorithm(AMA). On an Intel(R) Xeon 2.80GHz CPU running Linux the MATLAB version of AMA computes the rational expectations solution in 21 seconds while the the MATLAB version of gensys, a popular alternative procedure, requires 16,953 seconds. See (Anderson, 2006) for a systematic comparison of this algorithm with the alternative procedures. The code is available for download at http://www.federalreserve.gov/Pubs/oss/oss4/aimindex.html 16

References Anderson, Gary. 2006. Solving Linear Rational Expections Models: A Horse Race. Tech. rept. Federal Reserve System, Finance and Economics Discussion Series. Anderson, Gary, & Moore, George. 1985. A Linear Algebraic Procedure For Solving Linear Perfect Foresight Models. Economics Letters, 17(3). Blanchard, Olivier Jean, & Kahn, Charles M. 1980. The Solution of Linear Difference Models under Rational Expectations. Econometrica, 48(5), 1305–11. Bomfim,Antulio.1996. Forecasting the Forecasts of Others: Expectational Heterogeneity and Aggregate Dynamics. Tech. rept. 1996-41. Federal Reserve Board, Finance and Economics Discussion Series. Edge, Rochelle M., Laubach, Thomas, & Williams, John C. 2003. The responses of wages and prices to technology shocks. Tech. rept. 2003-21. Federal Reserve Bank of San Francisco, Working Papers in Applied Economic Theory. Fuhrer, Jeffrey C. 1994. Optimal Monetary Policy and the Sacrifice Ratio. Pages 43–69 of: Fuhrer, Jeffrey C. (ed), ‘Goals, Guidelines, and Constraints Facing Monetary Policymakers’. Federal Reserve Bank of Boston Conference Series No. 38. Fuhrer, Jeffrey C. 1996. Monetary Policy Shifts and Long-term Interest Rates. Quarterly Journal of Economics, 111(November), 1183–1209. Fuhrer, Jeffrey C. 1997a. Inflation/Output Variance Trade-offs and Optimal MonetaryPolicy. Journal of Money Credit and Banking,29, No. 2(May), 214–234. Fuhrer, Jeffrey C. 1997b. The (Un)Importance of Forward-Looking Behavior in Price Specifications. Journal of Money Credit and Banking, 29, No. 3(August), 338–350. Fuhrer, Jeffrey C. 1997c. Towards a Compact, Empirically-Verified Rational Expectations Model for Monetary Policy Analysis. Carnegie-Rochester Conference Series on Public Policy, forthcoming. Fuhrer, Jeffrey C., & Hooker, Mark W. 1993. Learning About Monetary Regime Shifts in an Overlapping Wage Contract Economy. Journal of Economic Dynamics and Control, 17, 531–553. Fuhrer, Jeffrey C., & Madigan, Brian. 1997. Monetary Policy When Interest Rates are Bounded at Zero. Review of Economics and Statistics, 79, No. 4(November), 573–585. Fuhrer, Jeffrey C., & Moore, George R. 1995. Forward-Looking Behavior and the Stability of a Conventional Monetary Policy Rule. Journal of Money Credit and Banking, 27, No. 4(November), 1060–70. Fuhrer,JeffreyC.,&Moore,GeorgeR.1995a. InflationPersistence. Quarterly Journal of Economics, 110(February), 127–159. Fuhrer, Jeffrey C., & Moore, George R. 1995b. Monetary Policy Trade-Offs 17

and the Correlation Between Nominal Interest Rates and Real Output. American Economic Review, 85(March), 219–239. Fuhrer, Jeffrey C., Moore, George, & Schuh, Scott. 1995. A Comparison of GMM and Maximum-Likelihood Estimates of Quadratic Cost Inventory Models. Journal of Monetary Economics, 35(February), 115–157. Hallet, Andrew Hughes, & McAdam, Peter (eds). 1999. Analyses in Macroeconomic Modelling. Kluwer Academic Publishers. Chap. Accelerating Non Linear Perfect Foresight Solution by Exploiting the Steady State Linearization. Lehoucq, R.B., Sorensen, D.C., & Yang, C. 1996. ARPACK USERS GUIDE: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods. Tech. rept. Department of Computational and Applied Mathematics Rice University. Levin, Andrew, Williams, John, & Wieland, Volker. 1998 (January). Are Simple Rules Robust under Model Uncertainty? Seminar Paper. Luenberger, David G. 1977. Dynamic Equations in Descriptor Form. IEEE Transactions on Automatic Control, 312–321. Luenberger, David G. 1978. Time-Invariant Descriptor Systems. Automatica, 14, 473–480. Noble, Ben. 1969. Applied Linear Algebra. Prentice-Hall, Inc. Numerical Analysis Group. 1995. The Harwell Sparse Matrix Library Manual. Tech. rept. The Numerical Analysis Group. Orphanides, Athanasios. 1998. Monetary Policy Evaluation with Noisy Information. Tech. rept. 98-50. Finance and Economics Discussion Series. Orphanides, Athanasios, & Wieland, Volker. 1998 (January). Price Stability and Monetary Policy Effectiveness when Nominal Interest Rates are Bounded at Zero. Working Paper. Orphanides, Athanasios, & Williams, John C. 2002. Robust Monetary Policy Rules with Unknown Natural Rates. Brookings Papers on Economic Activity. Orphanides, Athanasios, Small, David, Wilcox, David, & Wieland, Volker. 1997 (May). A Quantitative Exploration of the Opportunistic Approach to Disinflation. Tech. rept. 97-36. Federal Reserve Board, Finance and Economics Discussion Series. Saad, Youcef. 1994. SPARSKIT: A Basic Tool Kit for Sparse Matrix Computations. Tech. rept. Computer Science Department, University of Minnesota. Taylor, John. 1979. Staggered Wage Setting in a Macro Model. The American Economic Review, 69(2), 108–113. Zagaglia, Paolo. 2002. Matlab Implementation of the AIM Algorithm: A Beginner’s Guide. Tech. rept. 169. Universita’ Politecnica delle Marche (I), Dipartimento di Economia. 18

A Proofs A.1 Unconstrained Autoregression Theorem 3 Let   H −τ ... H θ         H −τ ... H θ     H =   ...    τ+θ+1 (A.1)        H ... H  −τ θ There are two cases: • When H is full rank the algorithm terminates with Z]∗ (Z[∗) and non-singular H]∗ (H[∗) θ τ • When H is not full rank the algorithm terminates when some row (cid:20) (cid:21) of Hk ...Hk is zero. −τ θ Proof Consider the case when H is full rank. Each step of the algorithm applies a rank preserving pre-multiplication by a non singular matrix. Each step of the algorithm where H],k is singular, increases the row rank of Z],k θ by at least one. At each step Z],k are full rank. The rank of Z],k cannot exceed L(τ +θ). (cid:4) Thefollowingcorollaryindicatesthatmodelswithuniquesteady-statesalways terminate with non singular H],∗. θ Corollary 1 If (Pθ H ) is non singular then i=−τ i (1) H is full rank. (2) The origin is the unique steady state of equation 1. (3) there exists a sequence of elementary row operations that transforms H into H∗ Proof Suppose H is not full rank. Then there is a non zero vector V such that VH = 0. Consequently,   I ... I   (cid:20) (cid:21) . .  V ...V H. . ... . .  = 0 (A.2) −τ θ       I ... I 19

and θ X V ( H ) = 0∀i (A.3) i j j=−τ So that (Pθ H ) must be singular. (cid:4) i=−τ i A.2 Asymptotic Constraints Theorem 4 Consider a left invariant space and a right invariant space with no eigenvalues in common. Suppose V spans the left 1 invariant space and W spans the right invariant space. 2 V A = T V (A.4) 1 1 1 AW = W T (A.5) 2 2 2 With eigenvalues of T 6= T . Then V W = 0 1 2 1 2 Proof Arighteigenvectorx andaleft-eigenvectory correspondingtodisi j tinct eigenvalues λ and λ are orthogonal.(Noble, 1969) Finite dimensional i j matrices have finite dimensional Jordan blocks. Raising a given matrix to a power produces a matrix with smaller Jordan blocks. Raising the matrix to a high enough power ultimately eliminates all nontrivial Jordan Blocks. Consequently, the left invariant space vectors are linear combination of the left eigenvectors and the right invariant space vectors are a linear combination of the right eigenvectors of the transition matrix raised to some finite power.   y 1    .  V = β  . .  (A.6) 1 1      y J (cid:20) (cid:21) W = x ... x α (A.7) 2 1 K 2   y 1    . (cid:20) (cid:21) V W = β  . .  x ... x α = 0 (A.8) 1 2 1  1 K 2     y J (cid:4) Theorem 5 Let {xconv}, t = −τ,...,∞ be a non explosive solution t satisfying equation 1. Let A be the state space transition matrix for 20

equation 1 and V be a set of invariant space vectors spanning the invariant space associated with roots of A of magnitude bigger than 1. Then for t = 0,...,∞   xconv  t−τ   .  V  . .  = 0 (A.9)       xconv t+θ−1 Proof Using W, the left generalized eigenvectors of A, one can employ the Jordan Canonical Form of A to write WHA = JWH (A.10) so that At = (WH)−1JtWH (A.11) y = Aty (A.12) t 0 WHy = JtWHy (A.13) t 0 lim y = 0 ⇒ lim WHy = 0 (A.14) t t t→∞ t→∞ Consequently, WHy = 0∀i such that |J | > 1. (A.15) i 0 ii so that Vy = αWHy = 0 (A.16) 0 0 (cid:4) Corollary 2 Let {x }, t = −τ,...,∞ be a solution satisfying equat tion 8. If A has no roots with magnitude 1 then the path converges to the origin if and only if   x t−τ    .  V  . .  = 0 (A.17)       x t+θ−1 for some t. Proof WHy = 0∀i such that |J | > 1. (A.18) i 0 ii 21

means y t J 6= 1. (A.19) ii y = Aty (A.20) t 0 (cid:4) A.3 Existence and Uniqueness Theorem 6 Identify Q ,Q from L R   Z] (cid:20) (cid:21) Q =   = Q Q (A.21)   L R V with Q an (η×Lθ) matrix, Q an (η×Lτ) matrix, where η represent R L the number of rows in the matrix Q. The existence of convergent solutions depends on the magnitudes of the ranks of the augmented matrix (cid:18)(cid:20) (cid:21)(cid:19) r = rank Q −Q x (A.22) 1 R L data and r = rank(Q ). (A.23) 2 R By construction, r ≥ r and r ≤ Lθ. There are three cases. 1 2 2 (1) If r > r there is no nontrivial convergent solution 1 2 (2) If r = r = Lθ there is a unique convergent solution 1 2 (3) If r = r < Lθ the system has an infinity of convergent solutions 1 2 Corollary 3 When η = Lθ, Q is square. If Q is non-singular, the R R system has a unique solution   B     B   2 = Q−1Q (A.24) . . R L .       B θ 22

and solutions are of the form       x x x t−τ t−τ t−τ        .   .   .  x = B . . , x = B . . , ... x = B . .  (A.25) t   t+1   t+θ     2   θ         x x x t−1 t−1 t−1 If Q is singular, the system has an infinity of solutions. R When η < Lθ, The system has an infinity of solutions. When Q has more than Lθ rows, The system has a unique nontrivial solution only for specific values of x data Proof of rank of Q Proof Thetheoremapplieswellknownresultsonexistenceanduniqueness   x data of solutions to linear equation systems(Noble, 1969). If M =   does 2   0   I 0 not lie in the column space of M =  , then there is no solution. If 1   Q Q L R M lies in the column space of M and the latter matrix is full rank, then 2 1 there is a unique solution. If M lies in the column space of M and the 2 1 latter matrix is not full rank, then there are multiple solutions. (cid:4) A.4 Proof of Theorem 1 Proof Construct the Lτ ×Lτ matrix:   0 I B ˜ =   (A.26)   B Define ˜ B = B B (A.27) k+1 k 23

Applying equation 8 to the unique convergent solution, it follows that   I     (cid:20) (cid:21) B    H − H 0 H +   . .   = 0 (A.28)  .      B θ+1 where (cid:20) (cid:21) (cid:20) (cid:21) H = H ...H H = H ...H (A.29) − −τ −1 + 1 θ Which can also be written as:   I       0 I    (cid:20) (cid:21)      H H H  B   = 0 (A.30) − 0 +     .  B ˜   .    .         B θ So that:   B   (cid:20) (cid:21)  .  H +( 0 H +H  . . )B ˜ = 0 (A.31) − 0 +      B θ   B B L R   (cid:20) (cid:21)  . .  H +( 0 H +H  . . . . )B ˜ = 0 (A.32) − 0 +      B B Lθ Rθ where the Bk matrices are L×L. R Note that       B B 0 B B L R L R        . .  . .   .   . . . . B ˜ = . . . . + . . B (A.33)                   B B 0 B B Lθ Rθ Lθ Rθ 24

Now     0 B B L R     . .   .  H +H . . . . +(H +H  . . )B = 0 (A.34) − +  0 +          0 B B Lθ Rθ Define   B R    .  φ = (H +H  . . )−1 (A.35) 0 +      B Rθ So that   0 B L   . .  φH +φH . . . . +B = 0 (A.36) − +      0 B Lθ Consider the impact that the time t + s value z has on the value of t+s x We can write t+s   I 0     ... . . .           I 0  x t+s−τ    (cid:20) (cid:21)  .  H H H 0 ... 0 I  . .  = Ψz (A.37) − 0 +    t+s        B B  x  L R  t+s  . .    0 . . . .       Bθ Bθ L R or equivalently, φΨz = t+s (cid:20) (cid:21) (cid:20) (cid:21) φ( H 0 + 0 H + − 0        0 B B x (A.38) L R t+s−τ         . .   .   .  H . . . .  H  . . ) . .   +  +                  0 B B x Lθ Rθ t+s 25

or φΨz = t+s (cid:20) (cid:21) φ( H 0 + −       0 B x L t+s−τ        . .    .  H . . . .  0) . . +  +                 0 B Lθ x t+s (A.39) (cid:20) (cid:21) φ( 0 H + 0      B x R t+s−τ        .   .  0 H  . . ) . .   +              B x Rθ t+s Which by equations A.35 and A.36 can be written φΨz = t+s   x t+s−τ   (cid:20) (cid:21) .  −B 0  . . +       x t+s (A.40)   x t+s−τ   (cid:20) (cid:21) .  0 I  . .        x t+s So we have   x t+s−τ   (cid:20) (cid:21) .  −B I  . .  = φΨz (A.41)   t+s     x t+s   x t+s−τ    .  x = B . . +φΨz (A.42) t+s   t+s     x t+s−1 26

Now consider the impact of z on x . We can write t+s t+s−1 (cid:20) (cid:21) (cid:20) (cid:21) φ( H 0 + 0 H + − 0        0 B B x L R t+s−τ−1         . .   .   .  H . . . .  H  . . ) . . +  +  +                  0 B B x Lθ θ t+s−1   I      B   R  φH + . φΨz t+s = 0  .   .      B θ−1 where the last term captures the impact z has on values of x t+s and t+s later. Using equations A.35 and A.36 we can write     I x   (cid:20) −B I (cid:21)      t+s . . . −τ−1      +φH +       B . . . R       φΨz t+s = 0 (A.43)  x    t+s−1   B θ−1     I x   x t+s−1 = B      t+s . . . −τ−1      −φH +       B . . . R       φΨz t+s (A.44)  x    t+s−2   B θ−1 and more generally     I x   x t+s−i = B      t+s . . . −τ−1      +(−φH +       B . . . R       )iφΨz t+s (A.45)  x    t+s−2   B Rθ−1 To accommodate lagged expectations, suppose that information on all the endogenous variables becomes available with the same lag (D∗) in time: ∃K∗ such that x ∈ I ,∀K ≥ K∗ t−k t 27

then set,       E[x |I ] B xdata  t+1 t     t−τ+1−k∗  .   .   .   . .  =  . . B ˜k∗ . . + (A.46)                   E[x |I ] B xdata t+θ t θ t−k∗     I           B      P∞ s=0 ((−φH +    . . R    )sφΨz t+s+1 )     .            B   Rθ−1   .    . .   (A.47)        I              B      P∞ s=0 ((−φH +    . . R    )sφΨz t+s+θ )      .           B Rθ−1 So that x = t   x t−τ   (cid:20) (cid:21) .  B B  . . + L R       x t−1   (cid:20) (cid:21) ∞ 0 0 ...0 I X (Fs )   s=0 φΨz t+s Where   B B L R    . .  B =  . . . .  = Q−1Q (A.48)   R L     Bθ Bθ L R 28

  B R    .  φ = (H +H  . . )−1 (A.49) 0 +      B Rθ   0 I     . . . ...          0 I           F =  0 0 I  (A.50)           . .  . .     .  .   B    −φH +     −φH +     ... −φH +   . R           .   0  I   .                I B B R Rθ−1     I x   x t = B      t . . . −τ      + s X ∞ =0 ((−φH +       B . . . R       )sφΨz t+s ) (A.51)  x    t−1   B Rθ−1 (cid:4) 29

A.5 Proof of Theorem 2 Proof     ∞ 0 ∞ 0 X (Fs ) = X (Fs Υs)z = ϑz (A.52)     t t s=0 φΨz t+s s=0 φΨ where   0    . .  (cid:20) (cid:21)  .  vec(ϑ) = 0 ...0 I (I −ΥT ⊗F)−1vec     (A.53)    0      φΨ (cid:4) A.6 Observable Structure Since one can write     xdata E[x |I ]  t−τ   t+1 t  (cid:20) (cid:21) .  (cid:20) (cid:21) .  (cid:15) = H ...H  . . + H ...H  . . −Ψz (A.54) t −τ 0   1 θ   t         xdata E[x |I ] t t+θ t We find that   xdata  t−τ+1−max(1,k∗)  .  (cid:15) = S . . −Ψz (A.55) t   t     xdata t 30

where (cid:20) (cid:21) S = 0 H ... H + (A.56) L×Lmax(0,k∗−1) −τ 0     B     (cid:20) (cid:21) .    H ...H  . . B ˜k∗ 0  (A.57)  1 θ   L×Lmax(0,k∗−1)         B θ Note that for k∗ ≥ 1 ∂(cid:15) t = H (A.58) ∂xdata 0 t and for k∗ = 0   B   ∂(cid:15) (cid:20) (cid:21) .  t = H + H ...H  . .  = φ−1 (A.59) ∂xdata 0 1 θ   t     B θ A.7 Stochastic Transition Matrices One can write     x x  t−τ−max(k∗−1,0)+1  t−τ−max(k∗−1,0)  .   .   . .  = A . . +B(cid:15) (A.60)     t         x x t t−1 where   0 I     . . . ...     A =   (A.61)    0 I      S−1S ... ... S−1S 0 t−τ−max(k∗−1,0)+1 0 −1   0    . .   .    B =   (A.62)    0      S−1 0 31

B An Example This section applies algorithms 1-3 to the following model with N=2.(Taylor, 1979) 1 N−1 X w = E [ W ]−α(u −u )+ν (B.1) t t t+i t n t N i=0 1 N−1 X W = w (B.2) t t−i N i=0 u = ϑu +γW +µ+(cid:15) (B.3) t t−1 t t E [ν ] = E [(cid:15) ] = 0∀i ≥ 0 (B.4) t t+i t t+i TheinitialmatrixH],0,fortheexamplemodelwithvariableorder{(cid:15),ν,u,w,W} is   0 0 0 0 0 0 −1 α 1 −1 0 0 0 0 −1  2 2   0 0 0 −1 0 0 0 0 −1 1 0 0 0 0 0   2 2    H],0 =  0 0 −ϑ 0 0 −1 0 1 0 −γ 0 0 0 0 0       0 0 0 0 0 1 0 0 0 0 0 0 0 0 0        0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 We would like to express the time t+1 variables in terms of the time t and t−1 variables. If the sub-matrix corresponding to the t+1 variables were non 32

singular we could immediately write:   (cid:15) t−1     ν   t−1          u t−1   (cid:15) t+1     w     t−1 ν     t+1      (cid:20) (cid:21)W   u   = (H],0)−1 H],0 H],0   t−1   t+1 θ −1 0      (cid:15) t  w     t+1      ν     t  W   t+1    u   t     w   t      W t Since (H],0) is singular, we use equations dated subsequent to the present time θ period to construct a set of linear constraints where the leading block is non singular.   0 0 0 0 0 0 0 0 −1 0 0 0 0 −1 1  2 2    0 0 0 0 0 0 0 −ϑ 0 0 −1 0 1 0 −γ     H],? =  0 0 0 0 0 0 0 0 0 0 1 0 0 0 0       0 0 0 0 0 0 0 0 0 0 0 1 0 0 0        0 0 0 0 0 0 −1 α 1 −1 0 0 0 0 −1 2 2 So that 33

  0 0 0 0 0 1 0 0 0 0     0 0 0 0 0 0 1 0 0 0        0 0 0 0 0 0 0 1 0 0      0 0 0 0 0 0 0 0 1 0        0 0 0 0 0 0 0 0 0 1    A =   (B.5)   0 0 0 0 0 0 0 0 0 0      0 0 0 0 0 0 0 0 0 0        0 0 0 0 0 0 −2γ ρ 2γ −γ     0 0 0 0 0 0 −4 4α 3 −2       0 0 0 0 0 0 −2 2α 2 −1 where ρ = θ+2γα (B.6) For the example model   0 0 0 −1 0 0 0 0 −1 1  2 2    0 0 −ϑ 0 0 −1 0 1 0 −γ Z] =       0 0 0 0 0 1 0 0 0 0      0 0 0 0 0 0 1 0 0 0   0 0 0 −1 0 0 0 0 −1 1  2 2    0 0 −ϑ 0 0 −1 0 1 0 −γ       Q = 0 0 0 0 0 1 0 0 0 0      0 0 0 0 0 0 1 0 0 0        0 0 φ 0 0 0 2 φ φ 1 4 5 6 34

where the φ’s are φ = 1+γφ +2φ 7 5 6 −2 (cid:16) 18αγ +(−1+ρ)2 +(2+ρ) φ 2 1 3 +φ 2 2 3 (cid:17) φ = 6 φ 3 −1−3αγ+3ρ+21αγρ−3ρ2+ρ3+φ1+ 108α2γ2+(−1+ρ)4+6αγ(− 1 1+ρ)(1+5ρ)+2(−1+ρ)φ1+ ( 12αγ+(−1+ρ)2) φ23 1 φ5=− γφ3 φ23 (cid:16) (cid:17) 3αγ 216α2γ2 −4(−1+ρ)3 +9αγ (1+(−14+ρ) ρ) +φ 2 1 φ = − 4 2 γφ 23 φ 3 φ 3 = 18αγ +(−1+ρ)2 +(5+ρ) φ 2 1 3 +φ 2 2 3 φ = (−1+ρ)3 +27αγ (1+ρ)+3φ 2 1 r (cid:16) (cid:16) (cid:17)(cid:17) φ = − αγ 216α2γ2 −4(−1+ρ)3 +9αγ (1−14ρ+ρ2) 1   0 0 0 0 0     0 0 0 0 0   1     B = 0 0 ϑ−γφ +2ϑφ γφ 0 φ  4 6 6  7   0 0 −2(φ +ϑφ ) −1−γφ 0  4 5 5      0 0 −φ −ϑφ φ 0 4 5 6 35