Compare revisions

Sébastien Villemot · Houtan Bastani · Houtan Bastani · Houtan Bastani · Houtan Bastani · Sébastien Villemot
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -15,7 +15,7 @@ build_linux_32:
  stage: build
  script:
    - autoreconf -si
-    - './configure LDFLAGS="-m32 -static -static-libgcc -static-libstdc++" CXXFLAGS=-m32'
+    - './configure --host=i686-linux-gnu LDFLAGS="-static -static-libgcc -static-libstdc++"'
    - make -j$(nproc)
    - strip src/dynare_m
    - mkdir -p bin
@@ -32,7 +32,7 @@ build_linux_64:
  stage: build
  script:
    - autoreconf -si
-    - './configure LDFLAGS="-static -static-libgcc -static-libstdc++"'
+    - './configure --host=x86_64-linux-gnu LDFLAGS="-static -static-libgcc -static-libstdc++"'
    - make -j$(nproc)
    - strip src/dynare_m
    - mkdir -p bin
@@ -93,7 +93,7 @@ build_macOS:
    - macOS
  script:
    - autoreconf -si
-    - './configure CXX=g++-9 CXXFLAGS=-static-libgcc'
+    - './configure CXX=g++-10 CXXFLAGS=-static-libgcc LEX=/usr/local/opt/flex/bin/flex YACC=/usr/local/opt/bison/bin/bison'
    - make -j$(nproc)
    - strip src/dynare_m
    - mkdir -p bin

--- a/Makefile.am
+++ b/Makefile.am
@@ -2,9 +2,7 @@ SUBDIRS = src doc
 ACLOCAL_AMFLAGS = -I m4
-EXTRA_DIST = \
+EXTRA_DIST = COPYING
-	license.txt \
-	COPYING
 dist-hook:
 	rm -rf `find $(distdir) -name '.git*'`

--- a/configure.ac
+++ b/configure.ac
@@ -18,7 +18,7 @@ dnl You should have received a copy of the GNU General Public License
 dnl along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 AC_PREREQ([2.62])
-AC_INIT([dynare-preprocessor], [4.6-unstable])
+AC_INIT([dynare-preprocessor], [4.6.5])
 AC_CONFIG_AUX_DIR([build-aux])
 AC_CONFIG_SRCDIR([src/DynareMain.cc])
 AM_INIT_AUTOMAKE([1.11 -Wall -Wno-portability foreign no-dist-gzip dist-xz tar-pax])
@@ -55,8 +55,7 @@ AM_PROG_AR
 AM_PROG_LEX
 # Hack to get lex include dir, ticket #575
-AC_PATH_PROG([LEXPATH], [$LEX])
+AC_SUBST([LEXINC], [$(echo "$(dirname "$(which $LEX)")"/../include)])
-AC_SUBST([LEXINC], [$(eval "echo $LEXPATH | sed 's|\(.*\)$LEX$|\1../include|'")])
 AC_CHECK_PROG([YACC], [bison], [bison])
 if test -z "$YACC"; then

--- a/license.txt
+++ b/license.txt
-Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
-Upstream-Name: Dynare-preprocessor
-Upstream-Contact: Dynare Team, whose members in 2019 are:
-                  Stéphane Adjemian <stephane.adjemian@univ-lemans.fr>
-                  Houtan Bastani <houtan@dynare.org>
-                  Michel Juillard <michel.juillard@mjui.fr>
-                  Frédéric Karamé <frederic.karame@univ-lemans.fr>
-                  Junior Maih <junior.maih@gmail.com>
-                  Ferhat Mihoubi <fmihoubi@univ-evry.fr>
-                  Johannes Pfeifer <jpfeifer@gmx.de>
-                  Marco Ratto <marco.ratto@ec.europa.eu>
-                  Sébastien Villemot <sebastien@dynare.org>
-Source: https://www.dynare.org
-Files: *
-Copyright: 1996-2019 Dynare Team
-License: GPL-3+
-Files: m4/ax_boost_base.m4
-       m4/ax_boost_system.m4
-       m4/ax_boost_filesystem.m4
-Copyright: 2008-2009 Thomas Porschberg <thomas@randspringer.de>
-           2009 Peter Adolphs
-           2008-2009 Michael Tindal
-           2008 Daniel Casimiro <dan.casimiro@gmail.com>
-           2009 Roman Rybalko <libtorrent@romanr.info>
-License: FSFAP
-Files: m4/ax_cxx_compile_stdcxx.m4
-       m4/ax_cxx_compile_stdcxx_17.m4
-Copyright: 2008 Benjamin Kosnik <bkoz@redhat.com>
-           2012 Zack Weinberg <zackw@panix.com>
-           2013 Roy Stogner <roystgnr@ices.utexas.edu>
-           2014, 2015 Google Inc.; contributed by Alexey Sokolov <sokolov@google.com>
-           2015 Paul Norman <penorman@mac.com>
-           2015 Moritz Klammler <moritz@klammler.eu>
-           2016, 2018 Krzesimir Nowak <qdlacz@gmail.com>
-           2019 Enji Cooper <yaneurabeya@gmail.com>
-License: FSFAP
-Files: m4/ax_latex_class.m4 m4/ax_tex_test.m4
-Copyright: 2008 Boretti Mathieu <boretti@eig.unige.ch>
-           2009 Dynare Team
-License: LGPL-2.1+
-License: GPL-3+
- Dynare is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
- .
- Dynare is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU General Public License for more details.
- .
- You should have received a copy of the GNU General Public License
- along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
-License: GPL-3+ with Autoconf exception
- This program is free software: you can redistribute it and/or modify it
- under the terms of the GNU General Public License as published by the
- Free Software Foundation, either version 3 of the License, or (at your
- option) any later version.
- .
- This program is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
- Public License for more details.
- .
- You should have received a copy of the GNU General Public License along
- with this program. If not, see <http://www.gnu.org/licenses/>.
- .
- As a special exception, the respective Autoconf Macro's copyright owner
- gives unlimited permission to copy, distribute and modify the configure
- scripts that are the output of Autoconf when processing the Macro. You
- need not follow the terms of the GNU General Public License when using
- or distributing such scripts, even though portions of the text of the
- Macro appear in them. The GNU General Public License (GPL) does govern
- all other use of the material that constitutes the Autoconf Macro.
- .
- This special exception to the GPL applies to versions of the Autoconf
- Macro released by the Autoconf Archive. When you make and distribute a
- modified version of the Autoconf Macro, you may extend this special
- exception to the GPL to apply to your modified version as well.
-License: LGPL-2.1+
- This library is free software; you can redistribute it and/or modify it
- under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation; either version 2.1 of the License, or (at
- your option) any later version.
- .
- This library is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
- General Public License for more details.
- .
- You should have received a copy of the GNU Lesser General Public License
- along with this library. If not, see <http://www.gnu.org/licenses/>.
-License: LGPL-3+
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU Lesser General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
- .
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- GNU Lesser General Public License for more details.
- .
- You should have received a copy of the GNU Lesser General Public License
- along with this program.  If not, see <http://www.gnu.org/licenses/>.
-License: FSFAP
- Copying and distribution of this file, with or without modification, are
- permitted in any medium without royalty provided the copyright notice
- and this notice are preserved. This file is offered as-is, without any
- warranty.
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -34,6 +34,7 @@ using namespace std;
 #pragma GCC diagnostic pop
 #include <utility>
+#include <algorithm>
 SteadyStatement::SteadyStatement(OptionsList options_list_arg) :
  options_list{move(options_list_arg)}
@@ -477,7 +478,7 @@ VarRestrictionsStatement::writeOutput(ostream &output, const string &basename, b
  if (itvs == var_map.end())
    {
      cerr << "ERROR: you are imposing restrictions on a VAR named " << var_model_name
-           << " but this VAR has not been declared via thevar_model statement." << endl;
+           << " but this VAR has not been declared via the var_model statement." << endl;
      exit(EXIT_FAILURE);
    }
  vector<string> vars = itvs->second;
@@ -646,7 +647,7 @@ StochSimulStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -667,7 +668,7 @@ StochSimulStatement::writeOutput(ostream &output, const string &basename, bool m
  options_list.writeOutput(output);
  symbol_list.writeOutput("var_list_", output);
-  output << "[info, oo_, options_] = stoch_simul(M_, options_, oo_, var_list_);" << endl;
+  output << "[info, oo_, options_, M_] = stoch_simul(M_, options_, oo_, var_list_);" << endl;
 }
 void
@@ -699,7 +700,7 @@ ForecastStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolida
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -733,56 +734,6 @@ ForecastStatement::writeJsonOutput(ostream &output) const
  output << "}";
 }
-DetCondForecast::DetCondForecast(SymbolList symbol_list_arg,
-                                 OptionsList options_list_arg,
-                                 const bool linear_decomposition_arg) :
-  options_list{move(options_list_arg)},
-  symbol_list{move(symbol_list_arg)},
-  linear_decomposition{linear_decomposition_arg}
-{
-}
-void
-DetCondForecast::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
-{
-  options_list.writeOutput(output);
-  if (linear_decomposition)
-    {
-      output << "first_order_solution_to_compute = 1;" << endl
-             << "if isfield(oo_, 'dr')" << endl
-             << "    if isfield(oo_.dr, 'ghx') && isfield(oo_.dr, 'ghu') && isfield(oo_.dr, 'state_var') && isfield(oo_.dr, 'order_var')" << endl
-             << "        first_order_solution_to_compute = 0;" << endl
-             << "    end" << endl
-             << "end" << endl
-             << "if first_order_solution_to_compute" << endl
-             << "  fprintf('%s','Computing the first order solution ...');" << endl
-             << "  options_.nograph = true;" << endl
-             << "  options_.order = 1;" << endl
-             << "  options_.noprint = true;" << endl
-             << "  options_.nocorr = true;" << endl
-             << "  options_.nomoments = true;" << endl
-             << "  options_.nodecomposition = true;" << endl
-             << "  options_.nofunctions = true;" << endl
-             << "  options_.irf = 0;" << endl
-             << "  tmp_periods = options_.periods;" << endl
-             << "  options_.periods = 0;" << endl
-             << "  var_list_ = char();" << endl
-             << "  info = stoch_simul(var_list_);" << endl
-             << R"(  fprintf('%s\n','done');)" << endl
-             << "  options_.periods = tmp_periods;" << endl
-             << "end" << endl;
-    }
-  vector<string> symbols = symbol_list.get_symbols();
-  if (symbols.size() > 0)
-    output << symbols[1] << " = det_cond_forecast(";
-  for (unsigned int i = 0; i < symbols.size() - 1; i++)
-    output << symbols[i] << ", ";
-  if (symbols.size() > 0)
-    output << symbols[symbols.size() - 1];
-  output << ");" << endl;
-}
 RamseyModelStatement::RamseyModelStatement(OptionsList options_list_arg) :
  options_list{move(options_list_arg)}
 {
@@ -989,7 +940,7 @@ RamseyPolicyStatement::checkPass(ModFileStructure &mod_file_struct, WarningConso
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -1093,7 +1044,7 @@ DiscretionaryPolicyStatement::checkPass(ModFileStructure &mod_file_struct, Warni
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -1114,7 +1065,7 @@ DiscretionaryPolicyStatement::writeOutput(ostream &output, const string &basenam
  options_list.writeOutput(output);
  symbol_list.writeOutput("var_list_", output);
-  output << "[info, oo_, options_] = discretionary_policy(M_, options_, oo_, var_list_);" << endl;
+  output << "[info, oo_, options_, M_] = discretionary_policy(M_, options_, oo_, var_list_);" << endl;
 }
 void
@@ -1134,8 +1085,10 @@ DiscretionaryPolicyStatement::writeJsonOutput(ostream &output) const
  output << "}";
 }
-EstimationStatement::EstimationStatement(SymbolList symbol_list_arg,
+EstimationStatement::EstimationStatement(const SymbolTable &symbol_table_arg,
+                                         SymbolList symbol_list_arg,
                                         OptionsList options_list_arg) :
+  symbol_table{symbol_table_arg},
  symbol_list{move(symbol_list_arg)},
  options_list{move(options_list_arg)}
 {
@@ -1212,9 +1165,25 @@ EstimationStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
      exit(EXIT_FAILURE);
    }
+  /* Check that we are not trying to estimate a parameter appearing in the
+     planner discount factor (see dynare#1173) */
+  vector<int> estimated_params_in_planner_discount;
+  set_intersection(mod_file_struct.estimated_parameters.begin(),
+                   mod_file_struct.estimated_parameters.end(),
+                   mod_file_struct.parameters_in_planner_discount.begin(),
+                   mod_file_struct.parameters_in_planner_discount.end(),
+                   back_inserter(estimated_params_in_planner_discount));
+  if (!estimated_params_in_planner_discount.empty())
+    {
+      cerr << "ERROR: It is not possible to estimate a parameter ("
+           << symbol_table.getName(estimated_params_in_planner_discount[0])
+           << ") that appears in the discount factor of the planner (i.e. in the 'planner_discount' option)." << endl;
+      exit(EXIT_FAILURE);
+    }
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -1275,7 +1244,11 @@ DynareSensitivityStatement::checkPass(ModFileStructure &mod_file_struct, Warning
 {
  if (auto it = options_list.num_options.find("identification");
      it != options_list.num_options.end() && it->second == "1")
-    mod_file_struct.identification_present = true;
+    {
+      mod_file_struct.identification_present = true;
+      // The following triggers 3rd order derivatives, see preprocessor#40
+      mod_file_struct.identification_order = max(mod_file_struct.identification_order, 2);
+    }
  mod_file_struct.sensitivity_present = true;
 }
@@ -1318,6 +1291,20 @@ RplotStatement::RplotStatement(SymbolList symbol_list_arg) :
 {
 }
+void
+RplotStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
+{
+  try
+    {
+      symbol_list.checkPass(warnings, { SymbolType::endogenous, SymbolType::exogenous});
+    }
+  catch (SymbolList::SymbolListException &e)
+    {
+      cerr << "ERROR: rplot: " << e.message << endl;
+      exit(EXIT_FAILURE);
+    }
+}
 void
 RplotStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
@@ -1886,6 +1873,16 @@ OsrParamsStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolid
  if (mod_file_struct.osr_params_present)
    cerr << "WARNING: You have more than one osr_params statement in the .mod file." << endl;
  mod_file_struct.osr_params_present = true;
+  try
+    {
+      symbol_list.checkPass(warnings, { SymbolType::parameter });
+    }
+  catch (SymbolList::SymbolListException &e)
+    {
+      cerr << "ERROR: osr: " << e.message << endl;
+      exit(EXIT_FAILURE);
+    }
 }
 void
@@ -1994,7 +1991,7 @@ OsrStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -2120,7 +2117,7 @@ DynaSaveStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolida
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous, SymbolType::exogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -2162,7 +2159,7 @@ DynaTypeStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolida
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous, SymbolType::exogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -2561,7 +2558,7 @@ MSSBVARIrfStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsoli
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -2707,7 +2704,7 @@ IdentificationStatement::checkPass(ModFileStructure &mod_file_struct, WarningCon
      mod_file_struct.identification_order = max(mod_file_struct.identification_order, order);
    }
  else
-    // The default value for order is 1 (which triggers 2nd order dynamic derivatives)
+    // The default value for order is 1 (which triggers 2nd order dynamic derivatives, see preprocessor#40)
    mod_file_struct.identification_order = max(mod_file_struct.identification_order, 1);
 }
@@ -2838,7 +2835,7 @@ ShockDecompositionStatement::checkPass(ModFileStructure &mod_file_struct, Warnin
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -2888,7 +2885,7 @@ RealtimeShockDecompositionStatement::checkPass(ModFileStructure &mod_file_struct
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -2929,6 +2926,20 @@ PlotShockDecompositionStatement::PlotShockDecompositionStatement(SymbolList symb
 {
 }
+void
+PlotShockDecompositionStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
+{
+  try
+    {
+      symbol_list.checkPass(warnings, { SymbolType::endogenous, SymbolType::epilogue });
+    }
+  catch (SymbolList::SymbolListException &e)
+    {
+      cerr << "ERROR: plot_shock_decomposition: " << e.message << endl;
+      exit(EXIT_FAILURE);
+    }
+}
 void
 PlotShockDecompositionStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
@@ -2971,7 +2982,7 @@ InitialConditionDecompositionStatement::checkPass(ModFileStructure &mod_file_str
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -3017,7 +3028,7 @@ SqueezeShockDecompositionStatement::checkPass(ModFileStructure &mod_file_struct,
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -3096,7 +3107,7 @@ PlotConditionalForecastStatement::checkPass(ModFileStructure &mod_file_struct,
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -4670,7 +4681,7 @@ CalibSmootherStatement::checkPass(ModFileStructure &mod_file_struct, WarningCons
  mod_file_struct.calib_smoother_present = true;
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -4802,7 +4813,7 @@ GMMEstimationStatement::checkPass(ModFileStructure &mod_file_struct, WarningCons
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {
@@ -4849,7 +4860,7 @@ SMMEstimationStatement::checkPass(ModFileStructure &mod_file_struct, WarningCons
 {
  try
    {
-      symbol_list.checkPass(warnings);
+      symbol_list.checkPass(warnings, { SymbolType::endogenous });
    }
  catch (SymbolList::SymbolListException &e)
    {

--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -95,19 +95,6 @@ public:
  void writeJsonOutput(ostream &output) const override;
 };
-class DetCondForecast : public Statement
-{
-private:
-  const OptionsList options_list;
-  const SymbolList symbol_list;
-  const bool linear_decomposition;
-public:
-  DetCondForecast(SymbolList symbol_list_arg,
-                  OptionsList options_list_arg,
-                  const bool linear_decompositiontion_arg);
-  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
-};
 class ModelInfoStatement : public Statement
 {
 private:
@@ -248,7 +235,6 @@ public:
                        SymbolList symbol_list_arg,
                        OptionsList options_list_arg);
  void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
-  void checkRamseyPolicyList();
  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
  void writeJsonOutput(ostream &output) const override;
 };
@@ -279,6 +265,7 @@ private:
  const SymbolList symbol_list;
 public:
  explicit RplotStatement(SymbolList symbol_list_arg);
+  void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
  void writeJsonOutput(ostream &output) const override;
 };
@@ -314,10 +301,12 @@ public:
 class EstimationStatement : public Statement
 {
 private:
+  const SymbolTable &symbol_table;
  const SymbolList symbol_list;
  const OptionsList options_list;
 public:
-  EstimationStatement(SymbolList symbol_list_arg,
+  EstimationStatement(const SymbolTable &symbol_table_arg,
+                      SymbolList symbol_list_arg,
                      OptionsList options_list_arg);
  void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
@@ -747,6 +736,7 @@ private:
 public:
  PlotShockDecompositionStatement(SymbolList symbol_list_arg,
                                  OptionsList options_list_arg);
+  void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
  void writeJsonOutput(ostream &output) const override;
 };

--- a/src/DataTree.cc
+++ b/src/DataTree.cc
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2021 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -98,6 +98,14 @@ DataTree::operator=(const DataTree &d)
  // Constants must be initialized first because they are used in some Add* methods
  initConstants();
+  /* Model local variables must be next, because they can be evaluated in Add*
+     methods when the model equations are added. They need to be cloned in
+     order of appearance in the model block (hence with
+     local_variables_vector), because if there is a model_local_variable statement
+     the symbol IDs ordering may not be the right one (see dynare#1782) */
+  for (const auto &id : d.local_variables_vector)
+    local_variables_table[id] = d.local_variables_table.at(id)->clone(*this);
  for (const auto &it : d.node_list)
    it->clone(*this);
@@ -105,9 +113,6 @@ DataTree::operator=(const DataTree &d)
  local_variables_vector = d.local_variables_vector;
-  for (const auto &it : d.local_variables_table)
-    local_variables_table[it.first] = it.second->clone(*this);
  return *this;
 }

--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -130,7 +130,9 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
  xref_exo{m.xref_exo},
  xref_exo_det{m.xref_exo_det},
  nonzero_hessian_eqs{m.nonzero_hessian_eqs},
+  dynJacobianColsNbr{m.dynJacobianColsNbr},
  v_temporary_terms_inuse{m.v_temporary_terms_inuse},
+  variableMapping{m.variableMapping},
  map_idx{m.map_idx},
  global_temporary_terms{m.global_temporary_terms},
  block_type_firstequation_size_mfs{m.block_type_firstequation_size_mfs},
@@ -194,10 +196,14 @@ DynamicModel::operator=(const DynamicModel &m)
  xref_exo_det = m.xref_exo_det;
  nonzero_hessian_eqs = m.nonzero_hessian_eqs;
+  dynJacobianColsNbr = m.dynJacobianColsNbr;
  v_temporary_terms.clear();
  v_temporary_terms_inuse = m.v_temporary_terms_inuse;
+  variableMapping = m.variableMapping;
  first_chain_rule_derivatives.clear();
  map_idx = m.map_idx;
@@ -306,8 +312,6 @@ DynamicModel::computeTemporaryTermsOrdered()
              expr_t id = get<3>(it);
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
            }
-          for (const auto &it : derivative_endo[block])
-            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
          for (const auto &it : derivative_other_endo[block])
            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
          v_temporary_terms_inuse[block] = {};
@@ -337,8 +341,6 @@ DynamicModel::computeTemporaryTermsOrdered()
              expr_t id = get<3>(it);
              id->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
            }
-          for (const auto &it : derivative_endo[block])
-            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
          for (const auto &it : derivative_other_endo[block])
            it.second->computeTemporaryTerms(reference_count, temporary_terms, first_occurence, block, v_temporary_terms, block_size-1);
        }
@@ -364,8 +366,6 @@ DynamicModel::computeTemporaryTermsOrdered()
              expr_t id = get<3>(it);
              id->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
            }
-          for (const auto &it : derivative_endo[block])
-            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          for (const auto &it : derivative_other_endo[block])
            it.second->collectTemporary_terms(temporary_terms, temporary_terms_in_use, block);
          for (const auto &it : derivative_exo[block])
@@ -415,7 +415,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
      //recursive_variables.clear();
      feedback_variables.clear();
      //For a block composed of a single equation determines wether we have to evaluate or to solve the equation
-      nze = derivative_endo[block].size();
+      nze = blocks_derivatives[block].size();
      nze_other_endo = derivative_other_endo[block].size();
      nze_exo = derivative_exo[block].size();
      nze_exo_det = derivative_exo_det[block].size();
@@ -599,7 +599,7 @@ DynamicModel::writeModelEquationsOrdered_M(const string &basename) const
      else if (simulation_type == SOLVE_TWO_BOUNDARIES_COMPLETE || simulation_type == SOLVE_TWO_BOUNDARIES_SIMPLE)
        output << "  residual=zeros(" << block_mfs << ",y_kmin+periods);" << endl;
      if (simulation_type == EVALUATE_BACKWARD)
-        output << "  for it_ = (y_kmin+periods):y_kmin+1" << endl;
+        output << "  for it_ = (y_kmin+periods):-1:y_kmin+1" << endl;
      if (simulation_type == EVALUATE_FORWARD)
        output << "  for it_ = y_kmin+1:(y_kmin+periods)" << endl;
@@ -1691,7 +1691,7 @@ DynamicModel::writeDynamicJuliaFile(const string &basename) const
 }
 void
-DynamicModel::writeDynamicCFile(const string &basename, int order) const
+DynamicModel::writeDynamicCFile(const string &basename) const
 {
  filesystem::create_directories(basename + "/model/src");
  string filename = basename + "/model/src/dynamic.c";
@@ -1780,7 +1780,7 @@ DynamicModel::writeDynamicCFile(const string &basename, int order) const
                  << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
                  << "{" << endl
                  << "  /* Check that no derivatives of higher order than computed are being requested */" << endl
-                  << "  if (nlhs > " << order + 1 << ")" << endl
+                  << "  if (nlhs > " << computed_derivs_order + 1 << ")" << endl
                  << R"(    mexErrMsgTxt("Derivatives of higher order than computed have been requested");)" << endl
                  << "  /* Create a pointer to the input matrix y. */" << endl
                  << "  double *y = mxGetPr(prhs[0]);" << endl
@@ -2148,7 +2148,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                            << "  oo_.deterministic_simulation.block(blck_num).error = 0;" << endl
                            << "  oo_.deterministic_simulation.block(blck_num).iterations = 0;" << endl
                            << "  g1=[];g2=[];g3=[];" << endl
-                            << "  " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);" << endl
+                            << "  y = " << basename << ".block.dynamic_" << block + 1 << "(y, x, params, steady_state, 0, y_kmin, periods);" << endl
                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
                            << "    disp(['Inf or Nan value during the evaluation of block " << block <<"']);" << endl
@@ -2178,7 +2178,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                            << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
                            << ", y, x, params, steady_state, y_index, " << nze
                            << ", options_.periods, " << blocks_linear[block]
-                            << ", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);" << endl
+                            << ", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0, M_, options_, oo_);" << endl
                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
                            << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);" << endl
@@ -2209,7 +2209,7 @@ DynamicModel::writeSparseDynamicMFile(const string &basename) const
                            << "  y = solve_one_boundary('" << basename << ".block.dynamic_" <<  block + 1 << "'"
                            <<", y, x, params, steady_state, y_index, " << nze
                            <<", options_.periods, " << blocks_linear[block]
-                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0);" << endl
+                            <<", blck_num, y_kmin, options_.simul.maxit, options_.solve_tolf, options_.slowc, " << cutoff << ", options_.stack_solve_algo, 1, 1, 0, M_, options_, oo_);" << endl
                            << "  tmp = y(:,M_.block_structure.block(" << block + 1 << ").variable);" << endl
                            << "  if any(isnan(tmp) | isinf(tmp))" << endl
                            << "    disp(['Inf or Nan value during the resolution of block " << block <<"']);" << endl
@@ -3080,7 +3080,7 @@ DynamicModel::includeExcludeEquations(const string &eqs, bool exclude_eqs)
 }
 void
-DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const
+DynamicModel::writeOutput(ostream &output, const string &basename, bool block_decomposition, bool linear_decomposition, bool byte_code, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const
 {
  /* Writing initialisation for M_.lead_lag_incidence matrix
     M_.lead_lag_incidence is a matrix with as many columns as there are
@@ -3177,6 +3177,14 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
      for (size_t i = 0; i < temporary_terms_derivatives.size(); i++)
        output << temporary_terms_derivatives[i].size() + (i == 0 ? temporary_terms_mlv.size() : 0) << "; ";
      output << "];" << endl;
+      /* Write mapping between model local variables and indices in the temporary
+         terms vector (dynare#1722) */
+      output << modstruct << "model_local_variables_dynamic_tt_idxs = {" << endl;
+      for (auto [mlv, value] : temporary_terms_mlv)
+        output << "  '" << symbol_table.getName(mlv->symb_id) << "', "
+               << temporary_terms_idxs.at(mlv)+1 << ';' << endl;
+      output << "};" << endl;
    }
  // Write equation tags
@@ -3689,7 +3697,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
  // Use -1 if the derivatives have not been computed
  output << modstruct << (julia ? "nnzderivatives" : "NNZDerivatives") << " = [";
  for (int i = 1; i < static_cast<int>(NNZDerivatives.size()); i++)
-    output << (i > order ? -1 : NNZDerivatives[i]) << "; ";
+    output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
  output << "];" << endl;
  // Write Pac Model Consistent Expectation parameter info
@@ -4634,7 +4642,7 @@ DynamicModel::walkPacParameters(const string &name)
            }
          if (lhs.first == -1)
            {
-              cerr << "walkPacParameters: error obtaining LHS varibale." << endl;
+              cerr << "walkPacParameters: error obtaining LHS variable." << endl;
              exit(EXIT_FAILURE);
            }
          if (ec_params_and_vars.second.empty() || ar_params_and_vars.empty())
@@ -4959,11 +4967,8 @@ DynamicModel::computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsO
  computeDerivatives(derivsOrder, vars);
  if (derivsOrder > 1)
-    {
+    for (const auto &[indices, d2] : derivatives[2])
-      hessian_computed = true;
+      nonzero_hessian_eqs.insert(indices[0]);
-      for (const auto &[indices, d2] : derivatives[2])
-        nonzero_hessian_eqs.insert(indices[0]);
-    }
  if (paramsDerivsOrder > 0)
    {
@@ -5442,7 +5447,7 @@ DynamicModel::collectBlockVariables()
 }
 void
-DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, int order, bool julia) const
+DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const
 {
  if (block && bytecode)
    writeModelEquationsCode_Block(basename, map_idx, linear_decomposition);
@@ -5456,7 +5461,7 @@ DynamicModel::writeDynamicFile(const string &basename, bool block, bool linear_d
    writeSparseDynamicMFile(basename);
  else if (use_dll)
    {
-      writeDynamicCFile(basename, order);
+      writeDynamicCFile(basename);
      compileDll(basename, "dynamic", mexext, matlabroot, dynareroot);
    }
  else if (julia)
@@ -5516,13 +5521,25 @@ DynamicModel::writeAuxVarRecursiveDefinitions(ostream &output, ExprNodeOutputTyp
    }
 }
+void
+DynamicModel::clearEquations()
+{
+  equations.clear();
+  equations_lineno.clear();
+  equation_tags.clear();
+  equation_tags_xref.clear();
+}
 void
 DynamicModel::replaceMyEquations(DynamicModel &dynamic_model) const
 {
-  dynamic_model.equations.clear();
+  dynamic_model.clearEquations();
  for (size_t i = 0; i < equations.size(); i++)
-    dynamic_model.addEquation(equations[i]->clone(dynamic_model),
+    dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i]);
-                              equations_lineno[i]);
+  dynamic_model.equation_tags = equation_tags;
+  dynamic_model.equation_tags_xref = equation_tags_xref;
 }
 void
@@ -5539,9 +5556,9 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
    }
  cout << "Ramsey Problem: added " << i << " Multipliers." << endl;
-  // Add Planner Objective to equations to include in computeDerivIDs
+  // Add Planner Objective to equations so that it appears in Lagrangian
  assert(static_model.equations.size() == 1);
-  addEquation(static_model.equations[0]->clone(*this), static_model.equations_lineno[0]);
+  addEquation(static_model.equations[0]->clone(*this), -1);
  // Get max endo lead and max endo lag
  set<pair<int, int>> dynvars;
@@ -5550,9 +5567,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
  for (auto &equation : equations)
    equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
-  for (const auto &dynvar : dynvars)
+  for (const auto &[symb_id, lag] : dynvars)
    {
-      int lag = dynvar.second;
      if (max_eq_lead < lag)
        max_eq_lead = lag;
      else if (-max_eq_lag > lag)
@@ -5584,21 +5600,47 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model)
                                      equations[i]->getNonZeroPartofEquation()->decreaseLeadsLags(lag)), lagrangian);
      }
-  equations.clear();
+  // Save line numbers and tags, see below
+  auto old_equations_lineno = equations_lineno;
+  auto old_equation_tags = equation_tags;
+  // Prepare derivation of the Lagrangian
+  clearEquations();
  addEquation(AddEqual(lagrangian, Zero), -1);
  computeDerivIDs();
-  //Compute derivatives and overwrite equations
+  /* Compute Lagrangian derivatives.
+     Also restore line numbers and tags for FOCs w.r.t. a Lagrange multiplier
+     (i.e. a FOC identical to an equation of the original model) */
  vector<expr_t> neweqs;
-  for (auto &it : deriv_id_table)
+  vector<int> neweqs_lineno;
-    // For all endogenous variables with zero lag
+  map<int, vector<pair<string, string>>> neweqs_tags;
-    if (symbol_table.getType(it.first.first) == SymbolType::endogenous && it.first.second == 0)
+  for (auto &[symb_id_and_lag, deriv_id] : deriv_id_table)
-      neweqs.push_back(AddEqual(equations[0]->getNonZeroPartofEquation()->getDerivative(it.second), Zero));
+    {
+      auto &[symb_id, lag] = symb_id_and_lag;
+      if (symbol_table.getType(symb_id) == SymbolType::endogenous && lag == 0)
+        {
+          neweqs.push_back(AddEqual(equations[0]->getNonZeroPartofEquation()->getDerivative(deriv_id), Zero));
+          if (int i = symbol_table.getEquationNumberForMultiplier(symb_id);
+              i != -1)
+            {
+              // This is a derivative w.r.t. a Lagrange multiplier
+              neweqs_lineno.push_back(old_equations_lineno[i]);
+              vector<pair<string, string>> tags;
+              for (auto &[j, tagpair] : old_equation_tags)
+                if (j == i)
+                  tags.emplace_back(tagpair);
+              neweqs_tags[neweqs.size()-1] = tags;
+            }
+          else
+            neweqs_lineno.push_back(-1);
+        }
+    }
-  // Add new equations
+  // Overwrite equations with the Lagrangian derivatives
-  equations.clear();
+  clearEquations();
-  for (auto &neweq : neweqs)
+  for (size_t i = 0; i < neweqs.size(); i++)
-    addEquation(neweq, -1);
+    addEquation(neweqs[i], neweqs_lineno[i], neweqs_tags[i]);
 }
 void
@@ -6404,6 +6446,13 @@ DynamicModel::substituteLeadLagInternal(AuxVarType type, bool deterministic_mode
 void
 DynamicModel::substituteAdl()
 {
+  /* Contrary to other substitution methods, we do the substitution in MLV
+     definitions here, instead of doing it at the ExprNode method level,
+     because otherwise this would substitute MLV in the original model (see
+     #65). */
+  for (auto &[id, definition] : local_variables_table)
+    definition = definition->substituteAdl();
  for (auto &equation : equations)
    equation = dynamic_cast<BinaryOpNode *>(equation->substituteAdl());
 }
@@ -6798,6 +6847,9 @@ DynamicModel::isChecksumMatching(const string &basename, bool block) const
 void
 DynamicModel::writeJsonOutput(ostream &output) const
 {
+  deriv_node_temp_terms_t tef_terms;
+  writeJsonModelLocalVariables(output, false, tef_terms);
+  output << ", ";
  writeJsonModelEquations(output, false);
  output << ", ";
  writeJsonXrefs(output);
@@ -6850,18 +6902,24 @@ void
 DynamicModel::writeJsonVariableMapping(ostream &output) const
 {
  output << R"("variable_mapping":[)" << endl;
-  int ii = 0;
+  for (auto it = variableMapping.begin(); it != variableMapping.end(); ++it)
-  int end_idx_map = static_cast<int>(variableMapping.size()-1);
-  for (const auto &variable : variableMapping)
    {
-      output << R"({"name": ")" << symbol_table.getName(variable.first) << R"(", "equations":[)";
+      if (it != variableMapping.begin())
-      int it = 0;
+        output << ", ";
-      int end_idx_eq = static_cast<int>(variable.second.size())-1;
+      auto [var, eqs] = *it;
-      for (const auto &equation : variable.second)
+      output << R"({"name": ")" << symbol_table.getName(var) << R"(", "equations":[)";
+      bool first_eq = true;
+      for (auto it2 = eqs.begin(); it2 != eqs.end(); ++it2)
        for (const auto &equation_tag : equation_tags)
-          if (equation_tag.first == equation && equation_tag.second.first == "name")
+          if (equation_tag.first == *it2 && equation_tag.second.first == "name")
-            output << R"(")" << equation_tag.second.second << (it++ == end_idx_eq ? R"("])" : R"(", )");
+            {
-      output << (ii++ == end_idx_map ? R"(})" : R"(},)") << endl;
+              if (first_eq)
+                first_eq = false;
+              else
+                output << ", ";
+              output << '"' << equation_tag.second.second << '"';
+            }
+      output << "]}" << endl;
    }
  output << "]";
 }
@@ -6991,7 +7049,7 @@ DynamicModel::writeJsonComputingPassOutput(ostream &output, bool writeDetails) c
  deriv_node_temp_terms_t tef_terms;
  temporary_terms_t temp_term_union;
-  writeJsonModelLocalVariables(model_local_vars_output, tef_terms);
+  writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
  writeJsonTemporaryTerms(temporary_terms_derivatives[0], temp_term_union, d_output[0], tef_terms, "");
  d_output[0] << ", ";
@@ -7079,7 +7137,7 @@ DynamicModel::writeJsonParamsDerivativesFile(ostream &output, bool writeDetails)
  ostringstream g3p_output; // 1st deriv. of 3rd deriv. matrix w.r.t. parameters
  deriv_node_temp_terms_t tef_terms;
-  writeJsonModelLocalVariables(model_local_vars_output, tef_terms);
+  writeJsonModelLocalVariables(model_local_vars_output, true, tef_terms);
  temporary_terms_t temp_term_union;
  for (const auto &it : params_derivs_temporary_terms)

--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -99,8 +99,6 @@ private:
  //! Nonzero equations in the Hessian
  set<int> nonzero_hessian_eqs;
-  //! Whether the hessian has actually been computed
-  bool hessian_computed{false};
  //! Number of columns of dynamic jacobian
  /*! Set by computeDerivID()s and computeDynJacobianCols() */
@@ -123,7 +121,7 @@ private:
  void writeDynamicJuliaFile(const string &dynamic_basename) const;
  //! Writes dynamic model file (C version)
  /*! \todo add third derivatives handling */
-  void writeDynamicCFile(const string &basename, int order) const;
+  void writeDynamicCFile(const string &basename) const;
  //! Writes dynamic model file when SparseDLL option is on
  void writeSparseDynamicMFile(const string &basename) const;
  //! Writes the dynamic model equations and its derivatives
@@ -327,7 +325,7 @@ public:
  void computingPass(bool jacobianExo, int derivsOrder, int paramsDerivsOrder,
                     const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll, bool bytecode, bool linear_decomposition);
  //! Writes model initialization and lead/lag incidence matrix to output
-  void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, int order, bool estimation_present, bool compute_xrefs, bool julia) const;
+  void writeOutput(ostream &output, const string &basename, bool block, bool linear_decomposition, bool byte_code, bool use_dll, bool estimation_present, bool compute_xrefs, bool julia) const;
  //! Write JSON AST
  void writeJsonAST(ostream &output) const;
@@ -363,7 +361,7 @@ public:
  inline bool
  isHessianComputed() const
  {
-    return hessian_computed;
+    return computed_derivs_order >= 2;
  }
  //! Returns equations that have non-zero second derivatives
  inline set<int>
@@ -415,7 +413,7 @@ public:
  void Write_Inf_To_Bin_File_Block(const string &basename,
                                   int num, int &u_count_int, bool &file_open, bool is_two_boundaries, bool linear_decomposition) const;
  //! Writes dynamic model file
-  void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, int order, bool julia) const;
+  void writeDynamicFile(const string &basename, bool block, bool linear_decomposition, bool bytecode, bool use_dll, const string &mexext, const filesystem::path &matlabroot, const filesystem::path &dynareroot, bool julia) const;
  //! Writes file containing parameters derivatives
  void writeParamsDerivativesFile(const string &basename, bool julia) const;
@@ -446,6 +444,10 @@ public:
  //! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
  void computeRamseyPolicyFOCs(const StaticModel &static_model);
+  //! Clears all equations
+  void clearEquations();
  //! Replaces the model equations in dynamic_model with those in this model
  void replaceMyEquations(DynamicModel &dynamic_model) const;

--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -80,14 +80,14 @@ class ParsingDriver;
 %token BVAR_REPLIC BYTECODE ALL_VALUES_REQUIRED PROPOSAL_DISTRIBUTION REALTIME VINTAGE
 %token CALIB_SMOOTHER CHANGE_TYPE CHECK CONDITIONAL_FORECAST CONDITIONAL_FORECAST_PATHS CONF_SIG CONSTANT CONTROLLED_VAREXO CORR CUTOFF CYCLE_REDUCTION LOGARITHMIC_REDUCTION
 %token COMMA CONSIDER_ALL_ENDOGENOUS CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION
-%token DATAFILE FILE SERIES DET_COND_FORECAST DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
+%token DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_TOL DR_LOGARITHMIC_REDUCTION_MAXITER DR_ALGO DROP DSAMPLE DYNASAVE DYNATYPE CALIBRATION DIFFERENTIATE_FORWARD_VARS
 %token END ENDVAL EQUAL ESTIMATION ESTIMATED_PARAMS ESTIMATED_PARAMS_BOUNDS ESTIMATED_PARAMS_INIT EXTENDED_PATH ENDOGENOUS_PRIOR EXPRESSION
 %token FILENAME DIRNAME FILTER_STEP_AHEAD FILTERED_VARS FIRST_OBS LAST_OBS SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED
-%token <string> FLOAT_NUMBER DATES
+%token <string> FALSE FLOAT_NUMBER DATES
 %token DEFAULT FIXED_POINT FLIP OPT_ALGO COMPILATION_SETUP COMPILER ADD_FLAGS SUBSTITUTE_FLAGS ADD_LIBS SUBSTITUTE_LIBS
 %token FORECAST K_ORDER_SOLVER INSTRUMENTS SHIFT MEAN STDEV VARIANCE MODE INTERVAL SHAPE DOMAINN
 %token GAMMA_PDF GRAPH GRAPH_FORMAT CONDITIONAL_VARIANCE_DECOMPOSITION NOCHECK STD
-%token HISTVAL HISTVAL_FILE HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HOMOTOPY_FORCE_CONTINUE HP_FILTER HP_NGRID HYBRID ONE_SIDED_HP_FILTER
+%token HISTVAL HISTVAL_FILE HOMOTOPY_SETUP HOMOTOPY_MODE HOMOTOPY_STEPS HOMOTOPY_FORCE_CONTINUE HP_FILTER HP_NGRID FILTERED_THEORETICAL_MOMENTS_GRID HYBRID ONE_SIDED_HP_FILTER
 %token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE BOUNDS JSCALE INIT INFILE INVARS
 %token <string> INT_NUMBER
 %token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS IRF_PLOT_THRESHOLD IRF_CALIBRATION
@@ -116,7 +116,7 @@ class ParsingDriver;
 %token STDERR STEADY STOCH_SIMUL SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME REALTIME_SHOCK_DECOMPOSITION CONDITIONAL UNCONDITIONAL
 %token TEX RAMSEY_MODEL RAMSEY_POLICY RAMSEY_CONSTRAINTS PLANNER_DISCOUNT PLANNER_DISCOUNT_LATEX_NAME
 %token DISCRETIONARY_POLICY DISCRETIONARY_TOL EVALUATE_PLANNER_OBJECTIVE
-%token <string> TEX_NAME
+%token <string> TEX_NAME TRUE
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED
 %token VALUES VAR VAREXO VAREXO_DET VARIABLE VAROBS VAREXOBS PREDETERMINED_VARIABLES VAR_EXPECTATION VAR_EXPECTATION_MODEL PLOT_SHOCK_DECOMPOSITION MODEL_LOCAL_VARIABLE
 %token WRITE_LATEX_DYNAMIC_MODEL WRITE_LATEX_STATIC_MODEL WRITE_LATEX_ORIGINAL_MODEL CROSSEQUATIONS COVARIANCE WRITE_LATEX_STEADY_STATE_MODEL
@@ -180,7 +180,7 @@ class ParsingDriver;
 %type <string> vec_of_vec_value vec_value_list date_expr number
 %type <string> vec_value_1 vec_value signed_inf signed_number_w_inf
 %type <string> range vec_value_w_inf vec_value_1_w_inf
-%type <string> integer_range signed_integer_range
+%type <string> integer_range signed_integer_range boolean
 %type <string> sub_sampling_options list_sub_sampling_option
 %type <SymbolType> change_type_arg
 %type <vector<string>> vec_str vec_str_1
@@ -310,7 +310,6 @@ statement : parameters
          | smm_estimation
          | shock_groups
          | init2shocks
-          | det_cond_forecast
          | var_expectation_model
          | compilation_setup
          ;
@@ -1436,16 +1435,6 @@ prior_posterior_function_options : o_function
                                 | o_sampling_draws
                                 ;
-det_cond_forecast : DET_COND_FORECAST '(' symbol ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3); }
-                  | DET_COND_FORECAST '(' symbol COMMA symbol ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3,$5); }
-                  | DET_COND_FORECAST '(' symbol COMMA simul_options_list ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3); }
-                  | DET_COND_FORECAST '(' symbol COMMA symbol COMMA simul_options_list ')' ';'
-                    { driver.det_cond_forecast_linear_decomposition($3,$5); }
-                  ;
 simul : SIMUL ';'
        { driver.simul(); }
      | SIMUL '(' simul_options_list ')' ';'
@@ -1515,6 +1504,7 @@ stoch_simul_primary_options : o_dr_algo
                            | o_diagonal_only
                            | o_hp_filter
                            | o_hp_ngrid
+                            | o_filtered_theoretical_moments_grid
                            | o_periods
                            | o_simul
                            | o_simul_seed
@@ -1596,6 +1586,10 @@ signed_number_w_inf : signed_inf
                    | signed_number
                    ;
+boolean : TRUE
+        | FALSE
+        ;
 estimated_params : ESTIMATED_PARAMS ';' estimated_list END ';' { driver.estimated_params(); };
 estimated_list : estimated_list estimated_elem
@@ -2105,6 +2099,8 @@ list_optim_option : QUOTED_STRING COMMA QUOTED_STRING
                    { driver.optim_options_string($1, $3); }
                  | QUOTED_STRING COMMA signed_number
                    { driver.optim_options_num($1, $3); }
+                  | QUOTED_STRING COMMA boolean
+                    { driver.optim_options_num($1, $3); }
                  ;
 optim_options : list_optim_option
@@ -2822,8 +2818,14 @@ shock_decomposition_option : o_parameter_set
                           | o_first_obs
                           | o_nobs
                           | o_init_state
-			   | o_forecast_type
+                           | o_forecast_type
                           | o_shock_decomposition_with_epilogue
+                           | o_prefilter
+                           | o_loglinear
+                           | o_diffuse_kalman_tol
+                           | o_diffuse_filter
+                           | o_xls_sheet
+                           | o_xls_range
                           ;
 realtime_shock_decomposition_options_list : realtime_shock_decomposition_option COMMA realtime_shock_decomposition_options_list
@@ -2992,6 +2994,8 @@ calib_smoother_option : o_filtered_vars
                      | o_diffuse_filter
                      | o_smoothed_state_uncertainty
                      | o_parameter_set
+                      | o_xls_sheet
+                      | o_xls_range
                      ;
 generate_irfs : GENERATE_IRFS ';' END ';'
@@ -3167,7 +3171,11 @@ o_nomoments : NOMOMENTS { driver.option_num("nomoments", "true"); };
 o_irf : IRF EQUAL INT_NUMBER { driver.option_num("irf", $3); };
 o_irf_shocks : IRF_SHOCKS EQUAL '(' symbol_list ')' { driver.option_symbol_list("irf_shocks"); };
 o_hp_filter : HP_FILTER EQUAL non_negative_number { driver.option_num("hp_filter", $3); };
-o_hp_ngrid : HP_NGRID EQUAL INT_NUMBER { driver.option_num("hp_ngrid", $3); };
+o_hp_ngrid : HP_NGRID EQUAL INT_NUMBER {
+                                         driver.warning("The 'hp_ngrid' option is deprecated. It has been superseded by the 'filtered_theoretical_moments_grid' option.");
+                                         driver.option_num("filtered_theoretical_moments_grid", $3);
+                                       };
+o_filtered_theoretical_moments_grid : FILTERED_THEORETICAL_MOMENTS_GRID EQUAL INT_NUMBER { driver.option_num("filtered_theoretical_moments_grid", $3); };
 o_one_sided_hp_filter : ONE_SIDED_HP_FILTER EQUAL non_negative_number { driver.option_num("one_sided_hp_filter", $3); };
 o_periods : PERIODS EQUAL INT_NUMBER { driver.option_num("periods", $3); };
 o_solver_periods : SOLVER_PERIODS EQUAL INT_NUMBER { driver.option_num("ep.periods", $3); };
@@ -3955,6 +3963,8 @@ symbol : NAME
       | NONE
       | DR
       | PRIOR
+       | TRUE
+       | FALSE
       ;

--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -71,7 +71,7 @@ string eofbuff;
 #define YY_USER_ACTION location_increment(yylloc, yytext);
 %}
-DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|w([1-9]{1}|[1-4][0-9]|5[0-2]))
+DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 %%
 /* Code put at the beginning of yylex() */
@@ -187,7 +187,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|w([1-9]{1}|[1-4][0-9]|5[0-2]))
 <INITIAL>smoother2histval {BEGIN DYNARE_STATEMENT; return token::SMOOTHER2HISTVAL;}
 <INITIAL>perfect_foresight_setup {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SETUP;}
 <INITIAL>perfect_foresight_solver {BEGIN DYNARE_STATEMENT; return token::PERFECT_FORESIGHT_SOLVER;}
-<INITIAL>det_cond_forecast {BEGIN DYNARE_STATEMENT; return token::DET_COND_FORECAST;}
 <INITIAL>compilation_setup {BEGIN DYNARE_STATEMENT; return token::COMPILATION_SETUP;}
 <DYNARE_STATEMENT>; {
@@ -765,6 +764,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|w([1-9]{1}|[1-4][0-9]|5[0-2]))
 <DYNARE_STATEMENT>one_sided_hp_filter {return token::ONE_SIDED_HP_FILTER;}
 <DYNARE_STATEMENT>bandpass_filter {return token::BANDPASS_FILTER;}
 <DYNARE_STATEMENT>hp_ngrid {return token::HP_NGRID;}
+<DYNARE_STATEMENT>filtered_theoretical_moments_grid {return token::FILTERED_THEORETICAL_MOMENTS_GRID;}
 <DYNARE_STATEMENT>simul_seed {return token::SIMUL_SEED;}
 <DYNARE_STATEMENT>qz_criterium {return token::QZ_CRITERIUM;}
 <DYNARE_STATEMENT>qz_zero_threshold {return token::QZ_ZERO_THRESHOLD;}
@@ -796,6 +796,14 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|w([1-9]{1}|[1-4][0-9]|5[0-2]))
 <DYNARE_STATEMENT>coefficients {return token::COEFFICIENTS;}
 <DYNARE_STATEMENT>variances {return token::VARIANCES;}
 <DYNARE_STATEMENT>equations {return token::EQUATIONS;}
+<DYNARE_STATEMENT>true {
+  yylval->build<string>(yytext);
+  return token::TRUE;
+}
+<DYNARE_STATEMENT>false {
+    yylval->build<string>(yytext);
+    return token::FALSE;
+}
 <DYNARE_STATEMENT,DYNARE_BLOCK>\. {return Dynare::parser::token_type (yytext[0]);}
 <DYNARE_STATEMENT>\\ {return Dynare::parser::token_type (yytext[0]);}

--- a/src/DynareMain.cc
+++ b/src/DynareMain.cc
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -49,16 +49,15 @@ void main2(stringstream &in, const string &basename, bool debug, bool clear_all,
           const string &exclude_eqs, const string &include_eqs,
           JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, bool jsonderivsimple,
           const string &mexext, const filesystem::path &matlabroot,
-           const filesystem::path &dynareroot, bool onlymodel);
+           const filesystem::path &dynareroot, bool onlymodel, bool gui);
 void main1(const string &filename, const string &basename, istream &modfile, bool debug, bool save_macro, string &save_macro_file,
-           bool no_line_macro, bool no_empty_line_macro, const vector<pair<string, string>> &defines,
+           bool line_macro, const vector<pair<string, string>> &defines, vector<filesystem::path> &paths, stringstream &macro_output);
-           vector<filesystem::path> &paths, stringstream &macro_output);
 void
 usage()
 {
-  cerr << "Dynare usage: dynare mod_file [debug] [noclearall] [onlyclearglobals] [savemacro[=macro_file]] [onlymacro] [nolinemacro] [noemptylinemacro] [notmpterms] [nolog] [warn_uninit]"
+  cerr << "Dynare usage: dynare mod_file [debug] [noclearall] [onlyclearglobals] [savemacro[=macro_file]] [onlymacro] [linemacro] [notmpterms] [nolog] [warn_uninit]"
       << " [console] [nograph] [nointeractive] [parallel[=cluster_name]] [conffile=parallel_config_path_and_filename] [parallel_slave_open_mode] [parallel_test]"
       << " [-D<variable>[=<value>]] [-I/path] [nostrict] [stochastic] [fast] [minimal_workspace] [compute_xrefs] [output=dynamic|first|second|third] [language=matlab|julia]"
       << " [params_derivs_order=0|1|2] [transform_unary_ops] [exclude_eqs=<equation_tag_list_or_file>] [include_eqs=<equation_tag_list_or_file>]"
@@ -135,8 +134,7 @@ main(int argc, char **argv)
  bool debug = false;
  bool no_tmp_terms = false;
  bool only_macro = false;
-  bool no_line_macro = false;
+  bool line_macro = false;
-  bool no_empty_line_macro = false;
  bool no_log = false;
  bool no_warn = false;
  int params_derivs_order = 2;
@@ -155,6 +153,7 @@ main(int argc, char **argv)
  bool minimal_workspace = false;
  bool compute_xrefs = false;
  bool transform_unary_ops = false;
+  bool gui = false;
  string exclude_eqs, include_eqs;
  vector<pair<string, string>> defines;
  vector<filesystem::path> paths;
@@ -207,10 +206,8 @@ main(int argc, char **argv)
              save_macro_file = s.substr(10);
            }
        }
-      else if (s == "nolinemacro")
+      else if (s == "linemacro")
-        no_line_macro = true;
+        line_macro = true;
-      else if (s == "noemptylinemacro")
-        no_empty_line_macro = true;
      else if (s == "notmpterms")
        no_tmp_terms = true;
      else if (s == "nolog")
@@ -402,6 +399,8 @@ main(int argc, char **argv)
        }
      else if (s == "onlymodel")
        onlymodel = true;
+      else if (s == "gui")
+        gui = true;
      else
        {
          cerr << "Unknown option: " << s << endl;
@@ -432,7 +431,7 @@ main(int argc, char **argv)
  // Do macro processing
  stringstream macro_output;
-  main1(filename, basename, modfile, debug, save_macro, save_macro_file, no_line_macro, no_empty_line_macro,
+  main1(filename, basename, modfile, debug, save_macro, save_macro_file, line_macro,
        defines, paths, macro_output);
  if (only_macro)
@@ -450,7 +449,7 @@ main(int argc, char **argv)
        parallel, config_file, warnings, nostrict, stochastic, check_model_changes, minimal_workspace,
        compute_xrefs, output_mode, language, params_derivs_order, transform_unary_ops, exclude_eqs, include_eqs,
        json, json_output_mode, onlyjson, jsonderivsimple,
-        mexext, matlabroot, dynareroot, onlymodel);
+        mexext, matlabroot, dynareroot, onlymodel, gui);
  return EXIT_SUCCESS;
 }
--- a/src/DynareMain1.cc
+++ b/src/DynareMain1.cc
 /*
- * Copyright © 2015-2019 Dynare Team
+ * Copyright © 2015-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -21,17 +21,17 @@
 #include <fstream>
 #include <filesystem>
 #include <algorithm>
+#include <regex>
 #include "macro/Driver.hh"
 void
 main1(const string &filename, const string &basename, istream &modfile, bool debug, bool save_macro, string &save_macro_file,
-      bool no_line_macro_arg, bool no_empty_line_macro, const vector<pair<string, string>> &defines,
+      bool line_macro, const vector<pair<string, string>> &defines, vector<filesystem::path> &paths, stringstream &macro_output)
-      vector<filesystem::path> &paths, stringstream &macro_output)
 {
  // Do macro processing
  macro::Environment env = macro::Environment();
-  macro::Driver m(env, no_line_macro_arg);
+  macro::Driver m(env);
  m.parse(filename, basename, modfile, macro_output, debug, defines, paths);
  if (save_macro)
    {
@@ -45,8 +45,9 @@ main1(const string &filename, const string &basename, istream &modfile, bool deb
        }
      string str(macro_output.str());
-      if (no_empty_line_macro)
+      if (!line_macro)
        {
+          str = regex_replace(str, regex(R"((^|\n)\s*@#line.*)"), "");
          auto compareNewline = [](char i, char j) { return i == '\n' && j == '\n'; };
          str.erase(0, str.find_first_not_of('\n'));
          str.erase(unique(str.begin(), str.end(), compareNewline), str.end());

--- a/src/DynareMain2.cc
+++ b/src/DynareMain2.cc
 /*
- * Copyright © 2008-2019 Dynare Team
+ * Copyright © 2008-2020  Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -35,7 +35,7 @@ main2(stringstream &in, const string &basename, bool debug, bool clear_all, bool
      const string &exclude_eqs, const string &include_eqs,
      JsonOutputPointType json, JsonFileOutputType json_output_mode, bool onlyjson, bool jsonderivsimple,
      const string &mexext, const filesystem::path &matlabroot,
-      const filesystem::path &dynareroot, bool onlymodel)
+      const filesystem::path &dynareroot, bool onlymodel, bool gui)
 {
  ParsingDriver p(warnings, nostrict);
@@ -71,7 +71,7 @@ main2(stringstream &in, const string &basename, bool debug, bool clear_all, bool
  else
    mod_file->writeOutputFiles(basename, clear_all, clear_global, no_log, no_warn, console, nograph,
                               nointeractive, config_file, check_model_changes, minimal_workspace, compute_xrefs,
-                               mexext, matlabroot, dynareroot, onlymodel);
+                               mexext, matlabroot, dynareroot, onlymodel, gui);
  cout << "Preprocessing completed." << endl;
 }
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
 /*
- * Copyright © 2007-2019 Dynare Team
+ * Copyright © 2007-2021 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -451,9 +451,7 @@ void
 NumConstNode::writeJsonAST(ostream &output) const
 {
  output << R"({"node_type" : "NumConstNode", "value" : )";
-  if (double testval = datatree.num_constants.getDouble(id); testval < 1.0 && testval > -1.0 && testval != 0.0)
+  output << std::stof(datatree.num_constants.get(id)) << "}";
-    output << "0";
-  output << datatree.num_constants.get(id) << "}";
 }
 void
@@ -889,6 +887,9 @@ VariableNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, t
 bool
 VariableNode::containsExternalFunction() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->containsExternalFunction();
  return false;
 }
@@ -976,9 +977,9 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
    {
      if (output_type == ExprNodeOutputType::latexDynamicSteadyStateOperator)
        output << R"(\bar)";
-      output << "{" << datatree.symbol_table.getTeXName(symb_id);
+      output << "{" << datatree.symbol_table.getTeXName(symb_id) << "}";
      if (output_type == ExprNodeOutputType::latexDynamicModel
-          && (type == SymbolType::endogenous || type == SymbolType::exogenous || type == SymbolType::exogenousDet || type == SymbolType::modelLocalVariable || type == SymbolType::trend || type == SymbolType::logTrend))
+          && (type == SymbolType::endogenous || type == SymbolType::exogenous || type == SymbolType::exogenousDet || type == SymbolType::trend || type == SymbolType::logTrend))
        {
          output << "_{t";
          if (lag != 0)
@@ -989,7 +990,6 @@ VariableNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
            }
          output << "}";
        }
-      output << "}";
      return;
    }
@@ -1248,6 +1248,9 @@ VariableNode::substituteStaticAuxiliaryVariable() const
 double
 VariableNode::eval(const eval_context_t &eval_context) const noexcept(false)
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->eval(eval_context);
  auto it = eval_context.find(symb_id);
  if (it == eval_context.end())
    throw EvalException();
@@ -1408,29 +1411,17 @@ VariableNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recur
    case SymbolType::logTrend:
      if (deriv_id == datatree.getDerivID(symb_id, lag))
        return datatree.One;
-      else
+      // If there is in the equation a recursive variable we could use a chaine rule derivation
+      else if (auto it = recursive_variables.find(datatree.getDerivID(symb_id, lag));
+               it != recursive_variables.end())
        {
-          //if there is in the equation a recursive variable we could use a chaine rule derivation
+          map<int, expr_t> recursive_vars2(recursive_variables);
-          if (auto it = recursive_variables.find(datatree.getDerivID(symb_id, lag));
+          recursive_vars2.erase(it->first);
-              it != recursive_variables.end())
+          return datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id, recursive_vars2));
-            {
-              if (auto it2 = derivatives.find(deriv_id);
-                  it2 != derivatives.end())
-                return it2->second;
-              else
-                {
-                  map<int, expr_t> recursive_vars2(recursive_variables);
-                  recursive_vars2.erase(it->first);
-                  //expr_t c = datatree.AddNonNegativeConstant("1");
-                  expr_t d = datatree.AddUMinus(it->second->getChainRuleDerivative(deriv_id, recursive_vars2));
-                  //d = datatree.AddTimes(c, d);
-                  derivatives[deriv_id] = d;
-                  return d;
-                }
-            }
-          else
-            return datatree.Zero;
        }
+      else
+        return datatree.Zero;
    case SymbolType::modelLocalVariable:
      return datatree.getLocalVariable(symb_id)->getChainRuleDerivative(deriv_id, recursive_variables);
    case SymbolType::modFileLocalVariable:
@@ -1476,10 +1467,12 @@ VariableNode::computeXrefs(EquationInfo &ei) const
    case SymbolType::parameter:
      ei.param.emplace(symb_id, 0);
      break;
+    case SymbolType::modFileLocalVariable:
+      datatree.getLocalVariable(symb_id)->computeXrefs(ei);
+      break;
    case SymbolType::trend:
    case SymbolType::logTrend:
    case SymbolType::modelLocalVariable:
-    case SymbolType::modFileLocalVariable:
    case SymbolType::statementDeclaredVariable:
    case SymbolType::unusedEndogenous:
    case SymbolType::externalFunction:
@@ -1644,6 +1637,9 @@ VariableNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
 int
 VariableNode::PacMaxLag(int lhs_symb_id) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->PacMaxLag(lhs_symb_id);
  if (lhs_symb_id == symb_id)
    return -lag;
  return 0;
@@ -1658,28 +1654,40 @@ VariableNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
 expr_t
 VariableNode::substituteAdl() const
 {
+  /* Do not recurse into model-local variables definition, rather do it at the
+     DynamicModel method level (see the comment there) */
  return const_cast<VariableNode *>(this);
 }
 expr_t
 VariableNode::substituteVarExpectation(const map<string, expr_t> &subst_table) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->substituteVarExpectation(subst_table);
  return const_cast<VariableNode *>(this);
 }
 void
 VariableNode::findDiffNodes(lag_equivalence_table_t &nodes) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    datatree.getLocalVariable(symb_id)->findDiffNodes(nodes);
 }
 void
 VariableNode::findUnaryOpNodesForAuxVarCreation(lag_equivalence_table_t &nodes) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    datatree.getLocalVariable(symb_id)->findUnaryOpNodesForAuxVarCreation(nodes);
 }
 int
 VariableNode::findTargetVariable(int lhs_symb_id) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->findTargetVariable(lhs_symb_id);
  return -1;
 }
@@ -1687,18 +1695,27 @@ expr_t
 VariableNode::substituteDiff(const lag_equivalence_table_t &nodes, subst_table_t &subst_table,
                             vector<BinaryOpNode *> &neweqs) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->substituteDiff(nodes, subst_table, neweqs);
  return const_cast<VariableNode *>(this);
 }
 expr_t
 VariableNode::substituteUnaryOpNodes(const lag_equivalence_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->substituteUnaryOpNodes(nodes, subst_table, neweqs);
  return const_cast<VariableNode *>(this);
 }
 expr_t
 VariableNode::substitutePacExpectation(const string &name, expr_t subexpr)
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->substitutePacExpectation(name, subexpr);
  return const_cast<VariableNode *>(this);
 }
@@ -1723,6 +1740,9 @@ VariableNode::decreaseLeadsLags(int n) const
 expr_t
 VariableNode::decreaseLeadsLagsPredeterminedVariables() const
 {
+  /* Do not recurse into model-local variables definitions, since MLVs are
+     already handled by DynamicModel::transformPredeterminedVariables().
+     This is also necessary because of #65. */
  if (datatree.symbol_table.isPredetermined(symb_id))
    return decreaseLeadsLags(1);
  else
@@ -1864,6 +1884,9 @@ VariableNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *
 expr_t
 VariableNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->substituteExpectation(subst_table, neweqs, partial_information_model);
  return const_cast<VariableNode *>(this);
 }
@@ -1922,12 +1945,18 @@ VariableNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
 bool
 VariableNode::containsPacExpectation(const string &pac_model_name) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->containsPacExpectation(pac_model_name);
  return false;
 }
 bool
 VariableNode::containsEndogenous() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->containsEndogenous();
  if (get_type() == SymbolType::endogenous)
    return true;
  else
@@ -1937,12 +1966,18 @@ VariableNode::containsEndogenous() const
 bool
 VariableNode::containsExogenous() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->containsExogenous();
  return get_type() == SymbolType::exogenous || get_type() == SymbolType::exogenousDet;
 }
 expr_t
 VariableNode::replaceTrendVar() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->replaceTrendVar();
  if (get_type() == SymbolType::trend)
    return datatree.One;
  else if (get_type() == SymbolType::logTrend)
@@ -1954,6 +1989,9 @@ VariableNode::replaceTrendVar() const
 expr_t
 VariableNode::detrend(int symb_id, bool log_trend, expr_t trend) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->detrend(symb_id, log_trend, trend);
  if (this->symb_id != symb_id)
    return const_cast<VariableNode *>(this);
@@ -1976,12 +2014,18 @@ VariableNode::detrend(int symb_id, bool log_trend, expr_t trend) const
 int
 VariableNode::countDiffs() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->countDiffs();
  return 0;
 }
 expr_t
 VariableNode::removeTrendLeadLag(const map<int, expr_t> &trend_symbols_map) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->removeTrendLeadLag(trend_symbols_map);
  if ((get_type() != SymbolType::trend && get_type() != SymbolType::logTrend) || lag == 0)
    return const_cast<VariableNode *>(this);
@@ -2027,24 +2071,36 @@ VariableNode::removeTrendLeadLag(const map<int, expr_t> &trend_symbols_map) cons
 bool
 VariableNode::isInStaticForm() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->isInStaticForm();
  return lag == 0;
 }
 bool
 VariableNode::isParamTimesEndogExpr() const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->isParamTimesEndogExpr();
  return false;
 }
 bool
 VariableNode::isVarModelReferenced(const string &model_info_name) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->isVarModelReferenced(model_info_name);
  return false;
 }
 void
 VariableNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
 {
+  if (get_type() == SymbolType::modelLocalVariable)
+    return datatree.getLocalVariable(symb_id)->getEndosAndMaxLags(model_endos_and_lags);
  if (get_type() == SymbolType::endogenous)
    if (string varname = datatree.symbol_table.getName(symb_id);
        model_endos_and_lags.find(varname) == model_endos_and_lags.end())
@@ -2062,6 +2118,9 @@ VariableNode::findConstantEquations(map<VariableNode *, NumConstNode *> &table)
 expr_t
 VariableNode::replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table) const
 {
+  /* Do not recurse into model-local variables definitions, since MLVs are
+     already handled by ModelTree::simplifyEquations().
+     This is also necessary because of #65. */
  for (auto &it : table)
    if (it.first->symb_id == symb_id)
      return it.second;
@@ -3079,7 +3138,7 @@ UnaryOpNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int>> &r
 pair<int, expr_t>
 UnaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>> &List_of_Op_RHS) const
 {
-  pair<bool, expr_t> res = arg->normalizeEquation(var_endo, List_of_Op_RHS);
+  pair<int, expr_t> res = arg->normalizeEquation(var_endo, List_of_Op_RHS);
  int is_endogenous_present = res.first;
  expr_t New_expr_t = res.second;
@@ -4930,7 +4989,7 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
    {
      if (op_code == BinaryOpcode::equal) /* The end of the normalization process :
                                             All the operations needed to normalize the equation are applied. */
-        for (int i = 0; i < static_cast<int>(List_of_Op_RHS1.size()); i++)
+        while (!List_of_Op_RHS1.empty())
          {
            tuple<int, expr_t, expr_t> it = List_of_Op_RHS1.back();
            List_of_Op_RHS1.pop_back();
@@ -4949,7 +5008,7 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
  else if (is_endogenous_present_2)
    {
      if (op_code == BinaryOpcode::equal)
-        for (int i = 0; i < static_cast<int>(List_of_Op_RHS2.size()); i++)
+        while (!List_of_Op_RHS2.empty())
          {
            tuple<int, expr_t, expr_t> it = List_of_Op_RHS2.back();
            List_of_Op_RHS2.pop_back();
@@ -4969,12 +5028,9 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
    {
    case BinaryOpcode::plus:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        {
+        return { 0, datatree.AddPlus(expr_t_1, expr_t_2) };
-          List_of_Op_RHS.emplace_back(static_cast<int>(BinaryOpcode::minus), datatree.AddPlus(expr_t_1, expr_t_2), nullptr);
-          return { 0, datatree.AddPlus(expr_t_1, expr_t_2) };
-        }
      else if (is_endogenous_present_1 && is_endogenous_present_2)
-        return { 1, nullptr };
+        return { 2, nullptr };
      else if (!is_endogenous_present_1 && is_endogenous_present_2)
        {
          List_of_Op_RHS.emplace_back(static_cast<int>(BinaryOpcode::minus), expr_t_1, nullptr);
@@ -4988,12 +5044,9 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
      break;
    case BinaryOpcode::minus:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
-        {
+        return { 0, datatree.AddMinus(expr_t_1, expr_t_2) };
-          List_of_Op_RHS.emplace_back(static_cast<int>(BinaryOpcode::minus), datatree.AddMinus(expr_t_1, expr_t_2), nullptr);
-          return { 0, datatree.AddMinus(expr_t_1, expr_t_2) };
-        }
      else if (is_endogenous_present_1 && is_endogenous_present_2)
-        return { 1, nullptr };
+        return { 2, nullptr };
      else if (!is_endogenous_present_1 && is_endogenous_present_2)
        {
          List_of_Op_RHS.emplace_back(static_cast<int>(UnaryOpcode::uminus), nullptr, nullptr);
@@ -5020,7 +5073,7 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
          return { 1, expr_t_2 };
        }
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::divide:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
@@ -5036,7 +5089,7 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
          return { 1, expr_t_2 };
        }
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::power:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
@@ -5078,49 +5131,49 @@ BinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>>
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddMax(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::min:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddMin(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::less:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddLess(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::greater:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddGreater(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::lessEqual:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddLessEqual(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::greaterEqual:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddGreaterEqual(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::equalEqual:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddEqualEqual(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    case BinaryOpcode::different:
      if (!is_endogenous_present_1 && !is_endogenous_present_2)
        return { 0, datatree.AddDifferent(expr_t_1, expr_t_2) };
      else
-        return { 1, nullptr };
+        return { 2, nullptr };
      break;
    default:
      cerr << "Binary operator not handled during the normalization process" << endl;
@@ -6501,7 +6554,7 @@ TrinaryOpNode::normalizeEquation(int var_endo, vector<tuple<int, expr_t, expr_t>
  if (!is_endogenous_present_1 && !is_endogenous_present_2 && !is_endogenous_present_3)
    return { 0, datatree.AddNormcdf(expr_t_1, expr_t_2, expr_t_3) };
  else
-    return { 1, nullptr };
+    return { 2, nullptr };
 }
 expr_t
@@ -7429,7 +7482,7 @@ AbstractExternalFunctionNode::normalizeEquation(int var_endo, vector<tuple<int,
  if (!present)
    return { 0, datatree.AddExternalFunction(symb_id, V_expr_t) };
  else
-    return { 1, nullptr };
+    return { 2, nullptr };
 }
 void
@@ -7667,9 +7720,18 @@ ExternalFunctionNode::writeJsonOutput(ostream &output,
      return;
    }
-  output << datatree.symbol_table.getName(symb_id) << "(";
+  try
-  writeJsonExternalFunctionArguments(output, temporary_terms, tef_terms, isdynamic);
+    {
-  output << ")";
+      int tef_idx = getIndxInTefTerms(symb_id, tef_terms);
+      output << "TEF_" << tef_idx;
+    }
+  catch (UnknownFunctionNameAndArgs &)
+    {
+      // When writing the JSON output at parsing pass, we don’t use TEF terms
+      output << datatree.symbol_table.getName(symb_id) << "(";
+      writeJsonExternalFunctionArguments(output, temporary_terms, tef_terms, isdynamic);
+      output << ")";
+    }
 }
 void

--- a/src/ModFile.cc
+++ b/src/ModFile.cc
 /*
- * Copyright © 2006-2019 Dynare Team
+ * Copyright © 2006-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -386,6 +386,8 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
     — except substituting out variables which we know are constant (they
       appear in an equation of the form: X = constant)
     — except adl operators which we always want expanded
+     — except predetermined variables (which must be handled before the call to
+       setLeadsLagsOrig(), see preprocessor#47)
     — except diff operators with a lead which have been expanded by
       DataTree:AddDiff()
  */
@@ -393,6 +395,8 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
  dynamic_model.includeExcludeEquations(include_eqs, false);
  dynamic_model.simplifyEquations();
  dynamic_model.substituteAdl();
+  if (symbol_table.predeterminedNbr() > 0)
+    dynamic_model.transformPredeterminedVariables();
  dynamic_model.setLeadsLagsOrig();
  original_model = dynamic_model;
  dynamic_model.expandEqTags();
@@ -501,9 +505,6 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
  dynamic_model.addEquationsForVar();
-  if (symbol_table.predeterminedNbr() > 0)
-    dynamic_model.transformPredeterminedVariables();
  // Create auxiliary vars for Expectation operator
  dynamic_model.substituteExpectation(mod_file_struct.partial_information);
@@ -609,7 +610,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
    }
  // And finally perform the substitutions
  dynamic_model.substituteVarExpectation(var_expectation_subst_table);
-  dynamic_model.createVariableMapping(original_model.equation_number());
+  dynamic_model.createVariableMapping(mod_file_struct.orig_eq_nbr +
+                                      (mod_file_struct.ramsey_model_present ?
+                                       mod_file_struct.ramsey_eq_nbr : 0));
  if (mod_file_struct.stoch_simul_present
      || mod_file_struct.estimation_present
@@ -779,7 +782,7 @@ ModFile::computingPass(bool no_tmp_terms, FileOutputType output, int params_deri
                  exit(EXIT_FAILURE);
                }
              int derivsOrder = max(mod_file_struct.order_option,
-                                    mod_file_struct.identification_order + 1);
+                                    mod_file_struct.identification_order + 1); // See preprocessor#40
              if (mod_file_struct.sensitivity_present || linear || output == FileOutputType::second)
                derivsOrder = max(derivsOrder, 2);
              if (mod_file_struct.estimation_analytic_derivation || output == FileOutputType::third)
@@ -839,7 +842,7 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
                          bool check_model_changes, bool minimal_workspace, bool compute_xrefs,
                          const string &mexext,
                          const filesystem::path &matlabroot,
-                          const filesystem::path &dynareroot, bool onlymodel) const
+                          const filesystem::path &dynareroot, bool onlymodel, bool gui) const
 {
  bool hasModelChanged = !dynamic_model.isChecksumMatching(basename, block) || !check_model_changes;
  if (hasModelChanged)
@@ -1018,13 +1021,13 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
  if (dynamic_model.equation_number() > 0)
    {
      if (linear_decomposition)
-        non_linear_equations_dynamic_model.writeOutput(mOutputFile, basename, block, true, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
+        non_linear_equations_dynamic_model.writeOutput(mOutputFile, basename, block, true, byte_code, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
-      dynamic_model.writeOutput(mOutputFile, basename, block, false, byte_code, use_dll, mod_file_struct.order_option, mod_file_struct.estimation_present, compute_xrefs, false);
+      dynamic_model.writeOutput(mOutputFile, basename, block, false, byte_code, use_dll, mod_file_struct.estimation_present, compute_xrefs, false);
      if (!no_static)
        static_model.writeOutput(mOutputFile, block);
    }
-  if (onlymodel)
+  if (onlymodel || gui)
    for (auto &statement : statements)
      {
        /* Special treatment for initval block: insert initial values for the
@@ -1054,6 +1057,10 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
        if (auto sgs = dynamic_cast<ShockGroupsStatement *>(statement.get()); sgs)
          sgs->writeOutput(mOutputFile, basename, minimal_workspace);
+        if (gui)
+          if (auto it = dynamic_cast<NativeStatement *>(statement.get()); it)
+            it->writeOutput(mOutputFile, basename, minimal_workspace);
      }
  else
    {
@@ -1126,11 +1133,11 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
          if (linear_decomposition)
            {
-              non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mexext, matlabroot, dynareroot, mod_file_struct.order_option, false);
+              non_linear_equations_dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll, mexext, matlabroot, dynareroot, false);
              non_linear_equations_dynamic_model.writeParamsDerivativesFile(basename, false);
            }
-          dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mexext, matlabroot, dynareroot, mod_file_struct.order_option, false);
+          dynamic_model.writeDynamicFile(basename, block, false, byte_code, use_dll, mexext, matlabroot, dynareroot, false);
          dynamic_model.writeParamsDerivativesFile(basename, false);
@@ -1245,7 +1252,6 @@ ModFile::writeExternalFilesJulia(const string &basename) const
  if (dynamic_model.equation_number() > 0)
    {
      dynamic_model.writeOutput(jlOutputFile, basename, false, false, false, false,
-                                mod_file_struct.order_option,
                                mod_file_struct.estimation_present, false, true);
      if (!no_static)
        {
@@ -1253,7 +1259,7 @@ ModFile::writeExternalFilesJulia(const string &basename) const
          static_model.writeParamsDerivativesFile(basename, true);
        }
      dynamic_model.writeDynamicFile(basename, block, linear_decomposition, byte_code, use_dll,
-                                     "", {}, {}, mod_file_struct.order_option, true);
+                                     "", {}, {}, true);
      dynamic_model.writeParamsDerivativesFile(basename, true);
    }
  steady_state_model.writeSteadyStateFile(basename, mod_file_struct.ramsey_model_present, true);

--- a/src/ModFile.hh
+++ b/src/ModFile.hh
@@ -168,7 +168,7 @@ public:
                        bool console, bool nograph, bool nointeractive, const ConfigFile &config_file,
                        bool check_model_changes, bool minimal_workspace, bool compute_xrefs,
                        const string &mexext, const filesystem::path &matlabroot,
-                        const filesystem::path &dynareroot, bool onlymodel) const;
+                        const filesystem::path &dynareroot, bool onlymodel, bool gui) const;
  void writeExternalFiles(const string &basename, LanguageOutputType language) const;
  void writeExternalFilesJulia(const string &basename) const;

--- a/src/ModelTree.cc
+++ b/src/ModelTree.cc
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -33,6 +33,8 @@
 # include <mach-o/dyld.h>
 #endif
+#include <regex>
 using namespace MFS;
 void
@@ -72,7 +74,7 @@ ModelTree::copyHelper(const ModelTree &m)
  for (const auto &it : m.temporary_terms)
    temporary_terms.insert(f(it));
  for (const auto &it : m.temporary_terms_mlv)
-    temporary_terms_mlv[f(it.first)] = f(it.second);
+    temporary_terms_mlv[dynamic_cast<VariableNode *>(f(it.first))] = f(it.second);
  for (const auto &it : m.temporary_terms_derivatives)
    temporary_terms_derivatives.push_back(convert_temporary_terms_t(it));
  for (const auto &it : m.temporary_terms_idxs)
@@ -110,6 +112,7 @@ ModelTree::ModelTree(const ModelTree &m) :
  equations_lineno{m.equations_lineno},
  equation_tags{m.equation_tags},
  equation_tags_xref{m.equation_tags_xref},
+  computed_derivs_order{m.computed_derivs_order},
  NNZDerivatives{m.NNZDerivatives},
  equation_reordered{m.equation_reordered},
  variable_reordered{m.variable_reordered},
@@ -136,6 +139,7 @@ ModelTree::operator=(const ModelTree &m)
  aux_equations.clear();
  equation_tags = m.equation_tags;
  equation_tags_xref = m.equation_tags_xref;
+  computed_derivs_order = m.computed_derivs_order;
  NNZDerivatives = m.NNZDerivatives;
  derivatives.clear();
@@ -1243,6 +1247,8 @@ ModelTree::computeDerivatives(int order, const set<int> &vars)
 {
  assert(order >= 1);
+  computed_derivs_order = order;
  // Do not shrink the vectors, since they have a minimal size of 4 (see constructor)
  derivatives.resize(max(static_cast<size_t>(order+1), derivatives.size()));
  NNZDerivatives.resize(max(static_cast<size_t>(order+1), NNZDerivatives.size()), 0);
@@ -1361,6 +1367,8 @@ ModelTree::writeModelLocalVariableTemporaryTerms(temporary_terms_t &temp_term_un
  for (auto &it : temporary_terms_mlv)
    {
+      it.second->writeExternalFunctionOutput(output, output_type, temp_term_union, tt_idxs, tef_terms);
      if (isJuliaOutput(output_type))
        output << "    @inbounds const ";
@@ -1608,7 +1616,7 @@ ModelTree::compileTemporaryTerms(ostream &code_file, unsigned int &instruction_n
 }
 void
-ModelTree::writeJsonModelLocalVariables(ostream &output, deriv_node_temp_terms_t &tef_terms) const
+ModelTree::writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, deriv_node_temp_terms_t &tef_terms) const
 {
  /* Collect all model local variables appearing in equations, and print only
     them. Printing unused model local variables can lead to a crash (see
@@ -1632,23 +1640,24 @@ ModelTree::writeJsonModelLocalVariables(ostream &output, deriv_node_temp_terms_t
          printed = true;
        int id = it;
-        vector<string> efout;
        expr_t value = local_variables_table.find(id)->second;
-        value->writeJsonExternalFunctionOutput(efout, tt, tef_terms);
+        if (write_tef_terms)
-        for (auto it1 = efout.begin(); it1 != efout.end(); ++it1)
          {
-            if (it1 != efout.begin())
+            vector<string> efout;
+            value->writeJsonExternalFunctionOutput(efout, tt, tef_terms);
+            for (auto it1 = efout.begin(); it1 != efout.end(); ++it1)
+              {
+                if (it1 != efout.begin())
+                  output << ", ";
+                output << *it1;
+              }
+            if (!efout.empty())
              output << ", ";
-            output << *it1;
          }
-        if (!efout.empty())
+        output << R"({"variable": ")" << symbol_table.getName(id)
-          output << ", ";
+               << R"(", "value": ")";
-        /* We append underscores to avoid name clashes with "g1" or "oo_" (see
-           also VariableNode::writeOutput) */
-        output << R"({"variable": ")" << symbol_table.getName(id) << R"(__")"
-               << R"(, "value": ")";
        value->writeJsonOutput(output, tt, tef_terms);
        output << R"("})" << endl;
      }
@@ -1852,6 +1861,20 @@ ModelTree::writeLatexModelFile(const string &mod_basename, const string &latex_b
      content_output << "% Equation " << eq + 1 << endl;
      if (write_equation_tags)
        {
+          auto escape_special_latex_symbols
+            = [](string str)
+              {
+                const regex special_latex_chars (R"([&%$#_{}])");
+                const regex backslash (R"(\\)");
+                const regex tilde (R"(~)");
+                const regex carrot (R"(\^)");
+                const regex textbackslash (R"(\\textbackslash)");
+                str = regex_replace(str, backslash, R"(\textbackslash)");
+                str = regex_replace(str, special_latex_chars, R"(\$&)");
+                str = regex_replace(str, carrot, R"(\^{})");
+                str = regex_replace(str, tilde, R"(\textasciitilde{})");
+                return regex_replace(str, textbackslash, R"(\textbackslash{})");
+              };
          bool wrote_eq_tag = false;
          for (const auto & [tagged_eq, tag_pair] : equation_tags)
            if (tagged_eq == eq)
@@ -1861,16 +1884,16 @@ ModelTree::writeLatexModelFile(const string &mod_basename, const string &latex_b
                else
                  content_output << ", ";
-                content_output << tag_pair.first;
+                content_output << escape_special_latex_symbols(tag_pair.first);
                if (!(tag_pair.second.empty()))
-                  content_output << "= `" << tag_pair.second << "'";
+                  content_output << "= `" << escape_special_latex_symbols(tag_pair.second) << "'";
                wrote_eq_tag = true;
              }
          if (wrote_eq_tag)
-            content_output << "]";
+            content_output << "]" << endl;
        }
      content_output << R"(\begin{dmath})" << endl;
@@ -1999,22 +2022,26 @@ ModelTree::simplifyEquations()
 {
  size_t last_subst_table_size = 0;
  map<VariableNode *, NumConstNode *> subst_table;
-  findConstantEquations(subst_table);
+  // Equations with tags are excluded, in particular because of MCPs, see dynare#1697
+  findConstantEquationsWithoutTags(subst_table);
  while (subst_table.size() != last_subst_table_size)
    {
      last_subst_table_size = subst_table.size();
+      for (auto &[id, definition] : local_variables_table)
+        definition = definition->replaceVarsInEquation(subst_table);
      for (auto &equation : equations)
        equation = dynamic_cast<BinaryOpNode *>(equation->replaceVarsInEquation(subst_table));
      subst_table.clear();
-      findConstantEquations(subst_table);
+      findConstantEquationsWithoutTags(subst_table);
    }
 }
 void
-ModelTree::findConstantEquations(map<VariableNode *, NumConstNode *> &subst_table) const
+ModelTree::findConstantEquationsWithoutTags(map<VariableNode *, NumConstNode *> &subst_table) const
 {
-  for (auto &equation : equations)
+  for (size_t i = 0; i < equations.size(); i++)
-    equation->findConstantEquations(subst_table);
+    if (getEquationTags(i).empty())
+      equations[i]->findConstantEquations(subst_table);
 }
 void

--- a/src/ModelTree.hh
+++ b/src/ModelTree.hh
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -77,6 +77,8 @@ protected:
   * The following structures keep track of the model equations and must all be updated
   * when adding or removing an equation. Hence, if a new parallel structure is added
   * in the future, it must be maintained whereever these structures are updated
+   * See in particular methods clearEquations(), replaceMyEquations() and
+   * computeRamseyPolicyFOCs() of DynamicModel class.
   * NB: This message added with the introduction of the `exclude_eqs` option, hence
   *     that's a place to update future structures.
   */
@@ -104,6 +106,9 @@ protected:
      tests/optimal_policy/nk_ramsey_expectation.mod */
  vector<BinaryOpNode *> aux_equations;
+  //! Maximum order at which (endogenous) derivatives have been computed
+  int computed_derivs_order{0};
  //! Stores derivatives
  /*! Index 0 is not used, index 1 contains first derivatives, ...
     For each derivation order, stores a map whose key is a vector of integer: the
@@ -131,7 +136,7 @@ protected:
  //! Used model local variables, that will be treated as temporary terms
  /*! See the comments in ModelTree::computeTemporaryTerms() */
-  map<expr_t, expr_t, ExprNodeLess> temporary_terms_mlv;
+  map<VariableNode *, expr_t, ExprNodeLess> temporary_terms_mlv;
  //! Temporary terms for residuals and derivatives
  /*! Index 0 is temp. terms of residuals, index 1 for first derivatives, ... */
@@ -200,7 +205,9 @@ protected:
  //! if residuals = true, we are writing the dynamic/static model.
  //! Otherwise, just the model equations (with line numbers, no tmp terms)
  void writeJsonModelEquations(ostream &output, bool residuals) const;
-  void writeJsonModelLocalVariables(ostream &output, deriv_node_temp_terms_t &tef_terms) const;
+  /* Writes JSON model local variables.
+     Optionally put the external function variable calls into TEF terms */
+  void writeJsonModelLocalVariables(ostream &output, bool write_tef_terms, deriv_node_temp_terms_t &tef_terms) const;
  //! Compiles model equations
  void compileModelEquations(ostream &code_file, unsigned int &instruction_number, const temporary_terms_t &tt, const map_idx_t &map_idx, bool dynamic, bool steady_dynamic) const;
@@ -369,12 +376,14 @@ public:
  bool isNonstationary(int symb_id) const;
  void set_cutoff_to_zero();
  //! Simplify model equations: if a variable is equal to a constant, replace that variable elsewhere in the model
+  /*! Equations with tags are excluded, in particular because of MCPs, see
+      dynare#1697 */
  void simplifyEquations();
  /*! Reorder auxiliary variables so that they appear in recursive order in
      set_auxiliary_variables.m and dynamic_set_auxiliary_series.m */
  void reorderAuxiliaryEquations();
-  //! Find equations where variable is equal to a constant
+  //! Find equations of the form “variable=constant”, excluding equations with tags
-  void findConstantEquations(map<VariableNode *, NumConstNode *> &subst_table) const;
+  void findConstantEquationsWithoutTags(map<VariableNode *, NumConstNode *> &subst_table) const;
  //! Helper for writing the Jacobian elements in MATLAB and C
  /*! Writes either (i+1,j+1) or [i+j*no_eq] */
  void jacobianHelper(ostream &output, int eq_nb, int col_nb, ExprNodeOutputType output_type) const;

--- a/src/NumericalInitialization.cc
+++ b/src/NumericalInitialization.cc
 /*
- * Copyright © 2003-2019 Dynare Team
+ * Copyright © 2003-2020 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -39,6 +39,10 @@ InitParamStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolid
 {
  if (symbol_table.getName(symb_id) == "dsge_prior_weight")
    mod_file_struct.dsge_prior_weight_initialized = true;
+  // Needed for the workaround discussed in dynare#1173
+  if (symbol_table.getName(symb_id) == "optimal_policy_discount_factor")
+    param_value->collectVariables(SymbolType::parameter, mod_file_struct.parameters_in_planner_discount);
 }
 void
@@ -593,7 +597,7 @@ LoadParamsAndSteadyStateStatement::writeOutput(ostream &output, const string &ba
 void
 LoadParamsAndSteadyStateStatement::writeJsonOutput(ostream &output) const
 {
-  output << R"({"statementName": "load_params_and_steady_state")"
+  output << R"({"statementName": "load_params_and_steady_state",)"
         << R"("values": [)";
  for (auto it = content.begin(); it != content.end(); ++it)
    {
No results found