From 784dd4122afe0f68ac2d0f903c0625b9c7e50b4f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 31 Aug 2021 17:54:20 +0200
Subject: [PATCH] =?UTF-8?q?PAC:=20generate=20the=20=E2=80=9Cgrowth=5Fneutr?=
 =?UTF-8?q?ality=5Fparam=5Findex=E2=80=9D=20MATLAB=20field=20for=20MCE=20m?=
 =?UTF-8?q?odels?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Previously it was only generated in the case of backward PAC models.
---
 src/DynamicModel.cc | 42 ++++++++++++++++++++++++------------------
 src/DynamicModel.hh | 16 +++++++++++++---
 src/ModFile.cc      |  3 +++
 3 files changed, 40 insertions(+), 21 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 81d3612c..78e102e3 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -114,6 +114,7 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
   pac_h1_indices{m.pac_h1_indices},
   pac_mce_z1_symb_ids{m.pac_mce_z1_symb_ids},
   pac_eqtag_and_lag{m.pac_eqtag_and_lag},
+  pac_growth_neutrality_params{m.pac_growth_neutrality_params},
   pac_model_info{m.pac_model_info},
   pac_equation_info{m.pac_equation_info}
 {
@@ -178,6 +179,7 @@ DynamicModel::operator=(const DynamicModel &m)
   pac_mce_z1_symb_ids = m.pac_mce_z1_symb_ids;
   pac_eqtag_and_lag = m.pac_eqtag_and_lag;
   pac_expectation_substitution.clear();
+  pac_growth_neutrality_params = m.pac_growth_neutrality_params;
   pac_model_info = m.pac_model_info;
   pac_equation_info = m.pac_equation_info;
 
@@ -3033,20 +3035,19 @@ DynamicModel::writeDriverOutput(ostream &output, const string &basename, bool bl
       output << "];" << endl;
     }
 
-  for (auto &it : pac_model_info)
+  for (auto &[model, growth_neutrality_param_index] : pac_growth_neutrality_params)
+    output << "M_.pac." << model << ".growth_neutrality_param_index = "
+           << symbol_table.getTypeSpecificID(growth_neutrality_param_index) + 1 << ";" << endl;
+
+  for (auto &[model, values] : pac_model_info)
     {
-      vector<int> lhs = get<0>(it.second);
-      output << "M_.pac." << it.first << ".lhs = [";
+      auto &[lhs, aux_model_type] = values;
+      output << "M_.pac." << model << ".lhs = [";
       for (auto it : lhs)
         output << it + 1 << " ";
       output << "];" << endl;
 
-      if (int growth_param_index = get<1>(it.second);
-          growth_param_index >= 0)
-        output << "M_.pac." << it.first << ".growth_neutrality_param_index = "
-               << symbol_table.getTypeSpecificID(growth_param_index) + 1 << ";" << endl;
-
-      output << "M_.pac." << it.first << ".auxiliary_model_type = '" << get<2>(it.second) << "';" << endl;
+      output << "M_.pac." << model << ".auxiliary_model_type = '" << aux_model_type << "';" << endl;
     }
 
   for (auto &pit : pac_equation_info)
@@ -4194,6 +4195,14 @@ DynamicModel::addPacModelConsistentExpectationEquation(const string &name, int d
   cout << "Pac Model Consistent Expectation: added " << neqs << " auxiliary variables and equations for model " << name << "." << endl;
 }
 
+void
+DynamicModel::createPacGrowthNeutralityParameter(const string &pac_model_name)
+{
+  int param_idx = symbol_table.addSymbol(pac_model_name +"_pac_growth_neutrality_correction",
+                                         SymbolType::parameter);
+  pac_growth_neutrality_params[pac_model_name] = param_idx;
+}
+
 void
 DynamicModel::fillPacModelInfo(const string &pac_model_name,
                                vector<int> lhs,
@@ -4208,12 +4217,6 @@ DynamicModel::fillPacModelInfo(const string &pac_model_name,
   bool stationary_vars_present = any_of(nonstationary.begin(), nonstationary.end(), logical_not<bool>());
   bool nonstationary_vars_present = any_of(nonstationary.begin(), nonstationary.end(), [](bool b) { return b; }); // FIXME: use std::identity instead of an anonymous function when we upgrade to C++20
 
-  int growth_param_index = -1;
-  if (growth)
-    growth_param_index = symbol_table.addSymbol(pac_model_name
-                                                +"_pac_growth_neutrality_correction",
-                                                SymbolType::parameter);
-
   for (auto pac_models_and_eqtags : pac_eqtag_and_lag)
     {
       if (pac_models_and_eqtags.first.first != pac_model_name)
@@ -4276,12 +4279,15 @@ DynamicModel::fillPacModelInfo(const string &pac_model_name,
         }
 
       if (growth)
-        subExpr = AddPlus(subExpr,
-                          AddTimes(AddVariable(growth_param_index), growth));
+        {
+          int growth_param_index = pac_growth_neutrality_params.at(pac_model_name);
+          subExpr = AddPlus(subExpr,
+                            AddTimes(AddVariable(growth_param_index), growth));
+        }
 
       pac_expectation_substitution[{pac_model_name, eqtag}] = subExpr;
     }
-  pac_model_info[pac_model_name] = {move(lhs), growth_param_index, move(aux_model_type)};
+  pac_model_info[pac_model_name] = {move(lhs), move(aux_model_type)};
 }
 
 void
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index dead4814..7b97e8d4 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -134,13 +134,16 @@ private:
   //! Store lag info for pac equations
   //! (pac_model_name, equation_tag) -> (standardized_eqtag, lag)
   map<pair<string, string>, pair<string, int>> pac_eqtag_and_lag;
+  // Store indices for growth neutrality parameters
+  // pac_model_name -> growth_param_index
+  map<string, int> pac_growth_neutrality_params;
 
   //! (pac_model_name, equation_tag) -> expr_t
   map<pair<string, string>, expr_t> pac_expectation_substitution;
 
-  //! Store info about pac models:
-  //! pac_model_name -> (lhsvars, growth_param_index, aux_model_type)
-  map<string, tuple<vector<int>, int, string>> pac_model_info;
+  //! Store info about backward PAC models:
+  //! pac_model_name -> (lhsvars, aux_model_type)
+  map<string, pair<vector<int>, string>> pac_model_info;
 
   //! Store info about pac models specific to the equation they appear in
   //! (pac_model_name, standardized_eqtag) ->
@@ -414,7 +417,14 @@ public:
 
   //! Get Pac equation parameter info
   map<pair<string, string>, pair<string, int>> walkPacParameters(const string &name);
+
+  // Create the growth neutrality parameter of a given PAC model (when the
+  // “growth” option has been passed)
+  void createPacGrowthNeutralityParameter(const string &pac_model_name);
+
   //! Add var_model info to pac_expectation nodes
+  // Must be called after createPacGrowthNeutralityParameter() has been called
+  // on the relevant models
   void fillPacModelInfo(const string &pac_model_name,
                         vector<int> lhs,
                         int max_lag,
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 68081881..4cd14be5 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -484,6 +484,9 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, bool
             cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl;
             exit(EXIT_FAILURE);
           }
+        if (pms->growth)
+          dynamic_model.createPacGrowthNeutralityParameter(pms->name);
+
         auto eqtag_and_lag = dynamic_model.walkPacParameters(pms->name);
         if (pms->aux_model_name.empty())
           dynamic_model.addPacModelConsistentExpectationEquation(pms->name, symbol_table.getID(pms->discount),
-- 
GitLab