Add new “matched_moments” block

Closes: #51

src/ComputingTasks.cc

@@ -5146,3 +5146,65 @@ VarExpectationModelStatement::writeJsonOutput(ostream &output) const
   discount->writeOutput(output);
   output << R"("})";
 }
+
+MatchedMomentsStatement::MatchedMomentsStatement(const SymbolTable &symbol_table_arg,
+                                                 vector<tuple<vector<int>, vector<int>, vector<int>>> moments_arg) :
+  symbol_table{symbol_table_arg}, moments{move(moments_arg)}
+{
+}
+
+void
+MatchedMomentsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  output << "M_.matched_moments = {" << endl;
+  for (const auto &[symb_ids, lags, powers] : moments)
+    {
+      output << "  [";
+      for (int s : symb_ids)
+        output << symbol_table.getTypeSpecificID(s)+1 << ',';
+      output << "], [";
+      for (int l : lags)
+        output << l << ',';
+      output << "], [";
+      for (int p : powers)
+        output << p << ',';
+      output << "]," << endl;
+    }
+  output << "};" << endl;
+}
+
+void
+MatchedMomentsStatement::writeJsonOutput(ostream &output) const
+{
+  output << R"({"statementName": "matched_moments", "moments": [)" << endl;
+  for (auto it = moments.begin(); it != moments.end(); ++it)
+    {
+      const auto &[symb_ids, lags, powers] = *it;
+      output << R"(  { "endos": [)";
+      for (auto it2 = symb_ids.begin(); it2 != symb_ids.end(); ++it2)
+        {
+          if (it2 != symb_ids.begin())
+            output << ',';
+          output << symbol_table.getTypeSpecificID(*it2)+1;
+        }
+      output << R"(], "lags": [)";
+      for (auto it2 = lags.begin(); it2 != lags.end(); ++it2)
+        {
+          if (it2 != lags.begin())
+            output << ',';
+          output << *it2;
+        }
+      output << R"(], "powers": [)";
+      for (auto it2 = powers.begin(); it2 != powers.end(); ++it2)
+        {
+          if (it2 != powers.begin())
+            output << ',';
+          output << *it2;
+        }
+      output << "]}";
+      if (next(it) != moments.end())
+        output << ',';
+      output << endl;
+    }
+  output << "]}" << endl;
+}

src/ComputingTasks.hh

@@ -1198,4 +1198,18 @@ public:
   void writeJsonOutput(ostream &output) const override;
 };
 
+class MatchedMomentsStatement : public Statement
+{
+private:
+  const SymbolTable &symbol_table;
+public:
+  /* Each moment is represented by a three vectors: symb_ids, lags, powers.
+     See the definition of ExprNode::matchMatchedMoment() for more details */
+  const vector<tuple<vector<int>, vector<int>, vector<int>>> moments;
+  MatchedMomentsStatement(const SymbolTable &symbol_table_arg,
+                          vector<tuple<vector<int>, vector<int>, vector<int>>> moments_arg);
+  void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
+  void writeJsonOutput(ostream &output) const override;
+};
+
 #endif

src/DynareBison.yy

@@ -94,7 +94,8 @@ class ParsingDriver;
 %token INV_GAMMA_PDF INV_GAMMA1_PDF INV_GAMMA2_PDF IRF IRF_SHOCKS IRF_PLOT_THRESHOLD IRF_CALIBRATION
 %token FAST_KALMAN_FILTER KALMAN_ALGO KALMAN_TOL DIFFUSE_KALMAN_TOL SUBSAMPLES OPTIONS TOLF TOLX PLOT_INIT_DATE PLOT_END_DATE
 %token LAPLACE LIK_ALGO LIK_INIT LINEAR LINEAR_DECOMPOSITION LOAD_IDENT_FILES LOAD_MH_FILE LOAD_RESULTS_AFTER_LOAD_MH LOAD_PARAMS_AND_STEADY_STATE LOGLINEAR LOGDATA LYAPUNOV LINEAR_APPROXIMATION
-%token LYAPUNOV_COMPLEX_THRESHOLD LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR MARKOWITZ MARGINAL_DENSITY MAX MAXIT
+%token LYAPUNOV_COMPLEX_THRESHOLD LYAPUNOV_FIXED_POINT_TOL LYAPUNOV_DOUBLING_TOL LOG_DEFLATOR LOG_TREND_VAR LOG_GROWTH_FACTOR
+%token MATCHED_MOMENTS MARKOWITZ MARGINAL_DENSITY MAX MAXIT
 %token MFS MH_CONF_SIG MH_DROP MH_INIT_SCALE MH_JSCALE MH_TUNE_JSCALE MH_TUNE_GUESS MH_MODE MH_NBLOCKS MH_REPLIC MH_RECOVER 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 MCMC_JUMPING_COVARIANCE MOMENT_CALIBRATION
@@ -195,6 +196,7 @@ class ParsingDriver;
 %type <pair<string,string>> named_var_elem subsamples_eq_opt integer_range_w_inf
 %type <vector<pair<string,string>>> named_var named_var_1
 %type <tuple<string,string,string,string>> prior_eq_opt options_eq_opt
+%type <vector<expr_t>> matched_moments_list
 %%
 
 %start statement_list;
@@ -316,6 +318,7 @@ statement : parameters
           | det_cond_forecast
           | var_expectation_model
           | compilation_setup
+          | matched_moments
           ;
 
 dsample : DSAMPLE INT_NUMBER ';'
@@ -930,6 +933,19 @@ compilation_setup_option : SUBSTITUTE_FLAGS EQUAL QUOTED_STRING
                            { driver.compilation_setup_compiler($3); }
                          ;
 
+matched_moments : MATCHED_MOMENTS ';' { driver.begin_matched_moments(); }
+                  matched_moments_list END ';' { driver.end_matched_moments($4); }
+                ;
+
+matched_moments_list : hand_side ';'
+                       { $$ = { $1 }; }
+                     | matched_moments_list hand_side ';'
+                       {
+                         $$ = $1;
+                         $$.push_back($2);
+                       }
+                     ;
+
 model_options : BLOCK { driver.block(); }
               | o_cutoff
               | o_mfs

src/DynareFlex.ll

@@ -223,6 +223,7 @@ DATE -?[0-9]+([ya]|m([1-9]|1[0-2])|q[1-4])
 <INITIAL>ramsey_constraints {BEGIN DYNARE_BLOCK; return token::RAMSEY_CONSTRAINTS;}
 <INITIAL>restrictions {BEGIN DYNARE_BLOCK; return token::RESTRICTIONS;}
 <INITIAL>generate_irfs {BEGIN DYNARE_BLOCK; return token::GENERATE_IRFS;}
+<INITIAL>matched_moments {BEGIN DYNARE_BLOCK; return token::MATCHED_MOMENTS;}
 
  /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}

src/ExprNode.cc

@@ -387,6 +387,13 @@ ExprNode::fillErrorCorrectionRow(int eqn,
     }
 }
 
+void
+ExprNode::matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const
+{
+  throw MatchFailureException{"Unsupported expression"};
+}
+
+
 NumConstNode::NumConstNode(DataTree &datatree_arg, int idx_arg, int id_arg) :
   ExprNode{datatree_arg, idx_arg},
   id{id_arg}
@@ -1968,6 +1975,17 @@ VariableNode::replaceVarsInEquation(map<VariableNode *, NumConstNode *> &table)
   return const_cast<VariableNode *>(this);
 }
 
+void
+VariableNode::matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const
+{
+  if (get_type() != SymbolType::endogenous)
+    throw MatchFailureException{"Variable " + datatree.symbol_table.getName(symb_id) + " is not an endogenous"};
+
+  symb_ids.push_back(symb_id);
+  lags.push_back(lag);
+  powers.push_back(1);
+}
+
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, int idx_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, int param1_symb_id_arg, int param2_symb_id_arg, string adl_param_name_arg, vector<int> adl_lags_arg) :
   ExprNode{datatree_arg, idx_arg},
   arg{arg_arg},
@@ -5566,6 +5584,32 @@ BinaryOpNode::substituteStaticAuxiliaryDefinition() const
   return buildSimilarBinaryOpNode(arg1, arg2subst, datatree);
 }
 
+void
+BinaryOpNode::matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const
+{
+  if (op_code == BinaryOpcode::times)
+    {
+      arg1->matchMatchedMoment(symb_ids, lags, powers);
+      arg2->matchMatchedMoment(symb_ids, lags, powers);
+    }
+  else if (op_code == BinaryOpcode::power)
+    {
+      if (!dynamic_cast<const VariableNode *>(arg1))
+        throw MatchFailureException("First argument of power expression must be a variable");
+      auto ncn = dynamic_cast<const NumConstNode *>(arg2);
+      if (!ncn)
+        throw MatchFailureException("Second argument of power expression must be a positive integer");
+      double c = datatree.num_constants.getDouble(ncn->id);
+      if (c <= 0 || round(c) != c)
+        throw MatchFailureException("Second argument of power expression must be a positive integer");
+      arg1->matchMatchedMoment(symb_ids, lags, powers);
+      powers.back() = static_cast<int>(c);
+    }
+  else
+    throw MatchFailureException("Unsupported binary operator");
+}
+
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, int idx_arg, const expr_t arg1_arg,
                              TrinaryOpcode op_code_arg, const expr_t arg2_arg, const expr_t arg3_arg) :
   ExprNode{datatree_arg, idx_arg},

src/ExprNode.hh

@@ -716,6 +716,13 @@ public:
     {
     };
   };
+
+  /* Match an expression of the form ∏ x(l)ᵏ, where x are endogenous, as used
+     in the match_moments block.
+     For each factor, adds an integer in the 3 vectors in argument (symb_id in
+     the first, lag in the second, exponent in the third).
+     Throws a MatchFailureException if not of the right form. */
+  virtual void matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -879,6 +886,7 @@ public:
   void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const override;
   //! Substitute auxiliary variables by their expression in static model
   expr_t substituteStaticAuxiliaryVariable() const override;
+  void matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const override;
 };
 
 //! Unary operator node
@@ -1126,6 +1134,7 @@ public:
   //! Substitute auxiliary variables by their expression in static model auxiliary variable definition
   expr_t substituteStaticAuxiliaryDefinition() const;
   void decomposeAdditiveTerms(vector<pair<expr_t, int>> &terms, int current_sign) const override;
+  void matchMatchedMoment(vector<int> &symb_ids, vector<int> &lags, vector<int> &powers) const override;
 };
 
 //! Trinary operator node

src/ParsingDriver.cc

@@ -3543,3 +3543,29 @@ ParsingDriver::var_expectation_model()
   var_expectation_model_discount = nullptr;
   var_expectation_model_expression = nullptr;
 }
+
+void
+ParsingDriver::begin_matched_moments()
+{
+  set_current_data_tree(&mod_file->dynamic_model);
+}
+
+void
+ParsingDriver::end_matched_moments(const vector<expr_t> &moments)
+{
+  vector<tuple<vector<int>, vector<int>, vector<int>>> parsed_moments;
+  for (auto m : moments)
+    try
+      {
+        vector<int> symb_ids, lags, powers;
+        m->matchMatchedMoment(symb_ids, lags, powers);
+        parsed_moments.emplace_back(symb_ids, lags, powers);
+      }
+    catch (ExprNode::MatchFailureException &e)
+      {
+        error("Matched moment expression has incorrect format: " + e.message);
+      }
+  mod_file->addStatement(make_unique<MatchedMomentsStatement>(mod_file->symbol_table, parsed_moments));
+
+  reset_data_tree();
+}

src/ParsingDriver.hh

@@ -913,6 +913,10 @@ public:
   void method_of_moments();
   //! Add a var_expectation_model statement
   void var_expectation_model();
+  //! Start parsing a matched_moments block
+  void begin_matched_moments();
+  //! Add a matched_moments block
+  void end_matched_moments(const vector<expr_t> &moments);
 };
 
 #endif // ! PARSING_DRIVER_HH