Implement explicit writing of first order conditions of Ramsey problem (ticket #5)

preprocessor/ComputingTasks.cc

@@ -216,6 +216,15 @@ RamseyPolicyStatement::writeOutput(ostream &output, const string &basename) cons
   output << "ramsey_policy(var_list_);\n";
 }
 
+string
+RamseyPolicyStatement::getPlannerDiscount() const
+{
+  OptionsList::num_options_t::const_iterator it = options_list.num_options.find("planner_discount");
+  if (it != options_list.num_options.end())
+    return it->second;
+  return "1.0";
+}
+
 DiscretionaryPolicyStatement::DiscretionaryPolicyStatement(const SymbolList &symbol_list_arg,
 							   const OptionsList &options_list_arg) :
   symbol_list(symbol_list_arg),
@@ -847,6 +856,12 @@ PlannerObjectiveStatement::checkPass(ModFileStructure &mod_file_struct)
   mod_file_struct.planner_objective_present = true;
 }
 
+StaticModel *
+PlannerObjectiveStatement::getPlannerObjective() const
+{
+  return model_tree;
+}
+
 void
 PlannerObjectiveStatement::computingPass()
 {

preprocessor/ComputingTasks.hh

@@ -101,6 +101,7 @@ public:
                         const OptionsList &options_list_arg);
   virtual void checkPass(ModFileStructure &mod_file_struct);
   virtual void writeOutput(ostream &output, const string &basename) const;
+  virtual string getPlannerDiscount() const;
 };
 
 class DiscretionaryPolicyStatement : public Statement
@@ -342,6 +343,8 @@ public:
   /*! \todo allow for the possibility of disabling temporary terms */
   virtual void computingPass();
   virtual void writeOutput(ostream &output, const string &basename) const;
+  //! Return the Planner Objective
+  StaticModel *getPlannerObjective() const;
 };
 
 class BVARDensityStatement : public Statement

preprocessor/DynamicModel.cc

@@ -18,6 +18,7 @@
  */
 
 #include <iostream>
+#include <sstream>
 #include <cmath>
 #include <cstdlib>
 #include <cassert>
@@ -3006,6 +3007,89 @@ DynamicModel::cloneDynamic(DynamicModel &dynamic_model) const
     dynamic_model.addAuxEquation((*it)->cloneDynamic(dynamic_model));
 }
 
+void
+DynamicModel::replaceMyEquations(DynamicModel &dynamic_model) const
+{
+  dynamic_model.equations.clear();
+  for (vector<BinaryOpNode *>::const_iterator it = equations.begin();
+       it != equations.end(); it++)
+    dynamic_model.addEquation((*it)->cloneDynamic(dynamic_model));
+}
+
+void
+DynamicModel::computeRamseyPolicyFOCs(const StaticModel &static_model, const string &discount_factor)
+{
+  // Add aux LM to constraints in equations
+  // equation[i]->lhs = rhs becomes equation[i]->AUX_LAMBDA_i*(lhs-rhs) = 0
+  int i;
+  for (i = 0; i < (int) equations.size(); i++)
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->addMultipliersToConstraints(i));
+      assert(substeq != NULL);
+      equations[i] = substeq;
+    }
+
+  cout << "Ramsey Problem: added " << i << " Multipliers." << endl;
+
+  // Add Planner Objective to equations to include in computeDerivIDs
+  assert(static_model.equations.size() == 1);
+  addEquation(static_model.equations[0]->cloneDynamic(*this));
+
+  // Get max endo lead and max endo lag
+  set<pair<int, int> > dynvars;
+  int max_eq_lead = 0;
+  int max_eq_lag = 0;
+  for (int i = 0; i < (int) equations.size(); i++)
+    equations[i]->collectVariables(eEndogenous, dynvars);
+
+  for (set<pair<int, int> >::const_iterator it = dynvars.begin();
+       it != dynvars.end(); it++)
+    {
+      int lag = it->second;
+      if (max_eq_lead < lag)
+        max_eq_lead = lag;
+      else if (-max_eq_lag > lag)
+        max_eq_lag = -lag;
+    }
+
+  // Create (modified) Lagrangian (so that we can take the derivative once at time t)
+  expr_t lagrangian = Zero;
+  expr_t discount_factor_node = AddNonNegativeConstant(discount_factor);
+  for (i = 0; i < (int) equations.size(); i++)
+    for (int lag = -max_eq_lag; lag <= max_eq_lead; lag++)
+      {
+        expr_t dfpower = NULL;
+        std::stringstream lagstream;
+        lagstream << abs(lag);
+        if (lag < 0)
+          dfpower = AddNonNegativeConstant(lagstream.str());
+        else if (lag == 0)
+          dfpower = Zero;
+        else
+          dfpower = AddMinus(Zero, AddNonNegativeConstant(lagstream.str()));
+
+        lagrangian = AddPlus(AddTimes(AddPower(discount_factor_node, dfpower),
+                                      equations[i]->getNonZeroPartofEquation()->decreaseLeadsLags(lag)), lagrangian);
+      }
+
+  equations.clear();
+  addEquation(AddEqual(lagrangian, Zero));
+  computeDerivIDs();
+
+  //Compute derivatives and overwrite equations
+  vector<expr_t> neweqs;
+  for (deriv_id_table_t::const_iterator it = deriv_id_table.begin();
+       it != deriv_id_table.end(); it++)
+    // For all endogenous variables with zero lag
+    if (symbol_table.getType(it->first.first)  == eEndogenous && it->first.second == 0)
+      neweqs.push_back(AddEqual(equations[0]->getNonZeroPartofEquation()->getDerivative(it->second), Zero));
+
+  // Add new equations
+  equations.clear();
+  for (int i = 0; i < (int) neweqs.size(); i++)
+    addEquation(neweqs[i]);
+}
+
 void
 DynamicModel::toStatic(StaticModel &static_model) const
 {
@@ -3634,6 +3718,9 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         case avExpectationRIS:
           cout << "expectation conditional on a restricted information set";
           break;
+        case avMultiplier:
+          cerr << "avMultiplier encountered: impossible case" << endl;
+          exit(EXIT_FAILURE);
         }
       cout << ": added " << neweqs.size() << " auxiliary variables and equations." << endl;
     }

preprocessor/DynamicModel.hh

@@ -284,6 +284,10 @@ public:
   /*! It assumes that the dynamic model given in argument has just been allocated */
   void cloneDynamic(DynamicModel &dynamic_model) const;
 
+  //! Replaces model equations with derivatives of Lagrangian w.r.t. endogenous
+  void computeRamseyPolicyFOCs(const StaticModel &static_model, const string &discount_factor);
+  void replaceMyEquations(DynamicModel &dynamic_model) const;
+
   //! Writes LaTeX file with the equations of the dynamic model
   void writeLatexFile(const string &basename) const;
 

preprocessor/ExprNode.cc

@@ -2238,6 +2238,15 @@ BinaryOpNode::prepareForDerivation()
             inserter(non_null_derivatives, non_null_derivatives.begin()));
 }
 
+expr_t
+BinaryOpNode::getNonZeroPartofEquation() const
+{
+  assert(arg1 == datatree.Zero || arg2 == datatree.Zero);
+  if (arg1 == datatree.Zero)
+    return arg2;
+  return arg1;
+}
+
 expr_t
 BinaryOpNode::composeDerivatives(expr_t darg1, expr_t darg2)
 {
@@ -3342,6 +3351,14 @@ BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+expr_t
+BinaryOpNode::addMultipliersToConstraints(int i)
+{
+  int symb_id = datatree.symbol_table.addMultiplierAuxiliaryVar(i);
+  expr_t newAuxLM = datatree.AddVariable(symb_id, 0);
+  return datatree.AddEqual(datatree.AddTimes(newAuxLM, datatree.AddMinus(arg1, arg2)), datatree.Zero);
+}
+
 bool
 BinaryOpNode::isNumConstNodeEqualTo(double value) const
 {

preprocessor/ExprNode.hh

@@ -657,6 +657,10 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   //! Function to write out the oPowerNode in expr_t terms as opposed to writing out the function itself
   expr_t unpackPowerDeriv() const;
+  //! Returns AUX_LAMBDA*(lhs-rhs) = 0, creating lagrange multiplier AUX_LAMBDA
+  expr_t addMultipliersToConstraints(int i);
+  //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
+  expr_t getNonZeroPartofEquation() const;
 };
 
 //! Trinary operator node

preprocessor/ModFile.cc

@@ -21,16 +21,20 @@
 #include <iostream>
 #include <fstream>
 #include <typeinfo>
+#include <cassert>
 #include "ModFile.hh"
 #include "ConfigFile.hh"
+#include "ComputingTasks.hh"
 
 ModFile::ModFile() : expressions_tree(symbol_table, num_constants, external_functions_table),
                      dynamic_model(symbol_table, num_constants, external_functions_table),
                      trend_dynamic_model(symbol_table, num_constants, external_functions_table),
+                     ramsey_FOC_equations_dynamic_model(symbol_table, num_constants, external_functions_table),
                      static_model(symbol_table, num_constants, external_functions_table),
                      steady_state_model(symbol_table, num_constants, external_functions_table, static_model),
                      linear(false), block(false), byte_code(false),
-                     use_dll(false), no_static(false), nonstationary_variables(false)
+                     use_dll(false), no_static(false), nonstationary_variables(false),
+                     ramsey_policy_orig_eqn_nbr(0)
 {
 }
 
@@ -236,6 +240,32 @@ ModFile::transformPass()
       dynamic_model.removeTrendVariableFromEquations();
     }
 
+  if (mod_file_struct.ramsey_policy_present)
+    {
+      StaticModel *planner_objective = NULL;
+      string planner_discount = "";
+      for (vector<Statement *>::iterator it = statements.begin(); it != statements.end(); it++)
+        {
+          PlannerObjectiveStatement *pos = dynamic_cast<PlannerObjectiveStatement *>(*it);
+          if (pos != NULL)
+            planner_objective = pos->getPlannerObjective();
+
+          RamseyPolicyStatement *rps = dynamic_cast<RamseyPolicyStatement *>(*it);
+          if (rps != NULL)
+            planner_discount = rps->getPlannerDiscount();
+        }
+      assert(planner_objective != NULL && !planner_discount.empty());
+      ramsey_policy_orig_eqn_nbr = dynamic_model.equation_number();
+
+      /*
+        clone the model then clone the new equations back to the original because
+        we have to call computeDerivIDs (in computeRamseyPolicyFOCs and computingPass)
+       */
+      dynamic_model.cloneDynamic(ramsey_FOC_equations_dynamic_model);
+      ramsey_FOC_equations_dynamic_model.computeRamseyPolicyFOCs(*planner_objective, planner_discount);
+      ramsey_FOC_equations_dynamic_model.replaceMyEquations(dynamic_model);
+    }
+
   if (mod_file_struct.stoch_simul_present
       || mod_file_struct.estimation_present
       || mod_file_struct.osr_present
@@ -286,14 +316,20 @@ ModFile::transformPass()
       cerr << "ERROR: There are " << dynamic_model.equation_number() << " equations but " << symbol_table.endo_nbr() << " endogenous variables!" << endl;
       exit(EXIT_FAILURE);
     }
-  
+
   if (symbol_table.exo_det_nbr() > 0 && mod_file_struct.simul_present)
     {
       cerr << "ERROR: A .mod file cannot contain both a simul command and varexo_det declaration (all exogenous variables are deterministic in this case)" << endl;
       exit(EXIT_FAILURE);
     }
-  
-  cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
+
+  if (!mod_file_struct.ramsey_policy_present)
+    cout << "Found " << dynamic_model.equation_number() << " equation(s)." << endl;
+  else
+    {
+      cout << "Found " << ramsey_policy_orig_eqn_nbr  << " equation(s)." << endl;
+      cout << "Found " << dynamic_model.equation_number() << " FOC equation(s) for Ramsey Problem." << endl;
+    }
 
   if (symbol_table.exists("dsge_prior_weight"))
     if (mod_file_struct.bayesian_irf_present)
@@ -518,14 +554,27 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
   if (block && !byte_code)
     mOutputFile << "addpath " << basename << ";" << endl;
 
+  if (mod_file_struct.ramsey_policy_present)
+    mOutputFile << "M_.orig_eq_nbr = " << ramsey_policy_orig_eqn_nbr << ";" << endl;
+
   if (dynamic_model.equation_number() > 0)
     {
       if (dynamic_model_needed)
+        /* to be removed before ramsey_policy commit
+        if (mod_file_struct.ramsey_policy_present)
+          ramsey_policy_FOC_dynamic_model.writeOutput(mOutputFile, basename, block, byte_code, use_dll, mod_file_struct.order_option);
+        else
+        */
         dynamic_model.writeOutput(mOutputFile, basename, block, byte_code, use_dll, mod_file_struct.order_option);
       else
         dynamic_model.writeOutput(mOutputFile, basename, false, false, false, mod_file_struct.order_option);
       if (!no_static)
-        static_model.writeOutput(mOutputFile, block);
+        /* to be removed before ramsey_policy commit
+        if (mod_file_struct.ramsey_policy_present)
+          ramsey_policy_FOC_static_model.writeOutput(mOutputFile, block);
+        else
+        */
+          static_model.writeOutput(mOutputFile, block);
     }
 
   // Print statements
@@ -565,13 +614,23 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool console,
   if (dynamic_model.equation_number() > 0)
     {
       if (!no_static)
-        static_model.writeStaticFile(basename, block, byte_code);
+        /* to be removed before ramsey_policy commit
+        if (mod_file_struct.ramsey_policy_present)
+          ramsey_policy_FOC_static_model.writeStaticFile(basename, block, byte_code);
+        else
+        */
+          static_model.writeStaticFile(basename, block, byte_code);
 
       if (dynamic_model_needed)
-        {
-          dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option);
-          dynamic_model.writeParamsDerivativesFile(basename);
-        }
+        /* to be removed before ramsey_policy commit
+        if (mod_file_struct.ramsey_policy_present)
+          ramsey_policy_FOC_dynamic_model.writeDynamicFile(basename, false, false, false, mod_file_struct.order_option);
+        else
+        */
+          {
+            dynamic_model.writeDynamicFile(basename, block, byte_code, use_dll, mod_file_struct.order_option);
+            dynamic_model.writeParamsDerivativesFile(basename);
+          }
       else
         {
           dynamic_model.writeDynamicFile(basename, false, false, false, mod_file_struct.order_option);

preprocessor/ModFile.hh

@@ -53,6 +53,8 @@ public:
   DynamicModel dynamic_model;
   //! A copy of Dynamic model, for testing trends declared by user
   DynamicModel trend_dynamic_model;
+  //! A model in which to create the FOC for the ramsey problem
+  DynamicModel ramsey_FOC_equations_dynamic_model;
   //! Static model, as derived from the "model" block when leads and lags have been removed
   StaticModel static_model;
   //! Static model, as declared in the "steady_state_model" block if present
@@ -79,6 +81,9 @@ public:
   /*! Filled using initval blocks and parameters initializations */
   eval_context_t global_eval_context;
 
+  //! Stores the original number of equations in the model_block
+  int ramsey_policy_orig_eqn_nbr;
+
 private:
   //! List of statements
   vector<Statement *> statements;

preprocessor/SymbolTable.cc

@@ -24,12 +24,14 @@
 
 #include "SymbolTable.hh"
 
-AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg, int orig_lead_lag_arg, string expectation_information_set_name_arg) :
+AuxVarInfo::AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id_arg, int orig_lead_lag_arg,
+                       string expectation_information_set_name_arg, int equation_number_for_multiplier_arg) :
   symb_id(symb_id_arg),
   type(type_arg),
   orig_symb_id(orig_symb_id_arg),
   orig_lead_lag(orig_lead_lag_arg),
-  expectation_information_set_name(expectation_information_set_name_arg)
+  expectation_information_set_name(expectation_information_set_name_arg),
+  equation_number_for_multiplier(equation_number_for_multiplier_arg)
 {
 }
 
@@ -235,6 +237,9 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
             output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl
                    << "M_.aux_vars(" << i+1 << ").orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
             break;
+          case avMultiplier:
+            output << "M_.aux_vars(" << i+1 << ").eq_nbr = '" << aux_vars[i].get_equation_number_for_multiplier() << "';" << endl;
+            break;
           }
       }
 
@@ -285,7 +290,7 @@ SymbolTable::addLeadAuxiliaryVarInternal(bool endo, int index) throw (FrozenExce
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLead : avExoLead), 0, 0, ""));
+  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLead : avExoLead), 0, 0, "", 0));
 
   return symb_id;
 }
@@ -311,7 +316,7 @@ SymbolTable::addLagAuxiliaryVarInternal(bool endo, int orig_symb_id, int orig_le
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLag : avExoLag), orig_symb_id, orig_lead_lag, ""));
+  aux_vars.push_back(AuxVarInfo(symb_id, (endo ? avEndoLag : avExoLag), orig_symb_id, orig_lead_lag, "", 0));
 
   return symb_id;
 }
@@ -362,8 +367,29 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, const st
       exit(EXIT_FAILURE);
     }
 
-  aux_vars.push_back(AuxVarInfo(symb_id, (information_set_name.empty() ? avExpectation : avExpectationRIS), 0, 0, information_set_name));
+  aux_vars.push_back(AuxVarInfo(symb_id, (information_set_name.empty() ? avExpectation : avExpectationRIS), 0, 0, information_set_name, 0));
+
+  return symb_id;
+}
+
+int
+SymbolTable::addMultiplierAuxiliaryVar(int index) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+  varname << "MULT_" << index;
+
+  try
+    {
+      symb_id = addSymbol(varname.str(), eEndogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
 
+  aux_vars.push_back(AuxVarInfo(symb_id, avMultiplier, 0, 0, "", index));
   return symb_id;
 }
 

preprocessor/SymbolTable.hh

@@ -33,12 +33,13 @@ using namespace std;
 //! Types of auxiliary variables
 enum aux_var_t
   {
-    avEndoLead = 0,      //!< Substitute for endo leads >= 2
-    avEndoLag = 1,       //!< Substitute for endo lags >= 2
-    avExoLead = 2,       //!< Substitute for exo leads >= 2
-    avExoLag = 3,        //!< Substitute for exo lags >= 2
-    avExpectation = 4,   //!< Substitute for Expectation Operator
-    avExpectationRIS = 5 //!< Substitute for Expectation Operator Conditional on Restricted Information Set
+    avEndoLead = 0,       //!< Substitute for endo leads >= 2
+    avEndoLag = 1,        //!< Substitute for endo lags >= 2
+    avExoLead = 2,        //!< Substitute for exo leads >= 2
+    avExoLag = 3,         //!< Substitute for exo lags >= 2
+    avExpectation = 4,    //!< Substitute for Expectation Operator
+    avExpectationRIS = 5, //!< Substitute for Expectation Operator Conditional on Restricted Information Set
+    avMultiplier = 6      //!< Multipliers for FOC of Ramsey Probelem
   };
 
 //! Information on some auxiliary variables
@@ -50,13 +51,15 @@ private:
   int orig_symb_id; //!< Symbol ID of the endo of the original model represented by this aux var. Only used for avEndoLag and avExoLag.
   int orig_lead_lag; //!< Lead/lag of the endo of the original model represented by this aux var. Only used for avEndoLag and avExoLag.
   string expectation_information_set_name; //!< Stores 'full' or 'varobs' for avExpectationRIS. Not used otherwise.
+  int equation_number_for_multiplier; //!< Stores the original constraint equation number associated with this aux var. Only used for avMultiplier.
 public:
-  AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id, int orig_lead_lag, string expectation_information_set_name_arg);
+  AuxVarInfo(int symb_id_arg, aux_var_t type_arg, int orig_symb_id, int orig_lead_lag, string expectation_information_set_name_arg, int equation_number_for_multiplier_arg);
   int get_symb_id() const { return symb_id; };
   aux_var_t get_type() const { return type; };
   int get_orig_symb_id() const { return orig_symb_id; };
   int get_orig_lead_lag() const { return orig_lead_lag; };
   string get_expectation_information_set_name() const { return expectation_information_set_name; };
+  int get_equation_number_for_multiplier() const { return equation_number_for_multiplier; };
 };
 
 //! Stores the symbol table
@@ -214,6 +217,12 @@ public:
     \return the symbol ID of the new symbol
   */
   int addExpectationAuxiliaryVar(int information_set, int index, const string &information_set_name) throw (FrozenException);
+  //! Adds an auxiliary variable for the multiplier for the FOCs of the Ramsey Problem
+  /*!
+    \param[in] index Used to construct the variable name
+    \return the symbol ID of the new symbol
+  */
+  int addMultiplierAuxiliaryVar(int index) throw (FrozenException);
   //! Searches auxiliary variables which are substitutes for a given symbol_id and lead/lag
   /*!
     The search is only performed among auxiliary variables of endo/exo lag.