From 38152c34a478e2ecc1094d3d210dd4cf12a7edd4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Tue, 11 Dec 2018 17:26:50 +0100
Subject: [PATCH] Make histval compatible with diff operator

The idea is to make use of the dynamic_set_auxiliary_dseries.m file to generate
the initial conditions for all auxiliary variables, including the diffs.

Also remove the check done by the preprocessor for the lags in histval, since
it does not work correctly with the diff operator.
---
 src/DynamicModel.cc            |  6 ++++
 src/DynamicModel.hh            |  2 +-
 src/ModFile.cc                 | 17 ----------
 src/NumericalInitialization.cc | 60 +++++++++++-----------------------
 src/NumericalInitialization.hh |  3 --
 src/ParsingDriver.cc           |  5 +--
 src/ParsingDriver.hh           |  2 --
 src/Statement.hh               |  2 --
 8 files changed, 27 insertions(+), 70 deletions(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index dbe07c4d..728030f4 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -121,6 +121,7 @@ DynamicModel::DynamicModel(const DynamicModel &m) :
   max_exo_det_lead {m.max_exo_det_lead},
   max_lag_orig {m.max_lag_orig},
   max_lead_orig {m.max_lead_orig},
+  max_lag_with_diffs_expanded_orig {m.max_lag_with_diffs_expanded_orig},
   max_endo_lag_orig {m.max_endo_lag_orig},
   max_endo_lead_orig {m.max_endo_lead_orig},
   max_exo_lag_orig {m.max_exo_lag_orig},
@@ -181,6 +182,7 @@ DynamicModel::operator=(const DynamicModel &m)
   max_exo_det_lead = m.max_exo_det_lead;
   max_lag_orig = m.max_lag_orig;
   max_lead_orig = m.max_lead_orig;
+  max_lag_with_diffs_expanded_orig = m.max_lag_with_diffs_expanded_orig;
   max_endo_lag_orig = m.max_endo_lag_orig;
   max_endo_lead_orig = m.max_endo_lead_orig;
   max_exo_lag_orig = m.max_exo_lag_orig;
@@ -3036,6 +3038,7 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
          << modstruct << "orig_maximum_exo_det_lead = " << max_exo_det_lead_orig << ";" << endl
          << modstruct << "orig_maximum_lag = " << max_lag_orig << ";" << endl
          << modstruct << "orig_maximum_lead = " << max_lead_orig << ";" << endl
+         << modstruct << "orig_maximum_lag_with_diffs_expanded = " << max_lag_with_diffs_expanded_orig << ";" << endl
          << modstruct << "lead_lag_incidence = [";
   // Loop on endogenous variables
   int nstatic = 0,
@@ -5086,6 +5089,9 @@ DynamicModel::setLeadsLagsOrig()
       equation->collectDynamicVariables(SymbolType::endogenous, dynvars);
       equation->collectDynamicVariables(SymbolType::exogenous, dynvars);
       equation->collectDynamicVariables(SymbolType::exogenousDet, dynvars);
+
+      max_lag_with_diffs_expanded_orig = max(equation->maxLagWithDiffsExpanded(),
+                                             max_lag_with_diffs_expanded_orig);
     }
 
     for (const auto & dynvar : dynvars)
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 2fd42318..02229a66 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -73,7 +73,7 @@ private:
   /*! Set by computeDerivIDs() */
   int max_exo_det_lag{0}, max_exo_det_lead{0};
   //! Maximum lag and lead over all types of variables (positive values) of original model
-  int max_lag_orig{0}, max_lead_orig{0};
+  int max_lag_orig{0}, max_lead_orig{0}, max_lag_with_diffs_expanded_orig{0};
   //! Maximum lag and lead over endogenous variables (positive values) of original model
   int max_endo_lag_orig{0}, max_endo_lead_orig{0};
   //! Maximum lag and lead over exogenous variables (positive values) of original model
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 43d89595..2247f557 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -498,23 +498,6 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
       mod_file_struct.ramsey_eq_nbr = dynamic_model.equation_number() - mod_file_struct.orig_eq_nbr;
     }
 
-  // Workaround for #1193
-  if (!mod_file_struct.hist_vals_wrong_lag.empty())
-    {
-      bool err = false;
-      for (map<int, int>::const_iterator it = mod_file_struct.hist_vals_wrong_lag.begin();
-           it != mod_file_struct.hist_vals_wrong_lag.end(); it++)
-          if (dynamic_model.minLagForSymbol(it->first) > it->second - 1)
-            {
-              cerr << "ERROR: histval: variable " << symbol_table.getName(it->first)
-                   << " does not appear in the model with the lag " << it->second - 1
-                   << " (see the reference manual for the timing convention in 'histval')" << endl;
-              err = true;
-            }
-      if (err)
-        exit(EXIT_FAILURE);
-    }
-
   /* Handle var_expectation_model statements: collect information about them,
      create the new corresponding parameters, and the expressions to replace
      the var_expectation statements.
diff --git a/src/NumericalInitialization.cc b/src/NumericalInitialization.cc
index e6878e09..88a0a9cb 100644
--- a/src/NumericalInitialization.cc
+++ b/src/NumericalInitialization.cc
@@ -304,11 +304,9 @@ EndValStatement::writeJsonOutput(ostream &output) const
 }
 
 HistValStatement::HistValStatement(hist_values_t hist_values_arg,
-                                   const hist_vals_wrong_lag_t hist_vals_wrong_lag_arg,
                                    const SymbolTable &symbol_table_arg,
                                    const bool &all_values_required_arg) :
   hist_values{move(hist_values_arg)},
-  hist_vals_wrong_lag{hist_vals_wrong_lag_arg},
   symbol_table{symbol_table_arg},
   all_values_required{all_values_required_arg}
 {
@@ -322,10 +320,9 @@ HistValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
       set<int> unused_endo = symbol_table.getEndogenous();
       set<int> unused_exo = symbol_table.getExogenous();
 
-      set<int>::iterator sit;
       for (const auto & hist_value : hist_values)
         {
-          sit = unused_endo.find(hist_value.first.first);
+          auto sit = unused_endo.find(hist_value.first.first);
           if (sit != unused_endo.end())
             unused_endo.erase(sit);
 
@@ -353,7 +350,6 @@ HistValStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidat
       if (unused_endo.size() > 0 || unused_exo.size() > 0)
         exit(EXIT_FAILURE);
     }
-  mod_file_struct.hist_vals_wrong_lag = hist_vals_wrong_lag;
 }
 
 void
@@ -362,9 +358,13 @@ HistValStatement::writeOutput(ostream &output, const string &basename, bool mini
   output << "%" << endl
          << "% HISTVAL instructions" << endl
          << "%" << endl
-         << "M_.endo_histval = zeros(M_.endo_nbr,M_.maximum_lag);" << endl
-         << "M_.exo_histval = zeros(M_.exo_nbr,M_.maximum_lag);" << endl
-         << "M_.exo_det_histval = zeros(M_.exo_det_nbr,M_.maximum_lag);" << endl;
+         << "M_.histval_dseries = dseries(zeros(M_.orig_maximum_lag_with_diffs_expanded, M_.orig_endo_nbr"
+         << (symbol_table.exo_nbr() > 0 ? "+M_.exo_nbr" : "")
+         << (symbol_table.exo_det_nbr() > 0 ? "+M_.exo_det_nbr" : "")
+         << "), dates(sprintf('%dY', -M_.orig_maximum_lag_with_diffs_expanded+1)), [ M_.endo_names(1:M_.orig_endo_nbr); "
+         << (symbol_table.exo_nbr() > 0 ? "M_.exo_names; " : "")
+         << (symbol_table.exo_det_nbr() > 0 ? "M_.exo_det_names; " : "")
+         << "]);" << endl;
 
   for (const auto & hist_value : hist_values)
     {
@@ -372,42 +372,20 @@ HistValStatement::writeOutput(ostream &output, const string &basename, bool mini
       int lag = hist_value.first.second;
       const expr_t expression = hist_value.second;
 
-      SymbolType type = symbol_table.getType(symb_id);
-
-      // For a lag greater than 1 on endo, or for any exo, lookup for auxiliary variable
-      if ((type == SymbolType::endogenous && lag < 0) || type == SymbolType::exogenous)
-        {
-          try
-            {
-              // This function call must remain the 1st statement in this block
-              symb_id = symbol_table.searchAuxiliaryVars(symb_id, lag);
-              lag = 0;
-              type = SymbolType::endogenous;
-            }
-          catch (SymbolTable::SearchFailedException &e)
-            {
-              if (type == SymbolType::endogenous)
-                {
-                  cerr << "HISTVAL: internal error of Dynare, please contact the developers";
-                  exit(EXIT_FAILURE);
-                }
-              // We don't fail for exogenous, because they are not replaced by
-              // auxiliary variables in deterministic mode.
-            }
-        }
-
-      int tsid = symbol_table.getTypeSpecificID(symb_id) + 1;
-
-      if (type == SymbolType::endogenous)
-        output << "M_.endo_histval( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
-      else if (type == SymbolType::exogenous)
-        output << "M_.exo_histval( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
-      else if (type == SymbolType::exogenousDet)
-        output << "M_.exo_det_histval( " << tsid << ", M_.maximum_lag + " << lag << ") = ";
-
+      output << "M_.histval_dseries{'" << symbol_table.getName(symb_id) << "'}(dates('" << lag << "Y'))=";
       expression->writeOutput(output);
       output << ";" << endl;
     }
+
+  output << "if exist(['+' M_.fname '/dynamic_set_auxiliary_series'])" << endl
+       << "  eval(['M_.histval_dseries = ' M_.fname '.dynamic_set_auxiliary_series(M_.histval_dseries, []);']);" << endl
+       << "end" << endl
+       << "M_.endo_histval = M_.histval_dseries{M_.endo_names{:}}(dates(sprintf('%dY', 1-M_.maximum_lag)):dates('0Y')).data';" << endl;
+
+  if (symbol_table.exo_nbr() > 0)
+    output << "M_.exo_histval = M_.histval_dseries{M_.exo_names{:}}(dates(sprintf('%dY', 1-M_.maximum_lag)):dates('0Y')).data';" << endl;
+  if (symbol_table.exo_det_nbr() > 0)
+    output << "M_.exo_det_histval = M_.histval_dseries{M_.exo_det_names{:}}(dates(sprintf('%dY', 1-M_.maximum_lag)):dates('0Y')).data';" << endl;
 }
 
 void
diff --git a/src/NumericalInitialization.hh b/src/NumericalInitialization.hh
index b0c589ba..d341ec82 100644
--- a/src/NumericalInitialization.hh
+++ b/src/NumericalInitialization.hh
@@ -107,15 +107,12 @@ public:
     Maps pairs (symbol_id, lag) to expr_t
   */
   using hist_values_t = map<pair<int, int>, expr_t>;
-  using hist_vals_wrong_lag_t = map<int, int>;
 private:
   const hist_values_t hist_values;
-  const hist_vals_wrong_lag_t hist_vals_wrong_lag;
   const SymbolTable &symbol_table;
   const bool all_values_required;
 public:
   HistValStatement(hist_values_t hist_values_arg,
-                   const hist_vals_wrong_lag_t hist_vals_wrong_lag_arg,
                    const SymbolTable &symbol_table_arg,
                    const bool &all_values_required_arg);
   //! Workaround for trac ticket #157
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index 5714c81b..deeb52ee 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -661,9 +661,6 @@ ParsingDriver::hist_val(const string &name, const string &lag, expr_t rhs)
 
   pair<int, int> key(symb_id, ilag);
 
-  if (mod_file->dynamic_model.minLagForSymbol(symb_id) > ilag - 1)
-    hist_vals_wrong_lag[symb_id] = ilag;
-
   if (hist_values.find(key) != hist_values.end())
     error("hist_val: (" + name + ", " + lag + ") declared twice");
 
@@ -807,7 +804,7 @@ ParsingDriver::end_endval(bool all_values_required)
 void
 ParsingDriver::end_histval(bool all_values_required)
 {
-  mod_file->addStatement(make_unique<HistValStatement>(hist_values, hist_vals_wrong_lag, mod_file->symbol_table, all_values_required));
+  mod_file->addStatement(make_unique<HistValStatement>(hist_values, mod_file->symbol_table, all_values_required));
   hist_values.clear();
 }
 
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index c6895454..b172b413 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -173,8 +173,6 @@ private:
   InitOrEndValStatement::init_values_t init_values;
   //! Temporary storage for histval blocks
   HistValStatement::hist_values_t hist_values;
-  //! Temporary storage for histval blocks
-  HistValStatement::hist_vals_wrong_lag_t hist_vals_wrong_lag;
   //! Temporary storage for homotopy_setup blocks
   HomotopyStatement::homotopy_values_t homotopy_values;
   //! Temporary storage for moment_calibration
diff --git a/src/Statement.hh b/src/Statement.hh
index 40adb8a3..fb672604 100644
--- a/src/Statement.hh
+++ b/src/Statement.hh
@@ -123,8 +123,6 @@ public:
   bool steady_state_model_present{false};
   //! Whether there is a write_latex_steady_state_model statement present
   bool write_latex_steady_state_model_present{false};
-  //! Histval values that do not have the appropriate lag
-  map<int, int> hist_vals_wrong_lag;
   //! Pac growth and discount
   set<int> pac_params;
 };
-- 
GitLab