From a6add6127bb51413c1e9ce6b866a6b78fb510614 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Villemot?= <sebastien@dynare.org>
Date: Fri, 16 Jul 2021 15:12:37 +0200
Subject: [PATCH] VAR models: output vector of constants in varmatrices.m

---
 src/DynamicModel.cc | 18 ++++++++++++++++++
 src/SubModel.cc     | 17 ++++++++++++++++-
 src/SubModel.hh     |  2 ++
 3 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 45853154..c5b2724d 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3537,6 +3537,7 @@ DynamicModel::fillVarModelTableMatrices()
 {
   map<string, map<tuple<int, int, int>, expr_t>> AR;
   map<string, map<tuple<int, int>, expr_t>> A0;
+  map<string, map<int, expr_t>> constants;
   for (const auto &[model_name, eqns] : var_model_table.getEqNums())
     {
       const vector<int> &lhs = var_model_table.getLhs(model_name);
@@ -3581,10 +3582,27 @@ DynamicModel::fillVarModelTableMatrices()
                   A0[model_name][{ i, lhs_symb_id }] = d;
                 }
             }
+
+          // Fill constants vector
+          // Constants are computed by replacing all (transformed) endos and exos by zero
+          constants[model_name] = {}; // Ensure that the map exists, even if constants are all zero
+          for (size_t i = 0; i < eqns.size(); i++)
+            {
+              auto rhs = equations[eqns[i]]->arg2;
+              map<VariableNode *, NumConstNode *> subst_table;
+              auto rhs_vars = var_model_table.getRhs(model_name)[i]; // All the (transformed) endogenous on RHS, as computed by updateVarAndTrendModel()
+              rhs->collectDynamicVariables(SymbolType::exogenous, rhs_vars); // Add exos
+              for (auto [symb_id, lag] : rhs_vars)
+                subst_table[AddVariable(symb_id, lag)] = Zero;
+              expr_t c = rhs->replaceVarsInEquation(subst_table);
+              if (c != Zero)
+                constants[model_name][i] = c;
+            }
         }
     }
   var_model_table.setAR(AR);
   var_model_table.setA0(A0);
+  var_model_table.setConstants(constants);
 }
 
 map<string, map<tuple<int, int, int>, expr_t>>
diff --git a/src/SubModel.cc b/src/SubModel.cc
index d0a827a0..c8c63e75 100644
--- a/src/SubModel.cc
+++ b/src/SubModel.cc
@@ -455,7 +455,7 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
       cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(EXIT_FAILURE);
     }
-  ar_output << "function [ar, a0] = varmatrices(model_name, params, reducedform)" << endl
+  ar_output << "function [ar, a0, constants] = varmatrices(model_name, params, reducedform)" << endl
             << "% File automatically generated by the Dynare preprocessor" << endl << endl
             << "if nargin<3" << endl
             << "    reducedform = false;" << endl
@@ -535,6 +535,15 @@ VarModelTable::writeOutput(const string &basename, ostream &output) const
                 << "            end" << endl
                 << "            a0 = eye(" << lhs.size() << ");" << endl
                 << "        end" << endl
+                << "        if nargout>2" << endl
+                << "            constants = zeros(" << lhs.size() << ");" << endl;
+      for (auto [eqn, expr] : constants.at(name))
+        {
+          ar_output << "            constants(" << eqn + 1 << ") = ";
+          expr->writeOutput(ar_output, ExprNodeOutputType::matlabDynamicModel);
+          ar_output << ";" << endl;
+        }
+      ar_output << "        end" << endl
                 << "    end" << endl
                 << "    return" << endl
                 << "end" << endl << endl;
@@ -687,6 +696,12 @@ VarModelTable::setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg)
   A0 = move(A0_arg);
 }
 
+void
+VarModelTable::setConstants(map<string, map<int, expr_t>> constants_arg)
+{
+  constants = move(constants_arg);
+}
+
 const vector<int> &
 VarModelTable::getMaxLags(const string &name_arg) const
 {
diff --git a/src/SubModel.hh b/src/SubModel.hh
index 37c35da7..06af8284 100644
--- a/src/SubModel.hh
+++ b/src/SubModel.hh
@@ -130,6 +130,7 @@ private:
      VAR context is not the same thing as in the trend-component model
      context. */
   map<string, map<tuple<int, int>, expr_t>> A0; // name -> (eqn, lhs_symb_id) -> param_expr_t
+  map<string, map<int, expr_t>> constants; // name -> eqn -> constant
 public:
   explicit VarModelTable(SymbolTable &symbol_table_arg);
 
@@ -161,6 +162,7 @@ public:
   void setOrigDiffVar(map<string, vector<int>> orig_diff_var_arg);
   void setAR(map<string, map<tuple<int, int, int>, expr_t>> AR_arg);
   void setA0(map<string, map<tuple<int, int>, expr_t>> A0_arg);
+  void setConstants(map<string, map<int, expr_t>> constants_arg);
 
   //! Write output of this class
   void writeOutput(const string &basename, ostream &output) const;