From c22d31330c7d3978f882003ffb999f65be70c4a5 Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@mjui.fr>
Date: Sun, 31 May 2015 12:18:06 +0200
Subject: [PATCH] Adding 'ramsey_constraints' block to declare contraints for
 Ramsey problems. Note that LMMCP (solve_algo=10) still doesn't work in all
 cases.

---
 doc/dynare.texi                               |  18 +++
 matlab/lmmcp/get_complementarity_conditions.m |  22 +++-
 .../private/simulation_core.m                 |   2 +-
 preprocessor/ComputingTasks.cc                | 103 ++++++++++++++++++
 preprocessor/ComputingTasks.hh                |  18 +++
 preprocessor/DynareBison.yy                   |  23 +++-
 preprocessor/DynareFlex.ll                    |   1 +
 preprocessor/ParsingDriver.cc                 |  51 +++++++++
 preprocessor/ParsingDriver.hh                 |  14 +++
 tests/Makefile.am                             |   1 +
 .../RamseyConstraints/test1.mod               |  38 +++++++
 11 files changed, 285 insertions(+), 6 deletions(-)
 create mode 100644 tests/optimal_policy/RamseyConstraints/test1.mod

diff --git a/doc/dynare.texi b/doc/dynare.texi
index 9ae87f06b..81216e4ae 100644
--- a/doc/dynare.texi
+++ b/doc/dynare.texi
@@ -6807,6 +6807,24 @@ instrument.
 
 @end deffn
 
+@deffn Block ramsey_constraints
+@anchor{ramsey_constraints}
+
+@descriptionhead
+
+This block lets you define constraints on the variables in the Ramsey
+problem. The constraints take the form of a variable, an inequality
+operator (@code{>} or @code{<}) and a constant.
+
+@examplehead
+
+@example
+ramsey_constraints;
+i > 0;
+end;
+@end example
+@end deffn
+ 
 @deffn Command ramsey_policy [@var{VARIABLE_NAME}@dots{}];
 @deffnx Command ramsey_policy (@var{OPTIONS}@dots{}) [@var{VARIABLE_NAME}@dots{}];
 @anchor{ramsey_policy}
diff --git a/matlab/lmmcp/get_complementarity_conditions.m b/matlab/lmmcp/get_complementarity_conditions.m
index 95ef094f9..cf5a40866 100644
--- a/matlab/lmmcp/get_complementarity_conditions.m
+++ b/matlab/lmmcp/get_complementarity_conditions.m
@@ -1,4 +1,4 @@
-function [lb,ub,eq_index] = get_complementarity_conditions(M)
+function [lb,ub,eq_index] = get_complementarity_conditions(M,ramsey_policy)
 
 % Copyright (C) 2014 Dynare Team
 %
@@ -17,11 +17,27 @@ function [lb,ub,eq_index] = get_complementarity_conditions(M)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-
-etags = M.equations_tags;
 ub = inf(M.endo_nbr,1);
 lb = -ub;
 eq_index = (1:M.endo_nbr)';
+if ramsey_policy
+    if isfield(M,'ramsey_model_constraints')
+        rc = M.ramsey_model_constraints;
+        for i = 1:length(rc)
+            switch rc{i}{2}
+              case {'>','>='}
+                lb(rc{i}{1}) = eval(rc{i}{3});
+              case {'<','<='}
+                ub(rc{i}{1}) = eval(rc{i}{3});
+              otherwise
+                error('Wrong operator in get_complementarity_conditions')
+            end
+            eq_index(i) = 1;
+        end
+    end
+end
+
+etags = M.equations_tags;
 for i=1:size(etags,1)
     if strcmp(etags{i,2},'mcp')
         str = etags{i,3};
diff --git a/matlab/perfect-foresight-models/private/simulation_core.m b/matlab/perfect-foresight-models/private/simulation_core.m
index cb5f6d540..fe22ce612 100644
--- a/matlab/perfect-foresight-models/private/simulation_core.m
+++ b/matlab/perfect-foresight-models/private/simulation_core.m
@@ -66,7 +66,7 @@ else
             elseif options_.stack_solve_algo == 7
                 periods = options_.periods;
                 if ~isfield(options_.lmmcp,'lb')
-                    [lb,ub,pfm.eq_index] = get_complementarity_conditions(M_);
+                    [lb,ub,pfm.eq_index] = get_complementarity_conditions(M_,options_.ramsey_policy);
                     options_.lmmcp.lb = repmat(lb,periods,1);
                     options_.lmmcp.ub = repmat(ub,periods,1);
                 end
diff --git a/preprocessor/ComputingTasks.cc b/preprocessor/ComputingTasks.cc
index 6057fad53..4e95890d2 100644
--- a/preprocessor/ComputingTasks.cc
+++ b/preprocessor/ComputingTasks.cc
@@ -397,6 +397,109 @@ RamseyModelStatement::writeOutput(ostream &output, const string &basename, bool
   options_list.writeOutput(output);
 }
 
+RamseyConstraintsStatement::RamseyConstraintsStatement(const constraints_t &constraints_arg) :
+  constraints(constraints_arg)
+{
+}
+
+void
+RamseyConstraintsStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
+{
+  if ((mod_file_struct.ramsey_model_present != true) || (  mod_file_struct.ramsey_policy_present != true))
+    cerr << "ramsey_constraints: can only be used with ramsey_model or ramsey_policy" << endl;
+}
+
+void
+RamseyConstraintsStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
+{
+  output << "M_.ramsey_model_constraints = {" << endl;
+  for (RamseyConstraintsStatement::constraints_t::const_iterator it = constraints.begin(); it != constraints.end(); ++it)
+    {
+      if (it != constraints.begin())
+	output << ", ";
+      output << "{" << it->endo + 1 << ", '";
+      switch(it->code)
+	{
+	case oLess:
+	  output << '<';
+	  break;
+	case oGreater:
+	  output << '>';
+	  break;
+	case oLessEqual:
+	  output << "<=";
+	  break;
+	case oGreaterEqual:
+	  output << ">=";
+	  break;
+	default:
+	  cerr << "Ramsey constraints: this shouldn't happen." << endl;
+	  exit(1);
+	}
+      output << "', '";
+      it->expression->writeOutput(output);
+      output << "'}" << endl;
+    }
+  output << "};" << endl;
+}
+
+// Statement *
+// RamseyConstraintsStatement::cloneAndReindexSymbIds(DataTree &dynamic_datatree, SymbolTable &orig_symbol_table)
+// {
+//   vector<string> errors;
+//   SymbolList new_symbol_list, new_options_symbol_list;
+//   OptionsList new_options_list = options_list;
+//   SymbolTable *new_symbol_table =  dynamic_datatree.getSymbolTable();
+//   vector<string> symbols = symbol_list.get_symbols();
+
+//   for (vector<string>::const_iterator it = symbols.begin(); it != symbols.end(); it++)
+//     try
+//       {
+//         new_symbol_table->getID(*it);
+//         new_symbol_list.addSymbol(*it);
+//       }
+//     catch (SymbolTable::UnknownSymbolIDException &e)
+//       {
+//         errors.push_back(orig_symbol_table.getName(e.id));
+//       }
+//     catch (SymbolTable::UnknownSymbolNameException &e)
+//       {
+//         errors.push_back(e.name);
+//       }
+
+//   OptionsList::symbol_list_options_t::const_iterator it = options_list.symbol_list_options.find("instruments");
+//   if (it != options_list.symbol_list_options.end())
+//     {
+//       symbols = it->second.get_symbols();
+//       for (vector<string>::const_iterator it1 = symbols.begin(); it1 != symbols.end(); it1++)
+//         try
+//           {
+//             new_symbol_table->getID(*it1);
+//             new_options_symbol_list.addSymbol(*it1);
+//           }
+//         catch (SymbolTable::UnknownSymbolIDException &e)
+//           {
+//             errors.push_back(orig_symbol_table.getName(e.id));
+//           }
+//         catch (SymbolTable::UnknownSymbolNameException &e)
+//           {
+//             errors.push_back(e.name);
+//           }
+//       new_options_list.symbol_list_options["instruments"] = new_options_symbol_list;
+//     }
+
+//   if (!errors.empty())
+//     {
+//       cerr << endl
+//            << "ERROR: The following vars were used in the ramsey_policy statement(s) but  were not declared." << endl
+//            << "       This likely means that you declared them as varexo and that they're not in the model" << endl;
+//       for (vector<string>::const_iterator it = errors.begin(); it != errors.end(); it++)
+//         cerr << *it << endl;
+//       exit(EXIT_FAILURE);
+//     }
+//   return new RamseyPolicyStatement(new_symbol_list, options_list);
+// }
+
 RamseyPolicyStatement::RamseyPolicyStatement(const SymbolList &symbol_list_arg,
                                              const OptionsList &options_list_arg) :
   symbol_list(symbol_list_arg),
diff --git a/preprocessor/ComputingTasks.hh b/preprocessor/ComputingTasks.hh
index 572d8016d..b9b65d465 100644
--- a/preprocessor/ComputingTasks.hh
+++ b/preprocessor/ComputingTasks.hh
@@ -125,6 +125,24 @@ public:
   virtual Statement *cloneAndReindexSymbIds(DataTree &dynamic_datatree, SymbolTable &orig_symbol_table);
 };
 
+class RamseyConstraintsStatement : public Statement
+{
+public:
+  struct Constraint {
+    int endo;
+    BinaryOpcode code;
+    expr_t expression;
+  }; 
+  typedef vector<Constraint> constraints_t;
+private:
+  const constraints_t constraints;
+public:
+  RamseyConstraintsStatement(const constraints_t &constraints_arg);
+  virtual void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings);
+  virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
+  //  virtual Statement *cloneAndReindexSymbIds(DataTree &dynamic_datatree, SymbolTable &orig_symbol_table);
+};
+
 class RamseyPolicyStatement : public Statement
 {
 private:
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 5893b6dc1..d886fb0f6 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -122,7 +122,7 @@ class ParsingDriver;
 %token SHOCKS SHOCK_DECOMPOSITION SIGMA_E SIMUL SIMUL_ALGO SIMUL_SEED ENDOGENOUS_TERMINAL_PERIOD
 %token SMOOTHER SMOOTHER2HISTVAL SQUARE_ROOT_SOLVER STACK_SOLVE_ALGO STEADY_STATE_MODEL SOLVE_ALGO SOLVER_PERIODS
 %token STDERR STEADY STOCH_SIMUL SURPRISE SYLVESTER SYLVESTER_FIXED_POINT_TOL REGIMES REGIME
-%token TEX RAMSEY_MODEL RAMSEY_POLICY PLANNER_DISCOUNT DISCRETIONARY_POLICY DISCRETIONARY_TOL
+%token TEX RAMSEY_MODEL RAMSEY_POLICY RAMSEY_CONSTRAINTS PLANNER_DISCOUNT DISCRETIONARY_POLICY DISCRETIONARY_TOL
 %token <string_val> TEX_NAME
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED
 %token VALUES VAR VAREXO VAREXO_DET VAROBS PREDETERMINED_VARIABLES
@@ -239,7 +239,8 @@ statement : parameters
           | planner_objective
           | ramsey_model
           | ramsey_policy
-          | discretionary_policy
+	  | ramsey_constraints
+	  | discretionary_policy
           | bvar_density
           | bvar_forecast
           | sbvar
@@ -1903,6 +1904,24 @@ ramsey_policy : RAMSEY_POLICY ';'
                 { driver.ramsey_policy(); }
               ;
 
+ramsey_constraints : RAMSEY_CONSTRAINTS ';' ramsey_constraints_list END ';'
+                     { driver.add_ramsey_constraints_statement(); }
+		   ;
+
+ramsey_constraints_list : ramsey_constraints_list ramsey_constraint 
+                 | ramsey_constraint
+		 ;
+
+ramsey_constraint : NAME  LESS expression ';'
+                    { driver.ramsey_constraint_add_less($1,$3); }
+		  | NAME  GREATER  expression ';'
+                    { driver.ramsey_constraint_add_greater($1,$3); }
+		  | NAME  LESS_EQUAL expression ';'		
+                    { driver.ramsey_constraint_add_less_equal($1,$3); }
+		  | NAME  GREATER  expression ';'
+                    { driver.ramsey_constraint_add_greater_equal($1,$3); }
+		  ;
+
 discretionary_policy : DISCRETIONARY_POLICY ';'
                        { driver.discretionary_policy(); }
                      | DISCRETIONARY_POLICY '(' discretionary_policy_options_list ')' ';'
diff --git a/preprocessor/DynareFlex.ll b/preprocessor/DynareFlex.ll
index 92d52ecfe..5dcd575d2 100644
--- a/preprocessor/DynareFlex.ll
+++ b/preprocessor/DynareFlex.ll
@@ -199,6 +199,7 @@ DATE -?[0-9]+([YyAa]|[Mm]([1-9]|1[0-2])|[Qq][1-4]|[Ww]([1-9]{1}|[1-4][0-9]|5[0-2
 <INITIAL>svar_identification {BEGIN DYNARE_BLOCK; return token::SVAR_IDENTIFICATION;}
 <INITIAL>moment_calibration {BEGIN DYNARE_BLOCK; return token::MOMENT_CALIBRATION;}
 <INITIAL>irf_calibration {BEGIN DYNARE_BLOCK; return token::IRF_CALIBRATION;}
+<INITIAL>ramsey_constraints {BEGIN DYNARE_BLOCK; return token::RAMSEY_CONSTRAINTS;}
 
  /* For the semicolon after an "end" keyword */
 <INITIAL>; {return Dynare::parser::token_type (yytext[0]);}
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index ba541ec81..94679f0aa 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -2772,3 +2772,54 @@ ParsingDriver::perfect_foresight_solver()
   mod_file->addStatement(new PerfectForesightSolverStatement(options_list));
   options_list.clear();
 }
+
+void
+ParsingDriver::add_ramsey_constraints_statement()
+{
+  mod_file->addStatement(new RamseyConstraintsStatement(ramsey_constraints));
+  ramsey_constraints.clear();
+}
+
+void
+ParsingDriver::ramsey_constraint_add_less(const string *name, const expr_t rhs)
+{
+  add_ramsey_constraint(name,oLess,rhs);
+}
+
+void
+ParsingDriver::ramsey_constraint_add_greater(const string *name, const expr_t rhs)
+{
+  add_ramsey_constraint(name,oGreater,rhs);
+}
+
+void
+ParsingDriver::ramsey_constraint_add_less_equal(const string *name, const expr_t rhs)
+{
+  add_ramsey_constraint(name,oLessEqual,rhs);
+}
+
+void
+ParsingDriver::ramsey_constraint_add_greater_equal(const string *name, const expr_t rhs)
+{
+  add_ramsey_constraint(name,oGreaterEqual,rhs);
+}
+
+void
+ParsingDriver::add_ramsey_constraint(const string *name, BinaryOpcode op_code, const expr_t rhs)
+{
+  check_symbol_existence(*name);
+  int symb_id = mod_file->symbol_table.getID(*name);
+  SymbolType type = mod_file->symbol_table.getType(symb_id);
+
+  if (type != eEndogenous)
+    error("ramsey_constraints: " + *name + " should be an endogenous variable");
+
+  RamseyConstraintsStatement::Constraint C;
+  C.endo = symb_id;
+  C.code = op_code;
+  C.expression = rhs;
+  ramsey_constraints.push_back(C);
+
+  delete name;
+}
+
diff --git a/preprocessor/ParsingDriver.hh b/preprocessor/ParsingDriver.hh
index dd93276c0..d7da536cd 100644
--- a/preprocessor/ParsingDriver.hh
+++ b/preprocessor/ParsingDriver.hh
@@ -157,6 +157,8 @@ private:
   MomentCalibration::constraints_t moment_calibration_constraints;
   //! Temporary storage for irf_calibration
   IrfCalibration::constraints_t irf_calibration_constraints;
+  //! Temporary storage for ramsey_constraints
+  RamseyConstraintsStatement::constraints_t ramsey_constraints;
   //! Temporary storage for svar_identification blocks
   SvarIdentificationStatement::svar_identification_restrictions_t svar_ident_restrictions;
   //! Temporary storage for mapping the equation number to the restrictions within an svar_identification block
@@ -506,6 +508,18 @@ public:
   void end_planner_objective(expr_t expr);
   //! Ramsey model statement
   void ramsey_model();
+  //! Ramsey constraints statement
+  void add_ramsey_constraints_statement();
+  //! Ramsey less constraint
+  void ramsey_constraint_add_less(const string *name, const expr_t rhs);
+  //! Ramsey greater constraint
+  void ramsey_constraint_add_greater(const string *name, const expr_t rhs);
+  //! Ramsey less or equal constraint
+  void ramsey_constraint_add_less_equal(const string *name, const expr_t rhs);
+  //! Ramsey greater or equal constraint
+  void ramsey_constraint_add_greater_equal(const string *name, const expr_t rhs);
+  //! Ramsey constraint helper function
+  void add_ramsey_constraint(const string *name, BinaryOpcode op_code, const expr_t rhs);
   //! Ramsey policy statement
   void ramsey_policy();
   //! Discretionary policy statement
diff --git a/tests/Makefile.am b/tests/Makefile.am
index c45755cd8..02f476feb 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -41,6 +41,7 @@ MODFILES = \
 	optimal_policy/Ramsey/ramsey_ex.mod \
 	optimal_policy/Ramsey/ramsey_ex_initval_AR2.mod \
 	optimal_policy/Ramsey/ramsey_ex_aux.mod \
+	optimal_policy/RamseyConstraints/test1.mod \
 	discretionary_policy/dennis_1.mod \
 	initval_file/ramst_initval_file.mod \
 	ramst_normcdf_and_friends.mod \
diff --git a/tests/optimal_policy/RamseyConstraints/test1.mod b/tests/optimal_policy/RamseyConstraints/test1.mod
new file mode 100644
index 000000000..2bb83b6d3
--- /dev/null
+++ b/tests/optimal_policy/RamseyConstraints/test1.mod
@@ -0,0 +1,38 @@
+var i y pi;
+varexo e_y e_pi;
+
+parameters beta1 beta2 beta3 lambda1 lambda2 pi_bar;
+beta1 = 0.6;
+beta2 = 0.25;
+beta3 = -0.2;
+lambda1 = 0.7;
+lambda2 = 0.1;
+pi_bar = 2.0;
+
+model;
+y = beta1*y(-1) + beta2*y(+1) + beta3*(i-pi(+1)) + e_y;
+pi = lambda1*pi(+1) + (1-lambda1)*pi(-1) + lambda2*y + e_pi;
+end;
+
+planner_objective (pi-pi_bar)^2 + y^2;
+
+ramsey_model(planner_discount=1.0);
+
+histval;
+y(0) = -2.0;
+pi(0) = 1.0;
+end;
+
+steady;
+
+ramsey_constraints;
+i > 0;
+end;
+
+perfect_foresight_setup(periods=50);
+options_.stack_solve_algo = 7;
+options_.solve_algo = 10;
+perfect_foresight_solver;
+
+rplot i;
+      
\ No newline at end of file
-- 
GitLab