From d2cbe2a19e55876e90d70a93d3008002810b1ae7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 11 Sep 2020 16:45:19 +0200
Subject: [PATCH] PAC model: fixes to M_.pac.MODEL.ar (information about
 autoregressive part)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

– Fix order of items in this structure. Previously, items were ordered
  according to the declaration order of parameters. Now, items are order
  according to lag order (first lag appears first)
– Gracefully handle the case where there is no autoregressive part
  (Closes: #52)
---
 src/DynamicModel.cc | 40 +++++++------------------
 src/DynamicModel.hh |  5 +---
 src/ExprNode.cc     | 73 ++++++---------------------------------------
 src/ExprNode.hh     | 18 ++++-------
 src/ModFile.cc      |  1 -
 5 files changed, 25 insertions(+), 112 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 27f5d3cf..3ae47128 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3144,16 +3144,16 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl
           }
       output << "];" << endl
              << modstruct << "pac." << substruct << "ar.params = [";
-      for (auto &it : ar_params_and_vars)
-        output << symbol_table.getTypeSpecificID(it.first) + 1 << " ";
+      for (auto &[pid, vid, vlag] : ar_params_and_vars)
+        output << (pid != -1 ? symbol_table.getTypeSpecificID(pid) + 1 : -1) << " ";
       output << "];" << endl
              << modstruct << "pac." << substruct << "ar.vars = [";
-      for (auto &it : ar_params_and_vars)
-        output << symbol_table.getTypeSpecificID(it.second.first) + 1 << " ";
+      for (auto &[pid, vid, vlag] : ar_params_and_vars)
+        output << (vid != -1 ? symbol_table.getTypeSpecificID(vid) + 1 : -1) << " ";
       output << "];" << endl
              << modstruct << "pac." << substruct << "ar.lags = [";
-      for (auto &it : ar_params_and_vars)
-        output << it.second.second << " ";
+      for (auto &[pid, vid, vlag] : ar_params_and_vars)
+        output << vlag << " ";
       output << "];" << endl;
       if (!non_optim_vars_params_and_constants.empty())
         {
@@ -3882,7 +3882,7 @@ DynamicModel::walkPacParameters(const string &name)
     {
       pair<int, int> lhs(-1, -1);
       pair<int, vector<tuple<int, bool, int>>> ec_params_and_vars;
-      set<pair<int, pair<int, int>>> ar_params_and_vars;
+      vector<tuple<int, int, int>> ar_params_and_vars;
       vector<tuple<int, int, int, double>> non_optim_vars_params_and_constants, optim_additive_vars_params_and_constants, additive_vars_params_and_constants;
 
       if (equation->containsPacExpectation(name))
@@ -3954,7 +3954,7 @@ DynamicModel::walkPacParameters(const string &name)
               cerr << "walkPacParameters: error obtaining LHS variable." << endl;
               exit(EXIT_FAILURE);
             }
-          if (ec_params_and_vars.second.empty() || ar_params_and_vars.empty())
+          if (ec_params_and_vars.second.empty())
             {
               cerr << "walkPacParameters: error obtaining RHS parameters." << endl;
               exit(EXIT_FAILURE);
@@ -3965,33 +3965,13 @@ DynamicModel::walkPacParameters(const string &name)
                                            non_optim_vars_params_and_constants,
                                            additive_vars_params_and_constants,
                                            optim_additive_vars_params_and_constants};
-          eqtag_and_lag[{name, eqtag}] = {eq, 0};
+          int max_lag = ar_params_and_vars.size();
+          eqtag_and_lag[{name, eqtag}] = {eq, max_lag};
         }
     }
   return eqtag_and_lag;
 }
 
-void
-DynamicModel::getPacMaxLag(const string &pac_model_name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag) const
-{
-  for (auto &equation : equations)
-    if (equation->containsPacExpectation(pac_model_name))
-      {
-        set<pair<int, int>> endogs;
-        equation->arg1->collectDynamicVariables(SymbolType::endogenous, endogs);
-        if (endogs.size() != 1)
-          {
-            cerr << "The LHS of the PAC equation may only be comprised of one endogenous variable"
-                 << endl;
-            exit(EXIT_FAILURE);
-          }
-
-        string eqtag = equation_tags.getTagValueByEqnAndKey(&equation - &equations[0], "name");
-        string eq = eqtag_and_lag[{pac_model_name, eqtag}].first;
-        eqtag_and_lag[{pac_model_name, eqtag}] = {eq, equation->PacMaxLag(endogs.begin()->first)};
-      }
-}
-
 int
 DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
 {
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 07c7b853..58a2f9b1 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -505,9 +505,6 @@ public:
   //! Substitute VarExpectation operators
   void substituteVarExpectation(const map<string, expr_t> &subst_table);
 
-  //! Return max lag of pac equation
-  void getPacMaxLag(const string &pac_model_name, map<pair<string, string>, pair<string, int>> &eqtag_and_lag) const;
-
   // Exception thrown by getPacTargetSymbId()
   class PacTargetNotIdentifiedException
   {
@@ -553,7 +550,7 @@ public:
   //! (pac_model_name, standardized_eqtag) ->
   //!     (lhs, optim_share_index, ar_params_and_vars, ec_params_and_vars, non_optim_vars_params_and_constants, additive_vars_params_and_constants, optim_additive_vars_params_and_constants)
   map<pair<string, string>,
-      tuple<pair<int, int>, int, set<pair<int, pair<int, int>>>, pair<int, vector<tuple<int, bool, int>>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>>> pac_equation_info;
+      tuple<pair<int, int>, int, vector<tuple<int, int, int>>, pair<int, vector<tuple<int, bool, int>>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>, vector<tuple<int, int, int, double>>>> pac_equation_info;
 
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index da5b3625..4bb67ba7 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -565,12 +565,6 @@ NumConstNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
   return 0;
 }
 
-int
-NumConstNode::PacMaxLag(int lhs_symb_id) const
-{
-  return 0;
-}
-
 expr_t
 NumConstNode::decreaseLeadsLags(int n) const
 {
@@ -1538,17 +1532,6 @@ VariableNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
   return maxLag();
 }
 
-int
-VariableNode::PacMaxLag(int lhs_symb_id) const
-{
-  if (get_type() == SymbolType::modelLocalVariable)
-    return datatree.getLocalVariable(symb_id)->PacMaxLag(lhs_symb_id);
-
-  if (lhs_symb_id == symb_id)
-    return -lag;
-  return 0;
-}
-
 expr_t
 VariableNode::substituteAdl() const
 {
@@ -3243,13 +3226,6 @@ UnaryOpNode::VarMinLag() const
   return arg->VarMinLag();
 }
 
-int
-UnaryOpNode::PacMaxLag(int lhs_symb_id) const
-{
-  //This will never be an UnaryOpcode::diff node
-  return arg->PacMaxLag(lhs_symb_id);
-}
-
 expr_t
 UnaryOpNode::substituteAdl() const
 {
@@ -4935,12 +4911,6 @@ BinaryOpNode::undiff() const
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
-int
-BinaryOpNode::PacMaxLag(int lhs_symb_id) const
-{
-  return max(arg1->PacMaxLag(lhs_symb_id), arg2->PacMaxLag(lhs_symb_id));
-}
-
 expr_t
 BinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -5248,7 +5218,7 @@ BinaryOpNode::findTargetVariable(int lhs_symb_id) const
 void
 BinaryOpNode::getPacAREC(int lhs_symb_id, int lhs_orig_symb_id,
                          pair<int, vector<tuple<int, bool, int>>> &ec_params_and_vars,
-                         set<pair<int, pair<int, int>>> &ar_params_and_vars,
+                         vector<tuple<int, int, int>> &ar_params_and_vars,
                          vector<tuple<int, int, int, double>> &additive_vars_params_and_constants) const
 {
   ec_params_and_vars.first = -1;
@@ -5303,7 +5273,7 @@ BinaryOpNode::getPacAREC(int lhs_symb_id, int lhs_orig_symb_id,
         get<3>(t) *= it.second; // Update sign of constants
 
       int pid = get<0>(m);
-      for (auto [vid, lag, pidtmp, constant] : m.second)
+      for (auto [vid, vlag, pidtmp, constant] : m.second)
         {
           if (pid == -1)
             pid = pidtmp;
@@ -5317,16 +5287,19 @@ BinaryOpNode::getPacAREC(int lhs_symb_id, int lhs_orig_symb_id,
               vidorig == lhs_symb_id || vidorig == lhs_orig_symb_id)
             {
               // This is an autoregressive term
-              if (constant != 1 || pid == -1)
+              if (constant != 1 || pid == -1 || !datatree.symbol_table.isDiffAuxiliaryVariable(vid))
                 {
-                  cerr << "BinaryOpNode::getPacAREC: autoregressive terms must be of the form 'parameter*lagged_variable" << endl;
+                  cerr << "BinaryOpNode::getPacAREC: autoregressive terms must be of the form 'parameter*diff_lagged_variable" << endl;
                   exit(EXIT_FAILURE);
                 }
-              ar_params_and_vars.insert({pid, { vid, lag }});
+              int ar_lag = datatree.symbol_table.getOrigLeadLagForDiffAuxVar(vid);
+              if (static_cast<int>(ar_params_and_vars.size()) < ar_lag)
+                ar_params_and_vars.resize(ar_lag, { -1, -1, 0 });
+              ar_params_and_vars[ar_lag-1] = { pid, vid, vlag };
             }
           else
             // This is a residual additive term
-            additive_vars_params_and_constants.emplace_back(vid, lag, pid, constant);
+            additive_vars_params_and_constants.emplace_back(vid, vlag, pid, constant);
         }
     }
 }
@@ -6198,12 +6171,6 @@ TrinaryOpNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
                  arg3->VarMaxLag(lhs_lag_equiv)));
 }
 
-int
-TrinaryOpNode::PacMaxLag(int lhs_symb_id) const
-{
-  return max(arg1->PacMaxLag(lhs_symb_id), max(arg2->PacMaxLag(lhs_symb_id), arg3->PacMaxLag(lhs_symb_id)));
-}
-
 expr_t
 TrinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -6634,15 +6601,6 @@ AbstractExternalFunctionNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
   return max_lag;
 }
 
-int
-AbstractExternalFunctionNode::PacMaxLag(int lhs_symb_id) const
-{
-  int val = 0;
-  for (auto argument : arguments)
-    val = max(val, argument->PacMaxLag(lhs_symb_id));
-  return val;
-}
-
 expr_t
 AbstractExternalFunctionNode::decreaseLeadsLags(int n) const
 {
@@ -8151,13 +8109,6 @@ VarExpectationNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
   exit(EXIT_FAILURE);
 }
 
-int
-VarExpectationNode::PacMaxLag(int lhs_symb_id) const
-{
-  cerr << "VarExpectationNode::PacMaxLag not implemented." << endl;
-  exit(EXIT_FAILURE);
-}
-
 expr_t
 VarExpectationNode::decreaseLeadsLags(int n) const
 {
@@ -8550,12 +8501,6 @@ PacExpectationNode::VarMaxLag(const set<expr_t> &lhs_lag_equiv) const
   return 0;
 }
 
-int
-PacExpectationNode::PacMaxLag(int lhs_symb_id) const
-{
-  return 0;
-}
-
 expr_t
 PacExpectationNode::decreaseLeadsLags(int n) const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index 4e3b2cd4..32d7d510 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -467,10 +467,6 @@ public:
     maxLagWithDiffsExpanded() returns 3. */
   virtual int maxLagWithDiffsExpanded() const = 0;
 
-  //! Get Max lag of var associated with Pac model
-  //! Takes account of undiffed LHS variables in calculating the max lag
-  virtual int PacMaxLag(int lhs_symb_id) const = 0;
-
   virtual expr_t undiff() const = 0;
 
   //! Returns a new expression where all the leads/lags have been shifted backwards by the same amount
@@ -757,7 +753,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -830,7 +825,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -931,7 +925,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1038,7 +1031,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1086,10 +1078,14 @@ public:
       is target ?, multiplicative scalar); this form theoretically allows for a
       linear combination in the cointegration, though for the time being we allow
       less than that
+    ar_params_and_vars: elements are indexed according to lag (index 0 is lag
+      1); each tuple is (parameter_id, variable_id, variable_lag) where
+      variable_lag is *not* the lag order in the AR
+      (because variable is an AUX_DIFF_LAG)
    */
   void getPacAREC(int lhs_symb_id, int lhs_orig_symb_id,
                   pair<int, vector<tuple<int, bool, int>>> &ec_params_and_vars,
-                  set<pair<int, pair<int, int>>> &ar_params_and_vars,
+                  vector<tuple<int, int, int>> &ar_params_and_vars,
                   vector<tuple<int, int, int, double>> &additive_vars_params_and_constants) const;
 
   //! Finds the share of optimizing agents in the PAC equation,
@@ -1170,7 +1166,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1282,7 +1277,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1453,7 +1447,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   void prepareForDerivation() override;
@@ -1526,7 +1519,6 @@ public:
   int maxLagWithDiffsExpanded() const override;
   int VarMinLag() const override;
   int VarMaxLag(const set<expr_t> &lhs_lag_equiv) const override;
-  int PacMaxLag(int lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   void prepareForDerivation() override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 4d528ebb..b60a35eb 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -494,7 +494,6 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
             exit(EXIT_FAILURE);
           }
         auto eqtag_and_lag = dynamic_model.walkPacParameters(pms->name);
-        original_model.getPacMaxLag(pms->name, eqtag_and_lag);
         if (pms->aux_model_name.empty())
           dynamic_model.addPacModelConsistentExpectationEquation(pms->name, symbol_table.getID(pms->discount),
                                                                  eqtag_and_lag, diff_subst_table);
-- 
GitLab