Compare revisions

d4bb0dce · bdb8ef99 · 329bdec3 · 8bfe3ef4 · 48863996 · 92a85273
--- a/.clang-format
+++ b/.clang-format
@@ -11,9 +11,12 @@ BreakInheritanceList: AfterColon
 Cpp11BracedListStyle: true
 DeriveLineEnding: false
 IndentPPDirectives: AfterHash
+InsertNewlineAtEOF: true
 PackConstructorInitializers: NextLine
 PPIndentWidth: 1
 PointerAlignment: Left
+# RemoveParentheses: ReturnStatement
+# RemoveSemicolon: true
 SpaceAfterTemplateKeyword: false
 SpaceBeforeParens: ControlStatements
 SpaceBeforeCpp11BracedList: true
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
 variables:
  TERM: linux
-  MINGW32_BOOST_VERSION: 1.84.0-1
+  MINGW64_BOOST_VERSION: 1.88.0-2
-  MINGW64_BOOST_VERSION: 1.84.0-1
  WGET_OPTIONS: '--no-verbose --no-use-server-timestamps --retry-connrefused --retry-on-host-error'
  # To ensure that "false && true" fails, see https://gitlab.com/gitlab-org/gitlab-runner/-/issues/25394#note_412609647
  FF_ENABLE_BASH_EXIT_CODE_CHECK: 'true'
@@ -72,3 +71,16 @@ test_clang_format:
    - meson setup build-clang-format
    - ninja -C build-clang-format clang-format-check
  needs: []
+test_clang_tidy:
+  stage: test
+  script:
+    # Hack needed for meson < 1.6.0 which only looks for unversioned clang-tidy
+    - mkdir -p ~/.local/bin && ln -s /usr/bin/clang-tidy-19 ~/.local/bin/clang-tidy
+    - export PATH="$HOME/.local/bin:$PATH"
+    - meson setup build-clang-tidy
+    # Generate Flex and Bison files
+    - meson compile -C build-clang-tidy
+    - ninja -C build-clang-tidy clang-tidy
+  needs: []
+  when: manual
--- a/doc/macroprocessor/macroprocessor.tex
+++ b/doc/macroprocessor/macroprocessor.tex
@@ -15,7 +15,7 @@
 \pgfdeclareimage[height=0.8cm]{logo}{dlogo}
 \institute[Dynare Team]{\pgfuseimage{logo}}
-\date{23 May 2023}
+\date{31 May 2024}
 \AtBeginSection[]
 {
@@ -156,9 +156,9 @@
    \item comparison operators: \texttt{< > <= >= == !=}
    \item logical operators: \verb+&& || !+
    \item range with unit increment: \texttt{1:4} is equivalent to
-      real array \texttt{[1, 2, 3, 4]}. (NB: \texttt{[1:4]} is equivalent to an
+      real array \texttt{[1, 2, 3, 4]} \\ (NB: \texttt{[1:4]} is equivalent to an
      array containing an array of reals, \textit{i.e.} \texttt{[[1, 2, 3, 4]]})
-    \item range with user-defined increment: \texttt{4:-1.1:-1} is equivalent to real array \texttt{[4, 2.9, 1.8, 0.7, -0.4]}.
+    \item range with user-defined increment: \\ \texttt{4:-1.1:-1} is equivalent to real array \texttt{[4, 2.9, 1.8, 0.7, -0.4]}
    \end{itemize}
  \end{block}
@@ -314,7 +314,7 @@ Then \texttt{distance(3, 4)} will be equivalent to \texttt{5}.
 \end{frame}
 \begin{frame}[fragile=singleslide]
-  \frametitle{Defining macro-variables}
+  \frametitle{Defining macro-variables (1/2)}
  The value of a macro-variable can be defined with the \verb+@#define+
  directive.
@@ -335,9 +335,31 @@ Then \texttt{distance(3, 4)} will be equivalent to \texttt{5}.
 @#define t = ("US" in w)    // Equals true
 \end{verbatim}
  \end{block}
-  NB: You can define macro variables on the Dynare command line by using the \texttt{-D} option
 \end{frame}
+\begin{frame}[fragile=singleslide]
+  \frametitle{Defining macro-variables (2/2)}
+  Macro variables can also be defined on the Dynare command line by using the
+  \texttt{-D} option, for easily switching between different flavours of a model.
+  \begin{block}{Example 1}
+\begin{verbatim}
+dynare myfile.mod -Dx=5
+\end{verbatim}
+The macro-variable \texttt{x} will be equal to \texttt{5} when running \texttt{myfile.mod}.
+\end{block}
+\begin{block}{Example 2}
+  Use single quotes around the \texttt{-D} option when there are spaces or
+  special characters in the variable definition.
+\begin{verbatim}
+dynare myfile.mod '-DA=[ i in [1,2,3] when i > 1 ]'
+\end{verbatim}
+The macro-variable \texttt{A} will be equal to \texttt{[2,3]} when running \texttt{myfile.mod}.
+\end{block}
+\end{frame}
 \begin{frame}[fragile=singleslide]
  \frametitle{Expression substitution}
  \framesubtitle{Dummy example}
@@ -611,7 +633,7 @@ end;
 \begin{frame}
  \frametitle{Macro-related command line options}
  \begin{itemize}
-  \item \texttt{savemacro}: Useful for debugging or learning purposes, saves the output of the macro processor. If your \texttt{.mod} file is called \texttt{file.mod}, the output is saved to \texttt{file-macroexp.mod}.
+  \item \texttt{savemacro}: Useful for debugging or learning purposes, saves the output of the macro processor. If your \texttt{.mod} file is called \texttt{file.mod}, the output is saved to \texttt{file\_macroexp.mod}.
  \item NB: \texttt{savemacro=filename} allows a user-defined file name
  \item \texttt{linemacro}: In the output of \texttt{savemacro}, print line numbers where the macro directives were placed.
  \item \texttt{onlymacro}: Stops processing after the macro processing step.
@@ -866,7 +888,7 @@ rhos = [ 0.8, 0.9, 1];
    \ccbysa
    \column{0.71\textwidth}
    \tiny
-    Copyright © 2008-2023 Dynare Team \\
+    Copyright © 2008-2024 Dynare Team \\
    License: \href{http://creativecommons.org/licenses/by-sa/4.0/}{Creative
      Commons Attribution-ShareAlike 4.0}
  \end{columns}

--- a/scripts/homebrew-native-arm64.ini
+++ b/scripts/homebrew-native-arm64.ini
 # Meson native file for compiling under Homebrew / arm64
 [binaries]
-cpp = '/opt/homebrew/bin/g++-13'
+cpp = '/opt/homebrew/bin/g++-15'
 flex = '/opt/homebrew/opt/flex/bin/flex'
 bison = '/opt/homebrew/opt/bison/bin/bison'

--- a/scripts/homebrew-native-x86_64.ini
+++ b/scripts/homebrew-native-x86_64.ini
 # Meson native file for compiling under Homebrew / x86-64
 [binaries]
-cpp = '/usr/local/bin/g++-13'
+cpp = '/usr/local/bin/g++-15'
 flex = '/usr/local/opt/flex/bin/flex'
 bison = '/usr/local/opt/bison/bin/bison'

--- a/scripts/wasm.ini
+++ b/scripts/wasm.ini
+# Meson cross file for creating a WebAssembly version of the preprocessor.
+#
+# Requires emscripten to be installed.
+# Was successfully tested with emscripten 3.1.69 installed through emsdk
+# tool, as described on: https://emscripten.org/docs/getting_started/downloads.html
+# Don’t forget to source script snippet in current shell before running meson.
+#
+# Compilation creates a .wasm and .js wrapper under <builddir>/src/
+#
+# Can be run locally with node.js using:
+#  node dynare-preprocessor.js file.mod
+# NB: a version of node.js is shipped with emscripten (under the node/
+# subdirectory), but another version should also work.
+[binaries]
+cpp = 'em++'
+[host_machine]
+system = 'emscripten'
+# Could be changed to wasm64 if 4GB memory constraint is hit
+# Some background: https://v8.dev/blog/4gb-wasm-memory
+cpu_family = 'wasm32'
+cpu = 'wasm32'
+endian = 'little'
+[built-in options]
+# Never do a debug build, because otherwise the lack of optimisations can
+# overflow the memory capacities.
+buildtype = 'release'
+# The -fexceptions flag (for both compilation and linking) is needed for an
+# unknown reason (C++ compilers are supposed to always add exception support).
+# The -Wno-unqualified-std-cast-call flag removes many warnings about “move”
+# not being qualified with “std::” namespace.
+# The -fexperimental-library flag is needed to get std::jthread support (it was
+# supposed to no longer be necessary for LLVM 20, but for some reason we still
+# need it).
+cpp_args = [ '-fexceptions', '-Wno-unqualified-std-cast-call', '-fexperimental-library' ]
+# NODERAWFS=1 is needed for accessing the local filesystem
+cpp_link_args = [ '-s', 'NODERAWFS=1', '-fexceptions' ]
+[properties]
+# It’s necessary to use a different copy of Boost than the one under
+# /usr/include, because otherwise GCC headers confuse Clang
+boost_root = '/tmp/boost_1_86_0'
--- a/src/Bytecode.cc
+++ b/src/Bytecode.cc
 /*
- * Copyright © 2022-2023 Dynare Team
+ * Copyright © 2022-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -82,8 +82,8 @@ operator<<(Writer& code_file, const FBEGINBLOCK& instr)
  write_member(instr.type);
  for (int i = 0; i < instr.size; i++)
    {
-      write_member(instr.variable[i]);
+      write_member(instr.variables[i]);
-      write_member(instr.equation[i]);
+      write_member(instr.equations[i]);
    }
  if (instr.type == BlockSimulationType::solveTwoBoundariesSimple
      || instr.type == BlockSimulationType::solveTwoBoundariesComplete
@@ -97,8 +97,8 @@ operator<<(Writer& code_file, const FBEGINBLOCK& instr)
  write_member(instr.det_exo_size);
  write_member(instr.exo_size);
-  for_each_n(instr.det_exogenous.begin(), instr.det_exo_size, write_member);
+  ranges::for_each_n(instr.det_exogenous.begin(), instr.det_exo_size, write_member);
-  for_each_n(instr.exogenous.begin(), instr.exo_size, write_member);
+  ranges::for_each_n(instr.exogenous.begin(), instr.exo_size, write_member);
  return code_file;
 }

--- a/src/Bytecode.hh
+++ b/src/Bytecode.hh
@@ -20,6 +20,7 @@
 #ifndef BYTECODE_HH
 #define BYTECODE_HH
+#include <concepts>
 #include <filesystem>
 #include <fstream>
 #include <ios>
@@ -121,11 +122,6 @@ enum class ExternalFunctionCallType
  numericalSecondDerivative
 };
-struct Block_contain_type
-{
-  int Equation, Variable, Own_Derivative;
-};
 class Writer;
 struct Instruction
@@ -144,6 +140,9 @@ protected:
  ~Instruction() = default;
 };
+template<typename T>
+concept IsInstruction = derived_from<T, Instruction>;
 struct FLDZ final : public Instruction
 {
  FLDZ() : Instruction {Tag::FLDZ}
@@ -431,7 +430,7 @@ struct FSTPV final : public Instruction
 class FCALL final : public Instruction
 {
-  template<typename B>
+  template<IsInstruction B>
  friend Writer& operator<<(Writer& code_file, const B& instr);
 private:
@@ -607,18 +606,17 @@ public:
 class FBEGINBLOCK final : public Instruction
 {
-  template<typename B>
+  template<IsInstruction B>
  friend Writer& operator<<(Writer& code_file, const B& instr);
 private:
  int size {0};
  BlockSimulationType type;
-  vector<int> variable;
+  vector<int> variables;
-  vector<int> equation;
+  vector<int> equations;
  vector<int> exogenous;
  vector<int> det_exogenous;
  bool is_linear {false};
-  vector<Block_contain_type> Block_Contain_;
  int u_count_int {0};
  int nb_col_jacob {0};
  int det_exo_size, exo_size;
@@ -628,16 +626,16 @@ public:
     sense when there is no block-decomposition, since there is no provision for
     derivatives w.r.t. endogenous not belonging to the block) */
  FBEGINBLOCK(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
-              const vector<int>& variable_arg, const vector<int>& equation_arg, bool is_linear_arg,
+              const vector<int>& variables_arg, const vector<int>& equations_arg,
-              int u_count_int_arg, int nb_col_jacob_arg, int det_exo_size_arg, int exo_size_arg,
+              bool is_linear_arg, int u_count_int_arg, int nb_col_jacob_arg, int det_exo_size_arg,
-              vector<int> det_exogenous_arg, vector<int> exogenous_arg) :
+              int exo_size_arg, vector<int> det_exogenous_arg, vector<int> exogenous_arg) :
      Instruction {Tag::FBEGINBLOCK},
      size {size_arg},
      type {type_arg},
-      variable {variable_arg.begin() + first_element,
+      variables {variables_arg.begin() + first_element,
-                variable_arg.begin() + (first_element + block_size)},
+                 variables_arg.begin() + (first_element + block_size)},
-      equation {equation_arg.begin() + first_element,
+      equations {equations_arg.begin() + first_element,
-                equation_arg.begin() + (first_element + block_size)},
+                 equations_arg.begin() + (first_element + block_size)},
      exogenous {move(exogenous_arg)},
      det_exogenous {move(det_exogenous_arg)},
      is_linear {is_linear_arg},
@@ -649,15 +647,15 @@ public:
  }
  // Constructor when derivatives w.r.t. exogenous are absent
  FBEGINBLOCK(int size_arg, BlockSimulationType type_arg, int first_element, int block_size,
-              const vector<int>& variable_arg, const vector<int>& equation_arg, bool is_linear_arg,
+              const vector<int>& variables_arg, const vector<int>& equations_arg,
-              int u_count_int_arg, int nb_col_jacob_arg) :
+              bool is_linear_arg, int u_count_int_arg, int nb_col_jacob_arg) :
      Instruction {Tag::FBEGINBLOCK},
      size {size_arg},
      type {type_arg},
-      variable {variable_arg.begin() + first_element,
+      variables {variables_arg.begin() + first_element,
-                variable_arg.begin() + (first_element + block_size)},
+                 variables_arg.begin() + (first_element + block_size)},
-      equation {equation_arg.begin() + first_element,
+      equations {equations_arg.begin() + first_element,
-                equation_arg.begin() + (first_element + block_size)},
+                 equations_arg.begin() + (first_element + block_size)},
      is_linear {is_linear_arg},
      u_count_int {u_count_int_arg},
      nb_col_jacob {nb_col_jacob_arg},
@@ -678,12 +676,12 @@ public:
    read_member(size);
    read_member(type);
+    variables.resize(size);
+    equations.resize(size);
    for (int i {0}; i < size; i++)
      {
-        Block_contain_type bc;
+        read_member(variables[i]);
-        read_member(bc.Variable);
+        read_member(equations[i]);
-        read_member(bc.Equation);
-        Block_Contain_.push_back(bc);
      }
    if (type == BlockSimulationType::solveTwoBoundariesSimple
        || type == BlockSimulationType::solveTwoBoundariesComplete
@@ -731,10 +729,15 @@ public:
  {
    return u_count_int;
  }
-  vector<Block_contain_type>
+  vector<int>
-  get_Block_Contain()
+  get_variables()
  {
-    return Block_Contain_;
+    return variables;
+  }
+  vector<int>
+  get_equations()
+  {
+    return equations;
  }
  int
  get_nb_col_jacob()
@@ -752,11 +755,6 @@ public:
    return det_exo_size;
  }
  vector<int>
-  get_endogenous()
-  {
-    return variable;
-  }
-  vector<int>
  get_exogenous()
  {
    return exogenous;
@@ -766,7 +764,7 @@ public:
 // Superclass of std::ofstream for writing a sequence of bytecode instructions
 class Writer : private ofstream
 {
-  template<typename B>
+  template<IsInstruction B>
  friend Writer& operator<<(Writer& code_file, const B& instr);
 private:
@@ -785,7 +783,7 @@ public:
     It is the responsibility of the caller to ensure that the new instruction
     occupies exactly as many bytes as the former one. */
  void
-  overwriteInstruction(int instruction_number, const auto& new_instruction)
+  overwriteInstruction(int instruction_number, const IsInstruction auto& new_instruction)
  {
    seekp(instructions_positions.at(instruction_number));
    *this << new_instruction;
@@ -796,7 +794,7 @@ public:
 // Overloads of operator<< for writing bytecode instructions
-template<typename B>
+template<IsInstruction B>
 Writer&
 operator<<(Writer& code_file, const B& instr)
 {

--- a/src/CommonEnums.hh
+++ b/src/CommonEnums.hh
 /*
- * Copyright © 2007-2023 Dynare Team
+ * Copyright © 2007-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -25,14 +25,17 @@
 * command */
 enum class SymbolType
 {
-  endogenous = 0,            //!< Endogenous
+  endogenous = 0,              // Endogenous (non-heterogeneous)
-  exogenous = 1,             //!< Exogenous
+  exogenous = 1,               // Exogenous (non-heterogeneous)
-  exogenousDet = 2,          //!< Exogenous deterministic
+  exogenousDet = 2,            // Exogenous deterministic (non-heterogeneous)
-  parameter = 4,             //!< Parameter
+  parameter = 4,               // Parameter (non-heterogeneous)
-  modelLocalVariable = 10,   //!< Local variable whose scope is model (pound expression)
+  heterogeneousEndogenous = 5, // Endogenous that is heterogeneous across some dimension
-  modFileLocalVariable = 11, //!< Local variable whose scope is mod file (model excluded)
+  heterogeneousExogenous = 6,  // Exogenous that is heterogeneous across some dimension
-  externalFunction = 12,     //!< External (user-defined) function
+  heterogeneousParameter = 7,  // Parameter that is heterogeneous across some dimension
-  trend = 13,                //!< Trend variable
+  modelLocalVariable = 10,     // Local variable whose scope is model (pound expression)
+  modFileLocalVariable = 11,   // Local variable whose scope is mod file (model excluded)
+  externalFunction = 12,       // External (user-defined) function
+  trend = 13,                  // Trend variable
  statementDeclaredVariable
  = 14, //!< Local variable assigned within a Statement (see subsample statement for example)
  logTrend = 15, //!< Log-trend variable
@@ -45,6 +48,13 @@ enum class SymbolType
  excludedVariable = 19 //!< Variable excluded via model_remove/var_remove/include_eqs/exclude_eqs
 };
+constexpr bool
+isHeterogeneous(SymbolType type)
+{
+  return type == SymbolType::heterogeneousEndogenous || type == SymbolType::heterogeneousExogenous
+         || type == SymbolType::heterogeneousParameter;
+}
 enum class UnaryOpcode
 {
  uminus,
@@ -75,7 +85,8 @@ enum class UnaryOpcode
  erf,
  erfc,
  diff,
-  adl
+  adl,
+  sum
 };
 enum class BinaryOpcode

--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -283,7 +283,9 @@ PerfectForesightWithExpectationErrorsSolverStatement::writeOutput(
    [[maybe_unused]] bool minimal_workspace) const
 {
  options_list.writeOutput(output);
-  output << "oo_ = perfect_foresight_with_expectation_errors_solver(M_, options_, oo_);" << endl;
+  output << "[oo_, Simulated_time_series] = perfect_foresight_with_expectation_errors_solver(M_, "
+            "options_, oo_);"
+         << endl;
 }
 void
@@ -508,90 +510,6 @@ RamseyModelStatement::writeJsonOutput(ostream& output) const
  output << "}";
 }
-RamseyConstraintsStatement::RamseyConstraintsStatement(const SymbolTable& symbol_table_arg,
-                                                       constraints_t constraints_arg) :
-    symbol_table {symbol_table_arg}, constraints {move(constraints_arg)}
-{
-}
-void
-RamseyConstraintsStatement::checkPass(ModFileStructure& mod_file_struct,
-                                      [[maybe_unused]] WarningConsolidation& warnings)
-{
-  mod_file_struct.ramsey_constraints_present = true;
-}
-void
-RamseyConstraintsStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
-                                        [[maybe_unused]] bool minimal_workspace) const
-{
-  output << "M_.ramsey_model_constraints = {" << endl;
-  for (bool printed_something {false}; const auto& it : constraints)
-    {
-      if (exchange(printed_something, true))
-        output << ", ";
-      output << "{" << it.endo + 1 << ", '";
-      switch (it.code)
-        {
-        case BinaryOpcode::less:
-          output << '<';
-          break;
-        case BinaryOpcode::greater:
-          output << '>';
-          break;
-        case BinaryOpcode::lessEqual:
-          output << "<=";
-          break;
-        case BinaryOpcode::greaterEqual:
-          output << ">=";
-          break;
-        default:
-          cerr << "Ramsey constraints: this shouldn't happen." << endl;
-          exit(EXIT_FAILURE);
-        }
-      output << "', '";
-      it.expression->writeOutput(output);
-      output << "'}" << endl;
-    }
-  output << "};" << endl;
-}
-void
-RamseyConstraintsStatement::writeJsonOutput(ostream& output) const
-{
-  output << R"({"statementName": "ramsey_constraints")"
-         << R"(, "ramsey_model_constraints": [)" << endl;
-  for (bool printed_something {false}; const auto& it : constraints)
-    {
-      if (exchange(printed_something, true))
-        output << ", ";
-      output << R"({"constraint": ")" << symbol_table.getName(it.endo) << " ";
-      switch (it.code)
-        {
-        case BinaryOpcode::less:
-          output << '<';
-          break;
-        case BinaryOpcode::greater:
-          output << '>';
-          break;
-        case BinaryOpcode::lessEqual:
-          output << "<=";
-          break;
-        case BinaryOpcode::greaterEqual:
-          output << ">=";
-          break;
-        default:
-          cerr << "Ramsey constraints: this shouldn't happen." << endl;
-          exit(EXIT_FAILURE);
-        }
-      output << " ";
-      it.expression->writeJsonOutput(output, {}, {});
-      output << R"("})" << endl;
-    }
-  output << "]" << endl;
-  output << "}";
-}
 RamseyPolicyStatement::RamseyPolicyStatement(SymbolList symbol_list_arg,
                                             OptionsList options_list_arg,
                                             const SymbolTable& symbol_table_arg) :
@@ -1004,10 +922,8 @@ EstimationStatement::checkPass(ModFileStructure& mod_file_struct, WarningConsoli
  /* 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(),
+  ranges::set_intersection(mod_file_struct.estimated_parameters,
-                   mod_file_struct.estimated_parameters.end(),
+                           mod_file_struct.parameters_in_planner_discount,
-                   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())
    {
@@ -1267,9 +1183,8 @@ AbstractEstimatedParamsStatement::commonCheckPass() const
      it.p4->collectVariables(SymbolType::parameter, used_params);
      it.jscale->collectVariables(SymbolType::parameter, used_params);
      vector<int> intersect;
-      set_intersection(declared_params.begin(), declared_params.end(), used_params.begin(),
+      ranges::set_intersection(declared_params, used_params, back_inserter(intersect));
-                       used_params.end(), back_inserter(intersect));
+      if (!intersect.empty())
-      if (intersect.size() > 0)
        {
          cerr << "ERROR: in `" << blockName() << "' block, the value of estimated parameter "
               << symbol_table.getName(intersect[0]) << " is used in the declaration for ";
@@ -1850,21 +1765,16 @@ DeterministicTrendsStatement::writeOutput(ostream& output, [[maybe_unused]] cons
                                          [[maybe_unused]] bool minimal_workspace) const
 {
  output << "options_.trend_coeff = {};" << endl;
-  for (const auto& trend_element : trend_elements)
+  for (const auto& [name, val] : trend_elements)
-    {
+    if (symbol_table.getType(name) == SymbolType::endogenous)
-      SymbolType type = symbol_table.getType(trend_element.first);
-      if (type == SymbolType::endogenous)
      {
-          output << "tmp1 = strmatch('" << trend_element.first << "',M_.endogenous_names,'exact');"
+        output << "tmp1 = strmatch('" << name << "',M_.endogenous_names,'exact');" << endl
-                 << endl;
+               << "options_.deterministic_trend_coeffs{tmp1} = '";
-          output << "options_.deterministic_trend_coeffs{tmp1} = '";
+        val->writeOutput(output);
-          trend_element.second->writeOutput(output);
        output << "';" << endl;
      }
    else
-        cerr << "Warning : Non-variable symbol used in deterministic_trends: "
+      cerr << "Warning: Non-variable symbol used in deterministic_trends: " << name << endl;
-             << trend_element.first << endl;
-    }
 }
 void
@@ -1872,20 +1782,17 @@ DeterministicTrendsStatement::writeJsonOutput(ostream& output) const
 {
  output << R"({"statementName": "deterministic_trends", )"
         << R"("trends" : {)";
-  for (bool printed_something {false}; const auto& trend_element : trend_elements)
+  for (bool printed_something {false}; const auto& [name, val] : trend_elements)
-    {
+    if (symbol_table.getType(name) == SymbolType::endogenous)
-      if (symbol_table.getType(trend_element.first) == SymbolType::endogenous)
      {
        if (exchange(printed_something, true))
          output << ", ";
-          output << R"(")" << trend_element.first << R"(": ")";
+        output << R"(")" << name << R"(": ")";
-          trend_element.second->writeJsonOutput(output, {}, {});
+        val->writeJsonOutput(output, {}, {});
        output << R"(")" << endl;
      }
    else
-        cerr << "Warning : Non-variable symbol used in deterministic_trends: "
+      cerr << "Warning: Non-variable symbol used in deterministic_trends: " << name << endl;
-             << trend_element.first << endl;
-    }
  output << "}"
         << "}";
 }
@@ -1901,21 +1808,16 @@ ObservationTrendsStatement::writeOutput(ostream& output, [[maybe_unused]] const
                                        [[maybe_unused]] bool minimal_workspace) const
 {
  output << "options_.trend_coeff = {};" << endl;
-  for (const auto& trend_element : trend_elements)
+  for (const auto& [name, val] : trend_elements)
-    {
+    if (symbol_table.getType(name) == SymbolType::endogenous)
-      SymbolType type = symbol_table.getType(trend_element.first);
-      if (type == SymbolType::endogenous)
      {
-          output << "tmp1 = strmatch('" << trend_element.first << "',options_.varobs,'exact');"
+        output << "tmp1 = strmatch('" << name << "',options_.varobs,'exact');" << endl
-                 << endl;
+               << "options_.trend_coeffs{tmp1} = '";
-          output << "options_.trend_coeffs{tmp1} = '";
+        val->writeOutput(output);
-          trend_element.second->writeOutput(output);
        output << "';" << endl;
      }
    else
-        cerr << "Warning : Non-variable symbol used in observation_trends: " << trend_element.first
+      cerr << "Warning: Non-variable symbol used in observation_trends: " << name << endl;
-             << endl;
-    }
 }
 void
@@ -1923,20 +1825,17 @@ ObservationTrendsStatement::writeJsonOutput(ostream& output) const
 {
  output << R"({"statementName": "observation_trends", )"
         << R"("trends" : {)";
-  for (bool printed_something {false}; const auto& trend_element : trend_elements)
+  for (bool printed_something {false}; const auto& [name, val] : trend_elements)
-    {
+    if (symbol_table.getType(name) == SymbolType::endogenous)
-      if (symbol_table.getType(trend_element.first) == SymbolType::endogenous)
      {
        if (exchange(printed_something, true))
          output << ", ";
-          output << R"(")" << trend_element.first << R"(": ")";
+        output << R"(")" << name << R"(": ")";
-          trend_element.second->writeJsonOutput(output, {}, {});
+        val->writeJsonOutput(output, {}, {});
        output << R"(")" << endl;
      }
    else
-        cerr << "Warning : Non-variable symbol used in observation_trends: " << trend_element.first
+      cerr << "Warning: Non-variable symbol used in observation_trends: " << name << endl;
-             << endl;
-    }
  output << "}"
         << "}";
 }
@@ -2014,7 +1913,7 @@ void
 OsrParamsStatement::checkPass(ModFileStructure& mod_file_struct, WarningConsolidation& warnings)
 {
  if (mod_file_struct.osr_params_present)
-    cerr << "WARNING: You have more than one osr_params statement in the .mod file." << endl;
+    warnings << "WARNING: You have more than one osr_params statement in the .mod file." << endl;
  mod_file_struct.osr_params_present = true;
  try
@@ -3816,8 +3715,7 @@ SvarGlobalIdentificationCheckStatement::writeJsonOutput(ostream& output) const
  output << R"({"statementName": "svar_global_identification"})";
 }
-SetTimeStatement::SetTimeStatement(OptionsList options_list_arg) :
+SetTimeStatement::SetTimeStatement(string period_arg) : period {move(period_arg)}
-    options_list {move(options_list_arg)}
 {
 }
@@ -3825,19 +3723,13 @@ void
 SetTimeStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                              [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output);
+  output << "options_.initial_period = " << period << endl;
 }
 void
 SetTimeStatement::writeJsonOutput(ostream& output) const
 {
-  output << R"({"statementName": "set_time")";
+  output << R"({"statementName": "set_time", "period": ")" << period << R"("})";
-  if (!options_list.empty())
-    {
-      output << ", ";
-      options_list.writeJsonOutput(output);
-    }
-  output << "}";
 }
 EstimationDataStatement::EstimationDataStatement(OptionsList options_list_arg) :
@@ -4315,18 +4207,13 @@ BasicPriorStatement::checkPass([[maybe_unused]] ModFileStructure& mod_file_struc
 bool
 BasicPriorStatement::is_structural_innovation(const SymbolType symb_type) const
 {
-  if (symb_type == SymbolType::exogenous)
+  return symb_type == SymbolType::exogenous;
-    return true;
-  return false;
 }
 void
 BasicPriorStatement::get_base_name(const SymbolType symb_type, string& lhs_field) const
 {
-  if (symb_type == SymbolType::exogenous)
+  lhs_field = (symb_type == SymbolType::exogenous ? "structural_innovation" : "measurement_error");
-    lhs_field = "structural_innovation";
-  else
-    lhs_field = "measurement_error";
 }
 void
@@ -4685,10 +4572,7 @@ BasicOptionsStatement::is_structural_innovation(const SymbolType symb_type) cons
 void
 BasicOptionsStatement::get_base_name(const SymbolType symb_type, string& lhs_field) const
 {
-  if (symb_type == SymbolType::exogenous)
+  lhs_field = (symb_type == SymbolType::exogenous ? "structural_innovation" : "measurement_error");
-    lhs_field = "structural_innovation";
-  else
-    lhs_field = "measurement_error";
 }
 void
@@ -4896,10 +4780,7 @@ OptionsEqualStatement::writeJsonOutput(ostream& output) const
 void
 OptionsEqualStatement::get_base_name(const SymbolType symb_type, string& lhs_field) const
 {
-  if (symb_type == SymbolType::exogenous)
+  lhs_field = (symb_type == SymbolType::exogenous ? "structural_innovation" : "measurement_error");
-    lhs_field = "structural_innovation";
-  else
-    lhs_field = "measurement_error";
 }
 void
@@ -5461,8 +5342,8 @@ void
 ResidStatement::writeOutput(ostream& output, [[maybe_unused]] const string& basename,
                            [[maybe_unused]] bool minimal_workspace) const
 {
-  options_list.writeOutput(output, "options_resid_");
+  options_list.writeOutput(output);
-  output << "display_static_residuals(M_, options_, oo_, options_resid_);" << endl;
+  output << "display_static_residuals(M_, options_, oo_);" << endl;
 }
 void

--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -179,28 +179,6 @@ public:
  void writeJsonOutput(ostream& output) const override;
 };
-class RamseyConstraintsStatement : public Statement
-{
-public:
-  struct Constraint
-  {
-    int endo;
-    BinaryOpcode code;
-    expr_t expression;
-  };
-  using constraints_t = vector<Constraint>;
-private:
-  const SymbolTable& symbol_table;
-  const constraints_t constraints;
-public:
-  RamseyConstraintsStatement(const SymbolTable& symbol_table_arg, constraints_t constraints_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;
-};
 class RamseyPolicyStatement : public Statement
 {
 private:
@@ -546,7 +524,7 @@ public:
  blockName() const override
  {
    return "estimated_params";
-  };
+  }
  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;
@@ -564,7 +542,7 @@ public:
  blockName() const override
  {
    return "estimated_params_init";
-  };
+  }
  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;
@@ -579,7 +557,7 @@ public:
  blockName() const override
  {
    return "estimated_params_bounds";
-  };
+  }
  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;
@@ -1029,10 +1007,10 @@ public:
 class SetTimeStatement : public Statement
 {
 private:
-  const OptionsList options_list;
+  const string period;
 public:
-  explicit SetTimeStatement(OptionsList options_list_arg);
+  explicit SetTimeStatement(string period_arg);
  void writeOutput(ostream& output, const string& basename, bool minimal_workspace) const override;
  void writeJsonOutput(ostream& output) const override;
 };

--- a/src/Configuration.cc
+++ b/src/Configuration.cc
 /*
- * Copyright © 2010-2023 Dynare Team
+ * Copyright © 2010-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -666,15 +666,14 @@ Configuration::transformPass()
      }
 #endif
-  auto cluster_it
+  auto& cluster = cluster_name.empty() ? clusters.at(firstClusterName) : clusters.at(cluster_name);
-      = cluster_name.empty() ? clusters.find(firstClusterName) : clusters.find(cluster_name);
  double weight_denominator {0.0};
-  for (const auto& it : cluster_it->second.member_nodes)
+  for (const auto& [name, weight] : cluster.member_nodes)
-    weight_denominator += it.second;
+    weight_denominator += weight;
-  for (auto& member_node : cluster_it->second.member_nodes)
+  for (auto& [name, weight] : cluster.member_nodes)
-    member_node.second /= weight_denominator;
+    weight /= weight_denominator;
 }
 vector<filesystem::path>
@@ -701,51 +700,35 @@ Configuration::writeCluster(ostream& output) const
  if (!parallel && !parallel_test)
    return;
-  auto cluster_it
+  const auto& cluster
-      = cluster_name.empty() ? clusters.find(firstClusterName) : clusters.find(cluster_name);
+      = cluster_name.empty() ? clusters.at(firstClusterName) : clusters.at(cluster_name);
-  for (int i {1}; const auto& follower_node : follower_nodes)
+  for (int i {1}; const auto& [name, node] : follower_nodes)
    {
-      bool follower_node_in_member_nodes = false;
+      if (!cluster.member_nodes.contains(name))
-      for (const auto& itmn : cluster_it->second.member_nodes)
+        continue; // Skip nodes not in the selected cluster
-        if (follower_node.first == itmn.first)
-          follower_node_in_member_nodes = true;
-      if (!follower_node_in_member_nodes)
-        continue;
      output << "options_.parallel";
      if (i > 1)
        output << "(" << i << ")";
      i++;
-      output << " = struct('Local', ";
+      output << " = struct('Local', " << noboolalpha << (node.computerName == "localhost") << ", "
-      if (follower_node.second.computerName == "localhost")
+             << "'ComputerName', '" << node.computerName << "', "
-        output << "1, ";
+             << "'Port', '" << node.port << "', "
-      else
+             << "'CPUnbr', [" << node.minCpuNbr << ":" << node.maxCpuNbr << "], "
-        output << "0, ";
+             << "'UserName', '" << node.userName << "', "
+             << "'Password', '" << node.password << "', "
-      output << "'ComputerName', '" << follower_node.second.computerName << "', "
+             << "'RemoteDrive', '" << node.remoteDrive << "', "
-             << "'Port', '" << follower_node.second.port << "', "
+             << "'RemoteDirectory', '" << node.remoteDirectory
-             << "'CPUnbr', [" << follower_node.second.minCpuNbr << ":"
-             << follower_node.second.maxCpuNbr << "], "
-             << "'UserName', '" << follower_node.second.userName << "', "
-             << "'Password', '" << follower_node.second.password << "', "
-             << "'RemoteDrive', '" << follower_node.second.remoteDrive << "', "
-             << "'RemoteDirectory', '" << follower_node.second.remoteDirectory
             << "', "
             // The following should be switched back to “ProgramPath” once we move to Dragonfly
-             << "'DynarePath', '" << follower_node.second.programPath << "', "
+             << "'DynarePath', '" << node.programPath << "', "
-             << "'ProgramConfig', '" << follower_node.second.programConfig << "', "
+             << "'ProgramConfig', '" << node.programConfig << "', "
-             << "'MatlabOctavePath', '" << follower_node.second.matlabOctavePath << "', "
+             << "'MatlabOctavePath', '" << node.matlabOctavePath << "', "
-             << "'OperatingSystem', '" << follower_node.second.operatingSystem << "', "
+             << "'OperatingSystem', '" << node.operatingSystem << "', "
-             << "'NodeWeight', '" << cluster_it->second.member_nodes.at(follower_node.first)
+             << "'NodeWeight', '" << cluster.member_nodes.at(name) << "', "
-             << "', "
+             << "'NumberOfThreadsPerJob', " << node.numberOfThreadsPerJob << ", "
-             << "'NumberOfThreadsPerJob', " << follower_node.second.numberOfThreadsPerJob << ", ";
+             << "'SingleCompThread', '" << boolalpha << node.singleCompThread << "');" << endl;
-      if (follower_node.second.singleCompThread)
-        output << "'SingleCompThread', 'true');" << endl;
-      else
-        output << "'SingleCompThread', 'false');" << endl;
    }
  // Default values for the following two are both in DynareMain.cc and

--- a/src/Configuration.hh
+++ b/src/Configuration.hh
@@ -47,7 +47,7 @@ private:
    get_paths() const
    {
      return paths;
-    };
+    }
  private:
    map<string, vector<string>> paths;

--- a/src/DataTree.cc
+++ b/src/DataTree.cc
 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -24,6 +24,7 @@
 #include <fstream>
 #include <iostream>
 #include <iterator>
+#include <ranges>
 #include "DataTree.hh"
@@ -47,10 +48,12 @@ DataTree::initConstants()
 }
 DataTree::DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-                   ExternalFunctionsTable& external_functions_table_arg, bool is_dynamic_arg) :
+                   ExternalFunctionsTable& external_functions_table_arg,
+                   HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg) :
    symbol_table {symbol_table_arg},
    num_constants {num_constants_arg},
    external_functions_table {external_functions_table_arg},
+    heterogeneity_table {heterogeneity_table_arg},
    is_dynamic {is_dynamic_arg}
 {
  initConstants();
@@ -60,19 +63,21 @@ DataTree::DataTree(const DataTree& d) :
    symbol_table {d.symbol_table},
    num_constants {d.num_constants},
    external_functions_table {d.external_functions_table},
+    heterogeneity_table {d.heterogeneity_table},
    is_dynamic {d.is_dynamic},
    local_variables_vector {d.local_variables_vector}
 {
  // Constants must be initialized first because they are used in some Add* methods
  initConstants();
+  // See commment in DataTree::operator=() for the rationale
+  for (int symb_id : d.local_variables_vector)
+    local_variables_table[symb_id] = d.local_variables_table.at(symb_id)->clone(*this);
  for (const auto& it : d.node_list)
    it->clone(*this);
  assert(node_list.size() == d.node_list.size());
-  for (const auto& [symb_id, value] : d.local_variables_table)
-    local_variables_table[symb_id] = value->clone(*this);
 }
 DataTree&
@@ -81,6 +86,7 @@ DataTree::operator=(const DataTree& d)
  assert(&symbol_table == &d.symbol_table);
  assert(&num_constants == &d.num_constants);
  assert(&external_functions_table == &d.external_functions_table);
+  assert(&heterogeneity_table == &d.heterogeneity_table);
  assert(is_dynamic == d.is_dynamic);
  num_const_node_map.clear();
@@ -791,11 +797,17 @@ DataTree::AddSecondDerivExternalFunction(int top_level_symb_id, const vector<exp
  return p;
 }
+expr_t
+DataTree::AddSum(expr_t arg)
+{
+  return AddUnaryOp(UnaryOpcode::sum, arg);
+}
 bool
 DataTree::isSymbolUsed(int symb_id) const
 {
-  for (const auto& [symb_lag, expr] : variable_node_map)
+  if (ranges::any_of(views::keys(variable_node_map),
-    if (symb_lag.first == symb_id)
+                     [=](const auto& symb_lag) { return symb_lag.first == symb_id; }))
    return true;
  if (local_variables_table.contains(symb_id))
@@ -842,18 +854,18 @@ DataTree::addAllParamDerivId([[maybe_unused]] set<int>& deriv_id_set)
 bool
 DataTree::isUnaryOpUsed(UnaryOpcode opcode) const
 {
-  return any_of(unary_op_node_map.begin(), unary_op_node_map.end(),
+  return ranges::any_of(views::keys(unary_op_node_map),
-                [=](const auto& it) { return get<1>(it.first) == opcode; });
+                        [=](const auto& key) { return get<1>(key) == opcode; });
 }
 bool
 DataTree::isUnaryOpUsedOnType(SymbolType type, UnaryOpcode opcode) const
 {
  set<int> var;
-  for (const auto& it : unary_op_node_map)
+  for (const auto& [key, value] : unary_op_node_map)
-    if (get<1>(it.first) == opcode)
+    if (get<1>(key) == opcode)
      {
-        it.second->collectVariables(type, var);
+        value->collectVariables(type, var);
        if (!var.empty())
          return true;
      }
@@ -863,18 +875,18 @@ DataTree::isUnaryOpUsedOnType(SymbolType type, UnaryOpcode opcode) const
 bool
 DataTree::isBinaryOpUsed(BinaryOpcode opcode) const
 {
-  return any_of(binary_op_node_map.begin(), binary_op_node_map.end(),
+  return ranges::any_of(views::keys(binary_op_node_map),
-                [=](const auto& it) { return get<2>(it.first) == opcode; });
+                        [=](const auto& key) { return get<2>(key) == opcode; });
 }
 bool
 DataTree::isBinaryOpUsedOnType(SymbolType type, BinaryOpcode opcode) const
 {
  set<int> var;
-  for (const auto& it : binary_op_node_map)
+  for (const auto& [key, value] : binary_op_node_map)
-    if (get<2>(it.first) == opcode)
+    if (get<2>(key) == opcode)
      {
-        it.second->collectVariables(type, var);
+        value->collectVariables(type, var);
        if (!var.empty())
          return true;
      }
@@ -900,7 +912,7 @@ DataTree::writeCHelpersDefinition(ostream& output) const
           << "getPowerDeriv(double x, double p, int k)" << endl
           << "{" << endl
           << "  if (fabs(x) < " << power_deriv_near_zero
-           << " && p > 0 && k > p && fabs(p-nearbyint(p)) < " << power_deriv_near_zero << ')'
+           << " && p >= 0 && k > p && fabs(p-nearbyint(p)) < " << power_deriv_near_zero << ')'
           << endl
           << "    return 0.0;" << endl
           << "  else" << endl
@@ -959,7 +971,7 @@ DataTree::writeToFileIfModified(stringstream& new_contents, const filesystem::pa
 {
  ifstream old_file {filename, ios::in | ios::binary};
  if (old_file.is_open()
-      && equal(istreambuf_iterator<char> {old_file}, istreambuf_iterator<char> {},
+      && ranges::equal(istreambuf_iterator<char> {old_file}, istreambuf_iterator<char> {},
                       istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {}))
    return;
  old_file.close();
@@ -972,7 +984,7 @@ DataTree::writeToFileIfModified(stringstream& new_contents, const filesystem::pa
      cerr << "ERROR: Can't open file " << filename.string() << " for writing" << endl;
      exit(EXIT_FAILURE);
    }
-  copy(istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {},
+  ranges::copy(istreambuf_iterator<char> {new_contents}, istreambuf_iterator<char> {},
               ostreambuf_iterator<char> {new_file});
  new_file.close();
 }
--- a/src/DataTree.hh
+++ b/src/DataTree.hh
 /*
- * Copyright © 2003-2023 Dynare Team
+ * Copyright © 2003-2024 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -33,6 +33,7 @@
 #include "ExprNode.hh"
 #include "ExternalFunctionsTable.hh"
+#include "HeterogeneityTable.hh"
 #include "NumericalConstants.hh"
 #include "SubModel.hh"
 #include "SymbolTable.hh"
@@ -48,7 +49,11 @@ public:
  NumericalConstants& num_constants;
  //! A reference to the external functions table
  ExternalFunctionsTable& external_functions_table;
+  // A reference to the heterogeneity table
+  HeterogeneityTable& heterogeneity_table;
  //! Is it possible to use leads/lags on variable nodes?
+  /* NB: This data member cannot be replaced by a virtual method, because this information is needed
+     in AddVariable(), which itself can be called from the copy constructor. */
  const bool is_dynamic;
 private:
@@ -137,7 +142,8 @@ private:
 public:
  DataTree(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
-           ExternalFunctionsTable& external_functions_table_arg, bool is_static_args = false);
+           ExternalFunctionsTable& external_functions_table_arg,
+           HeterogeneityTable& heterogeneity_table_arg, bool is_dynamic_arg = false);
  virtual ~DataTree() = default;
@@ -273,6 +279,9 @@ public:
  //! Adds an external function node for the second derivative of an external function
  expr_t AddSecondDerivExternalFunction(int top_level_symb_id, const vector<expr_t>& arguments,
                                        int input_index1, int input_index2);
+  // Adds "SUM(arg)" to model tree
+  expr_t AddSum(expr_t arg);
  //! Checks if a given symbol is used somewhere in the data tree
  [[nodiscard]] bool isSymbolUsed(int symb_id) const;
  //! Checks if a given unary op is used somewhere in the data tree
@@ -345,14 +354,6 @@ public:
  //! Adds to the set all the deriv IDs corresponding to parameters
  virtual void addAllParamDerivId(set<int>& deriv_id_set);
-  //! Returns bool indicating whether DataTree represents a Dynamic Model (returns true in
-  //! DynamicModel.hh)
-  [[nodiscard]] virtual bool
-  isDynamic() const
-  {
-    return false;
-  };
  struct UnknownLocalVariableException
  {
    //! Symbol ID
@@ -360,13 +361,19 @@ public:
  };
  [[nodiscard]] expr_t
-  getLocalVariable(int symb_id) const
+  getLocalVariable(int symb_id, int lead_lag) const
  {
    auto it = local_variables_table.find(symb_id);
    if (it == local_variables_table.end())
      throw UnknownLocalVariableException {symb_id};
+    /* In the following, the case without lead/lag is optimized. It makes a difference on models
+       with many nested model-local variables, see e.g.
+       https://forum.dynare.org/t/pre-processing-takes-very-long/26865 */
+    if (lead_lag == 0)
      return it->second;
+    else
+      return it->second->decreaseLeadsLags(-lead_lag);
  }
  static void
@@ -395,7 +402,7 @@ DataTree::AddPossiblyNegativeConstant(double v)
  if (isnan(v))
    return NaN;
  if (isinf(v))
-    return (v < 0 ? MinusInfinity : Infinity);
+    return v < 0 ? MinusInfinity : Infinity;
  bool neg = false;
  if (v < 0)

--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -34,17 +34,27 @@
 void
 DynamicModel::copyHelper(const DynamicModel& m)
 {
-  auto f = [this](const ExprNode* e) { return e->clone(*this); };
+  auto f = [this](const ExprNode* e) { return e ? e->clone(*this) : nullptr; };
  for (const auto& it : m.static_only_equations)
    static_only_equations.push_back(dynamic_cast<BinaryOpNode*>(f(it)));
+  for (const auto& it : m.static_only_complementarity_conditions)
+    if (it)
+      {
+        const auto& [symb_id, lb, ub] = *it;
+        static_only_complementarity_conditions.emplace_back(in_place, symb_id, f(lb), f(ub));
+      }
+    else
+      static_only_complementarity_conditions.emplace_back(nullopt);
 }
 DynamicModel::DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
                           ExternalFunctionsTable& external_functions_table_arg,
+                           HeterogeneityTable& heterogeneity_table_arg,
                           TrendComponentModelTable& trend_component_model_table_arg,
                           VarModelTable& var_model_table_arg) :
-    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg, true},
+    ModelTree {symbol_table_arg, num_constants_arg, external_functions_table_arg,
+               heterogeneity_table_arg, true},
    trend_component_model_table {trend_component_model_table_arg},
    var_model_table {var_model_table_arg}
 {
@@ -104,6 +114,7 @@ DynamicModel::operator=(const DynamicModel& m)
  static_only_equations_lineno = m.static_only_equations_lineno;
  static_only_equations_equation_tags = m.static_only_equations_equation_tags;
+  static_only_complementarity_conditions.clear();
  deriv_id_table = m.deriv_id_table;
  inv_deriv_id_table = m.inv_deriv_id_table;
  dyn_jacobian_cols_table = m.dyn_jacobian_cols_table;
@@ -171,31 +182,29 @@ DynamicModel::writeDynamicBytecode(const string& basename) const
                                                 + temporary_terms_derivatives[1].size())};
  // Declare the (single) block
-  vector<int> exo(symbol_table.exo_nbr()), exo_det(symbol_table.exo_det_nbr());
+  auto exo = views::iota(0, symbol_table.exo_nbr());
-  iota(exo.begin(), exo.end(), 0);
+  auto exo_det = views::iota(0, symbol_table.exo_det_nbr());
-  iota(exo_det.begin(), exo_det.end(), 0);
+  auto eq_idx = views::iota(0, static_cast<int>(equations.size()));
+  auto endo_idx = views::iota(0, symbol_table.endo_nbr());
-  int jacobian_ncols_endo {static_cast<int>(count_if(
-      dyn_jacobian_cols_table.begin(), dyn_jacobian_cols_table.end(),
+  int jacobian_ncols_endo {
-      [this](const auto& v) { return getTypeByDerivID(v.first) == SymbolType::endogenous; }))};
+      static_cast<int>(ranges::count_if(dyn_jacobian_cols_table, [this](const auto& v) {
-  vector<int> eq_idx(equations.size());
+        return getTypeByDerivID(v.first) == SymbolType::endogenous;
-  iota(eq_idx.begin(), eq_idx.end(), 0);
+      }))};
-  vector<int> endo_idx(symbol_table.endo_nbr());
-  iota(endo_idx.begin(), endo_idx.end(), 0);
  code_file << Bytecode::FBEGINBLOCK {symbol_table.endo_nbr(),
                                      simulation_type,
                                      0,
                                      symbol_table.endo_nbr(),
-                                      endo_idx,
+                                      {endo_idx.begin(), endo_idx.end()},
-                                      eq_idx,
+                                      {eq_idx.begin(), eq_idx.end()},
                                      false,
                                      u_count_int,
                                      jacobian_ncols_endo,
                                      symbol_table.exo_det_nbr(),
                                      symbol_table.exo_nbr(),
-                                      exo_det,
+                                      {exo_det.begin(), exo_det.end()},
-                                      exo};
+                                      {exo.begin(), exo.end()}};
  writeBytecodeHelper<true>(code_file);
 }
@@ -651,10 +660,10 @@ DynamicModel::parseIncludeExcludeEquations(const string& inc_exc_option_value, b
 }
 vector<int>
-DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
+DynamicModel::removeEquationsHelper(
-                                    bool excluded_vars_change_type,
+    set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
-                                    vector<BinaryOpNode*>& all_equations,
+    vector<BinaryOpNode*>& all_equations, vector<optional<int>>& all_equations_lineno,
-                                    vector<optional<int>>& all_equations_lineno,
+    vector<optional<tuple<int, expr_t, expr_t>>>& all_complementarity_conditions,
    EquationTags& all_equation_tags, bool static_equations) const
 {
  if (all_equations.empty())
@@ -686,6 +695,7 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
  // remove from equations, equations_lineno, equation_tags
  vector<BinaryOpNode*> new_equations;
  vector<optional<int>> new_equations_lineno;
+  vector<optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
  map<int, int> old_eqn_num_2_new;
  vector<int> excluded_vars;
  for (size_t i = 0; i < all_equations.size(); i++)
@@ -717,11 +727,13 @@ DynamicModel::removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag,
        new_equations.emplace_back(all_equations[i]);
        old_eqn_num_2_new[i] = new_equations.size() - 1;
        new_equations_lineno.emplace_back(all_equations_lineno[i]);
+        new_complementarity_conditions.emplace_back(all_complementarity_conditions[i]);
      }
  int n_excl = all_equations.size() - new_equations.size();
-  all_equations = new_equations;
+  all_equations = move(new_equations);
-  all_equations_lineno = new_equations_lineno;
+  all_equations_lineno = move(new_equations_lineno);
+  all_complementarity_conditions = move(new_complementarity_conditions);
  all_equation_tags.erase(eqs_to_delete_by_number, old_eqn_num_2_new);
@@ -754,12 +766,13 @@ DynamicModel::removeEquations(const vector<map<string, string>>& listed_eqs_by_t
  vector<int> excluded_vars
      = removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type, equations,
-                              equations_lineno, equation_tags, false);
+                              equations_lineno, complementarity_conditions, equation_tags, false);
  // Ignore output because variables are not excluded when equations marked 'static' are excluded
  removeEquationsHelper(listed_eqs_by_tag2, exclude_eqs, excluded_vars_change_type,
                        static_only_equations, static_only_equations_lineno,
-                        static_only_equations_equation_tags, true);
+                        static_only_complementarity_conditions, static_only_equations_equation_tags,
+                        true);
  if (!listed_eqs_by_tag2.empty())
    {
@@ -1030,14 +1043,8 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
         << (static_only_equations.size() > 0) << ";" << endl;
  // Say if model contains an external function call
-  bool has_external_function = false;
+  output << "M_.has_external_function = " << boolalpha
-  for (auto equation : equations)
+         << ranges::any_of(equations, &ExprNode::containsExternalFunction) << ';' << endl;
-    if (equation->containsExternalFunction())
-      {
-        has_external_function = true;
-        break;
-      }
-  output << "M_.has_external_function = " << boolalpha << has_external_function << ';' << endl;
  // Compute list of state variables, ordered in block-order
  vector<int> state_var;
@@ -1047,8 +1054,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
      try
        {
          getDerivID(symbol_table.getID(SymbolType::endogenous, endo_idx_block2orig[endoID]), lag);
-          if (find(state_var.begin(), state_var.end(), endo_idx_block2orig[endoID])
+          if (ranges::find(state_var, endo_idx_block2orig[endoID]) == state_var.end())
-              == state_var.end())
            state_var.push_back(endo_idx_block2orig[endoID]);
        }
      catch (UnknownDerivIDException& e)
@@ -1123,7 +1129,7 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
    output << (i > computed_derivs_order ? -1 : NNZDerivatives[i]) << "; ";
  output << "];" << endl;
-  writeDriverSparseIndicesHelper<true, false>(output);
+  writeDriverSparseIndicesHelper("dynamic", output);
  // Write LHS of each equation in text form
  output << "M_.lhs = {" << endl;
@@ -1134,6 +1140,11 @@ DynamicModel::writeDriverOutput(ostream& output, bool compute_xrefs) const
      output << "'; " << endl;
    }
  output << "};" << endl;
+  output << "M_.dynamic_mcp_equations_reordering = [";
+  for (auto i : mcp_equations_reordering)
+    output << i + 1 << "; ";
+  output << "];" << endl;
 }
 void
@@ -1203,8 +1214,7 @@ DynamicModel::updateVarAndTrendModel() const
                        catch (...)
                          {
                          }
-                      if (find(trend_lhs.begin(), trend_lhs.end(), *trend_var_symb_id)
+                      if (ranges::find(trend_lhs, *trend_var_symb_id) == trend_lhs.end())
-                          == trend_lhs.end())
                        {
                          cerr << "ERROR: trend found in trend_component equation #" << eqn << " ("
                               << symbol_table.getName(*trend_var_symb_id)
@@ -1343,8 +1353,7 @@ DynamicModel::fillVarModelTableFromOrigModel() const
          // save lhs variables
          equations[eqn]->arg1->collectVARLHSVariable(lhs);
-          equations[eqn]->arg1->countDiffs() > 0 ? diff_vec.push_back(true)
+          diff_vec.push_back(equations[eqn]->arg1->countDiffs() > 0);
-                                                 : diff_vec.push_back(false);
          if (diff_vec.back())
            {
              set<pair<int, int>> diff_set;
@@ -1632,13 +1641,13 @@ DynamicModel::computeErrorComponentMatrices(const ExprNode::subst_table_t& diff_
      for (int i {0}; auto eqn : eqns)
        {
-          if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
+          if (ranges::find(nontarget_eqnums, eqn) != nontarget_eqnums.end())
            parsed_undiff_nontarget_lhs.push_back(undiff_nontarget_lhs.at(i));
          i++;
        }
      for (int i {0}; auto eqn : eqns)
-        if (find(nontarget_eqnums.begin(), nontarget_eqnums.end(), eqn) != nontarget_eqnums.end())
+        if (ranges::find(nontarget_eqnums, eqn) != nontarget_eqnums.end())
          equations[eqn]->arg2->fillErrorCorrectionRow(i++, parsed_undiff_nontarget_lhs, target_lhs,
                                                       A0, A0star);
      A0r[model_name] = A0;
@@ -1688,10 +1697,7 @@ DynamicModel::fillTrendComponentModelTableFromOrigModel() const
          // save lhs variables
          equations[eqn]->arg1->collectVARLHSVariable(lhs);
-          if (equations[eqn]->arg1->countDiffs() > 0)
+          diff_vec.push_back(equations[eqn]->arg1->countDiffs() > 0);
-            diff_vec.push_back(true);
-          else
-            diff_vec.push_back(false);
          if (diff_vec.back())
            {
              set<pair<int, int>> diff_set;
@@ -1756,7 +1762,7 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
  for (auto eqn : nontrend_eqnums)
    {
-      auto i = distance(eqnumber.begin(), find(eqnumber.begin(), eqnumber.end(), eqn));
+      auto i = ranges::distance(eqnumber.begin(), ranges::find(eqnumber, eqn));
      if (eqnumber[i] != eqn)
        {
@@ -1764,7 +1770,7 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
          exit(EXIT_FAILURE);
        }
-      if (diff.at(i) != true)
+      if (!diff.at(i))
        {
          cerr << "ERROR: the variable on the LHS of equation #" << eqn
               << " does not have the diff operator applied to it yet you are trying to undiff it."
@@ -1772,7 +1778,6 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
          exit(EXIT_FAILURE);
        }
-      bool printerr = false;
      expr_t node = nullptr;
      expr_t aux_var = lhs_expr_t.at(i);
      for (const auto& it : diff_subst_table)
@@ -1789,11 +1794,8 @@ DynamicModel::getUndiffLHSForPac(const string& aux_model_name,
        }
      node = node->undiff();
-      auto it1 = diff_subst_table.find(node);
-      if (it1 == diff_subst_table.end())
-        printerr = true;
-      if (printerr)
+      if (auto it1 = diff_subst_table.find(node); it1 == diff_subst_table.end())
        { // we have undiffed something like diff(x), hence x is not in diff_subst_table
          lhs_expr_t.at(i) = node;
          lhs.at(i) = dynamic_cast<VariableNode*>(node)->symb_id;
@@ -1814,7 +1816,7 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
  for (auto& equation : equations)
    if (equation->containsPacExpectation(name))
      {
-        if (!pac_eq_name[name].empty())
+        if (pac_eq_name.contains(name))
          {
            cerr << "It is not possible to use 'pac_expectation(" << name
                 << ")' in several equations." << endl;
@@ -1913,6 +1915,13 @@ DynamicModel::analyzePacEquationStructure(const string& name, map<string, string
                                   move(additive_vars_params_and_constants),
                                   move(optim_additive_vars_params_and_constants)};
      }
+  if (!pac_eq_name.contains(name))
+    {
+      cerr << "ERROR: the model does not contain the 'pac_expectation(" << name << ")' operator."
+           << endl;
+      exit(EXIT_FAILURE);
+    }
 }
 int
@@ -2084,6 +2093,203 @@ DynamicModel::computePacModelConsistentExpectationSubstitution(
  pac_expectation_substitution[name] = AddPlus(AddVariable(mce_z1_symb_id), growth_correction_term);
 }
+void
+DynamicModel::computePacModelConsistentExpectationSubstitutionWithComponents(
+    const string& name, int discount_symb_id, int pac_eq_max_lag,
+    ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+    vector<PacModelTable::target_component_t>& pac_target_components,
+    map<string, expr_t>& pac_expectation_substitution)
+{
+  auto create_aux_param = [&](const string& param_name) {
+    try
+      {
+        return symbol_table.addSymbol(param_name, SymbolType::parameter);
+      }
+    catch (SymbolTable::AlreadyDeclaredException)
+      {
+        cerr << "ERROR: the variable/parameter '" << param_name
+             << "' conflicts with some auxiliary parameter that will be generated for the '" << name
+             << "' PAC model. Please rename that parameter." << endl;
+        exit(EXIT_FAILURE);
+      }
+  };
+  int neqs = 0;
+  // At the end of this loop:
+  //         ₘ                          ₘ
+  // A_1 = 1+∑ αᵢ = A(1) and A_beta = 1+∑ αᵢβⁱ
+  //        ᵢ₌₁                        ᵢ₌₁
+  expr_t A_1 = One;
+  expr_t A_beta = One;
+  expr_t beta = AddVariable(discount_symb_id);
+  for (int i = 1; i <= pac_eq_max_lag + 1; i++)
+    {
+      string param_name = "mce_alpha_" + name + "_" + to_string(i);
+      try
+        {
+          int alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+          pac_aux_param_symb_ids[name].push_back(alpha_i_symb_id);
+          A_1 = AddPlus(A_1, AddVariable(alpha_i_symb_id));
+          A_beta = AddPlus(A_beta, AddTimes(AddVariable(alpha_i_symb_id),
+                                            AddPower(beta, AddPossiblyNegativeConstant(i))));
+        }
+      catch (SymbolTable::AlreadyDeclaredException& e)
+        {
+          cerr << "The variable/parameter '" << param_name
+               << "' conflicts with a parameter that will be generated for the '" << name
+               << "' PAC model. Please rename it." << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
+  auto create_target_lag = [&](int variable, int lag) {
+    if (symbol_table.isAuxiliaryVariable(variable))
+      return AddVariable(symbol_table.getOrigSymbIdForAuxVar(variable), lag);
+    else
+      return AddVariable(variable, lag);
+  };
+  expr_t substexpr = Zero;
+  for (int component_idx {1};
+       auto& [component, growth, auxname, kind, coeff, growth_neutrality_param, h_indices,
+              original_growth, growth_info] : pac_target_components)
+    {
+      // Create the auxiliary variable for this component
+      int aux_id = symbol_table.addPacExpectationAuxiliaryVar(auxname, component);
+      // Get the component variable id
+      int component_id = dynamic_cast<VariableNode*>(component)->symb_id;
+      //      ₘ
+      // fp = ∑ αᵢβⁱZₜ₊ᵢ
+      //     ᵢ₌₁
+      expr_t fp = Zero;
+      for (int i = 1; i <= pac_eq_max_lag + 1; i++)
+        {
+          int alpha_i_symb_id = -1;
+          string param_name = "mce_alpha_" + name + "_" + to_string(i);
+          try
+            {
+              alpha_i_symb_id = symbol_table.getID(param_name);
+            }
+          catch (SymbolTable::UnknownSymbolNameException& e)
+            {
+              alpha_i_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+            }
+          fp = AddPlus(fp, AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
+                                             AddPower(beta, AddPossiblyNegativeConstant(i))),
+                                    AddVariable(aux_id, i)));
+        }
+      if (kind != PacTargetKind::ll) // non-stationary component y¹ₜ
+        {
+          // Add diff nodes and eqs for the non-stationary component
+          const VariableNode* aux_target_ns_var_diff_node;
+          expr_t diff_node_to_search = AddDiff(create_target_lag(component_id, 0));
+          if (auto sit = diff_subst_table.find(diff_node_to_search); sit != diff_subst_table.end())
+            aux_target_ns_var_diff_node = sit->second;
+          else
+            {
+              int symb_id
+                  = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
+              aux_target_ns_var_diff_node = AddVariable(symb_id);
+              auto neweq = AddEqual(const_cast<VariableNode*>(aux_target_ns_var_diff_node),
+                                    AddMinus(create_target_lag(component_id, 0),
+                                             create_target_lag(component_id, -1)));
+              addEquation(neweq, nullopt);
+              addAuxEquation(neweq);
+              neqs++;
+            }
+          map<int, VariableNode*> target_ns_aux_var_to_add;
+          const VariableNode* last_aux_var = aux_target_ns_var_diff_node;
+          for (int i = 1; i <= pac_eq_max_lag; i++, neqs++)
+            {
+              expr_t this_diff_node = AddDiff(create_target_lag(component_id, i));
+              int symb_id = symbol_table.addDiffLeadAuxiliaryVar(
+                  this_diff_node->idx, this_diff_node, last_aux_var->symb_id, 1);
+              VariableNode* current_aux_var = AddVariable(symb_id);
+              auto neweq = AddEqual(current_aux_var, AddVariable(last_aux_var->symb_id, 1));
+              addEquation(neweq, nullopt);
+              addAuxEquation(neweq);
+              last_aux_var = current_aux_var;
+              target_ns_aux_var_to_add[i] = current_aux_var;
+            }
+          // At the end of this loop,
+          //      ₘ₋₁ ⎛ ₘ₋₁     ⎞
+          //  fs = ∑  ⎢  ∑  αᵢβⁱ⎥ Δy¹ₜ₊ᵢ
+          //      ₖ₌₁ ⎝ ᵢ₌ₖ     ⎠
+          expr_t fs = Zero;
+          for (int k = 1; k <= pac_eq_max_lag; k++)
+            {
+              expr_t ssum = Zero;
+              for (int j = k + 1; j <= pac_eq_max_lag + 1; j++)
+                {
+                  int alpha_j_symb_id = -1;
+                  string param_name = "mce_alpha_" + name + "_" + to_string(j);
+                  try
+                    {
+                      alpha_j_symb_id = symbol_table.getID(param_name);
+                    }
+                  catch (SymbolTable::UnknownSymbolNameException& e)
+                    {
+                      alpha_j_symb_id = symbol_table.addSymbol(param_name, SymbolType::parameter);
+                    }
+                  ssum = AddPlus(ssum, AddTimes(AddVariable(alpha_j_symb_id),
+                                                AddPower(beta, AddPossiblyNegativeConstant(j))));
+                }
+              fs = AddPlus(fs, AddTimes(ssum, target_ns_aux_var_to_add[k]));
+            }
+          // Assembling the equation Zₜ¹ = A_1 ( Δyₜ¹ - fs ) - fp_1
+          auto neweq_1 = AddEqual(
+              AddVariable(aux_id),
+              AddMinus(
+                  AddTimes(A_1,
+                           AddMinus(const_cast<VariableNode*>(aux_target_ns_var_diff_node), fs)),
+                  fp));
+          addEquation(neweq_1, nullopt);
+          neqs++;
+          /* This equation is not added to the list of auxiliary equations, because it
+             is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
+        }
+      else // Stationary component yₜ⁰
+        {
+          // Assembling the equation Zₜ⁰ = A_beta A_1 yₜ⁰ - fp
+          auto neweq_0 = AddEqual(AddVariable(aux_id),
+                                  AddMinus(AddTimes(A_beta, AddTimes(A_1, component)), fp));
+          addEquation(neweq_0, nullopt);
+          neqs++;
+          /* This equation is not added to the list of auxiliary equations, because it
+             is recursive, and this would in particular break dynamic_set_auxiliary_series.m */
+        }
+      // If needed, add the growth neutrality correction for this component
+      expr_t growth_correction_term = Zero;
+      string name_component = name + "_component" + to_string(component_idx);
+      if (growth)
+        {
+          growth_neutrality_param
+              = create_aux_param(name_component + "_pac_growth_neutrality_correction");
+          growth_correction_term = AddTimes(growth, AddVariable(growth_neutrality_param));
+        }
+      else
+        growth_neutrality_param = -1;
+      substexpr = AddPlus(substexpr,
+                          AddTimes(coeff, AddPlus(AddVariable(aux_id), growth_correction_term)));
+      component_idx++;
+    }
+  cout << "PAC Model Consistent Expectation: added " << neqs
+       << " auxiliary variables and equations for model " << name << "." << endl;
+  /* The growth correction term is not added to the definition of Z₁
+     because the latter is recursive. Rather put it at the level of the
+     substition of pac_expectation operator. */
+  pac_expectation_substitution[name] = substexpr;
+}
 void
 DynamicModel::computePacBackwardExpectationSubstitution(
    const string& name, const vector<int>& lhs, int max_lag, const string& aux_model_type,
@@ -2310,6 +2516,8 @@ DynamicModel::computingPass(int derivsOrder, int paramsDerivsOrder,
      cerr << "ERROR: Block decomposition requested but failed." << endl;
      exit(EXIT_FAILURE);
    }
+  computeMCPEquationsReordering();
 }
 void
@@ -2613,10 +2821,13 @@ DynamicModel::writeDynamicFile(const string& basename, bool use_dll, const strin
  writeSetAuxiliaryVariablesFile<true>(basename, julia);
-  // Support for model debugging
  if (!julia)
+    {
+      writeComplementarityConditionsFile<true>(basename);
+      // Support for model debugging
      writeDebugModelMFiles<true>(basename);
    }
+}
 void
 DynamicModel::clearEquations()
@@ -2624,6 +2835,7 @@ DynamicModel::clearEquations()
  equations.clear();
  equations_lineno.clear();
  equation_tags.clear();
+  complementarity_conditions.clear();
 }
 void
@@ -2631,14 +2843,26 @@ DynamicModel::replaceMyEquations(DynamicModel& dynamic_model) const
 {
  dynamic_model.clearEquations();
+  auto clone_if_not_null = [&](expr_t e) { return e ? e->clone(dynamic_model) : nullptr; };
  for (size_t i = 0; i < equations.size(); i++)
-    dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i]);
+    {
+      optional<tuple<int, expr_t, expr_t>> cc_clone;
+      if (complementarity_conditions[i])
+        {
+          auto& [symb_id, lower_bound, upper_bound] = *complementarity_conditions[i];
+          cc_clone = {symb_id, clone_if_not_null(lower_bound), clone_if_not_null(upper_bound)};
+        }
+      dynamic_model.addEquation(equations[i]->clone(dynamic_model), equations_lineno[i],
+                                move(cc_clone));
+    }
  dynamic_model.equation_tags = equation_tags;
 }
 int
-DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
+DynamicModel::computeRamseyPolicyFOCs(const StaticModel& planner_objective,
+                                      map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints)
 {
  cout << "Ramsey Problem: added " << equations.size() << " multipliers." << endl;
@@ -2652,8 +2876,8 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
    }
  // Add Planner Objective to equations so that it appears in Lagrangian
-  assert(static_model.equations.size() == 1);
+  assert(planner_objective.equations.size() == 1);
-  addEquation(static_model.equations[0]->clone(*this), nullopt);
+  addEquation(planner_objective.equations[0]->clone(*this), nullopt);
  // Get max endo lead and max endo lag
  set<pair<int, int>> dynvars;
@@ -2695,9 +2919,10 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
                      lagrangian);
      }
-  // Save line numbers and tags, see below
+  // Save line numbers, tags and complementarity conditions, see below
  auto old_equations_lineno = equations_lineno;
  auto old_equation_tags = equation_tags;
+  auto old_complementarity_conditions = complementarity_conditions;
  // Prepare derivation of the Lagrangian
  clearEquations();
@@ -2714,6 +2939,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
  vector<expr_t> neweqs;
  vector<optional<int>> neweqs_lineno;
  map<int, map<string, string>> neweqs_tags;
+  map<int, optional<tuple<int, expr_t, expr_t>>> new_complementarity_conditions;
  int orig_endo_nbr {0};
  for (auto& [symb_id_and_lag, deriv_id] : deriv_id_table)
    {
@@ -2726,12 +2952,24 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
            {
              // This is a derivative w.r.t. a Lagrange multiplier
              neweqs_lineno.push_back(old_equations_lineno[*i]);
-              neweqs_tags[neweqs.size() - 1] = old_equation_tags.getTagsByEqn(*i);
+              neweqs_tags.emplace(neweqs.size() - 1, old_equation_tags.getTagsByEqn(*i));
+              new_complementarity_conditions.emplace(neweqs.size() - 1,
+                                                     old_complementarity_conditions.at(*i));
            }
          else
            {
              orig_endo_nbr++;
              neweqs_lineno.emplace_back(nullopt);
+              if (string eqname {"Ramsey FOC w.r.t. "s + symbol_table.getName(symb_id)};
+                  !equation_tags.exists("name", eqname))
+                neweqs_tags.emplace(neweqs.size() - 1, map<string, string> {{"name", eqname}});
+              if (cloned_ramsey_constraints.contains(symb_id))
+                {
+                  auto& [lower_bound, upper_bound] = cloned_ramsey_constraints.at(symb_id);
+                  new_complementarity_conditions.emplace(neweqs.size() - 1,
+                                                         tuple {symb_id, lower_bound, upper_bound});
+                }
            }
        }
    }
@@ -2739,7 +2977,7 @@ DynamicModel::computeRamseyPolicyFOCs(const StaticModel& static_model)
  // Overwrite equations with the Lagrangian derivatives
  clearEquations();
  for (size_t i = 0; i < neweqs.size(); i++)
-    addEquation(neweqs[i], neweqs_lineno[i], neweqs_tags[i]);
+    addEquation(neweqs[i], neweqs_lineno[i], new_complementarity_conditions[i], neweqs_tags[i]);
  return orig_endo_nbr;
 }
@@ -2771,8 +3009,8 @@ DynamicModel::expandEqTags()
    if (!existing_tags.contains(eq))
      {
        if (auto lhs_expr = dynamic_cast<VariableNode*>(equations[eq]->arg1);
-            lhs_expr && !equation_tags.exists("name", symbol_table.getName(lhs_expr->symb_id)))
+            lhs_expr && !equation_tags.exists("name", lhs_expr->getName()))
-          equation_tags.add(eq, "name", symbol_table.getName(lhs_expr->symb_id));
+          equation_tags.add(eq, "name", lhs_expr->getName());
        else if (!equation_tags.exists("name", to_string(eq + 1)))
          equation_tags.add(eq, "name", to_string(eq + 1));
        else
@@ -2794,8 +3032,7 @@ DynamicModel::findUnusedEndogenous()
  for (auto& equation : static_only_equations)
    equation->collectVariables(SymbolType::endogenous, usedEndo);
  set<int> allEndo = symbol_table.getEndogenous();
-  set_difference(allEndo.begin(), allEndo.end(), usedEndo.begin(), usedEndo.end(),
+  ranges::set_difference(allEndo, usedEndo, inserter(unusedEndo, unusedEndo.begin()));
-                 inserter(unusedEndo, unusedEndo.begin()));
  return unusedEndo;
 }
@@ -2809,10 +3046,8 @@ DynamicModel::findUnusedExogenous()
    equation->collectVariables(SymbolType::exogenous, usedExo);
  set<int> observedExo = symbol_table.getObservedExogenous();
  set<int> allExo = symbol_table.getExogenous();
-  set_difference(allExo.begin(), allExo.end(), observedExo.begin(), observedExo.end(),
+  ranges::set_difference(allExo, observedExo, inserter(unobservedExo, unobservedExo.begin()));
-                 inserter(unobservedExo, unobservedExo.begin()));
+  ranges::set_difference(unobservedExo, usedExo, inserter(unusedExo, unusedExo.begin()));
-  set_difference(unobservedExo.begin(), unobservedExo.end(), usedExo.begin(), usedExo.end(),
-                 inserter(unusedExo, unusedExo.begin()));
  return unusedExo;
 }
@@ -3252,10 +3487,8 @@ pair<lag_equivalence_table_t, ExprNode::subst_table_t>
 DynamicModel::substituteUnaryOps(VarExpectationModelTable& var_expectation_model_table,
                                 PacModelTable& pac_model_table)
 {
-  vector<int> eqnumbers(equations.size());
+  auto v = views::iota(0, static_cast<int>(equations.size()));
-  iota(eqnumbers.begin(), eqnumbers.end(), 0);
+  return substituteUnaryOps({v.begin(), v.end()}, var_expectation_model_table, pac_model_table);
-  return substituteUnaryOps({eqnumbers.begin(), eqnumbers.end()}, var_expectation_model_table,
-                            pac_model_table);
 }
 pair<lag_equivalence_table_t, ExprNode::subst_table_t>
@@ -3459,6 +3692,36 @@ DynamicModel::substituteLogTransform()
    }
 }
+void
+DynamicModel::substituteAggregationOperators()
+{
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode*> neweqs;
+  for (auto& [symb_id, expr] : local_variables_table)
+    expr = expr->substituteAggregationOperators(subst_table, neweqs);
+  for (auto& equation : equations)
+    {
+      equation = dynamic_cast<BinaryOpNode*>(
+          equation->substituteAggregationOperators(subst_table, neweqs));
+      assert(equation);
+    }
+  for (auto& equation : static_only_equations)
+    {
+      equation = dynamic_cast<BinaryOpNode*>(
+          equation->substituteAggregationOperators(subst_table, neweqs));
+      assert(equation);
+    }
+  for (auto neweq : neweqs)
+    {
+      addEquation(neweq, nullopt);
+      addAuxEquation(neweq);
+    }
+}
 void
 DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
 {
@@ -3469,8 +3732,7 @@ DynamicModel::checkNoWithLogTransform(const set<int>& eqnumbers)
  const set<int>& with_log_transform = symbol_table.getVariablesWithLogTransform();
  vector<int> intersect;
-  set_intersection(endos.begin(), endos.end(), with_log_transform.begin(), with_log_transform.end(),
+  ranges::set_intersection(endos, with_log_transform, back_inserter(intersect));
-                   back_inserter(intersect));
  if (!intersect.empty())
    {
      cerr << "ERROR: the following variables are declared with var(log) and therefore cannot "
@@ -3486,18 +3748,18 @@ void
 DynamicModel::detrendEquations()
 {
  // We go backwards in the list of trend_vars, to deal correctly with I(2) processes
-  for (const auto& it : std::ranges::reverse_view(nonstationary_symbols_map))
+  for (const auto& [symb_id, deflator] : std::ranges::reverse_view(nonstationary_symbols_map))
    {
      for (auto& equation : equations)
        {
          equation = dynamic_cast<BinaryOpNode*>(
-              equation->detrend(it.first, it.second.first, it.second.second));
+              equation->detrend(symb_id, deflator.first, deflator.second));
          assert(equation);
        }
      for (auto& equation : static_only_equations)
        {
          equation = dynamic_cast<BinaryOpNode*>(
-              equation->detrend(it.first, it.second.first, it.second.second));
+              equation->detrend(symb_id, deflator.first, deflator.second));
          assert(equation);
        }
    }
@@ -3576,6 +3838,7 @@ DynamicModel::fillEvalContext(eval_context_t& eval_context) const
 void
 DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
+                                    optional<tuple<int, expr_t, expr_t>> complementarity_condition,
                                    map<string, string> eq_tags)
 {
  auto beq = dynamic_cast<BinaryOpNode*>(eq);
@@ -3584,6 +3847,7 @@ DynamicModel::addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
  static_only_equations_equation_tags.add(static_only_equations.size(), move(eq_tags));
  static_only_equations.push_back(beq);
  static_only_equations_lineno.push_back(lineno);
+  static_only_complementarity_conditions.push_back(move(complementarity_condition));
 }
 size_t
@@ -3637,7 +3901,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
    {
      auto eq_tags_dynamic = eq_tags;
      eq_tags_dynamic["dynamic"] = "";
-      addEquation(AddEqual(term, Zero), lineno, eq_tags_dynamic);
+      addEquation(AddEqual(term, Zero), lineno, nullopt, eq_tags_dynamic);
    }
  // Create or update the static equation (corresponding to the pure relax regime)
@@ -3659,7 +3923,7 @@ DynamicModel::addOccbinEquation(expr_t eq, const optional<int>& lineno, map<stri
      else
        {
          eq_tags["static"] = "";
-          addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, move(eq_tags));
+          addStaticOnlyEquation(AddEqual(basic_term, Zero), lineno, nullopt, move(eq_tags));
        }
    }
 }
@@ -3974,8 +4238,8 @@ DynamicModel::simplifyEquations()
 {
  size_t last_subst_table_size = 0;
  map<VariableNode*, NumConstNode*> subst_table;
-  // Equations with “mcp” tag are excluded, see dynare#1697
+  // Equations with a complementarity condition are excluded, see dynare#1697
-  findConstantEquationsWithoutMcpTag(subst_table);
+  findConstantEquationsWithoutComplementarityCondition(subst_table);
  while (subst_table.size() != last_subst_table_size)
    {
      last_subst_table_size = subst_table.size();
@@ -3986,7 +4250,7 @@ DynamicModel::simplifyEquations()
      for (auto& equation : static_only_equations)
        equation = dynamic_cast<BinaryOpNode*>(equation->replaceVarsInEquation(subst_table));
      subst_table.clear();
-      findConstantEquationsWithoutMcpTag(subst_table);
+      findConstantEquationsWithoutComplementarityCondition(subst_table);
    }
 }
@@ -4101,7 +4365,7 @@ DynamicModel::OccbinRegimeTracker::addRegime(const vector<string>& constraints_b
  vector<bool> new_regime_template(constraints.size(), false);
  for (const auto& c : constraints_bind)
    {
-      int i = distance(constraints.begin(), ranges::find(constraints, c));
+      int i = ranges::distance(constraints.begin(), ranges::find(constraints, c));
      new_regime_template[i] = true;
    }
  set<vector<bool>> new_regimes {new_regime_template};
@@ -4111,7 +4375,7 @@ DynamicModel::OccbinRegimeTracker::addRegime(const vector<string>& constraints_b
                         back_inserter(constraints_not_mentioned));
  for (const auto& c : constraints_not_mentioned)
    {
-      int i = distance(constraints.begin(), ranges::find(constraints, c));
+      int i = ranges::distance(constraints.begin(), ranges::find(constraints, c));
      auto new_regimes_copy = new_regimes;
      for (const auto& r : new_regimes_copy)
        {
@@ -4150,8 +4414,9 @@ DynamicModel::OccbinRegimeTracker::checkAllRegimesPresent() const
      if (it == r.end())
        break;
      *it = true;
-      if (it != r.begin())
+      /* NB: cannot use ranges::fill since vector<bool>::iterator does not satisfy
-        fill(r.begin(), prev(it), false);
+         indirectly_writable concept. Should be fixed in C++23: https://wg21.link/p2321 */
+      fill(r.begin(), it, false);
    }
  while (true);
 }

--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -73,7 +73,8 @@ public:
    void checkAllRegimesPresent() const noexcept(false);
  private:
-    pair<vector<string>, vector<string>> convertBitVectorToRegimes(const vector<bool>& r) const;
+    [[nodiscard]] pair<vector<string>, vector<string>>
+    convertBitVectorToRegimes(const vector<bool>& r) const;
  };
 private:
@@ -92,6 +93,9 @@ private:
  //! Stores the equation tags of equations declared as [static]
  EquationTags static_only_equations_equation_tags;
+  // Complementarity conditions of equations declared as [static]
+  vector<optional<tuple<int, expr_t, expr_t>>> static_only_complementarity_conditions;
  using deriv_id_table_t = map<pair<int, int>, int>;
  //! Maps a pair (symbol_id, lag) to a deriv ID
  deriv_id_table_t deriv_id_table;
@@ -279,10 +283,10 @@ private:
     Returns a list of excluded variables (empty if
     excluded_vars_change_type=false) */
-  vector<int> removeEquationsHelper(set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs,
+  vector<int> removeEquationsHelper(
-                                    bool excluded_vars_change_type,
+      set<map<string, string>>& listed_eqs_by_tag, bool exclude_eqs, bool excluded_vars_change_type,
-                                    vector<BinaryOpNode*>& all_equations,
+      vector<BinaryOpNode*>& all_equations, vector<optional<int>>& all_equations_lineno,
-                                    vector<optional<int>>& all_equations_lineno,
+      vector<optional<tuple<int, expr_t, expr_t>>>& all_complementarity_conditions,
      EquationTags& all_equation_tags, bool static_equations) const;
  //! Compute autoregressive matrices of trend component models
@@ -322,6 +326,7 @@ protected:
 public:
  DynamicModel(SymbolTable& symbol_table_arg, NumericalConstants& num_constants_arg,
               ExternalFunctionsTable& external_functions_table_arg,
+               HeterogeneityTable& heterogeneity_table_arg,
               TrendComponentModelTable& trend_component_model_table_arg,
               VarModelTable& var_model_table_arg);
@@ -455,7 +460,8 @@ public:
     Returns the number of optimality FOCs, which is by construction equal to
     the number of endogenous before adding the Lagrange multipliers
     (internally called ramsey_endo_nbr). */
-  int computeRamseyPolicyFOCs(const StaticModel& static_model);
+  int computeRamseyPolicyFOCs(const StaticModel& planner_objective,
+                              map<int, pair<expr_t, expr_t>> cloned_ramsey_constraints);
  //! Clears all equations
  void clearEquations();
@@ -464,7 +470,9 @@ public:
  void replaceMyEquations(DynamicModel& dynamic_model) const;
  //! Adds an equation marked as [static]
-  void addStaticOnlyEquation(expr_t eq, const optional<int>& lineno, map<string, string> eq_tags);
+  void addStaticOnlyEquation(expr_t eq, const optional<int>& lineno,
+                             optional<tuple<int, expr_t, expr_t>> complementarity_condition,
+                             map<string, string> eq_tags);
  //! Returns number of static only equations
  size_t staticOnlyEquationsNbr() const;
@@ -533,13 +541,6 @@ public:
  void addAllParamDerivId(set<int>& deriv_id_set) override;
-  //! Returns true indicating that this is a dynamic model
-  bool
-  isDynamic() const override
-  {
-    return true;
-  };
  //! Drive test of detrended equations
  void runTrendTest(const eval_context_t& eval_context);
@@ -567,6 +568,10 @@ public:
  // Performs the transformations associated to variables declared with “var(log)”
  void substituteLogTransform();
+  /* Performs the transformations associated to aggregation operators in heterogeneous models, such
+     as SUM(…) */
+  void substituteAggregationOperators();
  // Check that no variable was declared with “var(log)” in the given equations
  void checkNoWithLogTransform(const set<int>& eqnumbers);
@@ -638,6 +643,18 @@ public:
      map<string, int>& pac_aux_var_symb_ids, map<string, vector<int>>& pac_aux_param_symb_ids,
      map<string, expr_t>& pac_expectation_substitution);
+  /* For a PAC MCE model with an associated pac_target_info, fill pac_expectation_substitution with
+     the expression that will be substituted for the pac_expectation operator. In the process, add
+     the variables and the equations defining Z₁ and Z₀ for each component. The new auxiliary
+     parameters are added to pac_mce_alpha_symb_ids. The routine also creates the auxiliary
+     variables for the components, and adds the corresponding equations.
+  */
+  void computePacModelConsistentExpectationSubstitutionWithComponents(
+      const string& name, int discount_symb_id, int pac_eq_max_lag,
+      ExprNode::subst_table_t& diff_subst_table, map<string, vector<int>>& pac_aux_param_symb_ids,
+      vector<PacModelTable::target_component_t>& pac_target_components,
+      map<string, expr_t>& pac_expectation_substitution);
  /* For a PAC backward model, fill pac_expectation_substitution with the
     expression that will be substituted for the pac_expectation operator.
     The symbol IDs of the new parameters are also added to pac_aux_param_symb_ids.
@@ -696,8 +713,8 @@ public:
  getStaticOnlyEquationsInfo() const
  {
    return tuple {static_only_equations, static_only_equations_lineno,
-                  static_only_equations_equation_tags};
+                  static_only_complementarity_conditions, static_only_equations_equation_tags};
-  };
+  }
  //! Returns true if a parameter was used in the model block with a lead or lag
  bool ParamUsedWithLeadLag() const;

--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
 // -*- C++ -*-
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -28,20 +28,20 @@
 %define parse.error verbose
 %define parse.trace
-%code top {
-class ParsingDriver;
-}
 %code requires {
+// Only headers needed for the value and location types go here
+// Headers needed by the Bison file itself go in the unqualified %code section
 #include <string>
 #include <vector>
 #include <map>
-#include <utility>
 #include <tuple>
 #include <variant>
 #include "CommonEnums.hh"
 #include "ExprNode.hh"
+#include "Shocks.hh"
+class ParsingDriver;
 }
 %param { ParsingDriver &driver }
@@ -54,6 +54,9 @@ class ParsingDriver;
 }
 %code {
+#include <ranges>
+#include <utility>
 /* Little hack: we redefine the macro which computes the locations, because
   we need to access the location from within the parsing driver for error
   and warning messages. */
@@ -87,7 +90,7 @@ string
 str_tolower(string s)
 {
  // Converting to unsigned char is needed, see https://en.cppreference.com/w/cpp/string/byte/tolower
-  transform(s.begin(), s.end(), s.begin(), [](unsigned char c){ return std::tolower(c); });
+  ranges::transform(s, s.begin(), [](unsigned char c){ return std::tolower(c); });
  return s;
 }
 }
@@ -100,15 +103,15 @@ str_tolower(string s)
 %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_ALL_ENDOGENOUS_AND_AUXILIARY CONSIDER_ONLY_OBSERVED INITIAL_CONDITION_DECOMPOSITION
-%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 DATAFILE FILE SERIES DOUBLING DR_CYCLE_REDUCTION_TOL DR_CYCLE_REDUCTION_MAXITER 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 FIRST_SIMULATION_PERIOD LAST_SIMULATION_PERIOD LAST_OBS
 %token SET_TIME OSR_PARAMS_BOUNDS KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED
-%token <string> FALSE FLOAT_NUMBER DATES
+%token <string> FALSE FLOAT_NUMBER DATE
 %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 FILTERED_THEORETICAL_MOMENTS_GRID 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 USE_FIRST_ORDER_SOLUTION ONE_SIDED_HP_FILTER
 %token IDENTIFICATION INF_CONSTANT INITVAL INITVAL_FILE BOUNDS JSCALE INIT INFILE INVARS
 %token <string> INT_NUMBER
 %token CONDITIONAL_LIKELIHOOD
@@ -121,7 +124,7 @@ str_tolower(string s)
 %token POSTERIOR_MAX_SUBSAMPLE_DRAWS MIN MINIMAL_SOLVING_PERIODS
 %token MODE_CHECK MODE_CHECK_NEIGHBOURHOOD_SIZE MODE_CHECK_SYMMETRIC_PLOTS MODE_CHECK_NUMBER_OF_POINTS MODE_COMPUTE MODE_FILE MODEL MODEL_COMPARISON MODEL_INFO MSHOCKS ABS SIGN
 %token MODEL_DIAGNOSTICS MODIFIEDHARMONICMEAN MOMENTS_VARENDO CONTEMPORANEOUS_CORRELATION DIFFUSE_FILTER SUB_DRAWS TAPER_STEPS GEWEKE_INTERVAL RAFTERY_LEWIS_QRS RAFTERY_LEWIS_DIAGNOSTICS BROOKS_GELMAN_PLOTROWS MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
-%token NUMBER_OF_PARTICLES RESAMPLING SYSTEMATIC GENERIC RESAMPLING_THRESHOLD RESAMPLING_METHOD KITAGAWA STRATIFIED SMOOTH
+%token NUMBER_OF_PARTICLES RESAMPLING SYSTEMATIC GENERIC RESAMPLING_THRESHOLD RESAMPLING_METHOD KITAGAWA STRATIFIED SMOOTH MULTINOMIAL
 %token CPF_WEIGHTS AMISANOTRISTANI MURRAYJONESPARSLOW WRITE_EQUATION_TAGS FILTER_INITIAL_STATE
 %token NONLINEAR_FILTER_INITIALIZATION FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
 %token <string> NAME
@@ -170,7 +173,8 @@ str_tolower(string s)
 %token VLISTLOG VLISTPER SPECTRAL_DENSITY INIT2SHOCKS
 %token RESTRICTION RESTRICTION_FNAME CROSS_RESTRICTIONS NLAGS CONTEMP_REDUCED_FORM REAL_PSEUDO_FORECAST
 %token DUMMY_OBS NSTATES INDXSCALESSTATES NO_BAYESIAN_PRIOR SPECIFICATION SIMS_ZHA
-%token <string> ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX DIRICHLET DIAGONAL OPTIMAL MFS
+%token <string> ALPHA BETA ABAND NINV CMS NCMS CNUM GAMMA INV_GAMMA INV_GAMMA1 INV_GAMMA2 NORMAL UNIFORM EPS PDF FIG DR NONE PRIOR DIRICHLET MFS RESIDUAL
+%token PRIOR_VARIANCE HESSIAN IDENTITY_MATRIX
 %token GSIG2_LMDM Q_DIAG FLAT_PRIOR NCSK NSTD WEIBULL WEIBULL_PDF
 %token INDXPARR INDXOVR INDXAP APBAND INDXIMF INDXFORE FOREBAND INDXGFOREHAT INDXGIMFHAT
 %token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE UPDATED_COVARIANCE FILTER_DECOMPOSITION SMOOTHED_STATE_UNCERTAINTY SMOOTHER_REDUX
@@ -215,13 +219,16 @@ str_tolower(string s)
 %token ENDVAL_STEADY STEADY_SOLVE_ALGO STEADY_MAXIT STEADY_TOLF STEADY_TOLX STEADY_MARKOWITZ
 %token HOMOTOPY_MAX_COMPLETION_SHARE HOMOTOPY_MIN_STEP_SIZE HOMOTOPY_INITIAL_STEP_SIZE HOMOTOPY_STEP_SIZE_INCREASE_SUCCESS_COUNT
 %token HOMOTOPY_LINEARIZATION_FALLBACK HOMOTOPY_MARGINAL_LINEARIZATION_FALLBACK HOMOTOPY_EXCLUDE_VAREXO FROM_INITVAL_TO_ENDVAL
-%token STATIC_MFS RELATIVE_TO_INITVAL MATCHED_IRFS MATCHED_IRFS_WEIGHTS WEIGHTS
+%token STATIC_MFS RELATIVE_TO_INITVAL MATCHED_IRFS MATCHED_IRFS_WEIGHTS WEIGHTS PERPENDICULAR
+%token HETEROGENEITY HETEROGENEITY_DIMENSION SUM PERFECT_FORESIGHT_CONTROLLED_PATHS EXOGENIZE ENDOGENIZE
+%token PRECONDITIONER UMFITER ITERSTACK ILU ITER_TOL ITER_MAXIT GMRES_RESTART ITERSTACK_MAXLU ITERSTACK_NPERIODS ITERSTACK_NLU ITERSTACK_RELU
+CHECK_JACOBIAN_SINGULARITY
 %token <vector<string>> SYMBOL_VEC
 %type <expr_t> expression expression_or_empty
 %type <expr_t> equation hand_side
-%type <string> non_negative_number signed_number signed_integer date_str
+%type <string> non_negative_number signed_number signed_integer
 %type <string> filename symbol namespace_qualified_filename namespace_qualified_symbol
 %type <string> date_expr signed_inf signed_number_w_inf range
 %type <string> integer_range signed_integer_range boolean
@@ -241,22 +248,26 @@ str_tolower(string s)
 %type <vector<map<string, string>>> tag_pair_list_for_selection
 %type <map<string, string>> tag_pair_list
 %type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
-%type <vector<pair<int, int>>> period_list
+%type <AbstractShocksStatement::period_range_t> period_range
-%type <vector<expr_t>> matched_moments_list value_list
+%type <vector<AbstractShocksStatement::period_range_t>> period_list
+%type <vector<expr_t>> matched_moments_list value_list ramsey_constraints_list
 %type <tuple<string, BinaryOpNode*, BinaryOpNode*, expr_t, expr_t>> occbin_constraints_regime
 %type <vector<tuple<string, BinaryOpNode*, BinaryOpNode*, expr_t, expr_t>>> occbin_constraints_regimes_list
 %type <map<string, expr_t>> occbin_constraints_regime_options_list
 %type <pair<string, expr_t>> occbin_constraints_regime_option
 %type <PacTargetKind> pac_target_kind
 %type <vector<tuple<string, string, vector<pair<string, string>>>>> symbol_list_with_tex_and_partition
-%type <map<string, variant<bool, string>>> mshocks_options_list
+%type <variant<int, string>> integer_or_date
-%type <pair<string, variant<bool, string>>> mshocks_option
+%type <map<string, variant<bool, variant<int, string>>>> mshocks_options_list
+%type <pair<string, variant<bool, variant<int, string>>>> mshocks_option
 %type <pair<vector<expr_t>, vector<expr_t>>> matched_irfs_elem_values_weights
 %type <pair<pair<string, string>, vector<tuple<int, int, expr_t, expr_t>>>> matched_irfs_elem
 %type <map<pair<string, string>, vector<tuple<int, int, expr_t, expr_t>>>> matched_irfs_list
 %type <tuple<string, string, string>> matched_irfs_weights_elem_var_varexo
 %type <pair<tuple<string, string, string, string, string, string>, expr_t>> matched_irfs_weights_elem
 %type <map<tuple<string, string, string, string, string, string>, expr_t>> matched_irfs_weights_list
+%type <tuple<string, vector<AbstractShocksStatement::period_range_t>, vector<expr_t>, string>> perfect_foresight_controlled_paths_elem
+%type <vector<tuple<string, vector<AbstractShocksStatement::period_range_t>, vector<expr_t>, string>>> perfect_foresight_controlled_paths_list
 %%
 %start statement_list;
@@ -272,6 +283,7 @@ statement : parameters
          | predetermined_variables
          | model_local_variable
          | change_type
+          | heterogeneity_dimension
          | model
          | initval
          | initval_file
@@ -377,6 +389,7 @@ statement : parameters
          | perfect_foresight_solver
          | perfect_foresight_with_expectation_errors_setup
          | perfect_foresight_with_expectation_errors_solver
+          | perfect_foresight_controlled_paths
          | prior_function
          | posterior_function
          | method_of_moments
@@ -531,9 +544,9 @@ log_trend_var : LOG_TREND_VAR '(' LOG_GROWTH_FACTOR EQUAL { driver.begin_model()
              ;
 var : VAR symbol_list_with_tex_and_partition ';'
-      { driver.var($2, false); }
+      { driver.var($2, {}, false); }
    | VAR '(' LOG ')' symbol_list_with_tex_and_partition ';'
-      { driver.var($5, true); }
+      { driver.var($5, {}, true); }
    | VAR '(' DEFLATOR EQUAL { driver.begin_model(); } hand_side ')' symbol_list_with_tex_and_partition ';'
      { driver.end_nonstationary_var(false, $6, $8, false); }
    | VAR '(' LOG COMMA DEFLATOR EQUAL { driver.begin_model(); } hand_side ')' symbol_list_with_tex_and_partition ';'
@@ -542,6 +555,8 @@ var : VAR symbol_list_with_tex_and_partition ';'
      { driver.end_nonstationary_var(true, $6, $8, false); }
    /* The case LOG + LOG_DEFLATOR is omitted, because it does not make much sense
       from an economic point of view (amounts to taking the log two times) */
+    | VAR '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+      { driver.var($7, $5, false); }
    ;
 var_remove : VAR_REMOVE symbol_list ';' { driver.var_remove($2); };
@@ -615,7 +630,9 @@ var_expectation_model_option : VARIABLE EQUAL symbol
                             ;
 varexo : VAREXO symbol_list_with_tex_and_partition ';'
-         { driver.varexo($2); }
+         { driver.varexo($2, {}); }
+       | VAREXO '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+         { driver.varexo($7, $5); }
       ;
 varexo_det : VAREXO_DET symbol_list_with_tex_and_partition ';'
@@ -627,7 +644,9 @@ predetermined_variables : PREDETERMINED_VARIABLES symbol_list ';'
                        ;
 parameters : PARAMETERS symbol_list_with_tex_and_partition ';'
-             { driver.parameters($2); }
+             { driver.parameters($2, {}); }
+           | PARAMETERS '(' HETEROGENEITY EQUAL symbol ')' symbol_list_with_tex_and_partition ';'
+             { driver.parameters($7, $5); }
           ;
 model_local_variable : MODEL_LOCAL_VARIABLE symbol_list_with_tex ';'
@@ -648,6 +667,10 @@ change_type_arg : PARAMETERS
                  { $$ = SymbolType::exogenousDet; }
                ;
+heterogeneity_dimension : HETEROGENEITY_DIMENSION symbol_list ';'
+                          { driver.heterogeneity_dimension($2); }
+                        ;
 init_param : symbol EQUAL expression ';' { driver.init_param($1, $3); };
 expression : '(' expression ')'
@@ -784,20 +807,24 @@ h_options: o_filename
          | o_first_obs
          | o_data_first_obs
          | o_first_simulation_period
-          | o_date_first_simulation_period
          | o_last_simulation_period
-          | o_date_last_simulation_period
          | o_last_obs
          | o_data_last_obs
          | o_nobs
          | o_series
          ;
+integer_or_date : INT_NUMBER
+                  { $$.emplace<int>(stoi($1)); }
+                | date_expr
+                  { $$.emplace<string>($1); }
+                ;
 endval : ENDVAL ';' endval_list END ';'
         { driver.end_endval(false); }
       | ENDVAL '(' ALL_VALUES_REQUIRED ')' ';' endval_list END ';'
         { driver.end_endval(true); }
-       | ENDVAL '(' LEARNT_IN EQUAL INT_NUMBER ')' ';' endval_list END ';'
+       | ENDVAL '(' LEARNT_IN EQUAL integer_or_date ')' ';' endval_list END ';'
         { driver.end_endval_learnt_in($5); }
       ;
@@ -995,6 +1022,8 @@ model : MODEL ';' { driver.begin_model(); }
        equation_list END ';' { driver.end_model(); }
      | MODEL '(' model_options_list ')' ';' { driver.begin_model(); }
        equation_list END ';' { driver.end_model(); }
+      | MODEL '(' HETEROGENEITY EQUAL symbol ')' { driver.begin_heterogeneous_model($5); }';'
+        equation_list END ';' { driver.end_model(); }
      ;
 equation_list : equation_list equation
@@ -1011,6 +1040,14 @@ equation : hand_side EQUAL hand_side ';'
           { $$ = driver.add_model_equal($4, $6, $2); }
         | '[' tag_pair_list ']' hand_side ';'
           { $$ = driver.add_model_equal_with_zero_rhs($4, $2); }
+         | hand_side EQUAL hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal($1, $3, {}, $5); }
+         | hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal_with_zero_rhs($1, {}, $3); }
+         | '[' tag_pair_list ']' hand_side EQUAL hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal($4, $6, $2, $8); }
+         | '[' tag_pair_list ']' hand_side PERPENDICULAR hand_side ';'
+           { $$ = driver.add_model_equal_with_zero_rhs($4, $2, $6); }
         ;
 tag_pair_list : tag_pair_list COMMA tag_pair
@@ -1144,6 +1181,8 @@ hand_side : '(' hand_side ')'
            { $$ = driver.add_erfc($3); }
          | STEADY_STATE '(' hand_side ')'
            { $$ = driver.add_steady_state($3); }
+          | SUM '(' hand_side ')'
+            { $$ = driver.add_sum($3); }
          ;
 comma_hand_side : hand_side
@@ -1194,9 +1233,15 @@ shocks : SHOCKS ';' shock_list END ';' { driver.end_shocks(false); }
       | SHOCKS '(' SURPRISE ')' ';' det_shock_list END ';' { driver.end_shocks_surprise(false); }
       | SHOCKS '(' SURPRISE COMMA OVERWRITE ')' ';' det_shock_list END ';' { driver.end_shocks_surprise(true); }
       | SHOCKS '(' OVERWRITE COMMA SURPRISE ')' ';' det_shock_list END ';' { driver.end_shocks_surprise(true); }
-       | SHOCKS '(' LEARNT_IN EQUAL INT_NUMBER ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, false); }
+       | SHOCKS '(' LEARNT_IN EQUAL integer_or_date ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, false); }
-       | SHOCKS '(' LEARNT_IN EQUAL INT_NUMBER COMMA OVERWRITE ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, true); }
+       | SHOCKS '(' LEARNT_IN EQUAL integer_or_date COMMA OVERWRITE ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($5, true); }
-       | SHOCKS '(' OVERWRITE COMMA LEARNT_IN EQUAL INT_NUMBER ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($7, true); }
+       | SHOCKS '(' OVERWRITE COMMA LEARNT_IN EQUAL integer_or_date ')' ';' det_shock_list END ';' { driver.end_shocks_learnt_in($7, true); }
+       | SHOCKS '(' HETEROGENEITY EQUAL symbol ')' ';' stoch_shock_list END ';'
+         { driver.end_heterogeneous_shocks($5, false); }
+       | SHOCKS '(' HETEROGENEITY EQUAL symbol COMMA OVERWRITE ')' ';' stoch_shock_list END ';'
+         { driver.end_heterogeneous_shocks($5, true); }
+       | SHOCKS '(' OVERWRITE COMMA HETEROGENEITY EQUAL symbol ')' ';' stoch_shock_list END ';'
+         { driver.end_heterogeneous_shocks($7, true); }
       ;
 shock_list : shock_list shock_elem
@@ -1204,7 +1249,10 @@ shock_list : shock_list shock_elem
           ;
 shock_elem : det_shock_elem
-           | VAR symbol ';' STDERR expression ';'
+           | stoch_shock_elem
+           ;
+stoch_shock_elem : VAR symbol ';' STDERR expression ';'
                   { driver.add_stderr_shock($2, $5); }
                 | VAR symbol EQUAL expression ';'
                   { driver.add_var_shock($2, $4); }
@@ -1214,6 +1262,10 @@ shock_elem : det_shock_elem
                   { driver.add_correl_shock($2, $4, $6); }
                 ;
+stoch_shock_list : stoch_shock_list stoch_shock_elem
+                 | stoch_shock_elem
+                 ;
 det_shock_elem : VAR symbol ';' PERIODS period_list ';' VALUES value_list ';'
                 { driver.add_det_shock($2, $5, $8, ParsingDriver::DetShockType::standard); }
               | VAR symbol ';' PERIODS period_list ';' ADD value_list ';'
@@ -1329,11 +1381,11 @@ mshocks : MSHOCKS ';' mshock_list END ';'
          { driver.end_mshocks(false, false); }
        | MSHOCKS '(' mshocks_options_list ')' ';' mshock_list END ';'
          {
-            /* NB: the following relies of the fact that bool is the first
+            /* NB: the following relies on the fact that bool is the first
               alternative in the variant, so that default initialization of the
               variant by the [] operator will give false */
            if ($3.contains("learnt_in"))
-              driver.end_mshocks_learnt_in(get<string>($3.at("learnt_in")),
+              driver.end_mshocks_learnt_in(get<variant<int, string>>($3.at("learnt_in")),
                                           get<bool>($3["overwrite"]),
                                           get<bool>($3["relative_to_initval"]));
            else
@@ -1354,7 +1406,7 @@ mshocks_options_list : mshocks_option
 mshocks_option : OVERWRITE
                 { $$ = {"overwrite", true}; }
-               | LEARNT_IN EQUAL INT_NUMBER
+               | LEARNT_IN EQUAL integer_or_date
                 { $$ = {"learnt_in", $3}; }
               | RELATIVE_TO_INITVAL
                 { $$ = {"relative_to_initval", true}; }
@@ -1364,46 +1416,36 @@ mshock_list : mshock_list det_shock_elem
            | det_shock_elem
            ;
-period_list : period_list COMMA INT_NUMBER
+period_list : period_range
-              {
+              { $$ = { $1 }; }
-                $$ = $1;
+            | period_list period_range
-                int p = stoi($3);
-                $$.emplace_back(p, p);
-              }
-            | period_list INT_NUMBER
              {
                $$ = $1;
-                int p = stoi($2);
+                $$.emplace_back($2);
-                $$.emplace_back(p, p);
              }
-            | period_list COMMA INT_NUMBER ':' INT_NUMBER
+            | period_list COMMA period_range
              {
                $$ = $1;
-                int p1 = stoi($3), p2 = stoi($5);
+                $$.emplace_back($3);
-                if (p1 > p2)
-                  driver.error("Can't have first period index greater than second index in range specification");
-                $$.emplace_back(p1, p2);
              }
-            | period_list INT_NUMBER ':' INT_NUMBER
+            ;
+period_range : INT_NUMBER
               {
-                $$ = $1;
+                 int p = stoi($1);
-                int p1 = stoi($2), p2 = stoi($4);
+                 $$.emplace<pair<int, int>>(p, p);
-                if (p1 > p2)
-                  driver.error("Can't have first period index greater than second index in range specification");
-                $$.emplace_back(p1, p2);
               }
             | INT_NUMBER ':' INT_NUMBER
               {
                 int p1 = stoi($1), p2 = stoi($3);
                 if (p1 > p2)
                   driver.error("Can't have first period index greater than second index in range specification");
-                $$ = {{p1, p2}};
+                 $$.emplace<pair<int, int>>(p1, p2);
-              }
-            | INT_NUMBER
-              {
-                int p = stoi($1);
-                $$ = {{p, p}};
               }
+             | date_expr
+               { $$.emplace<pair<string, string>>($1, $1); }
+             | date_expr ':' date_expr
+               { $$.emplace<pair<string, string>>($1, $3); }
             ;
 value_list : value_list COMMA '(' expression ')'
@@ -1457,9 +1499,11 @@ steady_options : o_solve_algo
               | o_markowitz
               | o_steady_maxit
               | o_nocheck
+               | o_noprint
               | o_steady_tolf
               | o_steady_tolx
               | o_fsolve_options
+               | o_non_zero
               ;
 check : CHECK ';'
@@ -1504,6 +1548,8 @@ perfect_foresight_setup_options_list : perfect_foresight_setup_options_list COMM
 perfect_foresight_setup_options : o_periods
                                | o_datafile
                                | o_endval_steady
+                                | o_pf_first_simulation_period
+                                | o_pf_last_simulation_period
                                ;
 perfect_foresight_solver : PERFECT_FORESIGHT_SOLVER ';'
@@ -1542,6 +1588,15 @@ perfect_foresight_solver_options : o_stack_solve_algo
                                 | o_homotopy_linearization_fallback
                                 | o_homotopy_marginal_linearization_fallback
                                 | o_homotopy_exclude_varexo
+                                 | o_preconditioner
+                                 | o_iter_tol
+                                 | o_iter_maxit
+                                 | o_gmres_restart
+                                 | o_iterstack_maxlu
+                                 | o_iterstack_nperiods
+                                 | o_iterstack_nlu
+                                 | o_iterstack_relu
+                                 | o_check_jacobian_singularity
                                 ;
 perfect_foresight_with_expectation_errors_setup : PERFECT_FORESIGHT_WITH_EXPECTATION_ERRORS_SETUP ';'
@@ -1556,6 +1611,8 @@ perfect_foresight_with_expectation_errors_setup_options_list : perfect_foresight
 perfect_foresight_with_expectation_errors_setup_options : o_periods
                                                        | o_datafile
+                                                        | o_pf_first_simulation_period
+                                                        | o_pf_last_simulation_period
                                                        ;
 perfect_foresight_with_expectation_errors_solver : PERFECT_FORESIGHT_WITH_EXPECTATION_ERRORS_SOLVER ';'
@@ -1572,6 +1629,31 @@ perfect_foresight_with_expectation_errors_solver_options : o_pfwee_constant_simu
                                                         | perfect_foresight_solver_options
                                                         ;
+perfect_foresight_controlled_paths : PERFECT_FORESIGHT_CONTROLLED_PATHS ';' perfect_foresight_controlled_paths_list END ';'
+                                     { driver.perfect_foresight_controlled_paths($3, 1); }
+| PERFECT_FORESIGHT_CONTROLLED_PATHS '(' LEARNT_IN EQUAL integer_or_date ')' ';' perfect_foresight_controlled_paths_list END ';'
+                                     { driver.perfect_foresight_controlled_paths($8, $5); }
+                                   ;
+perfect_foresight_controlled_paths_list : perfect_foresight_controlled_paths_list perfect_foresight_controlled_paths_elem
+                                          {
+                                            $$ = $1;
+                                            $$.push_back($2);
+                                          }
+                                        | perfect_foresight_controlled_paths_elem
+                                          { $$ = { $1 }; }
+                                        ;
+perfect_foresight_controlled_paths_elem : EXOGENIZE symbol ';' PERIODS period_list ';' VALUES value_list ';' ENDOGENIZE symbol ';'
+                                          {
+                                            driver.check_symbol_is_endogenous($2);
+                                            driver.check_symbol_is_exogenous($11, false);
+                                            if ($5.size() != $8.size())
+                                              driver.error("The number of periods is different from the number of values");
+                                            $$ = { $2, $5, $8, $11};
+                                          }
+                                        ;
 method_of_moments : METHOD_OF_MOMENTS ';'
                    { driver.method_of_moments(); }
                  | METHOD_OF_MOMENTS '(' method_of_moments_options_list ')' ';'
@@ -1594,6 +1676,7 @@ method_of_moments_option : o_add_tiny_number_to_cholesky
                         | o_datafile
                         | o_dirname
                         | o_dr
+                         | o_dr_cycle_reduction_maxiter
                         | o_dr_cycle_reduction_tol
                         | o_dr_logarithmic_reduction_maxiter
                         | o_dr_logarithmic_reduction_tol
@@ -1771,6 +1854,7 @@ stoch_simul_primary_options : o_solve_algo
                            | o_pruning
                            | o_dr
                            | o_dr_cycle_reduction_tol
+                            | o_dr_cycle_reduction_maxiter
                            | o_dr_logarithmic_reduction_tol
                            | o_dr_logarithmic_reduction_maxiter
                            | o_irf_plot_threshold
@@ -2053,9 +2137,8 @@ prior_pdf : BETA_PDF
            { $$ = PriorDistributions::weibull; }
          ;
-date_str : DATES
+date_expr : DATE
+            { $$ = "dates('" + $1 + "')"; }
-date_expr : date_str
          | date_expr PLUS INT_NUMBER
            { $$ = $1 + '+' + $3; }
          ;
@@ -2307,6 +2390,7 @@ estimation_options : o_datafile
                   | o_lyapunov_doubling_tol
                   | o_dr
                   | o_dr_cycle_reduction_tol
+                   | o_dr_cycle_reduction_maxiter
                   | o_dr_logarithmic_reduction_tol
                   | o_dr_logarithmic_reduction_maxiter
                   | o_analytic_derivation
@@ -2608,22 +2692,19 @@ ramsey_policy : RAMSEY_POLICY ';'
                { driver.ramsey_policy($5); }
              ;
-ramsey_constraints : RAMSEY_CONSTRAINTS ';' ramsey_constraints_list END ';'
+ramsey_constraints : RAMSEY_CONSTRAINTS ';'
-                     { driver.add_ramsey_constraints_statement(); }
+                     { driver.begin_ramsey_constraints(); }
+                     ramsey_constraints_list END ';'
+                     { driver.end_ramsey_constraints($4); }
 		   ;
-ramsey_constraints_list : ramsey_constraints_list ramsey_constraint
+ramsey_constraints_list : ramsey_constraints_list hand_side ';'
-                 | ramsey_constraint
+                          {
-		 ;
+                            $$ = $1;
+                            $$.push_back($2);
-ramsey_constraint : NAME  LESS expression ';'
+                          }
-                    { driver.ramsey_constraint_add_less($1, $3); }
+                        | hand_side ';'
-		  | NAME  GREATER  expression ';'
+                          { $$ = { $1 }; }
-                    { driver.ramsey_constraint_add_greater($1, $3); }
-		  | NAME  LESS_EQUAL expression ';'
-                    { driver.ramsey_constraint_add_less_equal($1, $3); }
-		  | NAME  GREATER_EQUAL  expression ';'
-                    { driver.ramsey_constraint_add_greater_equal($1, $3); }
                        ;
 evaluate_planner_objective : EVALUATE_PLANNER_OBJECTIVE ';'
@@ -3457,6 +3538,7 @@ extended_path_option : o_periods
                     | o_solver_periods
                     | o_extended_path_order
                     | o_hybrid
+                     | o_use_first_order_solution
                     | o_lmmcp
                     ;
@@ -3585,7 +3667,15 @@ matched_irfs_elem : matched_irfs_elem_var_varexo
                      vector<tuple<int, int, expr_t, expr_t>> v;
                      v.reserve($3.size());
                      for (size_t i {0}; i < $3.size(); i++)
-                        v.emplace_back($3[i].first, $3[i].second, $5.first[i], $5.second[i]);
+                        try
+                          {
+                            auto [p1, p2] = get<pair<int, int>>($3[i]);
+                            v.emplace_back(p1, p2, $5.first[i], $5.second[i]);
+                          }
+                        catch (bad_variant_access&)
+                          {
+                            driver.error("matched_irfs: dates are not allowed in the 'periods' keyword");
+                          }
                      $$ = {$1, v};
                    }
                  ;
@@ -3684,6 +3774,7 @@ 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); };
 o_extended_path_order : ORDER EQUAL INT_NUMBER { driver.option_num("ep.stochastic.order", $3); };
 o_hybrid : HYBRID { driver.option_num("ep.stochastic.hybrid_order", "2"); };
+o_use_first_order_solution : USE_FIRST_ORDER_SOLUTION { driver.option_num("ep.use_first_order_solution_as_initial_guess", "true"); };
 o_steady_maxit : MAXIT EQUAL INT_NUMBER { driver.option_num("steady.maxit", $3); };
 o_simul_maxit : MAXIT EQUAL INT_NUMBER { driver.option_num("simul.maxit", $3); };
 o_bandpass_filter : BANDPASS_FILTER { driver.option_num("bandpass.indicator", "true"); }
@@ -3752,13 +3843,17 @@ o_est_first_obs : FIRST_OBS EQUAL vec_int
 o_posterior_sampling_method : POSTERIOR_SAMPLING_METHOD EQUAL QUOTED_STRING
                              { driver.option_str("posterior_sampler_options.posterior_sampling_method", $3); } ;
 o_first_obs : FIRST_OBS EQUAL INT_NUMBER { driver.option_num("first_obs", $3); };
-o_data_first_obs : FIRST_OBS EQUAL date_expr { driver.option_date("firstobs", $3); } ;
+o_data_first_obs : FIRST_OBS EQUAL date_expr { driver.option_date("first_obs", $3); } ;
-o_first_simulation_period : FIRST_SIMULATION_PERIOD EQUAL INT_NUMBER { driver.option_num("first_simulation_period", $3); };
+o_first_simulation_period : FIRST_SIMULATION_PERIOD EQUAL INT_NUMBER { driver.option_num("first_simulation_period", $3); }
-o_date_first_simulation_period : FIRST_SIMULATION_PERIOD EQUAL date_expr { driver.option_date("firstsimulationperiod", $3); } ;
+                          | FIRST_SIMULATION_PERIOD EQUAL date_expr { driver.option_date("first_simulation_period", $3); }
-o_last_simulation_period : LAST_SIMULATION_PERIOD EQUAL INT_NUMBER { driver.option_num("last_simulation_period", $3); };
+                          ;
-o_date_last_simulation_period : LAST_SIMULATION_PERIOD EQUAL date_expr { driver.option_date("lastsimulationperiod", $3); } ;
+o_last_simulation_period : LAST_SIMULATION_PERIOD EQUAL INT_NUMBER { driver.option_num("last_simulation_period", $3); }
+                         | LAST_SIMULATION_PERIOD EQUAL date_expr { driver.option_date("last_simulation_period", $3); }
+                         ;
+o_pf_first_simulation_period : FIRST_SIMULATION_PERIOD EQUAL date_expr { driver.option_date("simul.first_simulation_period", $3); };
+o_pf_last_simulation_period : LAST_SIMULATION_PERIOD EQUAL date_expr { driver.option_date("simul.last_simulation_period", $3); };
 o_last_obs : LAST_OBS EQUAL INT_NUMBER { driver.option_num("last_obs", $3); };
-o_data_last_obs : LAST_OBS EQUAL date_expr { driver.option_date("lastobs", $3); } ;
+o_data_last_obs : LAST_OBS EQUAL date_expr { driver.option_date("last_obs", $3); } ;
 o_keep_kalman_algo_if_singularity_is_detected : KEEP_KALMAN_ALGO_IF_SINGULARITY_IS_DETECTED { driver.option_num("kalman.keep_kalman_algo_if_singularity_is_detected", "true"); } ;
 o_data_nobs : NOBS EQUAL INT_NUMBER { driver.option_num("nobs", $3); };
 o_shift : SHIFT EQUAL signed_number { driver.option_num("shift", $3); };
@@ -3944,6 +4039,7 @@ o_dr : DR EQUAL CYCLE_REDUCTION { driver.option_num("dr_cycle_reduction", "true"
       | DR EQUAL LOGARITHMIC_REDUCTION { driver.option_num("dr_logarithmic_reduction", "true"); }
       | DR EQUAL DEFAULT { driver.option_num("dr_cycle_reduction", "false"); driver.option_num("dr_logarithmic_reduction", "false"); };
 o_dr_cycle_reduction_tol : DR_CYCLE_REDUCTION_TOL EQUAL non_negative_number { driver.option_num("dr_cycle_reduction_tol", $3); };
+o_dr_cycle_reduction_maxiter : DR_CYCLE_REDUCTION_MAXITER EQUAL INT_NUMBER { driver.option_num("dr_cycle_reduction_maxiter", $3); };
 o_dr_logarithmic_reduction_tol : DR_LOGARITHMIC_REDUCTION_TOL EQUAL non_negative_number { driver.option_num("dr_logarithmic_reduction_tol", $3); };
 o_dr_logarithmic_reduction_maxiter : DR_LOGARITHMIC_REDUCTION_MAXITER EQUAL INT_NUMBER { driver.option_num("dr_logarithmic_reduction_maxiter", $3); };
 o_psd_detail_plot : DETAIL_PLOT { driver.option_num("plot_shock_decomp.detail_plot", "true"); };
@@ -3996,9 +4092,11 @@ o_resampling : RESAMPLING EQUAL SYSTEMATIC
              | RESAMPLING EQUAL NONE { driver.option_num("particle.resampling.status.systematic", "false"); driver.option_num("particle.resampling.status.none", "true"); }
              | RESAMPLING EQUAL GENERIC { driver.option_num("particle.resampling.status.systematic", "false"); driver.option_num("particle.resampling.status.generic", "true"); };
 o_resampling_threshold : RESAMPLING_THRESHOLD EQUAL non_negative_number { driver.option_num("particle.resampling.threshold", $3); };
-o_resampling_method : RESAMPLING_METHOD EQUAL KITAGAWA { driver.option_num("particle.resampling.method.kitagawa", "true"); driver.option_num("particle.resampling.method.smooth", "false"); driver.option_num("particle.resampling.smethod.stratified", "false"); }
+o_resampling_method : RESAMPLING_METHOD EQUAL KITAGAWA { driver.option_num("particle.resampling.method.kitagawa", "true"); driver.option_num("particle.resampling.method.smooth", "false"); driver.option_num("particle.resampling.method.stratified", "false"); driver.option_num("particle.resampling.method.residual", "false"); driver.option_num("particle.resampling.method.multinomial", "false") ;}
-              | RESAMPLING_METHOD EQUAL SMOOTH { driver.option_num("particle.resampling.method.kitagawa", "false"); driver.option_num("particle.resampling.method.smooth", "true"); driver.option_num("particle.resampling.smethod.stratified", "false"); }
+| RESAMPLING_METHOD EQUAL SMOOTH { driver.option_num("particle.resampling.method.kitagawa", "false"); driver.option_num("particle.resampling.method.smooth", "true"); driver.option_num("particle.resampling.method.stratified", "false"); driver.option_num("particle.resampling.method.residual", "false");  driver.option_num("particle.resampling.method.multinomial", "false"); };
-              | RESAMPLING_METHOD EQUAL STRATIFIED { driver.option_num("particle.resampling.method.kitagawa", "false"); driver.option_num("particle.resampling.method.smooth", "false"); driver.option_num("particle.resampling.method.stratified", "true"); };
+| RESAMPLING_METHOD EQUAL STRATIFIED { driver.option_num("particle.resampling.method.kitagawa", "false"); driver.option_num("particle.resampling.method.smooth", "false"); driver.option_num("particle.resampling.method.stratified", "true"); driver.option_num("particle.resampling.method.residual", "false");  driver.option_num("particle.resampling.method.multinomial", "false"); };
+| RESAMPLING_METHOD EQUAL RESIDUAL { driver.option_num("particle.resampling.method.kitagawa", "false"); driver.option_num("particle.resampling.method.smooth", "false"); driver.option_num("particle.resampling.method.stratified", "false"); driver.option_num("particle.resampling.method.residual", "true"); driver.option_num("particle.resampling.method.multinomial", "false"); };
+| RESAMPLING_METHOD EQUAL MULTINOMIAL { driver.option_num("particle.resampling.method.kitagawa", "false"); driver.option_num("particle.resampling.method.smooth", "false"); driver.option_num("particle.resampling.method.stratified", "false"); driver.option_num("particle.resampling.method.residual", "false"); driver.option_num("particle.resampling.method.multinomial", "true"); };
 o_cpf_weights : CPF_WEIGHTS EQUAL AMISANOTRISTANI { driver.option_num("particle.cpf_weights_method.amisanotristani", "true"); driver.option_num("particle.cpf_weights_method.murrayjonesparslow", "false"); }
 | CPF_WEIGHTS EQUAL MURRAYJONESPARSLOW { driver.option_num("particle.cpf_weights_method.amisanotristani", "false"); driver.option_num("particle.cpf_weights_method.murrayjonesparslow", "true"); };
 o_filter_algorithm : FILTER_ALGORITHM EQUAL symbol { driver.option_str("particle.filter_algorithm", $3); };
@@ -4055,6 +4153,21 @@ o_max_dim_cova_group : MAX_DIM_COVA_GROUP EQUAL INT_NUMBER { driver.option_num("
 o_homotopy_mode : HOMOTOPY_MODE EQUAL INT_NUMBER { driver.option_num("homotopy_mode", $3); };
 o_homotopy_steps : HOMOTOPY_STEPS EQUAL INT_NUMBER { driver.option_num("homotopy_steps", $3); };
 o_homotopy_force_continue: HOMOTOPY_FORCE_CONTINUE EQUAL INT_NUMBER { driver.option_num("homotopy_force_continue", $3); };
+o_preconditioner : PRECONDITIONER EQUAL UMFITER
+                   { driver.option_str("simul.preconditioner", "umfiter"); }
+                 | PRECONDITIONER EQUAL ITERSTACK
+                   { driver.option_str("simul.preconditioner", "iterstack"); }
+                 | PRECONDITIONER EQUAL ILU
+                   { driver.option_str("simul.preconditioner", "ilu"); };
+                 ;
+o_iter_tol : ITER_TOL EQUAL non_negative_number  { driver.option_num("simul.iter_tol", $3); };
+o_iter_maxit : ITER_MAXIT EQUAL INT_NUMBER  { driver.option_num("simul.iter_maxit", $3); };
+o_gmres_restart : GMRES_RESTART EQUAL INT_NUMBER  { driver.option_num("simul.gmres_restart", $3); };
+o_iterstack_maxlu : ITERSTACK_MAXLU EQUAL INT_NUMBER  { driver.option_num("simul.iterstack_maxlu", $3); };
+o_iterstack_nperiods : ITERSTACK_NPERIODS EQUAL INT_NUMBER  { driver.option_num("simul.iterstack_nperiods", $3); };
+o_iterstack_nlu : ITERSTACK_NLU EQUAL INT_NUMBER  { driver.option_num("simul.iterstack_nlu", $3); };
+o_iterstack_relu : ITERSTACK_RELU EQUAL non_negative_number  { driver.option_num("simul.iterstack_relu", $3); };
+o_check_jacobian_singularity : CHECK_JACOBIAN_SINGULARITY { driver.option_num("simul.check_jacobian_singularity", "true"); };
 o_nocheck : NOCHECK { driver.option_num("steadystate.nocheck", "true"); };
 o_controlled_varexo : CONTROLLED_VAREXO EQUAL '(' symbol_list ')' { driver.option_symbol_list("controlled_varexo", $4); };
@@ -4269,10 +4382,11 @@ o_analytic_derivation_mode : ANALYTIC_DERIVATION_MODE EQUAL signed_number { driv
 o_endogenous_prior : ENDOGENOUS_PRIOR { driver.option_num("endogenous_prior", "true"); }
 o_use_univariate_filters_if_singularity_is_detected : USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED EQUAL INT_NUMBER { driver.option_num("use_univariate_filters_if_singularity_is_detected", $3); }
 o_mcmc_jumping_covariance : MCMC_JUMPING_COVARIANCE EQUAL HESSIAN
-                            { driver.option_str("MCMC_jumping_covariance", $3); }                                  | MCMC_JUMPING_COVARIANCE EQUAL PRIOR_VARIANCE
+                            { driver.option_str("MCMC_jumping_covariance", "hessian"); }
-                            { driver.option_str("MCMC_jumping_covariance", $3); }
+                          | MCMC_JUMPING_COVARIANCE EQUAL PRIOR_VARIANCE
+                            { driver.option_str("MCMC_jumping_covariance", "prior_variance"); }
                          | MCMC_JUMPING_COVARIANCE EQUAL IDENTITY_MATRIX
-                            { driver.option_str("MCMC_jumping_covariance", $3); }
+                            { driver.option_str("MCMC_jumping_covariance", "identity_matrix"); }
                          | MCMC_JUMPING_COVARIANCE EQUAL filename
                            { driver.option_str("MCMC_jumping_covariance", $3); }
                          ;
@@ -4343,7 +4457,7 @@ o_emas_girf : EMAS_GIRF { driver.option_num("irf_opt.ergodic_mean_irf", "true");
 o_emas_drop : EMAS_DROP EQUAL INT_NUMBER { driver.option_num("irf_opt.EM.drop", $3); };
 o_emas_tolf : EMAS_TOLF EQUAL non_negative_number { driver.option_num("irf_opt.EM.tolf", $3); };
 o_emas_max_iter : EMAS_MAX_ITER EQUAL INT_NUMBER { driver.option_num("irf_opt.EM.iter", $3); };
-o_non_zero : NON_ZERO { driver.option_num("non_zero", "true"); };
+o_non_zero : NON_ZERO { driver.option_num("steady.non_zero", "true"); };
 // Some options to "identification"
 o_no_identification_strength : NO_IDENTIFICATION_STRENGTH { driver.option_num("no_identification_strength", "true"); };
@@ -4433,9 +4547,8 @@ vec_int_number : INT_NUMBER
 vec_int_elem : vec_int_number
             | INT_NUMBER ':' INT_NUMBER
               {
-                 $$ = {};
+                 auto v = views::iota(stoi($1), stoi($3) + 1);
-                 for (int i = stoi($1); i <= stoi($3); i++)
+                 $$ = {v.begin(), v.end()};
-                   $$.push_back(i);
               }
             ;
@@ -4563,6 +4676,7 @@ symbol : NAME
       | ADD
       | MULTIPLY
       | MFS
+       | RESIDUAL
       ;
 %%

--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
 /* -*- C++ -*- */
 /*
- * Copyright © 2003-2024 Dynare Team
+ * Copyright © 2003-2025 Dynare Team
 *
 * This file is part of Dynare.
 *
@@ -20,7 +20,6 @@
 %{
-#include <cstring>
 #include "ParsingDriver.hh"
 using namespace std;
@@ -60,7 +59,6 @@ string eofbuff;
 %x VERBATIM_BLOCK
 %x NATIVE
 %x NATIVE_COMMENT
-%x DATES_STATEMENT
 %x LINE1
 %x LINE2
 %x LINE3
@@ -70,7 +68,9 @@ string eofbuff;
 #define YY_USER_ACTION location_increment(yylloc, yytext);
 %}
-DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
+NAME [a-z_][a-z0-9_]*
+FLOAT_NUMBER ((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+)
+DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4]|[sh][12])
 %%
 /* Code put at the beginning of yylex() */
@@ -92,12 +92,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
                }
 /* spaces, tabs and carriage returns are ignored */
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,DATES_STATEMENT,LINE1,LINE2,LINE3>[[:space:]]+  { yylloc->step(); }
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,LINE1,LINE2,LINE3>[[:space:]]+  { yylloc->step(); }
 /* Comments */
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>%.*
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>%.*
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>"//".*
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>"//".*
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,DATES_STATEMENT>"/*"   {comment_caller = YY_START; BEGIN COMMENT;}
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK>"/*"   {comment_caller = YY_START; BEGIN COMMENT;}
 <COMMENT>"*/"        {BEGIN comment_caller;}
 <COMMENT>.
@@ -111,6 +111,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>predetermined_variables {BEGIN DYNARE_STATEMENT; return token::PREDETERMINED_VARIABLES;}
 <INITIAL>parameters {BEGIN DYNARE_STATEMENT; return token::PARAMETERS;}
 <INITIAL>model_local_variable {BEGIN DYNARE_STATEMENT; return token::MODEL_LOCAL_VARIABLE;}
+<INITIAL>heterogeneity_dimension {BEGIN DYNARE_STATEMENT; return token::HETEROGENEITY_DIMENSION;}
 <INITIAL>model_info {BEGIN DYNARE_STATEMENT; return token::MODEL_INFO;}
 <INITIAL>estimation {BEGIN DYNARE_STATEMENT; return token::ESTIMATION;}
 <INITIAL>set_time {BEGIN DYNARE_STATEMENT; return token::SET_TIME;}
@@ -235,6 +236,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>pac_target_info {BEGIN DYNARE_BLOCK; return token::PAC_TARGET_INFO;}
 <INITIAL>matched_irfs {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS;}
 <INITIAL>matched_irfs_weights {BEGIN DYNARE_BLOCK; return token::MATCHED_IRFS_WEIGHTS;}
+<INITIAL>perfect_foresight_controlled_paths {BEGIN DYNARE_BLOCK; return token::PERFECT_FORESIGHT_CONTROLLED_PATHS;}
 /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
@@ -245,7 +247,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>subsamples {return token::SUBSAMPLES;}
 <DYNARE_STATEMENT>options {return token::OPTIONS;}
 <DYNARE_STATEMENT>prior {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::PRIOR;
 }
 <INITIAL>std {BEGIN DYNARE_STATEMENT; return token::STD;}
@@ -255,31 +257,12 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>prior_function {BEGIN DYNARE_STATEMENT; return token::PRIOR_FUNCTION;}
 <INITIAL>posterior_function {BEGIN DYNARE_STATEMENT; return token::POSTERIOR_FUNCTION;}
- /* Inside  of a Dynare statement */
+<DYNARE_STATEMENT,DYNARE_BLOCK>{DATE} {
-<DYNARE_STATEMENT>{DATE} {
+  yylval->emplace<string>(yytext);
-                           char* yycopy = strdup(yytext);
+  return token::DATE;
-                           char* uput = yycopy + yyleng;
-                           unput(')');
-                           unput('\'');
-                           while (uput > yycopy)
-                             unput(*--uput);
-                           unput('\'');
-                           unput('(');
-                           unput('s');
-                           unput('e');
-                           unput('t');
-                           unput('a');
-                           unput('d');
-                           free( yycopy );
 }
-<DYNARE_STATEMENT>${DATE} { yylloc->step();
-#if (YY_FLEX_MAJOR_VERSION > 2) || (YY_FLEX_MAJOR_VERSION == 2 && YY_FLEX_MINOR_VERSION >= 6)
+ /* Inside a Dynare statement */
-                            yyout << yytext + 1;
-#else
-                            *yyout << yytext + 1;
-#endif
-                          }
-<DYNARE_STATEMENT>dates  {dates_parens_nb=0; BEGIN DATES_STATEMENT; yylval->build<string>("dates");}
 <DYNARE_STATEMENT>file                  {return token::FILE;}
 <DYNARE_STATEMENT>datafile 		{return token::DATAFILE;}
 <DYNARE_STATEMENT>dirname       {return token::DIRNAME;}
@@ -312,10 +295,10 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>posterior_nograph   		{return token::POSTERIOR_NOGRAPH;}
 <DYNARE_STATEMENT>nodisplay     {return token::NODISPLAY;}
 <DYNARE_STATEMENT>graph_format  {return token::GRAPH_FORMAT;}
-<DYNARE_STATEMENT>eps  {yylval->build<string>(yytext); return token::EPS;}
+<DYNARE_STATEMENT>eps  {yylval->emplace<string>(yytext); return token::EPS;}
-<DYNARE_STATEMENT>pdf  {yylval->build<string>(yytext); return token::PDF;}
+<DYNARE_STATEMENT>pdf  {yylval->emplace<string>(yytext); return token::PDF;}
-<DYNARE_STATEMENT>fig  {yylval->build<string>(yytext); return token::FIG;}
+<DYNARE_STATEMENT>fig  {yylval->emplace<string>(yytext); return token::FIG;}
-<DYNARE_STATEMENT>none  {yylval->build<string>(yytext); return token::NONE;}
+<DYNARE_STATEMENT>none  {yylval->emplace<string>(yytext); return token::NONE;}
 <DYNARE_STATEMENT>print   		{return token::PRINT;}
 <DYNARE_STATEMENT>noprint   		{return token::NOPRINT;}
 <DYNARE_STATEMENT>conf_sig  		{return token::CONF_SIG;}
@@ -423,6 +406,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>logarithmic_reduction {return token::LOGARITHMIC_REDUCTION;}
 <DYNARE_STATEMENT>use_univariate_filters_if_singularity_is_detected {return token::USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED;}
 <DYNARE_STATEMENT>hybrid {return token::HYBRID;}
+<DYNARE_STATEMENT>use_first_order_solution {return token::USE_FIRST_ORDER_SOLUTION;}
 <DYNARE_STATEMENT>default {return token::DEFAULT;}
 <DYNARE_STATEMENT>init2shocks {return token::INIT2SHOCKS;}
@@ -435,6 +419,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>kitagawa {return token::KITAGAWA;}
 <DYNARE_STATEMENT>smooth {return token::SMOOTH;}
 <DYNARE_STATEMENT>stratified {return token::STRATIFIED;}
+<DYNARE_STATEMENT>residual {
+  yylval->emplace<string>(yytext);
+  return token::RESIDUAL;
+}
+<DYNARE_STATEMENT>multinomial {return token::MULTINOMIAL;}
 <DYNARE_STATEMENT>cpf_weights {return token::CPF_WEIGHTS;}
 <DYNARE_STATEMENT>amisanotristani {return token::AMISANOTRISTANI;}
 <DYNARE_STATEMENT>murrayjonesparslow {return token::MURRAYJONESPARSLOW;}
@@ -456,43 +445,40 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>fsolve_options {return token::FSOLVE_OPTIONS;}
 <DYNARE_STATEMENT>alpha {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ALPHA;
 }
 <DYNARE_STATEMENT>beta {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::BETA;
 }
 <DYNARE_STATEMENT>gamma {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::GAMMA;
 }
 <DYNARE_STATEMENT>inv_gamma {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INV_GAMMA;
 }
 <DYNARE_STATEMENT>inv_gamma1 {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INV_GAMMA1;
 }
 <DYNARE_STATEMENT>inv_gamma2 {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INV_GAMMA2;
 }
 <DYNARE_STATEMENT>dirichlet {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DIRICHLET;
 }
-<DYNARE_STATEMENT>weibull {
+<DYNARE_STATEMENT>weibull {return token::WEIBULL;}
-  yylval->build<string>(yytext);
-  return token::WEIBULL;
-}
 <DYNARE_STATEMENT>normal {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NORMAL;
 }
 <DYNARE_STATEMENT>uniform {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::UNIFORM;
 }
 <DYNARE_STATEMENT>gsig2_lmdm {return token::GSIG2_LMDM;}
@@ -503,13 +489,13 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>ncsk {return token::NCSK;}
 <DYNARE_STATEMENT>nstd {return token::NSTD;}
 <DYNARE_STATEMENT>ninv {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NINV;
 }
 <DYNARE_STATEMENT>indxparr {return token::INDXPARR;}
 <DYNARE_STATEMENT>indxovr {return token::INDXOVR;}
 <DYNARE_STATEMENT>aband {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ABAND;
 }
 <DYNARE_STATEMENT>write_equation_tags {return token::WRITE_EQUATION_TAGS;}
@@ -526,18 +512,18 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>indxgdls {return token::INDXGDLS;}
 <DYNARE_STATEMENT>eq_ms {return token::EQ_MS;}
 <DYNARE_STATEMENT>cms {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::CMS;
 }
 <DYNARE_STATEMENT>ncms {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NCMS;
 }
 <DYNARE_STATEMENT>eq_cms {return token::EQ_CMS;}
 <DYNARE_STATEMENT>tlindx {return token::TLINDX;}
 <DYNARE_STATEMENT>tlnumber {return token::TLNUMBER;}
 <DYNARE_STATEMENT>cnum {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::CNUM;
 }
 <DYNARE_STATEMENT>nodecomposition {return token::NODECOMPOSITION;};
@@ -619,18 +605,9 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>substitute_libs {return token::SUBSTITUTE_LIBS;}
 <DYNARE_STATEMENT>compiler {return token::COMPILER;}
 <DYNARE_STATEMENT>instruments {return token::INSTRUMENTS;}
-<DYNARE_STATEMENT>hessian  {
+<DYNARE_STATEMENT>hessian  {return token::HESSIAN;}
-  yylval->build<string>(yytext);
+<DYNARE_STATEMENT>prior_variance  {return token::PRIOR_VARIANCE;}
-  return token::HESSIAN;
+<DYNARE_STATEMENT>identity_matrix  {return token::IDENTITY_MATRIX;}
-}
-<DYNARE_STATEMENT>prior_variance  {
-  yylval->build<string>(yytext);
-  return token::PRIOR_VARIANCE;
-}
-<DYNARE_STATEMENT>identity_matrix  {
-  yylval->build<string>(yytext);
-  return token::IDENTITY_MATRIX;
-}
 <DYNARE_STATEMENT>mcmc_jumping_covariance {return token::MCMC_JUMPING_COVARIANCE;}
 /* These four (var, varexo, varexo_det, parameters) are for change_type */
@@ -730,14 +707,6 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>lmmcp {return token::LMMCP;}
 <DYNARE_STATEMENT>additional_optimizer_steps	{return token::ADDITIONAL_OPTIMIZER_STEPS;}
 <DYNARE_STATEMENT>bartlett_kernel_lag {return token::BARTLETT_KERNEL_LAG; }
-<DYNARE_STATEMENT>optimal {
-  yylval->build<string>(yytext);
-  return token::OPTIMAL;
-}
-<DYNARE_STATEMENT>diagonal  {
-  yylval->build<string>(yytext);
-  return token::DIAGONAL;
-}
 <DYNARE_STATEMENT>gmm {return token::GMM;}
 <DYNARE_STATEMENT>smm {return token::SMM;}
 <DYNARE_STATEMENT>irf_matching {return token::IRF_MATCHING;}
@@ -781,10 +750,21 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>with_epilogue {return token::WITH_EPILOGUE;}
 <DYNARE_STATEMENT>heteroskedastic_filter {return token::HETEROSKEDASTIC_FILTER;}
 <DYNARE_STATEMENT>non_zero {return token::NON_ZERO;}
+<DYNARE_STATEMENT>preconditioner {return token::PRECONDITIONER;}
+<DYNARE_STATEMENT>umfiter {return token::UMFITER;}
+<DYNARE_STATEMENT>iterstack {return token::ITERSTACK;}
+<DYNARE_STATEMENT>ilu {return token::ILU;}
+<DYNARE_STATEMENT>iter_tol {return token::ITER_TOL;}
+<DYNARE_STATEMENT>iter_maxit {return token::ITER_MAXIT;}
+<DYNARE_STATEMENT>gmres_restart {return token::GMRES_RESTART;}
+<DYNARE_STATEMENT>iterstack_maxlu {return token::ITERSTACK_MAXLU;}
+<DYNARE_STATEMENT>iterstack_nperiods {return token::ITERSTACK_NPERIODS;}
+<DYNARE_STATEMENT>iterstack_nlu {return token::ITERSTACK_NLU;}
+<DYNARE_STATEMENT>iterstack_relu {return token::ITERSTACK_RELU;}
+<DYNARE_STATEMENT>check_jacobian_singularity {return token::CHECK_JACOBIAN_SINGULARITY;}
 <DYNARE_STATEMENT>\$[^$]*\$ {
-  strtok(yytext + 1, "$");
+  yylval->emplace<string>(yytext + 1).pop_back();
-  yylval->build<string>(yytext + 1);
  return token::TEX_NAME;
 }
@@ -797,20 +777,21 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>periods {return token::PERIODS;}
 <DYNARE_BLOCK>scales {return token::SCALES;}
 <DYNARE_BLOCK>add {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ADD;
 }
 <DYNARE_BLOCK>multiply {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::MULTIPLY;
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>cutoff {return token::CUTOFF;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>mfs {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::MFS;
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>static_mfs {return token::STATIC_MFS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>balanced_growth_test_tol {return token::BALANCED_GROWTH_TEST_TOL;}
+<DYNARE_STATEMENT,DYNARE_BLOCK>heterogeneity {return token::HETEROGENEITY;}
 <DYNARE_BLOCK>gamma_pdf {return token::GAMMA_PDF;}
 <DYNARE_BLOCK>beta_pdf {return token::BETA_PDF;}
 <DYNARE_BLOCK>normal_pdf {return token::NORMAL_PDF;}
@@ -822,19 +803,19 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>dsge_prior_weight {return token::DSGE_PRIOR_WEIGHT;}
 <DYNARE_BLOCK>surprise {return token::SURPRISE;}
 <DYNARE_BLOCK>bind {
-    yylval->build<string>(yytext);
+    yylval->emplace<string>(yytext);
    return token::BIND;
 }
 <DYNARE_BLOCK>relax {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::RELAX;
 }
 <DYNARE_BLOCK>error_bind {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ERROR_BIND;
 }
 <DYNARE_BLOCK>error_relax {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::ERROR_RELAX;
 }
 <DYNARE_BLOCK>relative_to_initval {return token::RELATIVE_TO_INITVAL;}
@@ -848,22 +829,24 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK,DYNARE_STATEMENT>auxname {return token::AUXNAME;}
 <DYNARE_BLOCK>auxname_target_nonstationary {return token::AUXNAME_TARGET_NONSTATIONARY;}
 <DYNARE_BLOCK,DYNARE_STATEMENT>kind {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::KIND;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>ll {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::LL;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>dl {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DL;
 }
 <DYNARE_BLOCK,DYNARE_STATEMENT>dd {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DD;
 }
 <DYNARE_BLOCK>weights {return token::WEIGHTS;}
+<DYNARE_BLOCK>exogenize {return token::EXOGENIZE;}
+<DYNARE_BLOCK>endogenize {return token::ENDOGENIZE;}
 /* Inside Dynare statement */
 <DYNARE_STATEMENT>solve_algo {return token::SOLVE_ALGO;}
@@ -875,13 +858,14 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>order {return token::ORDER;}
 <DYNARE_STATEMENT>lyapunov {return token::LYAPUNOV;}
 <DYNARE_STATEMENT>dr {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::DR;
 }
 <DYNARE_STATEMENT>lyapunov_complex_threshold {return token::LYAPUNOV_COMPLEX_THRESHOLD;}
 <DYNARE_STATEMENT>lyapunov_fixed_point_tol {return token::LYAPUNOV_FIXED_POINT_TOL;}
 <DYNARE_STATEMENT>lyapunov_doubling_tol {return token::LYAPUNOV_DOUBLING_TOL;}
 <DYNARE_STATEMENT>dr_cycle_reduction_tol {return token::DR_CYCLE_REDUCTION_TOL;}
+<DYNARE_STATEMENT>dr_cycle_reduction_maxiter {return token::DR_CYCLE_REDUCTION_MAXITER;}
 <DYNARE_STATEMENT>dr_logarithmic_reduction_tol {return token::DR_LOGARITHMIC_REDUCTION_TOL;}
 <DYNARE_STATEMENT>dr_logarithmic_reduction_maxiter {return token::DR_LOGARITHMIC_REDUCTION_MAXITER;}
 <DYNARE_STATEMENT>replic {return token::REPLIC;}
@@ -938,11 +922,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>time_shift {return token::TIME_SHIFT;}
 <DYNARE_STATEMENT>structural {return token::STRUCTURAL;}
 <DYNARE_STATEMENT>true {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::TRUE;
 }
 <DYNARE_STATEMENT>false {
-    yylval->build<string>(yytext);
+    yylval->emplace<string>(yytext);
    return token::FALSE;
 }
@@ -1010,12 +994,15 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
 <DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
 <DYNARE_BLOCK>pac_target_nonstationary {return token::PAC_TARGET_NONSTATIONARY;}
+<DYNARE_BLOCK>sum {return token::SUM;}
 <DYNARE_STATEMENT>discount {return token::DISCOUNT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>nan {return token::NAN_CONSTANT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>inf {return token::INF_CONSTANT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>constants {return token::CONSTANTS;}
+<DYNARE_BLOCK>⟂ {return token::PERPENDICULAR;}
+<DYNARE_BLOCK>_\|_ {return token::PERPENDICULAR;}
 /* options for GSA module by Marco Ratto */
 <DYNARE_STATEMENT>identification {return token::IDENTIFICATION;}
@@ -1063,38 +1050,27 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <DYNARE_STATEMENT>use_shock_groups {return token::USE_SHOCK_GROUPS;}
 <DYNARE_STATEMENT>colormap {return token::COLORMAP;}
-<DYNARE_STATEMENT,DYNARE_BLOCK>[a-z_][a-z0-9_]* {
+<DYNARE_STATEMENT,DYNARE_BLOCK>{NAME} {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::NAME;
 }
-<DYNARE_STATEMENT,DYNARE_BLOCK>((([0-9]*\.[0-9]+)|([0-9]+\.))([ed][-+]?[0-9]+)?)|([0-9]+[ed][-+]?[0-9]+) {
+<DYNARE_STATEMENT,DYNARE_BLOCK>{FLOAT_NUMBER} {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::FLOAT_NUMBER;
 }
 <DYNARE_STATEMENT,DYNARE_BLOCK>[0-9]+ {
-  yylval->build<string>(yytext);
+  yylval->emplace<string>(yytext);
  return token::INT_NUMBER;
 }
-<DATES_STATEMENT>\( { yylval->as<string>().append(yytext); dates_parens_nb++; }
-<DATES_STATEMENT>\) {
-                      yylval->as<string>().append(yytext);
-                      if (--dates_parens_nb == 0)
-                      {
-                        BEGIN DYNARE_STATEMENT;
-                        return token::DATES;
-                      }
-                    }
-<DATES_STATEMENT>.  { yylval->as<string>().append(yytext); }
 <DYNARE_BLOCK>\|e { return token::PIPE_E; }
 <DYNARE_BLOCK>\|x { return token::PIPE_X; }
 <DYNARE_BLOCK>\|p { return token::PIPE_P; }
 <DYNARE_STATEMENT,DYNARE_BLOCK>\'[^\']*\' {
-  yylval->build<string>(yytext + 1).pop_back();
+  yylval->emplace<string>(yytext + 1).pop_back();
  return token::QUOTED_STRING;
 }
@@ -1126,11 +1102,11 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
    element in initval (in which case Dynare recognizes the matrix name as an external
    function symbol), and may want to modify the matrix later with Matlab statements.
 */
-<INITIAL>[a-z_][a-z0-9_]* {
+<INITIAL>{NAME} {
  if (driver.symbol_exists_and_is_not_modfile_local_or_external_function(yytext))
    {
      BEGIN DYNARE_STATEMENT;
-      yylval->build<string>(yytext);
+      yylval->emplace<string>(yytext);
      return token::NAME;
    }
  else
@@ -1151,7 +1127,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
    be able to back out of the statement if we realize it's a native statement
    and move to the NATIVE context
 */
-<INITIAL>\[([[:space:]]*[a-z_][a-z0-9_]*[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*[a-z_][a-z0-9_]*[[:space:]]*){1}\] {
+<INITIAL>\[([[:space:]]*{NAME}[[:space:]]*,{1}[[:space:]]*)*([[:space:]]*{NAME}[[:space:]]*){1}\] {
  string yytextcpy{yytext};
  yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), '['), yytextcpy.end());
  yytextcpy.erase(remove(yytextcpy.begin(), yytextcpy.end(), ']'), yytextcpy.end());
@@ -1175,7 +1151,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
  if (dynare_statement)
    {
      BEGIN DYNARE_STATEMENT;
-      yylval->build<vector<string>>(val);
+      yylval->emplace<vector<string>>(val);
      return token::SYMBOL_VEC;
    }
 }
@@ -1220,7 +1196,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <NATIVE_COMMENT>"*/"[[:space:]]*\n   { BEGIN NATIVE; }
 <NATIVE_COMMENT>.
-<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,DATES_STATEMENT,LINE1,LINE2,LINE3,NATIVE_COMMENT><<EOF>> { yyterminate(); }
+<INITIAL,DYNARE_STATEMENT,DYNARE_BLOCK,COMMENT,LINE1,LINE2,LINE3,NATIVE_COMMENT><<EOF>> { yyterminate(); }
 <*>.      { driver.error(*yylloc, "character unrecognized by lexer"); }
 %%

--- a/src/DynareMain.cc
+++ b/src/DynareMain.cc
@@ -56,7 +56,7 @@ usage()
          "[conffile=path_to_config_file] [parallel_follower_open_mode] "
          "[parallel_test] [parallel_use_psexec=true|false]"
       << " [-D<variable>[=<value>]] [-I/path] [nostrict] [stochastic] [fast] [minimal_workspace] "
-          "[compute_xrefs] [output=second|third] [language=matlab|julia]"
+          "[compute_xrefs] [output=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>]"
       << " [json=parse|check|transform|compute] [jsonstdout] [onlyjson] [jsonderivsimple] "
@@ -320,7 +320,9 @@ main(int argc, char** argv)
          s.erase(0, 7);
-          if (s == "second")
+          if (s == "first")
+            output_mode = OutputType::first;
+          else if (s == "second")
            output_mode = OutputType::second;
          else if (s == "third")
            output_mode = OutputType::third;
No results found