Using the "Chandrasekhar Recursions" for Likelihood Evaluation of DSGE Models

Finance and Economics Discussion Series Divisions of Research & Statistics and Monetary Affairs Federal Reserve Board, Washington, D.C. Using the ”Chandrasekhar Recursions” for Likelihood Evaluation of DSGE Models Edward P. Herbst 2012-35 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.

Using the “Chandrasekhar Recursions” for Likelihood Evaluation of DSGE Models Ed Herbst Federal Reserve Board May 17, 2012 Abstract In likelihood-based estimation of linearized Dynamic Stochastic General Equilibrium (DSGE) models, the evaluation of the Kalman Filter dominates the running time of the entire algorithm. In this paper, we revisit a set of simple recursions known as the “Chandrasekhar Recursions” developed by Morf (1974) and Morf, Sidhu, and Kalaith (1974) for evaluating the likelihood of a Linear Gaussian State Space System. We show that DSGE models are ideally suited for the use of these recursions, which work best when the number of states is much greater than the number of observables. In several examples, we show that there are substantial benefits to using the recursions, withlikelihoodevaluationuptofivetimesfaster. Thisgainisespeciallypronouncedinlight ofthetrivialimplementationcosts–nomodelmodificationisrequired. Moreover,thealgorithmis complementarywithotherapproaches. JEL Classification: C18,C63,E20 Keywords: KalmanFilter,LikelihoodEstimation,ComputationalTechniques 1 Introduction Dynamic Stochastic General Equilibrium (DSGE) are increasingly estimated by Central Banks and academic economists. In estimation, the model equilibrium conditions can be linked to the data using a Linear Gaussian State Space (LGSS) representation. As model complexity increases, so too does computational time. In maximum likelihood and Bayesian contexts, the log likelihood has to be evaluated hundreds of thousands or millions of times, usually in sequential fashion. In this enviroment, likelihood computation time dominates the running time of the entire algorithm. So it becomes crucial to construct efficient algorithms for likelihood evaluation and state filtering. To wit, considerable effort Correspondence:BoardofGovernorsoftheFederalReserveSystem,20thStreetandConsitutionAvenueN.W.,Washington, D.C. 20551; edward.p.herbst@frb.gov. I thank, without implication, participants at the Research Computing Seminar at the Fed Board and John Roberts for comments. The views expressed in this paper are those of the author and do not necessarilyreflecttheviewsoftheFederalReserveBoardofGovernorsortheFederalReserveSystem. 1

has been expended constructing elaborate filters tailored to DSGE models–see, for example, Strid and Walentin(2009). The purpose of this paper is to report an old and simple algorithm for fast likelihood evaluation outlined in Morf (1974) and Morf, Sidhu, and Kalaith (1974) and show that it is ideally suited for DSGE models. The method, which we will call the “Chandrasekhar Recursions” (CR), is simple to implementandcanyieldconsiderablespeedimprovements. Thispaperiscloselyrelatedto StridandWalentin(2009), whodevelopa“BlockKalmanFilter”by exploitingtheaprioriknownstructureoftheDSGEmodeltoavoidsomelargematrixcalculations. The algorithmmustbeappliedonacase-by-casebasis. Wecomparethealgorithmsintheexamplesection. Moreover, in principle, one could use the CR and also exploit the block structure of the DSGE model. Finally, although the CR has not been used to aid the estimation of DSGE models, they have been employedinageneraltimeseriescontextbefore;see,forexample,Klein,Melard,andZahaf(1998). Thepaperisstructuredasfollows. Section2containsbackgroundinformationontheDSGEmodels and the Kalman Filter, Section 3 contains the derivation of the Chandrasekhar Recursions, Section 4 containsfourexamples,andSection5concludes. 2 The Kalman Filter and DSGE Models The starting point for our analysis is the LGSS representation of an economic model. The state vector s t isan n s ×1collectionofalltherelevantvariablesintheeconomicmodel,whileε t is nε ×1vectorof structuralshocksdrivingthemodel. Thestateequation1isderivedbyfirstlinearizingandthensolving themodel;seeAnandSchorfheide(2007)fordetails. s t = Ts t−1 +Rε t , ε t ∼N(0,Q) (1) y = D+Zs +η , η ∼N(0,H) (2) t t t t The vector y t is an n y ×1 vector of observables, and η t is an nη ×1 vector of measurement errors. This paper assumes that E[ε η(cid:48)] = 0 for ease of exposition. The objects (T,R,Q,D,Z,H) are matrix t t functionsofavectorθ ofstructureparameters. Wesuppressthedependenceforconvenience. Wemake twoassumptions,whichareveryoftenimposedinapplications. 1. Thesystemmatrices(T,R,Q,D,Z,H)aretime-invariant. 2. Theprocess{s }isstationary. Thismeansthatforourfilteringproblem,wecanusetheinvariant t distributiontoinitializeourstaterecursion. Thegoaloftheeconometricianistoevaluatetheloglikelihood p(Y|θ)whereY =[y (cid:48) ,...,y (cid:48)](cid:48) . Thisis 1 T accomplished by using the Kalman Filter. The Kalman Filter computes the log likelihood of the model usingthepredictiveerrordecomposition, (cid:76)(Y|θ)=− 1(cid:88) T (cid:128) n ln(2π)+ln(det(F ))+ν(cid:48) F −1ν (cid:138) (3) 2 y t t t t t=1 2

Where,theforecasterrorν andtheforecasterrorvariance F aregivenby, t t ν t = y t −D−Zˆs t|t−1 . (4) F t =ZP t|t−1 Z (cid:48)+H (5) Thequantitiesˆs t|t−1 and P t|t−1 arethepredictivestatemeanandvariance, ˆs t|t−1 =E t−1 [s t ], P t|t−1 =E t−1 [s t s (cid:48) t ]. Thepredictivemomentsofthestateequationevolveaccordingtotherecursions,1 ˆs t+1|t =Tˆs t|t−1 +K t F t −1ν t . (6) P t+1|t =TP t|t−1 T (cid:48)−K t F t −1K t (cid:48)+RQR (cid:48) . (7) Equation 7 is often called the matrix Riccati difference equation because of its resemblance to the univariate Riccati differential equation. The initial conditionals ˆs 1|0 and P 1|0 are, in principle, hyperparameters to be specified (or estimated) by the user, but this algorithm will use the unconditional distribution. Then ˆs 1|0 = 0 and P 1|0 will be the unconditional variance, i.e., the one that solves the discreteLyapunovEquation, TP 1|0 T (cid:48)−P 1|0 +RQR (cid:48)=0. Theformulaforthe K is t K t =TP t|t−1 Z (cid:48) . (8) This gain is usually written as K g,t = K t F t −1 = TP t|t−1 Z (cid:48) F t −1. Essentially it delivers the optimal extractionofinformationfromtheobservationsattime t. We are interested in computing the likelihood, which we could do by iterating (in suitable order) over the above equations. However, this filtering procedure slows down as the number of states, n , s grows. In particular, the performance is dominated by the updating of the state variance using Riccati equation(7). Inthatequation,wemustperformthematrixmultiplication, (cid:48) T P t|t−1 T . (cid:124)(cid:123)(cid:122)(cid:125) (cid:124)(cid:123)(cid:122)(cid:125) n s ×n s (cid:124) n (cid:123) × (cid:122) n (cid:125) n s ×n s s s ThisoperationisO(n3),dominatingallothermatrixoperationsinthefilter. StridandWalentin(2009) s (cid:48) reportthat,foralargeDSGEmodels,60%offilteringtimeisspentontheoperation TP t|t−1 T . 1Wecombinethestandardpredictionandupdatingequationsbecausewearetryingtoavoidunecessarycalculations. 3

3 The Chandrasekhar Recursions TheessenceoftheChandrasekharRecursionsistheavoidanceofdirectcomputationofP t|t−1 byinstead lookingatthedifferenceinthestatevariances, ∆P t|t−1 =P t|t−1 −P t−1|t−2 . ∆P t|t−1 is a real, symmetric, indefinite matrix.2 Suppose that rank(∆(P t ))=α. Then ∆P t|t−1 can be decomposedinto ∆P t|t−1 =W t−1 M t−1 W t (cid:48) −1 , where W is n ×α matrix and M is α×α matrix . There is room for improvement over the standard t s t algorithm if α<n , since the matrix multiplications will take place on lower dimensional matrices. In s general, however, this factorization is not unique. This seems like a difficult problem, because the factorization is not obvious, and, in our experience, even finding α consistently can be hard. Fortunately, theproblembecomesmucheasieronceitisknownthatthestatesarestationary. BeforediscussingtherecursionsforW and M ,it’sworthnotingthattheforecasterrorvariance F , t t t K t and K g,t canberewrittenasrecursionsdrivenby∆P t|t−1 . F t = F t−1 +Z∆P t|t−1 Z (cid:48) , (9) K t = K t−1 +T∆P t|t−1 Z (cid:48) , (10) K g,t = (K g,t−1 F t−1 +T∆P t|t−1 Z (cid:48))F t −1. (11) For verification of these recursions, see the Appendix. The heart of the algorithm lies in the recursions for M t and W t . Start with the initial state variance P 1|0 . Recall from Assumption 2 that the system is stationary,sowecanusetheunconditionalvarianceofthestates, E[s s (cid:48)]toinitializethestatevariance, t t P 1|0 =E t [s t s (cid:48) t ]. ThismatrixhasthepropertythatitsolvesthediscreteLyapunovequation, P 1|0 =TP 1|0 T (cid:48)+RQR (cid:48) . Consider now the next-period state variance forecast, P 2|1 . Using the Kalman Filter recursion (7), we canwritethisas, P 2|1 = TP 1|0 T (cid:48)+RQR (cid:48)−K 1 F 1 −1K 1 (cid:48) = P 1|0 −K 1 F 1 −1K 1 (cid:48) . (12) Thismeanswecanwritetheinitialdifferenceofstatevariancesas, ∆P 2|1 =−K 1 F 1 −1K 1 (cid:48) , 2Indeed,lookingatreducedrankdecompositionsof P t|t−1 ordifferencesthereof,hasawideliteratureknownasSquare RootKalmanFiltering.IpresenttheserelatedrecursionsbecauseIthinktheyaremucheasiertoimplement. 4

whichsuggestsnaturalchoicesfor M andW , 1 1 W t = K 1 =TP 1|0 Z (cid:48) , (13) M t = −F 1 −1=(ZP 1|0 Z (cid:48)+H)−1. (14) Note that thse equations already satistify the initial conditions, since F is positive definite. Moreover, 1 since K is n ×n and F is n ×n itiseasytoseethat 1 s y 1 y y rank(∆P )=α≤min{n ,n }. 2 y s Formanyeconomicmodels–inparticularlargescaleDSGEmodels–n <<n . Inthiscase,byusingthe y s Chandrasekhar recursions, very roughly speaking, the algorithm has “replaced” the matrix multiplication P t+1|t =TP t|t−1 T (cid:48)+...withW t M t W t (cid:48) ,whichinvolvesmatricesofmuchlowerdimension. Finally, we must derive the recursions for M and W . To do this, we utilize the following Lemma. t t Lemma. For∆P t =P t|t−1 −P t−1|t−2 , ∆P t+1|t = (T −K g,t Z)(∆P t|t−1 −∆P t|t−1 Z (cid:48) F t − − 1 1 ∆P t|t−1 Z (cid:48))(T −K g,t Z)(cid:48) (15) Proof. SeeAppendix. Thus,theRicatti-typedifferenceequationin P t+1 isreplacedbyanotherdifferenceequationintermsof ∆P t+1|t . In fact, this new difference equation gives the Chandrasekhar Recursion its name, because it evokesthe“so-calledChandrasekharalgorithmthroughwhichthematrix[RiccatiEquation]isreplaced by two coupled differential equations of lesser dimensionality [Aknouche and Hamdi (2007)].” Using theLemma,substituteourdecompositionW t−1 M t−1 W t (cid:48) −1 for∆P t|t−1 ,toobtain ∆P t+1|t =(T −K g,t Z)(W t−1 M t−1 W t (cid:48) −1 +W t−1 M t−1 W t (cid:48) −1 Z (cid:48) F t − − 1 1 Z (cid:48) W t−1 M t−1 W t (cid:48) −1 )(T −K g,t Z)(cid:48) . (16) RewritingwithW t−1 removedfromtheinnerproduct,wehave, ∆(P t+1 )=(T −K g,t Z)W t−1 (M t−1 +M t−1 W t (cid:48) −1 Z (cid:48) F t − − 1 1 Z (cid:48) W t−1 M t−1 )W t (cid:48) −1 (T −K g,t Z)(cid:48) . (17) Fromhereitiseasytoseethatwecanrewrite ∆P t+1|t =W t M t W t (cid:48) . (18) WhereW and M followtherecursions, t t W t =(T −K t F t −1Z)W t−1 , M t =M t−1 +M t−1 W t (cid:48) −1 Z (cid:48) F t − − 1 1 Z (cid:48) W t−1 M t−1 (19) Combiningtheequationsin19,withtherewrittenrecursionsfor F and M , t t 5

F t = F t−1 +ZW t−1 M t−1 W t (cid:48) −1 Z (cid:48) (20) K t = K t−1 +TW t−1 M t−1 W t (cid:48) −1 Z (cid:48) . (21) These equations, used in conjunction with with Equations 4 and 6, can be used to iteratively compute the log likelihood in equation 3. We have elimated the state variance prediction, P t|t−1 , from all the calculations in the algorithm and hence avoid the most computationally intensive calculation. Pseudo codeisreportedbelow. Initialization. 1. SolvetheLyapunovEquationfor P¯ 2. Set K =TP¯Z (cid:48) , F =ZP¯Z (cid:48)+H. 1 1 3. SetW =K , M =−F −1. 1 1 1 1 Iteration. Foreachtimeperiod t =1,...,T. 1. Computeν ,andevaluatethelikelihood. t 2. Compute aˆ t+1|t . 3. Compute F t+1 usingEquation 20,and F t − + 1 1 . 4. Compute K t+1 usingEquation 21, 5. ComputeW t+1 usingpartoneofEquation 19. 6. Compute M t+1 usingparttwoofEquation 19. Another advantage of this initialization is that it can shown that M will converge to a matrix as t t gets large. This is analgous to the steady state of the system expressed in usual form. With a general initialization, though, one cannot show that M will converge to anything. Finally, note that we can t recover P i|i−1 by i (cid:88) P i|i−1 =P 1|0 + ∆P j|j−1 , i>1. (22) j=1 3.1 Discussion ItisdifficulttocomputeanalyticallytheexactspeedgaingivenbytheChandrasekharRecursionsgiven the differences between highly optimized linear algebra routines across architectures. Still, one can perform a crude assessment of the differences in the algorithms without resorting to purely empirical studies. Welookedatallthematrixmultiplications,includingintermediatecalculations,intheKalman FilterandtheChandrasekharRecursions,takingcaretoavoidunnecessarycalculations,togaininsight intothedifferencesinthetwoalgorithms. Table 1 lists the remaining matrix multiplication operations after the “common” operations have been canceled out. The two algorithms appear almost the mirror image of one another, with n and s 6

n switched and a few additional operations for the Chandrasekhar recursions. Using naive linear y algebraic calculations, the running time of two additional distinct operations for the Kalman Filter is O(n3)andO(n2n ). FortheChandrasekharRecursions,therunningtimesareO(n2n )andO(n3). Itis s s y y s y clearthatifn >n ,theKalmanFilterwillhavegreateralgorithmiccomplexitythantheChandrasekhar s y recursions. If n <n ,thesituationwillbereversed. y s GiventhatDSGEmodelsfeaturemorestatesthanobservables,theChandrasekharrecursionsseema promisingalgorithmonthisbasis. However,thecalculationsinTable1arebasedon(1)acrudematrix multiplicationaccountingand(2)thenaivematrixmultiplicationalgorithmiccomplexity. Moreover,we haveabstractedfrommatrixadditionandtransposes. Weusefourexamplestogiveempiricalguidance ontherelativeperformance. 4 Four Examples Wecomparethealgorithmusingfourdifferentmodels: asmallRealBusinessCycleModel,theGeneric State Space model of Chib and Ramamurthy (2010), the Smets and Wouters (2007) model and the modelofSchmitt-GroheandUribe(2010). Foreachofthemodels,wecalculatethe“wall”(clock)time it takes to evaluate the likelihood at a particular point in the posterior 1000 times. We normalized these times, with the fastest algorithm being normalized to 100. This is a crude comparison, but it gives a sense of the actual user experience running the algorithms. We compare three algorithms, the standard Kalman Filter, the Block Kalman Filter of Strid and Walentin (2009), and the Chandrasekhar Recursions. WeimplementallthealgorithmsinIntelFortran11.1andMatlab2010a. WewrotecodeforthestandardKFandtheCRrecursions,andusedthecodeprovidedbyStridand Walentin (2009) for the Block Kalman Filter. This algorithm, specific for DSGE models, uses a priori knowledgeofthestructureof T –ithaslargematrixofzeroswheretheexogenousvariablesloadonto the endogenous variables – and Z, which is often quite sparse, to build a fast Kalman filter. It requires the user to prestructure the model in a particular way and apply it on a case-by-case basis. We did not benchmarkthisfortheGenericStateSpacemodel,sinceitisnotaDSGEmodel. TheFortrancodeutilizesIntel’sMathKernelLibraryimplementationofBLASandLAPACKfastlinear algebra routines, including dsymm, symmetric matrix multiplication. The Matlab code uses a standard distributionofBLASanddoesnotconsidersymmetricmatrixmultiplication,whichmightdisadvantage itsomewhat. BothprogramsutilizemultithreadedBLAS,usingfourthreadsformatrixoperations.3 All calculations were on a 12-core Intel Xeon x5670 CPU (2.93GHz) with L1, L2, and L3 cache sizes of 32K,256K,12288K,respectively. One other technical detail worth mentioning is method of computing the solution to the discrete Lyapunov Equation. For the Chandrasekhar Recursions, it is crucial that this solution be (at least approximately) correct. Repeating trials suggests that for large (n greater than 100), the Matlab s routine dlyap does not provide a good solution. Instead, the Matlab implementation of the CR uses lyap_symm, distributed as part of Dynare [Adjemian, Bastani, Juillard, Mihoubi, Ratto, and Villemot 3Theresultswerebroadlythesamewhenonlyserialcomputationwasconsidered. 7

(2011)],whichyieldsamuchbettersolution. The equilibriums for the RBC and SW models are computed with Sims (2002) GENSYS. GENSYS is widely used to compute equilibriums of Linear Rational Expectations Models in economics. The algorithm uses the Generalized (complex) Schur decomposition, which gives it the advantage that “controls”and“states”don’thavetobespecifiedaheadoftime;i.e.,redundantstatescanbeincluded. This means that for our examples, it is possible to reduce the size of the state vector and thus the efficacy of the CR. Given that our examples are of the small and medium size, we think that they are illustrativeofspeedgainsforlargermodels. 4.1 Real Business Cycle Model The first example is a simple Real Business Cycle model. There are n = 2 observables, labor and y output. In the GENSYS formulation of the model, there are n =12 states. The model is driven by two s shocks,neutraltechnologyanddemand. Table 2 reports the timings associated with evaluating the likelihood 1000 times, with the fastest normalized to 100. The language gain associated with Fortran is substantial, with the likelihood evaluation between four and 8 times faster than their Matlab counterparts. Within languages, that the Chandrasekharrecursionsarethefastestinbothlanguages. While,theBlockKFoutperformsthestandard KF in Matlab, that is not the case in Fortran. The slow performance in Fortran is consistent with the findings of Strid and Walentin (2009), who find that the additional overhead associated with the increasednumberofmatrixmultiplicationsoutweighsthesizegainforsmallmodels. 4.2 Generic State Space Model The next example is the Generic State Space model used in Chib and Ramamurthy (2010). This is not aDSGEmodel. Ithasnomeaningfuleconomicinterpretation. Itis,however,agoodexample,because there are fewer states than observables. In the model n = 5 while n = 10. Stochastic singularity is s y avoidedbymeasurementerror. Thedatausedissimulated,withlength200periods. Table 3 shows the timings associated with evaluating the likelihood 1000 times, with the fastest normalizedto100.4 Thelanguagegainisquitesubstantial,withFortranbeingroughlyfourtimesfaster for both algorithms. This large difference is mostly driven by the length of sample size. Concentrating on within-language performance, we see that there is a slight drop in speed when using the CR in Fortran. Relative to the standard KF, the CR about 17% percent slower. On the other hand, they are about 10% faster in Matlab. In both cases, the speed difference is small compared to the other examples. 4.3 Smets-Wouters Model TheSmetsandWouters(2007)isamedium-scalemacroeconomicmodelthathasbecomeabenchmark formanyissuesrelatedtotheestimationandidentificationofDSGEmodels. Inthisformulationofthe 4WeagainnotethatwedonotusetheBlockKFbecauseitisnotaDSGEmodel. 8

model, there are n =50 states and only n =7 observables.5 This example includes some redundant s y states, I but it is consistent with the conventional way of writing and estimating such models in economics. 50 states are not unusual for a medium-scale DSGE model, especially the kind used in central banks. Forthisreason,thisexampleservesasagoodbenchmark. Table 4 reports the timings associated with evaluating the likelihood 1000 times, with the fastest normalized to 100. It quite apparent that the language gain associated with Fortran relative Matlab is large,withtherunningtimeapproximatelydoubledforeachalgorithminMatlabcomparedtoFortran. Withinalanguage,weseethattheChandrasekharRecursionsdominateboththeBlockandthestandard KF. For Matlab, the CR offer algorithmic speed gains of 61% and 72%, for the Block KF and Standard KFrespectively. (cid:48) This speed gain is accomplished by eliminating of the matrix operation TP t|t−1 T . Indeed, inspection of the Matlab Profiler indicates that this operation is about 45% of all filtering time for the StandardKF.FortheBlockKF,about35%ofthefilteringtimeisspentonsimilar,slightlysmallermatrix multiplicationcorrespondingtothenon-exogenoussub-statevector(whichhas n −n =43states.) s y 4.4 Schmitt-Grohe and Uribe (2010) Model The final model comes from Schmitt-Grohe and Uribe (2010). The paper estimates a DSGE model augmentedwith“news”shocks. Specifically,theyconstructarealbusinesscyclemodelwithinvestment adjustment consts, capacity utilization, habits, and Jaimovich and Rebelo (2009) preferences. Each of the seven structural shocks is driven by three innovations. One of these innovations is unanticipated, while the other two are anticipated four and quarters ahead, respectively. The process for a given exogenousshock,z ,lookslike: t z t =ρ z z t−1 +εz t ,0+εz t− ,4 4 +εz t− ,8 8 . Writing this process recursively requires at least an additional eight states. In total, there are n = 98 s statesinthemodelandn =7observables,usingquarterlydatafrom1965-2006,yielding207obsery vations. Moreover,giventhenewsstructure,thismodelisnoteasilyconvertedintoaformusedbyStrid and Walentin. The comparison is thus restricted to the standard Kalman Filter and the Chandrasekhar Recursions. Table 4 reports the timings associated with evaluating the likelihood 1000 times, with the fastest normalized to 100. Once again, the gain from using Fortran is substantial, regardless of the algorithm used. Once again, the Chandrasekhar Recursions dominate the Kalman Filter irrespective of language. Much like the Smets-Wouters model, much of the gain comes from eliminating the matrix operation TP t|t−1 T. Overall, the CR posts algorithmic gains of about 38% and 68% for Matlab and Fortran, respectively. 5Wesetthemovingaveragecoefficientsonthemarkupshocksandautomaticstabilizeronthegovernmentspendingshock tozero.Thiswaywecanusethetwoblock(AR),sparseBlockKalmanFilter. 9

5 Conclusion ForDSGEmodels,thebiggestbottleneckincomputationistheevaluationoftheloglikelihood. Within this evaluation, the prediction of the state variance is the slowest operation, taking about half of all filtering time for large models. This paper has presented the Chandrasekhar recursions, a simple, fast algorithm for evaluating the likelihood of a Linear Gaussian State Space System. The CR algorithm works by eliminating the predicted state variance from the Kalman Filter equations altogether. In this way, it is ideally suited for DSGE models. The price paid to use this algorithm is relatively small. The systemmustbestationaryandtimeinvariant,assumptionsthataretypicallysatisfiedforDSGEmodels. Itshouldbenotedthatthismethodisentirelyconsistentwithotherfastfilteringtechniques,suchas DurbinandKoopman(2000). Theseadditionalspeedgainsarelargelyorthogonaltotheonespresented here. Given the ease of implementation and apparent speed gains, the Chandrasekhar Recursions shouldbecomepartofappliedmacroeconomists’computationaltoolkit. 10

References ADJEMIAN,S.,H.BASTANI,M.JUILLARD,F.MIHOUBI,M.RATTO, AND S.VILLEMOT (2011): “Dynare: ReferenceManual,Version4,”DynareWorkingPapers1,CEPREMAP. AKNOUCHE,A., AND F.HAMDI (2007): “PeriodicChandrasekharrecursions,”arXiv:0711.385v1. AN,S., AND F.SCHORFHEIDE (2007): “BayesianAnalysisofDSGEModels,”EconometricReviews,26(2-4), 113–172. CHIB,S., ANDS.RAMAMURTHY(2010): “TailoredRandomizedBlockMCMCMethodswithApplicationto DSGEModels,”JournalofEconometrics,155(1),19–38. DURBIN, J., AND S. KOOPMAN (2000): “Fast Filtering and Smoothing For Multivariate State Space Models,”JournalofTimeSeriesAnalysis,21,281–296. JAIMOVICH,N.,ANDS.REBELO(2009): “CanNewsabouttheFutureDrivetheBusinessCycle?,”American EconomicReview,9(4),1097–1118. KLEIN,A.,G.MELARD, AND T.ZAHAF (1998): “ComputationoftheExactInformationMatrixofGaussian DynamicRegressionTimeSeriesModels,”TheAnnalsofStatistics,26,1636–1650. MORF,M.(1974): “FastAlgorithmsforMultivariateSystems,”Ph.D.thesis,StanfordUniversity. MORF,M.,G.SIDHU,ANDT.KALAITH(1974): “Somenewalgorithmsforrecursiveestimationinconstant, linear,discrete-timesystems,”IEEETransactionsonAutomaticControl,19,315–323. SCHMITT-GROHE,S., AND M.URIBE (2010): “What’sNewsinBusinessCycles?,”WorkingPaper. SIMS,C.A.(2002): “SolvingLinearRationalExpectationsModels,”ComputationalEconomics,20,1–20. SMETS, F., AND R. WOUTERS (2007): “Shocks and Frictions in US Business Cycles: A Bayesian DSGE Approach,”AmericanEconomicReview,97,586–608. STRID, I., AND K. WALENTIN (2009): “Block Kalman Filtering for Large-Scale DSGE Models,” ComputationalEconomics,33,277–304. 11

6 Tables Table1: NoncommonMatrixMultiplications Algorithm StandardKF (2x)(n ×n )(n ×n ),(n ×n )(n ×n ) s s s s s s s y ChandrasekharRecursions (3x)(n ×n )(n ×n ),(2x)(n ×n )(n ×n ) s y y y y y y y Table2: RBCModel: WallTime,FastestNormalizedtoOneHundred Method Language WallTime StandardKF Matlab 982 BlockKF Matlab 616 ChandrasekharRecursions Matlab 410 StandardKF Fortran 123 BlockKF Fortran 151 ChandrasekharRecursions Fortran 100 Table3: GenericStateSpaceModel: WallTime,FastestNormalizedtoOneHundred. Method Language WallTime StandardKF Matlab 499 ChandrasekharRecursions Matlab 451 StandardKF Fortran 100 ChandrasekharRecursions Fortran 120 12

Table4: Smets-WoutersModel: WallTime,FastestNormalizedtoOneHundred Method Language WallTime StandardKF Matlab 643 BlockKF Matlab 451 ChandrasekharRecursions Matlab 174 StandardKF Fortran 249 BlockKF Fortran 214 ChandrasekharRecursions Fortran 100 Table5: NewsModel: WallTime,FastestNormalizedtoOneHundred Method Language WallTime StandardKF Matlab 1263 ChandrasekharRecursions Matlab 787 StandardKF Fortran 309 ChandrasekharRecursions Fortran 100 7 Appendix 7.1 Verification of recursions for F ,K ,and K . t t g,t For F : t F t = ZP t|t−1 Z (cid:48)+H = ZP t|t−1 Z (cid:48)+H+F t−1 −F t−1 = F t−1 +ZP t|t−1 Z (cid:48)−ZP t−1|t−2 Z (cid:48) = F t−1 +Z∆P t|t−1 Z (cid:48) . K and K aresimilar. t g,t 7.2 Verification of the difference equation for ∆P . t|t−1 To show that Chandrasekhar recursions work for the special case discussed above, we show how to writetherecursionsfor F and K intermsof∆(P ). t t t Proof. Using the definition of F t and adding and subtracting F t−1 , A similar algebraic manipulation 13

canbeusedfor K . Notethatwecanalsowrite K as t g,t K g,t =(K g,t−1 F t−1 +T∆(P t )Z (cid:48))F t −1. (23) Proof of Lemma. FromtheKalmanFilter,wehavethat P t+1|t =TP t|t−1 T (cid:48)+RQR−K g,t F t K g (cid:48) ,t . Subtracting P t|t−1 frombothsides,wehavethat ∆P t+1|t =T∆P t|t−1 T (cid:48)−K g,t F t K g (cid:48) ,t +K g,t−1 F t−1 K g,t−1 . Usingtherecursionsfor F showninthepreviouslemma,wehave t ∆P t+1|t =T∆P t|t−1 T (cid:48)−K g,t (F t−1 +Z∆P t−1|t Z (cid:48))K g (cid:48) ,t +K g,t−1 F t−1 K g,t−1 . Groupingliketermsandcompletingthesquarefor(T −K Z),wehave g,t ∆P t+1|t = (T −Kg,tZ)∆P t|t−1 (T −K g,t Z)(cid:48)+K g,t Z∆P t|t−1 T (cid:48)+T∆P t|t−1 Z (cid:48) K g (cid:48) ,t − 2K g,t ∆P t|t−1 K g (cid:48) ,t −K g,t F t−1 K g (cid:48) ,t +K g,t−1 F t−1 K g (cid:48) ,t−1 . (24) Notethatwecanwritethefinalproductintheequation,using23,20,andtediouslyexpandingterms, as, K g,t−1 F t−1 K g,t−1 = (K g,t F t −T∆P t|t−1 Z (cid:48))F t −1(K g,t F t −T∆P t|t−1 Z (cid:48)) = (K g,t (F t−1 +Z∆P t|t−1 Z (cid:48))−T∆P t|t−1 Z (cid:48))F t −1(K g,t (F t−1 +Z∆P t|t−1 Z (cid:48))T∆P t|t−1 Z (cid:48)) = (K g,t F t−1 +(K g,t Z−T)∆P t|t−1 Z (cid:48))F t − − 1 1 (K g,t F t−1 +(K g,t Z−T)∆P t|t−1 Z (cid:48)) = (T −K g,t Z)∆P t|t−1 Z (cid:48) F t − − 1 1 Z∆P t|t−1 (T −K g,t Z)+K g,t F t−1 K g,t + K g,t Z∆P t|t−1 (K g,t Z−T)(cid:48)+(K g,t Z−T)∆P t|t−1 Z (cid:48) K g (cid:48) ,t = (T −K g,t Z)∆P t|t−1 Z (cid:48) F t − − 1 1 Z∆P t|t−1 (T −K g,t Z)+K g,t F t−1 K g,t − K g,t Z∆P t|t−1 T (cid:48)−T∆P t|t−1 Z (cid:48) K g (cid:48) ,t +2K g,t ∆P t|t−1 K g (cid:48) ,t . (25) Combining24and25,wehaveverifiedthelemma. 14