Using JSON Output from the Dynare Preprocessor

We have recently added an option to produce JSON output from the Dynare Preprocessor. If you’re new to Dynare you should know that the Preprocessor is the part of Dynare that transforms your .mod file into a file usable by MATLAB, Octave, Julia, or the C compiler. Providing JSON output allows the preprocessor to communicate everything it knows about the model (e.g. the model equations, variables, static and dynamic derivatives, etc.) in a way that is easily parsed by many programming languages. This makes it possible to use the Dynare Modeling Language in any programming environment that can parse JSON.

In this post, I’d like to walk you through an example of using the JSON output of the Dynare Preprocessor. We will write a routine that parses the JSON output and estimates the parameters, equation by equation, via Ordinary Least Squares. We will then use this routine to estimate the Taylor rule parameters from Smets and Wouters (2007). However, before getting to the example, I’d like to briefly give you some background on the Dynare Preprocessor and the JSON output it produces.

On a final, practical, note, you should know that the OLS routine and the modified .mod file described herein work with Dynare 4.6 and later .

The Dynare Preprocessor

At the basic level, the Dynare Preprocessor takes as input a Dynare .mod file and outputs the derivatives of the static and dynamic versions of the model and their derivatives in addition to a driver file that guides the back-end actions to be taken. These outputs are provided for use with MATLAB, Octave, C (without the driver), and, as of Dynare 4.6, Julia.

In addition to the aforementioned outputs, Dynare 4.6 provides JSON output that represents the .mod file at every major preprocessing stage: Parsing, Check Pass, Transform Pass, and Computing Pass. To better understand the type of JSON output that can be obtained, it is helpful to see the Dynare Preprocessor Flow Chart and understand, in a general sense, what is done at each stage:

Dynare Preprocessor Flow Chart

As you can see from the flow chart above, there are 6 preprocessing stages:

  1. Macroprocessor: the Dynare Macroprocessing language is used to perform textual manipulations of the .mod file. The output from this stage is a .mod file that is ready to be parsed. You can read more about the Dynare Macroprocessing language here.
  2. Parsing: takes a potentially macro-expanded .mod file and parses it into an Abstract Syntax Tree (AST), comprising the internal representation of the .mod file. In doing so, among other cursory checks, it verifies that the .mod file has valid Dynare syntax, commands, and options.
  3. Check Pass: verifies the coherence of the .mod file. For example, this is where we ensure that the number of declared endogenous variables equals the number of equations in the model block.
  4. Transform Pass: among other transformations, adds auxiliary variables and equations for leaded and lagged variables, thereby transforming the model into t-1, t, t+1 form.
  5. Computing Pass: calculates the derivatives of the transformed static and dynamic models using the symbolic derivative engine.
  6. Write Output: writes MATLAB, Octave, C, Julia, or JSON files.

To read more about the Dynare Preprocessor, see this presentation.

More on JSON

JSON is a data interchange format that is easily understood by humans and easily parsed by many programming languages. In short, it associates keys with values like a dictionary. In JSON, keys are strings whereas values can be strings, numbers, arrays, objects, boolean, or null.

The easiest way to get a sense of what a JSON file looks like is to see it. This declaration of parameters in a .mod file

parameters beta $\beta$ (long_name='discount factor'), rho;

would produce the following lines in JSON

"parameters": [{"name":"beta", "texName":"\\beta", "longName":"discount factor"}
             , {"name":"rho", "texName":"rho", "longName":"rho"}]

This tells us that key "parameters" is associated with an array (enclosed by brackets) of objects (enclosed by braces). The array has two entries. The first entry in this array is an object where the key "name" is associated with the string "beta", the key "texName" is associated with the string "\\beta", and the key "longName" is associated with the string "discount factor". The second entry has similar keys but, for the case of rho, no specific \(\LaTeX\) name or long name was declared, so those keys take the default values. As you can see, understanding the contents of a JSON file and seeing how they correspond to the originating .mod file is straightforward. For more details on JSON visit https://www.json.org.

A JSON representation of the .mod file can be obtained after the Parsing, Check Pass, Transform Pass, and Computing Pass stages outlined above. To obtain JSON output from the Dynare Preprocessor, you must choose where you want that output to be produced by passing the command line option json=parse|check|transform|compute. Note that the output provided varies a bit, depending on where you want that output produced. For example, the JSON representation of the derivatives of the dynamic and static models will only be produced after the derivatives of the model have been calculated in the Computing Pass. Again, the details of what is produced after every stage are outlined in the Dynare manual.

An Example of Putting the JSON output to use: Ordinary Least Squares

As an example application of how one can use the Dynare JSON output, I will run OLS on the Taylor rule in Smets and Wouters (2007).

The original .mod file, Smets_Wouters_2007.mod, and data file, usmodel_data.mat, were downloaded from Johannes Pfeifer’s DSGE_mod repository.

Below, I show the .mod file and describe the necessary modifications to run OLS. After that, I describe the construction of the MATLAB routine that uses the Dynare Preprocessor JSON output to run OLS (this routine is general and would work with the JSON output provided for any .mod file). Finally, I run OLS on the monetary policy rule.

The .mod file

The following are the parts of Smets_Wouters_2007.mod that I modified for this post. The entire file used can be found here.

First Modification

@@ -1,3 +1,6 @@
+// --+ options: json=compute +--

The first line of the file tells the Dynare Preprocessor to produce JSON output after the Computing Pass. This creates the following files in the Smets_Wouters_2007/model/json directory:

  • modfile.json,
  • modfile-original.json,
  • static.json, and
  • dynamic.json.

The first file, modfile.json, is the equivalent of the standard .m driver file output by the Dynare Preprocessor only in JSON format. It contains lists of model variables, the model block (transformed into t-1, t, t+1 format), a list of Dynare statements, the list of equation cross references, and some general information about the model.

The second file, modfile-original.json contains a slightly modified version of the model as written in the model block. It contains no auxiliary variables or auxiliary equations, but it does expand adl nodes, if there are any. Here is what the Taylor rule looks like in JSON format:

{
"model":
[
  ...
  {
    "lhs": "r",
    "rhs": "pinf*crpiMcrpiXcrr+cryMcryXcrr*ygap+crdy*diff(ygap)+crr*r(-1)+ms",
    "line": 141,
    "tags": {
              "name": "taylor_rule"
            }
  },
  ...
]
}

This is the file of interest for the OLS routine as we want to maintain the lag information contained in the model block. This file is written when either json=compute or json=transform is passed as an option to the dynare command.

The final two files, static.json and dynamic.json, contain the derivatives of the dynamic and static models. These files are a byproduct of using json=compute. Our OLS routine doesn’t need them.

Second Modification

@@ -1,3 +1,6 @@
 // --+ options: json=compute +--
+path(['..' filesep 'ols'], path);

The second change I made to the original .mod file was to add the relative path to the ols folder containing the OLS routine and its helper functions. Since this routine is not part of the official Dynare release, Dynare will not actually take care of adding its folder to the MATLAB/Octave path. Hence, this line.

Third Modification

@@ -135,12 +138,13 @@ model(linear);
     y = cfc*( calfa*k+(1-calfa)*lab +a );
     pinf = (1/(1+cbetabar*cgamma*cindp)) * ( cbetabar*cgamma*pinf(1) +cindp*pinf(-1)+((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc-1)*curvp+1)*(mc)  )  + spinf;
     w = (1/(1+cbetabar*cgamma))*w(-1)+(cbetabar*cgamma/(1+cbetabar*cgamma))*w(1)+(cindw/(1+cbetabar*cgamma))*pinf(-1)-(1+cbetabar*cgamma*cindw)/(1+cbetabar*cgamma)*pinf+(cbetabar*cgamma)/(1+cbetabar*cgamma)*pinf(1)+(1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))*(csigl*lab + (1/(1-chabb/cgamma))*c - ((chabb/cgamma)/(1-chabb/cgamma))*c(-1) -w)+ 1*sw;
+    [name='taylor_rule']

On line 141, I add an equation tag to the monetary policy rule; this will be used to tell the OLS routine which equation to estimate. I make sure not to include any spaces in the equation tag name as the output of the OLS routine is stored in a substructure of oo_.ols with the name of the equation tag. Hence, here, the output will be stored in oo_.ols.taylor_rule.

Fourth Modification

@@ -34,15 +37,15 @@
  */

 var labobs robs pinfobs dy dc dinve dw ewma epinfma zcapf rkf kf pkf cf
-    invef yf labf wf rrf mc zcap rk k pk c inve y lab pinf w r a b g qs ms
-    spinf sw kpf kp;
+    invef yf labf wf rrf mc zcap rk k pk c inve y lab pinf w r a b g qs
+    spinf sw kpf kp ygap;

-varexo ea eb eg eqs em epinf ew;
+varexo ea eb eg eqs ms epinf ew;

 parameters curvw cgy curvp constelab constepinf constebeta cmaw cmap calfa
            czcap csadjcost ctou csigma chabb ccs cinvs cfc
            cindw cprobw cindp cprobp csigl clandaw
-           crdpi crpi crdy cry crr
+           crdpi crdy crr crpiMcrpiXcrr cryMcryXcrr
            crhoa crhoas crhob crhog crhols crhoqs crhoms crhopinf crhow
            ctrend cg;

@@ -68,10 +71,10 @@ cprobp=   0.6;
 cindw=    0.3243;
 cindp=    0.47;
 czcap=    0.2696;
-crpi=     1.488;
 crr=      0.8762;
-cry=      0.0593;
 crdy=     0.2347;
+crpiMcrpiXcrr = 0.1842;
+cryMcryXcrr   = 0.0073;

@@ -135,12 +138,13 @@ model(linear);
     y = cfc*( calfa*k+(1-calfa)*lab +a );
     pinf = (1/(1+cbetabar*cgamma*cindp)) * ( cbetabar*cgamma*pinf(1) +cindp*pinf(-1)+((1-cprobp)*(1-cbetabar*cgamma*cprobp)/cprobp)/((cfc-1)*curvp+1)*(mc)  )  + spinf;
     w = (1/(1+cbetabar*cgamma))*w(-1)+(cbetabar*cgamma/(1+cbetabar*cgamma))*w(1)+(cindw/(1+cbetabar*cgamma))*pinf(-1)-(1+cbetabar*cgamma*cindw)/(1+cbetabar*cgamma)*pinf+(cbetabar*cgamma)/(1+cbetabar*cgamma)*pinf(1)+(1-cprobw)*(1-cbetabar*cgamma*cprobw)/((1+cbetabar*cgamma)*cprobw)*(1/((clandaw-1)*curvw+1))*(csigl*lab + (1/(1-chabb/cgamma))*c - ((chabb/cgamma)/(1-chabb/cgamma))*c(-1) -w)+ 1*sw;
     [name='taylor_rule']
-    r = crpi*(1-crr)*pinf+cry*(1-crr)*(y-yf)+crdy*(y-yf-y(-1)+yf(-1))+crr*r(-1)+ms;
+    r = crpiMcrpiXcrr*pinf + cryMcryXcrr*ygap + crdy*diff(ygap) + crr*r(-1) + ms;
+    ygap = y - yf;

In Smets and Wouters (2007), potential output (yf) is calculated elsewhere in the model. But, if we are to run OLS on the Taylor rule equation, we need observations on this variable. Since that is not possible, I add a new equation on line 143, ygap = y - yf, replacing y-yf in the Taylor rule with ygap. This allows us to estimate the model as before while providing observed data on ygap for OLS estimation (more on that later). NB: I could have left y-yf in the equation and provided data for yf but this change makes it more clear that ygap is calculated differently for the call to the Dynare estimation command than for the call to dyn_ols.

I further create two new parameters, crpiMcrpiXcrr and cryMcryXcrr because the parsing algorithm implemented in the OLS routine only accounts for the situation where the linear combination of a single parameter multiplies one or more endogenous variables. This change implies the changes on lines 76 and 77 defining their initial values, and on lines 222 and 223 defining them as parameters to be estimated.

Finally, I redefine ms as an exogenous variable and remove the equation that defined it as well as the definition of the variable em.

Last Modification

@@ -227,6 +231,14 @@ end;

 varobs dy dc dinve labobs pinfobs dw robs;

-estimation(optim=('MaxIter',200),datafile=usmodel_data,mode_file=usmodel_shock_decomp_mode,mode_compute=0,first_obs=1, presample=4,lik_init=2,prefilter=0,mh_replic=0,mh_nblocks=2,mh_jscale=0.20,mh_drop=0.2, nograph, nodiagnostic, tex);
+ds = dseries('usmodel_dseries.csv');
+ds.ygap = ds.y.detrend(1);
+dyn_ols(ds, {}, {'taylor_rule'});
+crr           = 0.8762;
+crdy          = 0.2347;
+crpiMcrpiXcrr = 0.1842;
+cryMcryXcrr   = 0.0073;
+
+estimation(optim=('MaxIter',200),datafile=usmodel_data,mode_compute=4,first_obs=1, presample=4,lik_init=2,prefilter=0,mh_replic=0,mh_nblocks=2,mh_jscale=0.20,mh_drop=0.2, nograph, nodiagnostic, tex);

I first load the data into a dseries called ds on line 234. As mentioned above, I need data for ygap for the OLS estimation. This is obtained by detrending the output series. I then call the OLS routine, telling it to run OLS using the dataset ds on the equation specified by the equation tag taylor_rule.

As the OLS routine sets the parameter values it estimates in M_.params, I reset their initial values after the call to the routine on lines 238-241, in preparation for the call to the Dynare estimation routine.

The OLS routine in MATLAB: dyn_ols.m

The OLS routine outlined herein was written in MATLAB but could just as easily have been written in Julia, Python, C, or the language of your choice. There are three main steps involved in writing a routine that makes use of the Dynare JSON output:

  1. Parse the JSON file, loading it into a language-specific structure
  2. Parse this structure for your purposes
  3. Run your computational task, in our case estimation via OLS

The files described in this section can be found here.

Step 1: Parsing the JSON file

As JSON is widely supported, the first step is often straightforward, regardless of your choice of programming language. In our case, though MATLAB has only provided JSON support since R2016b (via the commands jsonencode/jsondecode), there’s a widely-used and well-tested toolbox called JSONLab that allows one to use JSON with older versions of MATLAB and with Octave. Using the JSONLab distributed with Dynare allows us to access the model block specified as an AST in just two lines:

45 tmp = loadjson(jsonfile);
46 ast = tmp.abstract_syntax_tree;

These lines reside in the function get_ast.m, a utility we wrote that is used by dyn_ols.m.

Line 45 reads in Smets_Wouters_2007/model/json/modfile-original.json (stored in the jsonfile variable) and loads it into a temporary MATLAB structure. Line 46 then selects the abstract_syntax_tree field as that is the only one we’re interested in and sets it to the output variable ast. When finished, ast is a cell array with 40 entries, one for each equation. Entry 23 of this cell array corresponds to the monetary policy equation and looks like:

>> ast{23}

ans =

  struct with fields:

    number: 22
      line: 141
      tags: [1x1 struct]
       AST: [1x1 struct]

>> ast{23}.tags

ans =

  struct with fields:

    name: 'taylor_rule'

>> ast{23}.AST

ans =

  struct with fields:

    node_type: 'BinaryOpNode'
           op: '='
         arg1: [1x1 struct]
         arg2: [1x1 struct]

As you can see, the JSON output contains a lot of information. First, we know the 0-indexed equation number (MATLAB indexing starts at 1, which explains the difference between the index in the structure and the index in the JSON), and the line it was on in the .mod file (line 141). Digging into the structure, we see that there is one tag associated with the equation that has key name and value 'taylor_rule'. We also see the root node of the AST.

As we only want to perform OLS on this equation, we don’t need the AST corresponding to the other equations in the .mod file. Hence, we overwrite it with the aid of a helper function, selecting only the necessary equation:

54     ast = getEquationsByTags(ast, 'name', eqtags);

Step 2: Parsing the model block

Below I describe the parsing of an equation as implemented in dyn_ols.m.

Given the pared-down ast variable returned by getEquationsByTags.m, we then call another helper function, common_parsing.m that handles the parsing for OLS-style estimation routines (e.g. SUR, pooled OLS, FGLS, …).

118 [Y, lhssub, X, fp, lp] = common_parsing(ds(ds_range), ast, true, param_names);

In turn, this function calls a helper function, parse_ols_style_equation.m that handles the parsing of one OLS-style equation at a time.

57 for i = 1:neqs
58     [Y{i}, lhssub{i}, X{i}, residnames{i}, startdates{i}, enddates{i}] = ...
59         parse_ols_style_equation(ds, ast{i});
60 end

Together, these functions return matrices corresponding to the observed variable, Y, the regressors, X, the constant term, lhssub, and the first and last periods of the observed data (fp and lp).

Parsing proceeds as follows in parse_ols_style_equation.m. First, the function arguments are verified (lines 38-66). After this, we know that the LHS contains either a VariableNode or a UnaryOpNode.

Note

Several different types of nodes that can appear in the AST, corresponding to the types of operations that are available in a .mod file equation. The nodes are:

  1. NumConstNode: non-negative integers or doubles
  2. VariableNode: endogenous, exogenous, or parameter
  3. UnaryOpNode: unary operation on a node, e.g. log, abs, unary minus, …
  4. BinaryOpNode: binary operation on a node, e.g. arithmetic operations, min, max, =, comparison operators, …
  5. TrinaryOpNode: trinary operation on a node, e.g. normcdf, normpdf
  6. ExternalFunctionNode: external function node

Assured that the type of the LHS is valid for OLS, the value of Y is set by evaluating the LHS of the specified equation. Given that all equations are represented as BinaryOpNode‘s (the two arguments being the LHS, arg1, and the RHS, arg2), the call to evaluate the LHS of the equation is:

69 Y = evalNode(ds, ast.AST.arg1, line, dseries());

evalNode, a local function in parse_ols_style_equation.m then evaluates this node given the dseries contained in ds.

309 function X = evalNode(ds, node, line, X)
310 global M_
311 if strcmp(node.node_type, 'NumConstNode')
312     X = dseries(node.value, ds.dates, 'const');
313 elseif strcmp(node.node_type, 'VariableNode')
314     if strcmp(node.type, 'endogenous') ...
315             || (strcmp(node.type, 'exogenous') && any(strcmp(ds.name, node.name)))
316         if ds.exist(node.name)
317             X = ds.(node.name)(node.lag);
318         else
319             error('Variable %s is not available in the database.', node.name)
320         end
321     elseif strcmp(node.type, 'parameter')
322         X = M_.params(not(cellfun('isempty', strfind(M_.param_names, node.name))));
323         if isnan(X) || isinf(X) || ~isreal(X)
324             parsing_error(['Value incorrectly set for parameter: ' node.name], line);
325         end
326     end
327 elseif strcmp(node.node_type, 'UnaryOpNode')
328     Xtmp = evalNode(ds, node.arg, line, X);
329     % Only works if dseries supports . notation for unary op (true for log/diff)
330     % Otherwise, use: X = eval([node.op '(Xtmp)']);
331     try
332         if strcmp(node.op, 'uminus')
333             X = -Xtmp;
334         else
335             X = Xtmp.(node.op);
336         end
337         if any(isinf(X)) || ~isreal(X)
338             parsing_error(['Error applying ' node.op], line, node);
339         end
340     catch
341         parsing_error(['Error applying ' node.op], line, node);
342     end
343 elseif strcmp(node.node_type, 'BinaryOpNode')
344     Xtmp1 = evalNode(ds, node.arg1, line, X);
345     Xtmp2 = evalNode(ds, node.arg2, line, X);
346     switch node.op
347         case '*'
348             Xtmp = Xtmp1 * Xtmp2;
349         case '/'
350             Xtmp = Xtmp1 / Xtmp2;
351         case '+'
352             Xtmp = Xtmp1 + Xtmp2;
353         case '-'
354             Xtmp = Xtmp1 - Xtmp2;
355         otherwise
356             parsing_error(['got unexpected binary op ' node.op], line, node);
357     end
358     if any(isinf(Xtmp)) || ~isreal(Xtmp)
359         parsing_error(['Error applying ' node.op], line, node);
360     end
361     X = X + Xtmp;
362 else
363     parsing_error(['got unexpected node type ' node.node_type], line, node);
364 end
365 end

In the evaluation of the observed variable (X in this local function, Y in the main routine), what’s important to us is lines 313-342 as we know it’s either a VariableNode or a UnaryOpNode. Regardless of the type of node it is, it’s evaluated using the dseries (ds) and the corresponding dseries vector is returned. In our case, the interest rate r is evaluated on line 317 by looking up its value in ds.

Though our work for Y was easy, computing the matrix of regressors and the vector of constants will be a bit more involved.

Back in the main function, before we can create the matrix of regressors X, we decompose the RHS (arg2) of the equation into additive terms, storing them in a cell array called terms. We do this by calling the locally defined function:

232 function terms = decomposeAdditiveTerms(terms, node, node_sign)
233 if strcmp(node.node_type, 'NumConstNode') || strcmp(node.node_type, 'VariableNode')
234     terms = [terms {{node node_sign}}];
235 elseif strcmp(node.node_type, 'UnaryOpNode')
236     if strcmp(node.op, 'uminus')
237         terms = decomposeAdditiveTerms(terms, node.arg, -node_sign);
238     else
239         terms = [terms {{node node_sign}}];
240     end
241 elseif strcmp(node.node_type, 'BinaryOpNode')
242     if strcmp(node.op, '+') || strcmp(node.op, '-')
243         terms = decomposeAdditiveTerms(terms, node.arg1, node_sign);
244         if strcmp(node.op, '+')
245             terms = decomposeAdditiveTerms(terms, node.arg2, node_sign);
246         else
247             terms = decomposeAdditiveTerms(terms, node.arg2, -node_sign);
248         end
249     else
250         terms = [terms {{node node_sign}}];
251     end
252 else
253     terms = [terms {{node node_sign}}];
254 end
255 end

As you can see, decomposeAdditiveTerms is a recursive, tree-traversal function that breaks down terms separated by + or -, storing them in the return value, terms. Each cell in the return value is comprised of a pair of elements: the root node of the sub-tree representing the additive node and the sign preceding this node (1 or -1).

Given the equation we want to estimate (shown here from Smets_Wouters_2007.mod as a reminder):

141     [name='taylor_rule']
142     r = crpiMcrpiXcrr*pinf + cryMcryXcrr*ygap + crdy*diff(ygap) + crr*r(-1) + ms;

the output of decomposeAdditiveTerms (stored in the variable terms) is a 5-element cell, corresponding to the RHS elements:

  1. crpiMcrpiXcrr*pinf
  2. cryMcryXcrr*ygap
  3. crdy*diff(ygap)
  4. crr*r(-1)
  5. ms

Hence, in this case, 4 nodes of terms contain BinaryOpNode‘s (corresponding to the binary operation *) and one term contains a VariableNode.

To understand how the terms are stored, one need only look at the following output for the first term, crpiMcrpiXcrr*pinf:

>> terms

terms =

  1x5 cell array

    {1x2 cell}    {1x2 cell}    {1x2 cell}    {1x2 cell}    {1x2 cell}

>> terms{1}

ans =

  1x2 cell array

    {1x1 struct}    {[1]}

>> terms{1}{1}

ans =

  struct with fields:

    node_type: 'BinaryOpNode'
           op: '*'
         arg1: [1x1 struct]
         arg2: [1x1 struct]

>> terms{1}{1}.arg1

ans =

  struct with fields:

    node_type: 'VariableNode'
         name: 'pinf'
         type: 'endogenous'
          lag: 0

>> terms{1}{1}.arg2

ans =

  struct with fields:

    node_type: 'VariableNode'
         name: 'crpiMcrpiXcrr'
         type: 'parameter'
          lag: 0

Here, the second element of terms{1} is 1 as it is not preceded by a minus sign. The node itself is a BinaryOpNode as it represents the multiplication of an endogenous variable (the first argument of the node) and a parameter (the second argument of the node).

Note

The elements of terms are not necessarily in the same order as written in the equation in the .mod file. This is because, for efficiency reasons (e.g. node sharing), the AST created by the preprocessor does not guarantee this ordering. By the same token, the elements of a node are not guaranteed to be in the same order as they appear in the .mod file, as you can see above where the first argument of the BinaryOpNode in terms{1}{1} is the endogenous variable pinf whereas pinf appears second in the multiplication as written in the .mod file: crpiMcrpiXcrr*pinf.

Now that we have terms set, we can construct the regressor matrix X by entering the loop on line 75 of parse_ols_style_equation.m:

 75 for i = 1:length(terms)
 76     Xtmp = dseries();
 77     node_sign = terms{i}{2};
 78     node_to_parse = terms{i}{1};
 79     if strcmp(node_to_parse.node_type, 'VariableNode')
 80         if strcmp(node_to_parse.type, 'parameter')
 81             % Intercept
 82             Xtmp = dseries(1, ds.dates, node_to_parse.name)*node_sign;
 83         elseif strcmp(node_to_parse.type, 'exogenous') && ~any(strcmp(ds.name, node_to_parse.name))
 84             % Residual if not contained in ds
 85             if isempty(residual)
 86                 residual = node_to_parse.name;
 87             else
 88                 parsing_error(['only one residual allowed per equation; encountered ' residual ' & ' node_to_parse.name], line);
 89             end
 90         elseif strcmp(node_to_parse.type, 'endogenous') ...
 91                 || (strcmp(node_to_parse.type, 'exogenous') && any(strcmp(ds.name, node_to_parse.name)))
 92             % Subtract VariableNode from LHS
 93             % NB: treat exogenous that exist in ds as endogenous
 94             lhssub = lhssub + evalNode(ds, node_to_parse, line, dseries())*node_sign;
 95         else
 96             parsing_error('unexpected variable type found', line, node_to_parse);
 97         end
 98     elseif strcmp(node_to_parse.node_type, 'UnaryOpNode')
 99         % Subtract UnaryOpNode from LHS
100         % NB: treat exogenous that exist in ds as endogenous
101         lhssub = lhssub + evalNode(ds, node_to_parse, line, dseries());
102     elseif strcmp(node_to_parse.node_type, 'BinaryOpNode') && strcmp(node_to_parse.op, '/')
103         % Subtract Node from LHS
104         % if either arg contains a parameter, it's a parsing error.
105         if containsParameter(node_to_parse, line)
106             parsing_error('unexpected node found', line, node_to_parse)
107         end
108         lhssub = lhssub + evalNode(ds, node_to_parse, line, dseries());
109     elseif strcmp(node_to_parse.node_type, 'BinaryOpNode') && strcmp(node_to_parse.op, '*')
110         % Parse param_expr * endog_expr
111         [Xtmp, names] = parseTimesNode(ds, node_to_parse, line);
112         Xtmp = Xtmp*node_sign;
113         if Xtmp.vobs > 1 || ...
114                 (Xtmp.vobs == 1 && ~isnan(str2double(Xtmp.name)))
115             % Handle constraints
116             % Look through Xtmp names for constant
117             % if found, subtract from LHS
118             for j = 1:length(names)
119                 if strcmp(names{j}{1}.node_type, 'NumConstNode')
120                     pname = num2str(names{j}{1}.value);
121                 elseif strcmp(names{j}{1}.node_type, 'VariableNode')
122                     pname = names{j}{1}.name;
123                 else
124                     parsing_error('unexpected node type', node_to_parse, line);
125                 end
126                 psign = names{j}{2};
127                 if ~isnan(str2double(pname))
128                     lhssub = lhssub + psign * str2double(pname) * Xtmp.(pname);
129                     Xtmp = Xtmp.remove(pname);
130                 else
131                     % Multiply by psign now so that it can be added together below
132                     % Otherwise, it would matter which was encountered first,
133                     % a parameter on its own or a linear constraint
134                     Xtmp.(pname) = psign * Xtmp.(pname);
135                 end
136             end
137         end
138     else
139         parsing_error('didn''t expect to arrive here', line, node_to_parse);
140     end
141 
142     names = Xtmp.name;
143     for j = length(names):-1:1
144         % Handle constraits
145         if any(strcmp(X.name, names{j}))
146             X.(names{j}) = X.(names{j}) + Xtmp.(names{j});
147             Xtmp = Xtmp.remove(names{j});
148         end
149     end
150     X = [X Xtmp];
151 end

Entering the loop, we set 3 variables: Xtmp, node_to_parse, and node_sign. Xtmp is where we will construct the regressor column to be appended to X at the end of each loop. node_to_parse and node_sign are simply the corresponding parts of the pair stored in each terms cell, as described above.

With these variables set, we take different actions, depending on the type of node encountered in node_to_parse. The following subsections explain those actions.

Condition 1: VariableNode (lines 79-97)

If node_to_parse is a VariableNode it is an additive variable in the equation. If it was declared as a parameter in the .mod file, then it’s the intercept of the equation. It is thus stored in Xtmp with the value of the parameter at every period. If it’s a lone exogenous variable that is not present in ds, we treat it as the residual. If there are more than one such variable, then it’s an error. If it’s an endogenous variable or an exogenous variable present in ds, we add it to the dseries lhssub, a dseries that will be subtracted from the LHS before returning. If the type of the node_to_parse doesn’t fall into any of these categories, parsing ends with an error.

Condition 2: UnaryOpNode (lines 98-101)

If node_to_parse is a UnaryOpNode, we evaluate it and add it to the lhssub variable to be subtracted from the LHS.

Condition 3: BinaryOpNode with division operator (lines 102-108)

In this case, if a parameter is found in this expression, we end with a parsing error. Otherwise, we evaluate the node and add it to lhssub.

Condition 4: BinaryOpNode with multiplication operator (lines 109-137)

In this case, we parse the node_to_parse, by calling the local function parseTimesNode:

257 function [X, pterms] = parseTimesNode(ds, node, line)
258 % Separate the parameter expression from the endogenous expression
259 assert(strcmp(node.node_type, 'BinaryOpNode') && strcmp(node.op, '*'))
260 if isOlsParamExpr(node.arg1, line)
261     pterms = decomposeAdditiveTerms([], node.arg1, 1);
262     X = evalNode(ds, node.arg2, line, dseries());
263 elseif isOlsParamExpr(node.arg2, line)
264     pterms = decomposeAdditiveTerms([], node.arg2, 1);
265     X = evalNode(ds, node.arg1, line, dseries());
266 else
267     parsing_error('expecting (param expr)*(var expr)', line, node);
268 end
269 if strcmp(pterms{1}{1}.node_type, 'NumConstNode')
270     X = X.rename(num2str(pterms{1}{1}.value));
271 elseif strcmp(pterms{1}{1}.node_type, 'VariableNode')
272     X = X.rename(pterms{1}{1}.name);
273 else
274     parsing_error('unexpected type', line, node)
275 end
276 for ii = 2:length(pterms)
277     if strcmp(pterms{ii}{1}.node_type, 'NumConstNode')
278         X = [X dseries(X{1}.data, X{1}.firstdate, num2str(pterms{ii}{1}.value))];
279     elseif strcmp(pterms{ii}{1}.node_type, 'VariableNode')
280         X = [X dseries(X{1}.data, X{1}.firstdate, pterms{ii}{1}.name)];
281     else
282         parsing_error('unexpected type', line, node)
283     end
284 end
285 end

This function returns the dseries vector/matrix X corresponding to the regressors found. It can handle parsing expressions of additively-separated parameters multiplied by an expression of variables. Each of the additively separated parameters in pterms corresponds to one column of the returned vector/matrix X.

Returning from this function, we have Xtmp corresponding to the evaluated variables and names corresponding to the parameter names. In the code that follows (lines 113-137), we add any columns of Xtmp that were multiplied by a constant to lhssub and then remove this column from Xtmp.

Condition 5: Otherwise (line 139)

If none of the previous 4 conditions were satisfied, we have encountered a parsing error.

End of loop (lines 141-150)

At the end of the loop, we combine the temporary vector/matrix of regressors, Xtmp with the matrix of regressors already created, X. The loop here handles the case where a parameter was used twice in two different terms. If any columns remain in Xtmp once we arrive at line 151, they are concatenated to X. The loop continues until the cell terms has been exhausted. When the loop ends on line 152, our matrix of regressors, X, has been constructed.

After the loop, the Y vector is adjusted, subtracting any constant terms on the RHS:

152 Y = Y - lhssub;

Following that, the first and last observed periods of the dseries are assigned to the variables fp and lp:

154 %% Set start and end dates
155 fp = Y.firstobservedperiod;
156 lp = Y.lastobservedperiod;
157 if ~isempty(X)
158     % X is empty when AR(1) without parameter is encountered
159     fp = max(fp, X.firstobservedperiod);
160     lp = min(lp, X.lastobservedperiod);
161 end
162 if ~isempty(lhssub)
163     fp = max(fp, lhssub.firstobservedperiod);
164     lp = min(lp, lhssub.lastobservedperiod);
165 end

Finally, the Y, X, and lhssub datasets are adjusted given fp and lp:

188 Y = Y(fp:lp);
189 if ~isempty(X)
190     X = X(fp:lp);
191     names = X.name;
192     for i = 1:length(names)
193         if all(X.(names{i}).data == 0)
194             X = X.remove(names{i});
195         end
196     end
197 end
198 if ~isempty(lhssub)
199     lhssub = lhssub(fp:lp);
200 end

Step 3: Estimation via OLS

Back in dyn_ols.m, after parsing the equations as discussed above, we have the cell arrays Y, X, and lhssub, where each entry in these cell arrays corresponds to the vector/matrix of observables and regressors in each equation to be estimated. Analogously, fp and lp are cells containing the first and last observed period of the estimation range for each equation.

118 [Y, lhssub, X, fp, lp] = common_parsing(ds(ds_range), ast, true, param_names);

We loop over these cell arrays, running our estimation for each equation: \(\hat{\beta} = (X'X)^{-1}X'Y\). [1] The output is set to the ols field of the standard Dynare MATLAB/Octave output structure, oo_. Each sub-field of oo_.ols corresponds to an equation tag that was either provided by the user or created by dyn_ols.m. In this case, they will be saved to oo_.ols.taylor_rule. Furthermore, the estimated parameter values are set in M_.

125     if ~isempty(model_name{i})
126         tag = model_name{i};
127     else
128         if isfield(ast{i}, 'tags') && isfield(ast{i}.tags, 'name')
129             tag = ast{i}.tags.('name');
130         else
131             tag = ['eq_line_no_' num2str(ast{i}.line)];
132         end
133     end
134 
135     %% Estimation
136     % From LeSage, James P. "Applied Econometrics using MATLAB"
137     oo_.ols.(tag).dof = nobs - nvars;
138 
139     % Estimated Parameters
140     [q, r] = qr(X{i}.data, 0);
141     xpxi = (r'*r)\eye(nvars);
142     oo_.ols.(tag).beta = r\(q'*Y{i}.data);
143     oo_.ols.(tag).param_idxs = zeros(length(pnames), 1);
144     for j = 1:length(pnames)
145         if ~strcmp(pnames{j}, 'intercept')
146             oo_.ols.(tag).param_idxs(j) = find(strcmp(M_.param_names, pnames{j}));
147             M_.params(oo_.ols.(tag).param_idxs(j)) = oo_.ols.(tag).beta(j);
148         end
149     end

And that’s it! The rest of the code simply takes care of calculating the various statistics and standard errors (all stored in oo_.ols) and displaying the estimated parameters in a table:

     OLS Estimation of equation 'taylor_rule' [name = 'taylor_rule']

Dependent Variable: r
No. Independent Variables: 4
Observations: 231 from 1947Q2 to 2004Q4

                    Estimates          t-statistic          Std. Error
                 ________________    ________________    ________________

crpiMcrpiXcrr            0.062335              2.6913            0.023161
cryMcryXcrr             0.0087422              2.7559           0.0031722
crdy                     0.056321              4.0914            0.013766
crr                       0.95707              59.384            0.016117

R^2: 0.942770
R^2 Adjusted: 0.942014
s^2: 0.043718
Durbin-Watson: 1.703493
_________________________________________________________________________

Conclusion

This was just one example of how Dynare’s new JSON output can be exploited to construct your own back-end routines in the language of your choosing. It essentially frees you from the Dynare back-end and allows you to build your own library routines while taking advantage of the Dynare modeling language.

We hope you find this development useful. If this has encouraged you to learn more about Dynare, please don’t hesitate to visit our web page where you can find guidelines for contributing. If you notice a bug in the JSON output, don’t hesitate to report it on the Dynare-preprocessor Issues page. Again, as a reminder, JSON output is provided in Dynare 4.6 and later.

[1]As matrix inversion is slow and numerically unstable for small values, we use the QR decomposition instead of the standard formula: \(\hat{\beta} = R^{-1}Q'Y\).

Comments

comments powered by Disqus