Retrieving matrices before solving model

Hello all,
How can I retrieve the matrices before Dynare solves a non-linear model, i.e, matrices f_{y+}, f_{y0} and f_{y-} of equation (4) in Villemot (2011), “Solving rational expectations models at first order:
what Dynare does”?

Dear Eduardo,

It is possible to retrieve these matrices, but you need to dig a bit into Dynare internals.

If your mod file is foo.mod, the matrices that you are looking for are generated by the file foo_dynamic.m, which is itself generated by the Dynare preprocessor. You will also need some information located in foo.m, which is also generated.

More precisely, the matrices can be retrieved from the second output argument (g1) of foo_dynamic.m, since it is the Jacobian of the whole model. You can retrieve g1 by calling the foo_dynamic function, giving as arguments the (column-)vector of endogenous (y), a row-vector for the exogenous (x), the vector of model parameters (params), the steady state (steady_state) and with it_=1.

The tricky part is that since we are considering the dynamic model, we need to disentangle the backwards part (f_{y-}), the contemporaneous one (f_{y0}), and the forward one (f_{y+}). The columns of g1 correspond to columns extracted from these 3 matrices.

To understand the correspondence, you need to look at matrix M_.lead_lag_incidence (in foo.m): this matrix has three columns and as many lines as there are endogenous variables. The first column corresponds to variables that appear with a lag, the second to contemporaneous variables, and the third to variables with a lead. The integers stored in this matrix correspond to columns of g1 (and incidentally also to indices in the y and steady_state vectors). For example if M_.lead_lag_incidence(5,1)=4, then this means that the fifth endogenous of your model (in the order of declaration), at the previous period, corresponds to column 4 of g1 (and to index 4 of y and steady_state).

Hope this helps,

Dear Sebastien
Thank you for the careful answer. I will work on it over the weekend and get back to you.
Best regards,

Hello Sebastien,
I was able to locate matrix g1 and interpret its lay-out. Now my challenge is to set up a code that rearranges it in the AB format of Blanchard-Kahn (with auxiliary variables, if necessary) to be able to apply the formulas that I developed for regime-switching. If you already know how to put g1 in the AB format, let me know.
Many thanks for your help.

Dear Eduardo,

I’m sorry but I don’t know what you mean by the “AB format”. You should be more specific.


Hello Sebastien,
I was able to do it on my own but I still have one pending question (elaborated below). The AB format of BK (I got this expression from Johannes Pfeiffer) is a linear(ized) DSGE represented as:


where x_{t+1} is a vector containing the predetermined variables as well as the expectation as of t of the jumpers as of t+1. I understand it is a transformation of the Jacobian g1 computed by Dynare:


that requires shifting equations back and creating more auxiliary variables. I was able to code Matlab to transform g1 into the AB format. The solution of the AB format yields a state-space model:


where x^p and x^j are the predetermined and jumper variables, respectively.
The pending question is the derivation of the formula of the matrix R (which pre-multiplies the shock vector) with the Blanchard-Kahn method. I would appreciate if you could help deriving the formula for this matrix R so that my IRFs and Dynare’s match.
The attachment “Question.pdf” elaborates the question in detail (with different notation, based on example 3.4 of Barthelemy and Marx, Solving RE models in: Oxford Univ Press, Handbook of Computational Ec.)
I look forward to your reply.
Eduardo.Question.pdf (462.0 KB)

I don’t understand what you are now trying to do. Before, you were working on the A and B matrix of the linearized equilibrium system. But now you seem to be working on the solution to this system. The matrix R shows up in your state transition equation.

Hi Johannes, thanks for joining. Sorry for the lack of clarity. Indeed, my first step was to collect the Jacobian of the system and convert it to the BK representation. This has been done. My second step is to solve the linearized system (reshaped as BK) with my own code. That’s where my stumbling block is now: finding the formula for shock-identifying matrix in the solution.

I need to solve the model with the BK method outside Dynare with my own code to be able to simulate a regime-switching process (more general than Markov); and then to replicate the simulation numerous times to search for optimal reaction function parameters (non-quadratic objective function). I’d be happy to post the paper here (work in progress) so that you have a complete view.

Regarding the matrix whose formula I’m looking for, did you have a look at the file Question.pdf attached to my previous post here? It elaborates the question (with different notation from my previous post) and is based on Barthelemy and Marx (2018), Solving RE models in: Oxford Univ Press, Handbook of Computational Ec., example 3.4.

You answered that R (the matrix whose formula I’m looking for, labeled \tilde{B} in the attachment) shows up in the state transition equation. Did you mean it shows up in the state transition equation of linearized unsolved system? Perhaphs referring to the attachment would make our dialogue easier.

I see. My problem is already the first line of your PDF. Generally, there must be a conformable matrix premultiplying the \varepsilon_{t+1}. Otherwise, you could only have as many shocks as variables and only one shock entering each equation.

Indeed, the first line of the attachment Question.pdf in a previous post is the special case of one shock per equation. In this new attachment Duvidas.pdf (78.2 KB)
, I present the general case. I derived the formula for \tilde{B} with the Blanchard-Kahn method: I plugged the solution for the jumpers twice in the upper rows of the system (the ones relative to predetermined variables). Can you let me know what is wrong with it that is preventing the IRF from matching Dynare’s?

For now I can only offer you my version of the codes doing this plus a Dynare mod-file that does the same:
MOUCandBK.m (3.5 KB)
myRBCmodel.mod (5.2 KB)
Chapter2_RBC.pdf (425.0 KB)

Thanks, Johannes.
I’m looking for the solution formula for the shock-identifying matrix for a general case. The example in the attachments you kindly posted (matrix \eta) applies only for the particular case where there is a single shock in an AR(1) equation (and therefore the solution is trivial as this original equation is part of the solution system).

I understand that my question is not specific to Dynare but to the BK method. But I would appreciate any further help that you or any reader could provide. My question is outlined in the attachment Duvidas.pdf in my more recent previous post. I have already searched many textbooks to no avail.

Have a look at
Dejong Dave 2011 - Structural Macroeconometrics 2nd Edition, Blanchard Kahn.pdf (48.4 KB)

Thanks for this, Johannes. I’m afraid Dejong&Dave only adds to the problem. As you can see with these codes, Dejong and Dave’s solution (equation 4.47) does not match Dynare’s. Duvidas.m (1.6 KB)
Toy2.mod (263 Bytes)
I’m still looking for the formula for shock-identifying matrix in a DSGE solution’s state equation.

Im looking forward to any suggestions on how to derive the shocks matrix with the Blanchard-Kahn method in a way that matches Dynare’s solution.

The solution consists of (4.46) and (4.47). Did you use 4.46 with \rho=0 as our shock is iid? And do you understand how (4.45) implies (4.46)? The D_2 is not raised to the power of i, but is nevertheless in the sum. That seems like another typo to add to the Errata list at In any case, as far as I can see, we have


to get the form in (4.37). That implies that


in your notation or
from (4.40).
Now we have from (4.46)
or in your code


which returns the correct result.


% Let A\B=P_inv*Lambda*P, ordered in increasing eigenvalues' modulus
[eVec,eVal]=eig(A\B); %Diagonalizing
P_inv= [eVec(:,2) eVec(:,1)]; %just reordering the eigenvectors
P=inv(P_inv); % finding P (the inverse of the eigenvector matrix)
Lambda=[ eVal(2,2) 0; 0   eVal(1,1)]; %just reordering the eigenvalues
H=-inv(P(2,2))*P(2,1) %signal equation coefficient
F=(P(2,2)-P(1,2)*H)\Lambda(1,1)*(P(2,2)-P(1,2)*H) %state equation coeff.


Thanks Johannes. Now I understand that current period shocks matter for the jumpers also directly, not only through their impact on predetermined variables. My mistake was that I was leaving out the term with the current period shocks from the jumpers’ formula.
Many thanks,

No problem, but can you confirm my suspicion that the formula in Dejong/Dave has a mistake?

I will check thoroughly the formula on Dejong/Dave for typos this weekend and get back to you, Johannes. Many thanks,

Indeed, Johannes, formula 4.46 of Dejong/Dave has a mistake. The term inside the parenthesis is not even square to be inverted. To correct, D_2 should be outside the parenthesis.
Sorry for the late response; Brazil moves slowly during World Cups until resolution with the Seleção.