Solving Linear DSGE Models with Bernoulli Iterations

Meyer-Gohde, Alexander Article — Published Version Solving Linear DSGE Models with Bernoulli Iterations Computational Economics Provided in Cooperation with: Springer Nature Suggested Citation: Meyer-Gohde, Alexander (2024) : Solving Linear DSGE Models with Bernoulli Iterations, Computational Economics, ISSN 1572-9974, Springer US, New York, NY, Vol. 66, Iss. 1, pp. 593-643, https://doi.org/10.1007/s10614-024-10708-z This Version is available at: https://hdl.handle.net/10419/330775 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Vol.:(0123456789) Computational Economics (2025) 66:593–643 https://doi.org/10.1007/s10614-024-10708-z Solving Linear DSGE Models withBernoulli Iterations AlexanderMeyer‑Gohde1 Accepted: 29 August 2024 / Published online: 26 September 2024 © The Author(s) 2024 Abstract This paper presents and compares Bernoulli iterative approaches for solving linear DSGE models. The methods are compared using 99 different models from the macroeconomic model data base (MMB) and different parameterizations of the monetary policy rule in the medium-scale new Keynesian model of Smets and Wouters (Am Econ Rev 97(3):586–606, 2007. https:// doi. org/ 10. 1257/ aer. 97.3. 586) iteratively. I find that Bernoulli methods providing similar accuracy as measured by the forward error of the solution at a somewhat higher computation burden to the standard method of Dynare when solving DSGE models. The method, however, has convergence properties useful when a specific solution, e.g., unique stable, is sought and can be combined with other iterative methods, such as the Newton method, lending themselves especially to refining solutions—either when standard methods fail or when one moves through a parameter space iteratively—as I show in applications of the methods. Keywords Functional iteration· Numerical accuracy· DSGE· Solution methods JEL Classification C61· C63· E17 1 Introduction Solving linear DSGE models can be reduced to solving a system of linear equations and a matrix quadratic equation with the latter being the more challenging computationally. Standard existing methods predominantly rely on a generalized Schur or QZ decomposition (Moler & Stewart, 1973; Golub & van Loan, 2013) to solve this matrix quadratic equation. Alternative methods from the applied mathematics literature to solve this matrix quadratic equation have yet to be systematically studied in a DSGE context. This paper fills part of that gap, using * Alexander Meyer-Gohde meyer[email protected] 1 Goethe-Universität Frankfurt andInstitute forMonetary andFinancial Stability (IMFS), Theodor-W.-Adorno-Platz 3, 60629FrankfurtamMain, Germany 594 A.Meyer-Gohde Bernoulli-based solution methods from the mathematics literature, introducing line-search and hybrid Newton based methods for matrix quadratic problems and applying them to the solution of linear DSGE models. The ability of Bernoulli methods to be formulated such that they converge to the solution with the smallest eigenvalues in magnitude (ideally the stable solution sought by the researcher) allow them to reliably compute the stable solution, generally with accuracy on the same order of magnitude but higher computational costs than QZ-based methods. When combined with Newton methods and their asymptotic quadratic rate of convergence, the iterative nature of the algorithms also enables the Bernoulli based methods I introduce to correct inaccurate solutions and to generate solutions more quickly than QZ-based methods when iteratively exploring a parameter space. Bernoulli’s method is originally a method for computing the largest root of a scalar polynomial1 and Bernoulli-based functional iterations, familiar to economists in such root-finding settings—see, e.g., Judd (1998), are an alternative to QZ-based methods. The Bernoulli algorithm for matrix polynomials, see Dennis etal. (1978) and explicitly formulated for minimal solvent to a matrix quadratic in Higham and Kim (2000), was introduced by Binder and Pesaran (1999) to linear DSGE models and Rendahl (2017) extended the theoretical convergence analysis of Higham and Kim (2000) to potentially singular (or non monic matrix polynomials) models that are pervasive in the DSGE literature. Yet this algorithm and recent extensions from the applied mathematics literature such as Bai and Gao (2007) have not yet been systematically studied for their numerical properties. This paper does precisely this and introduces a number of algorithms for solving matrix quadratic problems that incorporate insights from other studies and algorithms, including exact line searches from Higham and Kim (2001) and the solution of DSGE models using Newton methods from Meyer-Gohde and Saecker (2024). To build intuition I begin Sect.3 with the scalar case of the algorithm to build intuition before introducing the multidimensional algorithms. I show these extensions improve global convergence by making the Bernoulli method faster and more reliable, improving on the slow convergence (i.e., many iterations) of the baseline Bernoulli algorithm. Finally, I assess the accuracy of the different methods using the practical forward error bounds of Meyer-Gohde (2023), derived from a numerical stability approach incorporating conditioning number and backward error bounds. In this paper, I present eleven different Bernoulli-based solution algorithms using a unified notation and for the application to solving linear DSGE models as an alternative to QZ-based methods. I engage in a number of experiments to compare the algorithms to QZ-based methods2 and Meyer-Gohde and Saecker’s (2024) Newton algorithms, following exactly their experiments to ensure comparability. First I apply the different methods to the models in the 1 It can be cast as the power method, see Golub and van Loan (2013,pp.366–367) operating on the polynomial’s companion matrix. 2 I use Dynare’s (Adjemian etal., 2011) implementation of the QZ method, documented in Villemot (2011), for comparison. Varying implementations of the QZ or generalized Schur decomposition to solve linear DSGE models can be found in Uhlig (1999), Sims (2001), Klein (2000), Al-Sadoon (2018). 595 Solving Linear DSGE Models withBernoulli Iterations Macroeconomic Model Data Base (MMB) (see Wieland et al., 2012, 2016), comparing the performance to the QZ-based method of Dynare and the Newton method both unconditionally (i.e., replacing the QZ method) and then as a refinement (i.e., initializing the iterative methods with the solution generated from QZ). I find that the baseline Bernoulli method always converges to the unique stable solution, but provides a solution of the same order of magnitude of accuracy albeit at an order of magnitude higher computational cost than Newton methods. The different extended methods generally perform comparably—the exception being Bai and Gao’s (2007) modified Bernoulli iteration that is several orders of magnitude more slow and frequently faces convergence issues—trading the convergence of the baseline algorithm off against more accuracy and/or less computational costs. The algorithms are iterative in nature, enabling them to refine the solutions provided by the another method. Using the QZ solution from Dynare to initialize, the methods improve the accuracy of the solution at an additional computational cost generally between 0.1 and 2 times the original cost of Dynare’s QZ, with convergence of all of the algorithms to the unique stable solvent for all of the models in the MMB, except Bai and Gao’s (2007) modified Bernoulli and columnwise Newton–Bernoulli combinations. This iterative nature also lends itself to iterative parameter experiments or estimations and I compare the algorithms with Dynare’s QZ method and Newton algorithm in solving for different parameterizations of the monetary policy rule in the celebrated Smets and Wouters (2007) model of the US economy. Filling in a grid with different values of the reaction of the nominal interest rate rule to inflation and real activity, whereas Dynare’s QZ method starts anew at each parameterization, iterative methods can use the solution from the previous, nearby parameterization to initialize the algorithm. As the density of the grid increases, all of the methods eventually surpass QZ by roughly an order of magnitude in terms of computation cost. Finally, I show that when Dynare’s QZ provides a numerically unstable solution with high forward errors that leads the predicted moments of the model’s variables to be wrong in all digits, the Bernoulli method initialized at this solution provides a refined solution with improved predicted moments and forward errors. However, the Bernoulli algorithms are far from being uniform improvements over the existing standard in the literature and, in particular, outside of situations where an informative initialization is available, the standard QZ is to be preferred. The paper is organized as follows: Sect.2 lays out the general DSGE model and the unknown solution. In Sect.3, I begin by presenting the Bernoulli method for a scalar quadratic equation to build intuition before I then present the set of different Bernoulli methods from the applied mathematics literature in a unified notation as they apply to the class of DSGE models. Section4 examines practical and theoretical considerations such as the choice of initial value, accuracy, and convergence. In Sect.5, I compare the different Bernoulli methods to the standard QZ method and Newton method of Meyer-Gohde and Saecker (2024) in two applications, one using the MMB of 99 different models and the second over a range of parameterizations within the Smets and Wouters (2007) model. Finally, Sect.6 concludes. 596 A.Meyer-Gohde 2 Problem Statement Standard numerical DSGE solution packages available to economists and policy makers—e.g., Dynare (Adjemian etal., 2011), Gensys (Sims, 2001), (Perturbation) AIM (Anderson & Moore, 1985; Anderson et al., 2006), Uhlig’s Toolkit (Uhlig, 1999) and Solab (Klein, 2000)— analyze models that, when (log)linearized, can be expressed in the form of the system of ny linear expectational difference equations where A, B, and C are ny×ny real valued matrices, ℝn y →ℝny , that operate on yt∈ℝny the vector of ny endogenous variables; D is an ny×ne real valued matrix, ℝn e →ℝny , that operate on 𝜀t∈ℝne the vector of ne exogenous, serially uncorrelated shocks with a known distribution, assumed mean zero3; where ny and ne are positive integers ( ny,ne∈ℕ ). The solution to (1) is sought as the unknown linear solution in the form a recursive solution in the time domain–solutions that posit yt as a function of its own past, yt−1 , and exogenous innovations, 𝜀t , which reduces the problem to finding two real valued matrices, P and Q, ny×ny and ny×ne respectively.4 Inserting (2) into (1) and taking expectations ( Et[ 𝜀 t+1] = 0 ), yields the restrictions Generally, a unique P with eigenvalues inside the closed unit circle is sought. Lan and Meyer-Gohde (2014) prove the latter can be uniquely solved for Q if such a P can be found. Hence, the hurdle is the former, matrix quadratic equation. Most linear DSGE methods use a generalized Schur or QZ decomposition (Moler & Stewart, 1973; Golub & van Loan, 2013) of the companion linearization of (1) in some form or another. I will take a different route and instead now solve for P in (3) using Bernoulli iterations. (1) 0 =AE t[ y t+1] +By t +Cy t−1 +D𝜀 t (2) yt=Py t−1+Q𝜀t, ℝ n y +n e→ℝ ny (3) 0=AP2+BP +C,0=(AP +B)Q+D 3 This assumption follows Dynare for expediency at to facilitate the comparison and implementation in the applications of Sect.5. Any pattern of serial correlation than can be captured by a finite order ARMA process can brought into the above form by expanding yt appropriately. Meyer-Gohde and Neuhoff (2015) show that the matrix quadratic and its solution P is unchanged from the presentation here by this assumption—merely the mapping from exogenous shocks to endogenous variables becomes more involved. 4 As above this is a brief statement of the linear problem, Binder and Pesaran (1997), Uhlig (1999), Klein (2000), Sims (2001), and Anderson (2010) provide overviews of the multivariate approaches that have become standard and, more recently, Al-Sadoon (2018) and Al-Sadoon (2020) provide a much more rigorous approach that go beyond the problem of solving for two matrices that I address here. 597 Solving Linear DSGE Models withBernoulli Iterations 3 Bernoulli Iterations forLinear DSGE Models I will begin by analyzing a univariate equation, see, e.g., Higham (2002,p.508) for the functional iteration to find the unstable solution and Judd (1992, pp. 152–153) for functional iterations in general fix point problems for economists, to fix ideas and illustrate how Bernoulli functional iterations can be used to solve quadratic equations. The problem generated by (3) is a (matrix) quadratic problem. To fix ideas, consider its univariate equivalent where I consider (in accordance with the DSGE model), a, b, and c ∈ℝ1 . A functional iteration will reformulate this as giving an iterative procedure to generate a solution From the quadratic equation above, there are a number of possibilities, and so forth. I will focus on the first f(x)=− c ax+b as it is the univariate counterpart of the multivariate Bernoulli algorithm that will be introduced subsequently. Hence Characterize the two solutions via where the solutions are—in analogy to the generalized eigenvalue formulation with the generalized Schur of Klein (2000) (4) 0=ax2+bx +c (5) x=f(x) (6) x j+1 =f ( x j) , with some x 0 (7) f(x)=− c ax +b (8) f(x)=− b+ c x a (9) f(x)=− ax 2 + c b (10) x k+1=f ( xk ) =− c axk+b ,k=0, 1, 2..., x0 given (11) ( s1x−t1 )( s2x−t2 ) =s1s2 ⏟⏟⏟ a x 2 − ( s1t2+s2t1 ) ⏟⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏟ +b x+t1t2 ⏟⏟⏟ c (12) x (i)= {t i s i , if si≠0; �, if si=0; ℝ, if si=ti=0; i=1, 2 } 598 A.Meyer-Gohde Following, e.g., Judd (1998,pp.165–166), local convergence requires |f�(x)|<1 for an x that solves x=f(x) and as for, say, x( 1 ) = t 1 s 1 and hence the iteration is convergent if or there is local convergence to the minimal (smaller in modulus) root of the quadratic when the iteration is formulated via f(x)=− c ax+b . This highlights a significant advantage over Newton-based methods, namely that the functional iteration can be tailored to deliver convergence of the algorithm to a particular root: for example, f(x)=− c ax+b to the minimal (as above) and f(x)=− b+ c x a to the dominant (Higham & Kim, 2000) solution. As a particular solution or a solution with particular properties (namely the minimal solution in a saddle point stable problem) is sought in DSGE models, the Bernoulli iteration x k+1=f ( xk ) =− c axk+b is potentially very useful in solving DSGE models. Turning now to the matrix problem, I will formalize the matrix quadratic equation in (3). For A, B, and C ∈ℝn y ×ny , a matrix quadratic M(P)∶ ℂ n y ×n y→ℂ n y ×ny is defined as with its solutions, called solvents, given by P∈ℂn y ×ny if and only if M(P)=0 . The eigenvalues of the solvent, called latent roots of the associated lambda matrix,5 M(𝜆)∶ℂ → ℂn×n (here of degree two), are given via The latent roots are (i) values of 𝜆∈ℂ such that det M(𝜆)=0 and (ii) ny−rank(A) infinite roots. An explicit link between the quadratic matrix problem and the quadratic eigenvalue problem is given via (13) f�(x)= ac (ax +b)2= s 1 s 2 t 1 t 2 ( s 1 s 2 x k − ( s 1 t 2 +s 2 t 1)) 2 (14) f� ( x(1) ) = s 1 s 2 t 1 t 2 ( s1s2 t1 s1 − ( s1t2+s2t1 )) 2= s 1 s 2 t 1 t 2 ( s1t2 ) 2= s 2 t 1 s1t2 =x (1) x(2) (15) | | | f� ( x(1) )| | | = | | | | x (1) x (2) | | | | < 1 (16) M(P)≡AP 2+BP+C (17) M(𝜆) ≡ A𝜆2+B𝜆+C (18) 𝜆 ∈ℂ∶ ( A𝜆 2 +B𝜆+C ) x=0 for some x≠ 0 5 See, e.g., Dennis etal. (1976,p.835) or Gantmacher (1959, vol.I,p.228). 599 Solving Linear DSGE Models withBernoulli Iterations which has been reviewed extensively by Tisseur and Meerbergen (2001) and for which Hammarling etal. (2013) provide a comprehensive method to improve the accuracy of its solutions. The matrix quadratic (16) can be expanded following Higham and Kim (2001) as where DPM(ΔP) is the Fréchet derivative of M at P in the direction ΔP . 3.1 Baseline Bernoulli Iteration The multivariate counterpart to the algorithm x k+1=f ( xk ) =− c ax k + b above is the following baseline Bernoulli method for the minimal solvent. Beginning with (16) and defining the iteration via gives the baseline Bernoulli method This iteration was noted by Binder and Pesaran (1999) and is implemented in their RBCQDE Matlab and Gauss code,6 but without discussion of its convergence or computational properties. Rendahl (2017) studies the Bernoulli recursion in (25) following Higham and Kim (2000) by analyzing the theoretical properties and exploring a few examples. The baseline Bernoulli recursion is summarized in Algorithm1. (19) M(P+ΔP)=A(P+ΔP)2+B(P+ΔP)+C (20) =AP 2+BP+C+AΔP2+A(PΔP+ΔPP)+BΔP (21) =M(P)+(AΔPP +(AP+B)ΔP)+AΔP2 (22) =M(P)+ D P(ΔP)+AΔP2 (23) 0=(AP+B)P+C (24) 0 = ( AP j +B ) P j+1 +C (25) Pj+1 =− ( AP j +B )−1C 6 See https:// dge. repec. org/ codes/ binder/ handb ook/. 600 A.Meyer-Gohde Algorithm1: Baseline Bernoulli method Higham and Kim (2000) show that the above recursion converges asymptotically at a linear rate to the minimal solvent, but only under the assumption that this solvent is invertible (e.g., ruling out zero latent roots in (17)) and that a dominant solvent also exists (e.g., ruling out “infinite” latent roots in (17)), both of which abound in the DSGE literature. I return to this issue in Sect.4 and provide a proof that this recursion converges at least locally to the minimal solvent when this solvent is the unique stable solution. Both at such a solvent and for the first iteration if P0=0 and B is of full rank, see Binder and Pesaran (1997) and Uhlig (1999), the coefficient matrix AP j+B will be invertible.7 3.2 Bai andGao’s (2007) Modified Bernoulli Method The modified Bernoulli method of Bai and Gao (2007) is a straightforward extension of the Bernoulli method that updates the solution matrix column by column alongside the standard iteration above, incorporating the update of the previous column when solving for the current column. Bai and Gao (2007) note that the Bernoulli recursion solves n linear equations at each iteration with Pi,j+1 and Ci the i’th columns of Pj+1 and C respectively (26) Pj+1 =− ( AP j +B )−1C 7 In practice, see the applications later, singularity is not a general issue, but in the few iterations where it did occur, I set Pj + 1 to the minimum 2-norm solution of ( AP j +B ) P j +1=− C , see Higham and Kim (2000, p. 512) determination that the algorithm can be rerun or reinitialized if numerical difficulties should be encountered. 607 Solving Linear DSGE Models withBernoulli Iterations or the weighting above could be tilted towards one or the other increment with an exponent p<1 shifting the weight in favor of the Bernoulli increment (and hence sj(p) replacing sj in (54)), perhaps to reduce the likelihood of the algorithm converging to a solvent other than the unique stable one. This gives a combined Bernoulli–Newton procedure as Algorithm4. Algorithm4: Combined Bernoulli–Newton While this procedure will hopefully yield the best of both worlds: convergence to the desired solvent via Bernoulli and quadratic convergence via Newton, it might also do just the opposite—combining the linear convergence of Bernoulli to the unpredictable solvent of Newton. Tilting the weight towards one or the other procedure gives the user flexibility, but it is still not clear a priori how that tilt might be chosen to deliver a procedure with the desired properties. 3.5 Bernoulli andNewton Line Search As the weighting of the Bernoulli and Newton increments relied on the same intuition that guided their respective line search algorithms, the next logical step would be to combine both the line searches for the optimal magnitudes of the respective increments and then weight the two to produce a combined Bernoulli and Newton increment with line searches. That is, if combining Bernoulli and Newton might bring the slow convergence of Bernoulli together with the (62) cos 𝜃i,j= Δ NP � i,j Δ BPi,j ‖ ‖ ‖ ΔNPi,j ‖ ‖ ‖2 ‖ ‖ ‖ ΔBPi,j ‖ ‖ ‖2 (63)  s j(p)=s p j 608 A.Meyer-Gohde unpredictable convergence destination of Newton, performing line searches on the steps of both of these might attenuate this danger. From above, the step size of Bernoulli can be increased past 1 to improve the speed of convergence, and the step size of Newton can be adjusted optimally to balance the danger of taking too large of steps. This gives a combined Bernoulli–Newton procedure with line searches as Algorithm5. Algorithm5: Combined Bernoulli–Newton with line searches As above, the underlying steps in the procedure would tend to limit the drawbacks of the two procedures on their own. Yet, it is not a priori certain how they will perform when combined. 3.6 Optimal Bernoulli andNewton The final set of algorithms explicitly take the optimality approach when determining sj , minimizing the same merit function used in the line searches above. First for the initial increments and then for the line-search optimized increments where (64) s =argmin 0≤x≤1 ‖ ‖ M(xΔBP+(1−x)ΔNP+P) ‖ ‖ 2 F (65) s =argmin 0≤x≤1 ‖ ‖ M(xtBΔBP+(1−x)tNΔNP+P) ‖ ‖ 2 F 609 Solving Linear DSGE Models withBernoulli Iterations where XB restricts x to be greater or equal to one, see above, and XN restricts x to be between zero and two (Higham & Kim, 2001). This gives an optimally (in the sense of the merit function) combined Bernoulli–Newton procedure in algorithm6. Algorithm6: Optimal Bernoulli–Newton This algorithm has the advantage of weighting the Bernoulli and Newton increments in a non-arbitrary manner. Yet this comes at a cost, here of an additional optimization problem to be solved, and is likely biased towards the Newton increment as it—intuitively via quadratic convergence—generally takes larger steps, making each increment more likely to be favored over the timid Bernoulli one. This carries again the potential of losing the advantage of Bernoulli as formulated in the baseline algorithm that guarantees convergence to a particular solvent. 4 Theoretical andPractical Considerations 4.1 Initial Value All Bernoulli iterations need an initial value, P0 . In contrast to Newton methods, the baseline Bernoulli method has strong convergence results, see the next subsection, and hence this initial value is of lesser importance. Yet many of the algorithms presented above combine Bernoulli and Newton methods and are subject to this concern. As the goal is to obtain the minimal solvent P, I choose the initial value P0=0 . In the absence of any other guidance, this choice satisfies the requirement of having all eigenvalues inside the unit circle. (66) t i=argmin x∈Xi ‖ ‖ M(xΔiP+P) ‖ ‖ 2 F, for i=B, N 610 A.Meyer-Gohde Furthermore, Higham and Kim (2000,p.512) note that the Bernoulli algorithm can be rerun or restated with a different initialization should numerical difficulties be encountered. For the baseline Bernoulli iteration which solves the (near) singularity of AP j+B would pose such a difficulty. In this case, the rank deficiency of the leading coefficient matrix admits multiple solutions and I chose the norm ( min ‖ ‖ ‖( AP j+B ) Pj+1+C ‖ ‖ ‖F ) minimizing Pj+1 =− ( AP j +B )+C , where + indicates the Moore-Penrose inverse. 4.2 Convergence While Higham and Kim (2000) and Bai and Gao (2007) provide convergence results for Bernoulli iterations, their analyses assume that A is nonsingular and that both a minimal and a dominant solvent exist. This is untenable in DSGE models where singular A’s abound—associated with variables that arise only at time t and t−1 —and singular C’s—associated with variables that arise only at time t+1 and t—prevent the application of their results to the reverse quadratic. Rendahl (2017) adapts Higham and Kim’s (2000) approach to DSGE models and provides proofs of the global convergence of the Bernoulli iteration for arbitrarily small perturbations of an originally singular model that produce a nonsingular system. Complementing this, I provide a local convergence proof of the original unperturbed problem in the following. Recall the baseline Bernoulli method, the Fréchet derivative of F at Pj in the direction ΔPj is D P j F ( ΔPj ) and be defined implicity by differentiating ( AP j +B ) F ( P j) =− C using the Kronecker / vectorized representation ( vec (ABC)=(C�⊗A)vec (B) ) the algorithm converges (locally) to the minimal solvent. (67) ( AP j +B ) P j+1 =− C (68) Pj+1 =− ( AP j +B )−1 C=F ( P j) (69) A ΔPjF ( Pj ) + ( AP j+B ) DP j F ( ΔPj ) = 0 (70) D P j F ( ΔPj ) =− ( AP j+B )−1 AΔPjF ( Pj ) (71) D P j F ( ΔPj ) = ( AP j+B )−1 AΔPj ( AP j+B )−1 C (72) vec (DPjF ( ΔPj ) )= ([( AP j+B ) −1C ] � ⊗ [( AP j+B ) −1A ]) vec (ΔPj ) 611 Solving Linear DSGE Models withBernoulli Iterations Theorem 1 (Convergence to the unique stable solvent P) Assume there exists a unique solvent P of M(P)≡AP 2+BP+C in (16) such that the eigenvalues of P comprise all the latent roots, 𝜆 of M(𝜆) —see 17—stable with respect to the closed unit circle, i.e., less than or equal to one in absolute value. Then the Bernoulli iteration Pj+1 =− ( AP j +B )−1C is stable in the neighborhood of P. Proof See the “Appendix”. ◻ The conditions for the existence of the unique stable solvent P are Blanchard and Kahn’s (1980) celebrated rank and order conditions, see Lan and Meyer-Gohde (2014) and Meyer-Gohde (2023) for the conditions expressed in terms of the general class of multivariate singular leading A models pervasive in the literature today. So conditional on its existence, the Bernoulli method above will converge asymptotially to P. Equally interesting is the performance in the absence of a unique stable solvent P, either in the case of nonexistence or, especially, indeterminacy. In the latter case, the model permits an entire continuum of solutions and some studies, such as Lubik and Schorfheide (2003), analyze the implications of the model in such a case. Both the results here, theorem1, and Higham and Kim’s (2000) approach in the absence of singular coefficient matrices, show that the Bernoulli method converges to the minimal solvent. That is, that solvent P of 0=M(P)≡AP 2+BP+C in such that the eigenvalues of P the smallest latent roots, 𝜆 of M(𝜆) , that is, with the smallest possible spectral radius. In the case of indeterminacy, more latent roots than the dimension of the problem (i.e., more than ny the length of the vector yt and the number of equations in f), the Bernoulli method will return that stable solution with the smallest roots, leaving the largest stable root(s) out of the solution. For a non existence, the Bernoulli method will return that solution with the smallest (in absolute value) roots again, but some of them with be outside the unit circle and hence the solution returned will be unstable. This follows locally from theorem1 as the crossproducts of the ny smallest latest roots with inverses of the ny largest latest roots will be bounded above by one generically.8 While the convergence to a specific solvent (in this formulation, the unique stable one) is an advantage over Newton methods, which cannot guarantee convergence to a particular solvent (see Higham & Kim, 2001), Bernoulli methods converge only at a linear rate (given above by the ratio of the largest stable and smallest unstable eigenvalues) instead of Newton methods’ quadratic rate. In practice, I follow Higham and Kim (2001) and use the relative residual � � � M(Pj) � � �F ∕ �‖ A ‖ F � � � P2 j � � �F + ‖ B ‖ F � � � Pj � � �F + ‖ C ‖ F � <ny 𝜖 to assess whether convergence has occurred. 8 The exception being if eigenvalues of identical magnitude are included in both sets. Limiting the permissible cases to �𝜆 1 � ≤…≤ ‖𝜆 n y �<�𝜆 n y +1≤ ‖𝜆 2n y � rules out this exceptional case. 612 A.Meyer-Gohde 4.3 Accuracy The practical forward error bounds of Meyer-Gohde (2023) can be used to assess the accuracy of a computed solution  P where R P= A  P2 + B  P +C is the residual of the matrix quadratic and H  P=In y ⊗ ( A  P+B ) +  P�⊗ A . Stewart’s (1971) separation function, see also Kågström (1994), Kågström and Poromaa (1996), and Chen and Lv (2018), is where 𝜆( A,A  P+B ) is the spectrum or set of (generalized) eigenvalues of the pencil ( A,A  P+B ) (and, accordingly, 𝜆( P ) the set of eigenvalues of  P ) and the last line holds with equality for A=I and  P and  P+B regular—hence, the separation between the two pencils—the smallest singular value of H P —is generically smaller than the minimal separation between their spectra. Analogously to the generalized Sylvester and algebraic Riccati equations, the separation function provides the natural extension of the conditioning number from standard linear equations to these structured problems, and the a posteriori condition number for the matrix quadratic is given by sep −1 [( A,A P+B ) , ( I,− P )] = ‖ ‖ ‖ H−1  P‖ ‖ ‖2 =𝜎min ( H P )−1 , which—from above—can be arbitrarily larger than the inverse of the minimal distance between the spectra of the pencils ( A,A  P+B ) , ( I,−  P ) . This inverse of the separation relates an upper bound to the forward error directly to the residual, like the condition number for a standard linear system, and a tighter bound takes into account the structure more carefully and considers the linear operator H P and the residual R P jointly. 5 Applications I conduct two sets of experiments to assess the performance of the Bernoulli algorithms presented above. These two sets are chosen to assess the different methods in a specific, policy relevant model but also in as non-model specific an (73) � � �P− P � � �F ‖ P ‖ F ������� Forward Error ≤ � � �H−1  P vec � R P �� � �2 � � �  P � � � F ������������������� Forward Error Bound 1 ≤� � � H−1  P� � � 2 � �R P� �F � � �  P � � � F ��������������� Forward Error Bound 2 (74) sep�� A,A  P+B � , � I,−  P �� =min ‖X‖F=1� � � AX  P+ � A  P+B � X � � � F (75) = min ‖vec(X)‖2=1 � � H Pvec(X) � �2 (76) = 𝜎min ( H P ) ≤min | | | 𝜆 ( A,A P+B ) −𝜆 (  P )| | | 613 Solving Linear DSGE Models withBernoulli Iterations environment as possible. This is in contrast to much of the existing literature that compares alternate solution in a single model, frequently a stochastic growth or real business cycle model, such as Taylor and Uhlig (1990), Aruoba et al. (2006), or Caldara etal. (2012) or in a few stylistic models and in Rendahl (2017). To more systematically asses the different Bernoulli based algorithms from above, I compare these algorithms with Dynare’s QZ implementation9 and the baseline Newton algorithm from Meyer-Gohde and Saecker (2024).10 both in the model of Smets and Wouters (2007) and on the suite of models in the Macroeconomic Model Data Base (MMB) (see Wieland et al., 2012, 2016), a model comparison initiative at the Institute for Monetary and Financial Stability (IMFS),.11 The performance is measured in terms of accuracy, computational time, and convergence to the stable solvent, initializing both from zero matrix (an uninformed initialization of a stable solvent) and the output from the QZ algorithm. 5.1 Smets andWouters’s (2007) Model I begin with the medium scale, estimated model of Smets and Wouters (2007) that is arguably the benchmark for policy analysis. In their model they analyze and estimate a New Keynesian DSGE model with US data featuring the usual frictions, sticky prices and wages, inflation indexation, consumption habit formation as well as production frictions concerning investment, capital and fixed costs. Among the equations is the following log-linearized monetary policy rule that will be the focus of the final experiment, assessing the accuracy of the methods here when solving under alternate, but nearby parameterizations. This Taylor rule sets the interest rate rt according to inflation 𝜋t , the current output gap ( y t −y p t) and the change in the output gap, with the parameters r𝜋 , rY and rΔy describing the strength of each of these reactions and 𝜌 controlling the degree of interest rate smoothing. The monetary policy shock, 𝜀r t , follows an AR(1)-process with an iid normal error. The Bayesian estimation of the model employs seven macroeconomic time series from the US economy to estimate the model parameters and the authors show that the model matches the US macroeconomic data very closely and that out-of-sample forecasting performance is favorable compared to (B) VAR models. (77) rt = 𝜌 r t−1 +(1− 𝜌 )(r 𝜋𝜋t +r Y (y t −y p t )) + r Δy ((y t −y p t )−(y t−1 −y p t−1 )) + 𝜀r t, 9 A further development of Klein (2000), see Villemot (2011). Note that Dynare uses the real Schur decomposition as provided by LAPACK’s routine DGGES, see https:// git. dynare. org/ Dynare/ dynar e/-/ tree/ master/ mex/ sourc es/ mjdgg es. 10 Additionally, note that I follow Dynare and reduce the dimensionality of the problem by grouping variables and structuring the matrix quadratic according to the classification of “static”, “purely forward”, “purely backward looking”, and “mixed” variables. The details are in the appendix and Meyer-Gohde and Saecker (2024). 11 See http:// www. macro model base. com 614 A.Meyer-Gohde Table1 summarizes the results at the posterior mode calibration of the model of Smets and Wouters (2007). The baseline Bernoulli method takes 12 times longer than Dynare’s QZ, which would appear to put it at a disadvantage compared to the baseline Newton, however the large maximal absolute difference to the QZ solution of the Newton algorithm shows that it has converged to a different solution than the unique stable solution found by Dynare’s QZ. Indeed, this danger looms with Newton related algorithms as can be seen here for the Newton–Bernoulli and Newton–Bernoulli optimal algorithms, both of which also converged to a different solvent. The baseline Bernoulli required 440 iterations and line searches reduced this number to 420, but the reduction in iterations was outweighed by the costliness of the line search algorithm, resulting altogether in a longer computation time. The columnwise Newton–Bernoulli and the modified Bernoulli both took extraordinarily long times to solve the model, with the Newton–Bernoulli with line searches and optimal Newton–Bernoulli with line searches providing solutions within an order of magnitude of computation time relative to Dynare’s QZ, 33 and 19 iterations respectively to do so, and providing solutions that are an order of magnitude more accurate than QZ. Table 2 assesses the different methods as solution refinement techniques, by parameterizing the model of Smets and Wouters (2007) within the prior to Table 1 Results: model of Smets and Wouters (2007), posterior mode For Dynare, refer to Adjemian etal. (2011) Run Time for Dynare in seconds, for all others, run time relative to Dynare Max Abs. Diff. measures the largest absolute difference in the computed P of each method from the P produced by Dynare Forward error 1 and 2 are the upper bounds for the true forward error, see (73) Method Relative performance Forward errors Iterations Run time Max Abs. Diff Bound 1 Bound 2 Dynare (QZ) 0.00063 5.5e−14 2.4e−11 1 Baseline Newton 2.4 108 4.1e−14 1.6e−09 11 Baseline Bernoulli 12 7.7e−13 3.8e−14 2.6e−11 436 Modified Bernoulli (MBI) 132 8.8e−13 3.5e−14 2.5e−11 423 Bernoulli with Line Searches 20 6.9e−13 3.7e−14 2.4e−11 423 Newton–Bernoulli 8.9 108 2.3e−14 1.5e−09 39 Newton–Bernoulli Column 611 6.9e−13 0.49 3.8e+06 2385 Newton–Bernoulli 1/3 11 6.8e−13 4.7e−15 4.1e−12 47 Newton–Bernoulli LS 8.5 6.7e−13 1.1e−14 2.2e−12 33 Newton–Bernoulli Column LS 9.3 7.2e−13 5.1e−15 9.6e−12 32 Newton–Bernoulli LS 1/3 15 5.9e−13 1.1e−14 5.1e−12 57 Newton–Bernoulli Opt 5.6 108 2.3e−14 1e-09 18 Newton–Bernoulli Opt LS 6.5 7.7e−13 1.5e−15 2.2e−12 19 615 Solving Linear DSGE Models withBernoulli Iterations demonstrate a numerical instability whose consequences are economically relevant, with moments of variables predicted by Dynare’s QZ differing in all digits. The second column now displays the variance of inflation as predicted by the different solution methods,12 At this parameterization, Dynare’s QZ solution predicts an inflation variance of 0.28. However, even the lower of the two upper bounds on the forward error is consistent with a numerical instability being several orders of magnitude beyond machine precision. In the final experiment with the model of Smets and Wouters (2007), I use algorithms to solve iteratively for different parameterizations of the Taylor rule. The goal here is to explore whether solutions from previous, nearby parameterizations can be used to efficiently initialize the Bernoulli methods similarly to the experiment above with the QZ solution as the initial guess. For the parameters determining the Taylor rule reaction to inflation and the long run reaction to the output gap, the experiment iterates through a grid of 10 ×10 parameter values varying the size of the interval considered—setting r𝜋∈[1.5, 1.5 (1+10−x)] and rY∈[0.125, 0.125 (1+10−x)] , where x∈ [−1, 8] (Smets & Wouters (2007) calibrate them to r𝜋 =2.0443 and Table 2 Results: model of Smets and Wouters (2007), numerically problematic parameterization For Dynare, refer to Adjemian etal. (2011) Run time for Dynare in seconds, for all others, run time relative to Dynare Variance 𝜋t gives the associated value for the population or theoretical variance of inflation—note that two algorithms did not converge to a stable P and hence the variance could not be calculated for them Forward error 1 and 2 are the upper bounds for the true forward error, see (73) Method Relative performance Forward errors Iterations Run time Variance 𝜋t Bound 1 Bound 2 Dynare (QZ) 0.0046 0.28 1e-11 4.6 1 Baseline Newton 4.3 0.45 2.4e−14 0.00058 4 Baseline Bernoulli 3 0.39 3.8e−13 0.018 90 Modified Bernoulli (MBI) 2814 – 2.9 6.6e+06 5e+04 Bernoulli with Line Searches 504 – 0.99 1.9e+11 5e+04 Newton–Bernoulli 5.1 0.39 3.3e−15 0.00081 8 Newton–Bernoulli Column 3.7 0.4 3.7e−14 0.019 5 Newton–Bernoulli 1/3 1.2 0.46 1.4e−14 0.0027 8 Newton–Bernoulli LS 7.7 0.39 7.9e−15 0.0023 6 Newton–Bernoulli Column LS 5.7 0.44 1.7e−14 0.00079 6 Newton–Bernoulli LS 1/3 1 0.38 3.3e−14 0.0022 8 Newton–Bernoulli Opt 3.7 0.45 2.4e−14 0.00058 4 Newton–Bernoulli Opt LS 4.2 0.5 1.3e−15 0.0011 4 12 See Meyer-Gohde (2023) for more details on the parameterization. Smets and Wouters (2007) report a variance of inflation in the entire sample of 0.62 and 0.55 and 0.25 in two subsamples. 616 A.Meyer-Gohde rY=0.0882 ). The algorithm iterates through the two-dimensional grid taking the solution under the previous parameterization as the initialization for the next iteration. A decrease in the spacing between the 100 grid points thus increases the precision of the starting guess, the solution from the previous parameterization.13 Figure 1 summarizes the experiment graphically. Figure 1c confirms a decrease in run time per grid point with a narrower grid for the iterative algorithms here and an irrelevance of the grid spacing for QZ. As the grid becomes narrower, the iterative Bernoulli and Newton procedures increasingly benefit from starting from the solution of the previous iteration as it becomes closer to the unknown solution of the current iteration. The QZ algorithm does not operate Fig. 1 Forward errors and computation time per grid point for different parameterizations of the model by Smets and Wouters (2007). Figure1a, b plot the upper forward error bounds 1 and 2 against the grid size, log10 scale on both axes. Figure1c plots the computation per grid point against the number of grid points, log10 scale on both axes 13 Alternatively, one could fix the end points of the grid and increase the number of grid points—this would provide the same message but would possess a computationally prohibitive number of grid points at the narrowest spacing I examine here. 623 Solving Linear DSGE Models withBernoulli Iterations endogenous variables—and forward errors (bound 1 and bound 2) relative to QZ are plotted in log10. The most striking result is the difference of accuracy, with algorithms involving a Newton step more often below the x axis than the remaining methods, indicating that they are generally associated with more accuracy—an observation I return to shortly in a density comparison. There appears to be a downward trend, indicating that the methods become more accurate relative to Dynare’s QZ for larger models, which is confirmed by looking at Fig.4d, which shows the clear downward trend for the optimal Bernoulli–Newton method with line searches and is exemplary for many of the methods. Finally, Fig.4c shows that there appears to be a relationship between the size of the model and the computation time relative to Dynare for at least the very large models towards the right of the figure, indicating that the methods here are likely to be particularly competitive alternatives for larger scale applications, with Meyer-Gohde and Saecker’s (2024) Newton algorithm presenting the most convincing evidence in this regard. Fig. 6 Distribution of forward error bounds relative to dynare for the macroeconomic model data base (MMB). Figure6a, b plot the distribution of model solutions against the upper bounds of the forward error 1 and 2 for all algorithms, log10 scale on the x axis, 99 MMB models (starting guess: solution Dynare(QZ)) 624 A.Meyer-Gohde Figure 5 provides an overview of the entire distribution of forward errors, the upper row relative to those from Dynare’s QZ method and the lower in absolute terms. Forward errors left of the vertical line are thus smaller than Dynare for both figures in the upper row. For both the first, Fig.5a, and second, Fig.5b, upper bounds on the forward error, there is an obvious shift to the left of about one order of magnitude for all the methods involving a Newton step and less visually compelling evidence for Bernoulli algorithms not combined with Newton (again, the modified Bernoulli algorithm performs worst). From the lower row, this entails tightening the distributions as well as shifting them closer to machine precision—a lower convergence criterion would allow more iterations and likely bring yet more solutions below machine precision. To assess the potential for improving on solutions, I repeat the exercise, but now initialize with the solution provided by QZ, see table4. Here the baseline Bernoulli method is the top performer—with its low per iteration cost, it runs one iteration at a small fraction of the original Dynare QZ cost and provides a significant improvement in accuracy. The modified Bernoulli algorithm again performs unsatisfactorily and, interestingly, the combined Bernoulli and Newton methods all require more than one iteration to converge, which would seem to imply that the Newton and Bernoulli steps individually were generally pulling in different directions in the vicinity of the solution provided by Dynare’s QZ. Note that the modified Bernoulli algorithm along with the Newton Bernoulli Column algorithms did not always converge—all three of these algorithms operate column-wise on the problem which apparently can interfere with the convergence, even when initialized close to the solution. Figure6, like Fig.5 but now initialized at the Dynare’s QZsolution, provides an overview of the entire distribution of forward errors, the upper row relative to those from Dynare’s QZ method and the lower in absolute terms. Forward errors left of the vertical line are thus smaller than Dynare’s QZ for both figures in the upper row. For both the first, Fig.6a, and second, Fig.6b, upper bounds on the forward error, there is again an obvious shift to the left of about one order of magnitude for all the methods involving a Newton step and a marginal at best improvement using algorithms without a Newton step. This is consistent with the quadratic convergence properties of Newton methods, see Meyer-Gohde and Saecker (2024), when close to a solution. The modified Bernoulli algorithm now performs comparably to the other methods without a Newton step, highlighting some applicability of the conclusions of Bai and Gao (2007) to the DSGE context. 6 Conclusion I have applied and extended Bernoulli-based methods for solving the matrix quadratic equation underlying the solution of linear DSGE models. This adds a set of alternatives alongside Meyer-Gohde and Saecker’s (2024) Newton-based and Huber etal. (2023) doubling algorithms to the current standard of a generalized Schur or QZ decomposition (Moler & Stewart, 1973; Golub & van Loan, 2013). Applying the 625 Solving Linear DSGE Models withBernoulli Iterations methods to the suite of models in the Macroeconomic Model Data Base (MMB), I demonstrate that Bernoulli-based methods are a potential alternative, with a tradeoff between convergence to a particular solution (here the unique, stable solvent) and performance in terms of computational costs. Particularly in iterative environments or when a solution refinement is sought do these algorithms show potential for future application. In filling in an increasingly dense grid of parameterizations for the Taylor rule in the model of Smets and Wouters (2007), iterative methods like the Bernoulli-based methods here can initialize with the solution from the previous parameterization and significantly outperform the current generalized Schur or QZ method both in terms of computational costs and forward error. Taking the solution from QZ as the initialization, the methods provide roughly an order of magnitude improvement in the accuracy of the solution at a fraction of the original computational cost. This initialization and iteration makes applying the set of Bernoulli methods to improving the accuracy of solutions to linear DSGE models a potentially useful direction of application, which is demonstrated here with a problematic parameterization of the Smets and Wouters (2007) model as presented by Meyer-Gohde (2023). Starting with the QZ solution that underpredicts the variance of inflation by about 25% and demonstrates large forward errors, the Bernoulli method refines this solution, providing forward errors several orders of magnitude smaller and a predicted variance of inflation in line with Newton methods. That being said, the baseline Bernoulli algorithm is not a uniform improvement over the standard QZ method and, in particular, for most models when an informative initialization is not available, the standard QZ is to be preferred. Future research might explore the application of the methods here to reduce the computational burden associated with solving the model for iterative estimation procedures and might be adapted to more quickly and/or accurately perform likelihood calculations or solve heterogenous agent models. 626 A.Meyer-Gohde Appendix Overview The following flow chart visualizes the decision process that leads to each algorithm. The code to replicate the analysis can be found at https:// github. com/ AlexM eyerGohde/ LinearDSGEwithBerno ulli The directory algorithm contains the algorithms used here (with the exception of Dynare with can be downloaded at https:// www. dynare. org/). The file bernoulli_matrix_quadratic.m contains the Bernoulli algorithms introduced here and outlined by the flowchart above. The directory mmb_experiment contains the MMB experiment starting at the zero matrix in Table 3 and the Smets and Wouters (2007) comparison of Table1 (the model of Smets and Wouters (2007) is among the 99 MMB models used in the comparison, themselves contained in the directory mmb_experiment). 627 Solving Linear DSGE Models withBernoulli Iterations The directory mmb_experiment_2 contains the MMB experiment starting at the QZ solution of Table 4. The numerically problematic parameterization of Smets and Wouters (2007) displayed in Table2 can be found in the directory improvement_experiment. Finally the grid size experiment for Smets and Wouters (2007) of Fig.1 is in the directory policy_experiment_2. Detailed Dynare Variable Classification Here I summarize the details in the matrix quadratic that follows from the classification of variables from Dynare as laid out in Villemot (2011). See MeyerGohde and Saecker (2024) for details. Subdividing the system of equations in accordance with the QR decomposition yields where nd is the number of dynamic variables, the sum of number of purely backward-looking, n −− , mixed nm , and purely forward-looking variables, n++ . The number of forward-looking variables, n+ , is the sum of the number of mixed, nm , and purely forward-looking variables, n++ , and the number of backward-looking variables, n − , is the sum of the number of purely backward-looking, n −− and mixed variables nm . Hence, the number of endogenous variables is the sum of the number of static, ns , and dynamic variables, nd , or the sum of the number of static, ns , purely backward-looking, n −− , mixed nm , and purely forward-looking variables, n++ . The dimensions satisfy the following The transition matrix, P, from (2) that solves the matrix equation (16) can be subdivided in accordance to Dynare’s typology as n d =n−− +n m +n++ ,n+=n m +n++ , n −=n−− +nm , n=ns+nd=ns+n−− +nm+n ++ 628 A.Meyer-Gohde The matrix quadratic can be expressed as For a solvent P of the matrix quadratic, taking the structure of C from the Dynare typology above into account yields Following Meyer-Gohde and Saecker (2024), who apply corollary 4.5 of Lan and Meyer-Gohde (2014), if P is the unique solvent of M(P) stable with respect to the closed unit circle, 𝐆 has full rank and hence the columns of P associated with nonzero columns in C, the static and forward-looking variables are zero →𝐏 ∙,s =𝟎 n×n s ,𝐏 ∙,++ =𝟎 n×n++ , whence 𝐏 is and 𝐌 (𝐏)= [ 𝟎 n×ns &𝐌(𝐏)−− n×n −− 𝐌(𝐏)m n×n m 𝟎 n×n++ ] . Consequentially, the first ns rows of the matrix quadratic, taking as given, yield as 𝐌 (𝐏 n×n )= 𝐀 n×n 𝐏 2 +𝐁 n×n 𝐏+𝐂 n× n = ( 𝐀𝐏 +𝐁 ) ⏟⏞⏞⏟⏞⏞⏟ ≡𝐆 𝐏+𝐂 629 Solving Linear DSGE Models withBernoulli Iterations and the first ns rows of P are . The last nd columns and rows of P solve the reduced matrix quadratic equation Recalling that 𝐏 ∙,++ =𝟎 n×n++ ,  𝐏 can be reduced and two submatrices 𝐏 and  𝐏 defined via 630 A.Meyer-Gohde where and allow the matrix quadratic to be written as which can be reduced to This leads to the Bernoulli iteration 631 Solving Linear DSGE Models withBernoulli Iterations Detailed Dynare Topology‑Line Search The line search methods in the text require finding the minimum of the merit function The first order condition for an interior solution requires finding the zeros of the polynomial where 𝛼 = ‖ M(P) ‖2 F , 𝛽 = ‖ Ad𝐏⋅P ‖2 F , 𝛾 = ‖ ‖ ‖ A(d𝐏)2 ‖ ‖ ‖ 2 F , 𝜉 =tr ( (Ad𝐏⋅𝐏) ∗ A(d𝐏) 2+ ( A(d𝐏)2 ) ∗Ad𝐏⋅𝐏 ) , 𝜎 =tr ( M(P)∗A(d𝐏)2+ ( A(d𝐏)2 )∗ M(P) ) , and 𝛿= tr (M(P)∗ Ad 𝐏⋅𝐏+(Ad𝐏⋅𝐏) ∗ M(P) ) , which follows from ‖ M(P+xΔP) ‖2 F=(1−x) 2‖ M(P) ‖2 F+x 2‖ AΔP(xΔP+P) ‖2 F +(1−x)x∗tr(M(P)*AΔP(xΔP+P) +(xΔP*+xP*)ΔP*A*M(P)) =(1−x)2‖M(P)‖2 F+x2‖AΔP(xΔP+P)‖2 F +(1−x)x2∗tr(M(P)*AΔP2+ΔP*2A*M(P)) +(1−x)x∗tr(M(P)*AΔP⋅P+P*ΔP*A*M(P )) =(1−x)2‖M(P)‖2 F+x4� � �AΔP2� � � 2 F+x2‖AΔP⋅P‖2 F +x3∗tr(ΔP*2A*AΔP⋅P+P*ΔP*A*AΔP2) +(1−x)x2∗tr(M(P)*AΔP2+ΔP*2A*M(P)) +(1−x)x∗tr(M(P)*A𝛿P⋅P+P*ΔP*A*M(P)) ≡ g(x) (A1) g�(x)=4𝛾x3+3(𝜉−𝜎)x2+2(𝛼+𝛽−𝛿+𝜎)x+𝛿−2𝛼 632 A.Meyer-Gohde Using the typology from Dynare and the results above g� (x)=2(1−x)(−1)‖M(P)‖ 2 F +2x‖AΔP⋅P‖2 F+4x3� � �AΔP2� � � 2 F +3x2∗tr(ΔP*2A*AΔP⋅P+P*ΔP*A*AΔP2) + [(1−x)2x+ (−1)x2]∗tr(M(P)*AΔP2+ΔP*2A*M(P)) + [(1−x)x+ (−1)x]∗tr(M(P)*AΔP⋅P+P*ΔP*A*M(P)) =−2‖M(P)‖2 F+2‖M(P)‖2 Fx+2‖AΔP⋅P‖2 Fx+4� � �AΔP2� � � 2 Fx3 +3∗tr(ΔP*2A*AΔP⋅P+P*ΔP*A*AΔP2)x2 +tr(M(P)*AΔP⋅P+P*ΔP*A*M(P)) − 2∗tr(M(P)*AΔP⋅P+P*ΔP*A*M(P)) x +2∗tr(M(P)*AΔP2+ΔP*2A*M(P))x −3∗tr(M(P)*AΔP2+ΔP*2A*M(P))x2 =4� � �AΔP2� � � 2 Fx3+(3∗tr(ΔP*2A*AΔP⋅P+P*ΔP*A*AΔP2) −3tr(M(P)*AΔP2+ΔP*2A*M(P)))x2 +2‖M(P)‖2 Fx+2‖AΔP⋅P‖2 Fx −2∗tr(M(P)*AΔP⋅P+P*ΔP*A*M(P))x +2∗tr(M(P)*AΔP2+ΔP2*A*M(P))x +tr(M(P)*AΔP⋅P+P*ΔP*A*M(P)) − 2‖M(P)‖2 F 639 Solving Linear DSGE Models withBernoulli Iterations Unsing the Dynare typology above ‖ 𝐌(⋅)‖ 2 F=‖𝐚‖ 2 F+s(tr(𝐚 ∗ ⋅𝐛)+tr(𝐛 ∗ ⋅𝐚))+s 2 ⋅ � ‖𝐛‖ 2 F+tr(𝐚 ∗ 𝐜)+tr(𝐜 ∗ 𝐚) � +s3(tr(𝐛∗⋅𝐜)+tr(𝐜∗⋅𝐛))+s4 ⋅‖𝐜‖2 F =‖𝐚‖2 F+s⋅2⋅tr(𝐚∗𝐛)+s2 ⋅�‖𝐛‖2 F+2tr(𝐚∗𝐜)�+s3 ⋅2tr(𝐛∗⋅𝐜)+s4 ⋅‖𝐜‖2 F ≡𝐭(s) 𝐭�(s)=2⋅tr(𝐚∗𝐛)+2⋅( ‖ 𝐛 ‖ 2 F +2⋅tr(𝐚∗𝐜)⋅s+6⋅tr(𝐛∗⋅𝐜)⋅s2+4⋅ ‖ 𝐜 ‖ 2 F ⋅s3 640 A.Meyer-Gohde Proof ofTheorem1 Following Higham (2008,Section4.9.4.), local stability requires the Fréchet derivative of F at P to have bounded powers, which holds if its spectral radius is less than one. At P, (72) is Hence the spectral radius is equal to the eigenvalue of P� ⊗ [ (AP+B) −1 A ] with the largest magnitude. As the eigenvalues of the Kronecker product of two matrices is the set of all crossproducts of the eigenvalues of the respective matrices, the spectral radius is the product of the eigenvalues of P and (AP+B)−1A with the largest magnitudes. Following Lan and Meyer-Gohde (2014,Corollary 4.2.), (17) can be factored as Hence, the latent roots associated with A𝜆+AP+B are those roots of M(𝜆) not contained in set of eigenvalues of P. By assumption, the eigenvalues of P are inside the closed unit circle and the roots associated with A𝜆+AP+B are outside the open unit circle. By inspection, the eigenvalues of (AP+B)−1A are the inverses of the roots associated with A𝜆+AP+B and, therefore, are all inside the open unit circle. Hence the product of the largest magnitude eigenvalue of P and that of (AP+B)−1A is less than one in absolute value. That is, the spectral radius of the Fréchet derivative of F at P, the unique bounded solution, is less than one, completing the proof. Acknowledgements I am grateful to Johanna Saecker, Pablo Winant, and participants of the 29th International Conference on Computing in Economics and Finance (2023 CEF) for useful comments and suggestions and to Elena Schlipphack and Maximilian Thomin for invaluable research assistance. The code to replicate the analysis can be found at https:// github. com/ AlexM eyerGohde/ LinearDSGEwithBerno ulli. Any and all errors are entirely my own. Funding Open Access funding enabled and organized by Projekt DEAL. This research was supported by the DFG through grant no. 465469938 “Numerical diagnostics and improvements for the solution of linear dynamic macroeconomic models”. Declarations Competing Interests The authors have not disclosed any competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission (A2) vec (DPF(ΔP))= ([ (AP+B)−1C ]� ⊗ [ (AP+B)−1A ]) vec (ΔP ) (A3) =( P � ⊗ [ (AP+B) −1 A ]) vec (ΔP ) (A4) M(𝜆) ≡ A𝜆2+B𝜆+C=(A𝜆+AP+B)(𝜆−P) 641 Solving Linear DSGE Models withBernoulli Iterations directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/. References Adjemian, S., Bastani, H., Juillard, M., Mihoubi, F., Perendia, G., Ratto, M., & Villemot, S. (2011). Dynare: Reference manual, Version 4. Dynare Working Papers1. CEPREMAP. Al-Sadoon, M. M. (2018). The linear systems approach to linear rational expectations models. Econometric Theory, 34(3), 628–658. https:// doi. org/ 10. 1017/ S0266 46661 70001 60 Al-Sadoon, M. M. (2020). Regularized solutions to linear rational expectations models. Discussion paper. Anderson, G. S. (2010). A reliable and computationally efficient algorithm for imposing the saddle point property in dynamic models. Journal of Economic Dynamics and Control, 34(3), 472–489. https:// doi. org/ 10. 1016/ S14746670(17) 40506-4 Anderson, G.S., Levin, A., & Swanson, E. (2006). Higher-order perturbation solutions to dynamic discretetime rational expectations models. Discussion Paper 2006-01. Federal Reserve Bank of San Francisco Working Paper Series. Anderson, G. S., & Moore, G. (1985). A linear algebraic procedure for solving linear perfect foresight models. Economics Letters, 17(3), 247–252. https:// doi. org/ 10. 1016/ 01651765(85) 90211-3 Aruoba, S. B., Fernández-Villaverde, J., & Rubio-Ramírez, J. F. (2006). Comparing solution methods for dynamic equilibrium economies. Journal of Economic Dynamics and Control, 30(12), 2477–2508. https:// doi. org/ 10. 1016/j. jedc. 2005. 07. 008 Bai, Z.-Z., & Gao, Y.-H. (2007). Modified Bernoulli iteration methods for quadratic matrix equation. Journal of Computational Mathematics, 25(5), 498–511. Binder, M., & Pesaran, M. H. (1997). Multivariate linear rational expectations models: Characterization of the nature of the solutions and their fully recursive computation. Econometric Theory, 13(6), 877–88. https:// doi. org/ 10. 1017/ S0266 46660 00063 07 Binder, M., & Pesaran, M. H. (1999). Multivariate rational expectations models and macroeconometric modeling: A review and some new results. Macroeconomics. In M. H. Pesaran & M. Wickens (Eds.), Handbook of applied econometrics (3rd ed., Vol. 1, pp. 111–155). London: Wiley. https:// doi. org/ 10. 1111/b. 97806 31215 585. 1999. 00004.x Blanchard, O. J., & Kahn, C. M. (1980). The solution of linear difference models under rational expectations. Econometrica, 48(5), 1305–1311. https:// doi. org/ 10. 2307/ 19121 86 Caldara, D., Fernández-Villaverde, J., Rubio-Ramírez, J., & Yao, W. (2012). Computing DSGE models with recursive preferences and stochastic volatility. Review of Economic Dynamics, 15(2), 188–206. https:// doi. org/ 10. 1016/j. red. 2011. 10. 001 Chen, X. S., & Lv, P. (2018). On estimating the separation between (A, B) and (C, D) associated with the generalized Sylvester equation AXD −BXC =E . Journal of Computational and Applied Mathematics, 330, 128–140. https:// doi. org/ 10. 1016/j. cam. 2017. 07. 025 Dennis, J. E., Jr., Traub, J. F., & Weber, R. P. (1976). The algebraic theory of matrix polynomials. SIAM Journal on Numerical Analysis, 13(6), 831–845. https:// doi. org/ 10. 1137/ 071306 Dennis, J. E., Jr., Traub, J. F., & Weber, R. P. (1978). Algorithms for solvents of matrix polynomials. SIAM Journal on Numerical Analysis, 15(3), 523–533. https:// doi. org/ 10. 1137/ 07150 34 Gantmacher, F. R. (1959). The theory of matrices (Vol. I, II). Chelsea Publishing Company. Golub, G. H., & van Loan, C. F. (2013). Matrix computations (4th ed.). The Johns Hopkins University Press. Hammarling, S., Munro, C. J., & Tisseur, F. (2013). An algorithm for the complete solution of quadratic eigenvalue problems. ACM Transactions on Mathematical Software, 39(3), 18:1-18:19. https:// doi. org/ 10. 1145/ 24501 53. 24501 56 Higham, N. J. (2002). Accuracy and stability of numerical algorithms (2nd ed.). Society for Industrial and Applied Mathematics. Higham, N. J. (2008). Functions of matrices: Theory and computation. Society for Industrial and Applied Mathematics. Higham, N. J., & Kim, H.-M. (2000). Numerical analysis of a quadratic matrix equation. IMA Journal of Numerical Analysis, 20, 499–519. https:// doi. org/ 10. 1093/ imanum/ 20.4. 499 642 A.Meyer-Gohde Higham, N. J., & Kim, H.-M. (2001). Solving a quadratic matrix equation by Newton’s method with exact line searches. SIAM Journal on Matrix Analysis and Applications, 23(2), 499–519. https:// doi. org/ 10. 1137/ S0895 47989 93509 76 Horn, R. A., & Johnson, C. R. C. R. (2012). Matrix analysis (2nd ed.). Cambridge University Press. Huber, J., Meyer-Gohde, A., & Saecker, J. (2023). Solving linear DSGE models with structure-preserving doubling methods. IMFS Working Paper Series 195. Goethe University Frankfurt, Institute for Monetary and Financial Stability (IMFS). Judd, K. L. (1992). Projection methods for solving aggregate growth models. Journal of Economic Theory, 58(2), 410–452. https:// doi. org/ 10. 1016/ 00220531(92) 90061-L Judd, K. L. (1998). Numerical methods in economics. MIT Press. Klein, P. (2000). Using the generalized Schur form to solve a multivariate linear rational expectations model. Journal of Economic Dynamics and Control, 24(10), 1405–1423. https:// doi. org/ 10. 1016/ S01651889(99) 00045-7 Kågström, B. (1994). A perturbation analysis of the generalized Sylvester equation (AR −LB , DR −LE)=(C , F) . SIAM Journal on Matrix Analysis and Applications, 15(4), 1045–1060. https:// doi. org/ 10. 1137/ S0895 47989 32462 12 Kågström, B., & Poromaa, P. (1996). LAPACK-style algorithms and software for solving the generalized Sylvester equation and estimating the separation between regular matrix pairs. ACM Transactions on Mathematical Software (TOMS), 22(1), 78–103. https:// doi. org/ 10. 1145/ 225545. 225552 Lan, H., & Meyer-Gohde, A. (2014). Solvability of perturbation solutions in DSGE models. Journal of Economic Dynamics and Control, 45, 366–388. https:// doi. org/ 10. 1016/j. jedc. 2014. 06. 005 Lubik, T. A., & Schorfheide, F. (2003). Computing sunspot equilibria in linear rational expectations models. Journal of Economic Dynamics and Control, 28(2), 273–285. https:// doi. org/ 10. 1016/ S01651889(02) 00153-7 Meyer-Gohde, A. (2023). Numerical stability analysis of linear DSGE model—Backward errors, forward errors and condition numbers. IMFS Working Paper Series 193. Goethe University Frankfurt, Institute for Monetary and Financial Stability (IMFS). Meyer-Gohde, A., & Neuhoff, D. (2015). Solving and estimating linearized DSGE models with VARMA shock processes and filtered data. Economics Letters, 133, 89–91. https:// doi. org/ 10. 1016/j. econl et. 2015. 05. 024 Meyer-Gohde, A., & Saecker, J. (2024). Solving linear DSGE models with Newton methods. Economic Modelling, 133, 106670. https:// doi. org/ 10. 1016/j. econm od. 2024. 106670 Moler, C. B., & Stewart, G. W. (1973). An algorithm for generalized matrix eigenvalue problems. SIAM Journal on Numerical Analysis, 10(2), 241–256. https:// doi. org/ 10. 1137/ 07100 24 Polyak, B. (1964). Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5), 1–17. https:// doi. org/ 10. 1016/ 00415553(64) 90137-5 Rendahl, P. (2017). Linear time iteration. Discussion paper: IHS economics series. Sims, C. A. (2001). Solving linear rational expectations models. Computational Economics, 20(1–2), 1–20. https:// doi. org/ 10. 1023/A: 10205 17101 123 Smets, F., & Wouters, R. (2007). Shocks and frictions in US business cycles: A bayesian DSGE approach. The American Economic Review, 97(3), 586–606. https:// doi. org/ 10. 1257/ aer. 97.3. 586 Stewart, G. W. (1971). Error bounds for approximate invariant subspaces of closed linear operators. SIAM Journal on Numerical Analysis, 8(4), 796–808. https:// doi. org/ 10. 1137/ 07080 73 Taylor, J. B., & Uhlig, H. (1990). Solving nonlinear stochastic growth models: A comparison of alternative solution methods. Journal of Business & Economic Statistics, 8(1), 1–17. https:// doi. org/ 10. 1080/ 07350 015. 1990. 10509 766 Tisseur, F., & Meerbergen, K. (2001). The quadratic eigenvalue problem. SIAM Review, 43(2), 235–286. https:// doi. org/ 10. 1137/ S0036 14450 03819 88 Uhlig, H. (1999). A toolkit for analysing nonlinear dynamic stochastic models easily. In R. Marimon & A. Scott (Eds.), Computational methods for the study of dynamic economies (Vol. 3, pp. 30–61). Oxford University Press. Villemot, S. (2011). Solving rational expectations models at first order: What dynare does. Dynare Working Papers2. CEPREMAP. Wieland, V., Afanasyeva, E., Kuete, M., & Yoo, J. (2016). New methods for macro-financial model comparison and policy analysis. In J. B. Taylor & H. Uhlig (Eds.), Handbook of macroeconomics handbook of macroeconomics (Vol. 2, pp. 1241–1319). Elsevier. https:// doi. org/ 10. 1016/ bs. hesmac. 2016. 04. 004 643 Solving Linear DSGE Models withBernoulli Iterations Wieland, V., Cwik, T., Müller, G. J., Schmidt, S., & Wolters, M. (2012). A new comparative approach to macroeconomic modeling and policy analysis. Journal of Economic Behavior & Organization, 83(3), 523– 541. https:// doi. org/ 10. 1016/j. jebo. 2012. 01. 006 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.