Implement tags [static] and [dynamic] for equations

Closes #307

doc/dynare.texi

@@ -183,6 +183,7 @@ Steady state
 
 * Finding the steady state with Dynare nonlinear solver::  
 * Using a steady state file::   
+* Replace some equations during steady state computations::  
 
 Stochastic solution and simulation
 
@@ -2329,6 +2330,7 @@ providing it with a ``steady state file''.
 @menu
 * Finding the steady state with Dynare nonlinear solver::  
 * Using a steady state file::   
+* Replace some equations during steady state computations::  
 @end menu
 
 @node Finding the steady state with Dynare nonlinear solver
@@ -2677,6 +2679,50 @@ steady;
 
 @end deffn
 
+@node Replace some equations during steady state computations
+@subsection Replace some equations during steady state computations
+
+When there is no steady state file, Dynare computes the steady state
+by solving the static model, @i{i.e.} the model from the @file{.mod}
+file from which leads and lags have been removed.
+
+In some specific cases, one may want to have more control over the way
+this static model is created. Dynare therefore offers the possibility
+to explicitly give the form of equations that should be in the static
+model.
+
+More precisely, if an equation is prepended by a @code{[static]} tag,
+then it will appear in the static model used for steady state
+computation, but that equation will not be used for other
+computations. For every equation tagged in this way, you must tag
+another equation with @code{[dynamic]}: that equation will not be used
+for steady state computation, but will be used for other computations.
+
+This functionality can be useful on models with a unit root, where
+there is an infinity of steady states. An equation (tagged
+@code{[dynamic]}) would give the law of motion of the nonstationary
+variable (like a random walk). To pin down one specific steady state,
+an equation tagged @code{[static]} would affect a constant value to
+the nonstationary variable.
+
+@examplehead
+
+This is a trivial example with two endogenous variables. The second equation
+takes a different form in the static model.
+
+@example
+var c k;
+varexo x;
+
+@dots{}
+
+model;
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+[dynamic] c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+[static] k = ((delt+bet)/(x*aa*alph))^(1/(alph-1));
+end;
+@end example
+
 
 @node Getting information about the model
 @section Getting information about the model

matlab/evaluate_steady_state.m

@@ -22,7 +22,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
 % SPECIAL REQUIREMENTS
 %   none
 
-% Copyright (C) 2001-2012 Dynare Team
+% Copyright (C) 2001-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -97,6 +97,31 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
         return
     end
 
+    % If some equations are tagged [static] or [dynamic], verify consistency
+    if M.static_and_dynamic_models_differ
+        % Evaluate residual of *dynamic* model using the steady state
+        % computed on the *static* one
+        z = repmat(ys,1,M.maximum_lead + M.maximum_lag + 1);
+        zx = repmat([oo.exo_simul oo.exo_det_simul],M.maximum_lead + M.maximum_lag + 1, 1);
+        if options.bytecode
+            [chck, r, junk]= bytecode('dynamic','evaluate', z, zx, M.params, ys, 1);
+            mexErrCheck('bytecode', chck);
+        elseif options.block
+            [r, data] = feval([M.fname '_dynamic'], z', zx, M.params, ys, M.maximum_lag+1, data);
+        else
+            iyv = M.lead_lag_incidence';
+            iyr0 = find(iyv(:));
+            xys = z(iyr0);
+            r = feval([M.fname '_dynamic'], z(iyr0), zx, M.params, ys, M.maximum_lag + 1);
+        end
+
+        % Fail if residual greater than tolerance
+        if max(abs(r)) > options.solve_tolf
+            info(1) = 25;
+            return
+        end
+    end
+    
     if ~isreal(ys)
         info(1) = 21;
         info(2) = sum(imag(ys).^2);

matlab/print_info.m

@@ -10,7 +10,7 @@ function print_info(info,noprint)
 % SPECIAL REQUIREMENTS
 %    none
 
-% Copyright (C) 2005-2012 Dynare Team
+% Copyright (C) 2005-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -78,6 +78,8 @@ if ~noprint
         error('Some updated params are complex')
       case 24
         error('Some updated params contain NaN or Inf')
+      case 25
+        error('The solution to the static equations is not a steady state of the dynamic model: verify that the equations tagged by [static] and [dynamic] are consistent')
       case 30
         error('Variance can''t be computed')
       case 41

preprocessor/DynamicModel.cc

@@ -2354,6 +2354,12 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
            << equation_tags[i].second.second << "' ;" << endl;
   output << "};" << endl;
 
+  /* Say if static and dynamic models differ (because of [static] and [dynamic]
+     equation tags) */
+  output << "M_.static_and_dynamic_models_differ = "
+         << (static_only_equations.size() > 0 ? "1" : "0")
+         << ";" << endl;
+
   //In case of sparse model, writes the block_decomposition structure of the model
   if (block_decomposition)
     {
@@ -3403,6 +3409,11 @@ DynamicModel::cloneDynamic(DynamicModel &dynamic_model) const
   for (deque<BinaryOpNode *>::const_iterator it = aux_equations.begin();
        it != aux_equations.end(); it++)
     dynamic_model.addAuxEquation((*it)->cloneDynamic(dynamic_model));
+
+  // Convert static_only equations
+  for (vector<BinaryOpNode *>::const_iterator it = static_only_equations.begin();
+       it != static_only_equations.end(); it++)
+    dynamic_model.addStaticOnlyEquation((*it)->cloneDynamic(dynamic_model));
 }
 
 void
@@ -3506,9 +3517,28 @@ DynamicModel::toStatic(StaticModel &static_model) const
     static_model.AddLocalVariable(it->first, it->second->toStatic(static_model));
 
   // Convert equations
-  for (vector<BinaryOpNode *>::const_iterator it = equations.begin();
-       it != equations.end(); it++)
-    static_model.addEquation((*it)->toStatic(static_model));
+  int static_only_index = 0;
+  for (int i = 0; i < (int) equations.size(); i++)
+    {
+      // Detect if equation is marked [dynamic]
+      bool is_dynamic_only = false;
+      for (vector<pair<int, pair<string, string> > >::const_iterator it = equation_tags.begin();
+           it != equation_tags.end(); ++it)
+        if (it->first == i && it->second.first == "dynamic")
+          {
+            is_dynamic_only = true;
+            break;
+          }
+
+      // If yes, replace it by an equation marked [static]
+      if (is_dynamic_only)
+        {
+          static_model.addEquation(static_only_equations[static_only_index]->toStatic(static_model));
+          static_only_index++;
+        }
+      else
+        static_model.addEquation(equations[i]->toStatic(static_model));
+    }
 
   // Convert auxiliary equations
   for (deque<BinaryOpNode *>::const_iterator it = aux_equations.begin();
@@ -4150,3 +4180,32 @@ DynamicModel::isModelLocalVariableUsed() const
     }
   return used_local_vars.size() > 0;
 }
+
+void
+DynamicModel::addStaticOnlyEquation(expr_t eq)
+{
+  BinaryOpNode *beq = dynamic_cast<BinaryOpNode *>(eq);
+  assert(beq != NULL && beq->get_op_code() == oEqual);
+
+  static_only_equations.push_back(beq);
+}
+
+size_t
+DynamicModel::staticOnlyEquationsNbr() const
+{
+  return static_only_equations.size();
+}
+
+size_t
+DynamicModel::dynamicOnlyEquationsNbr() const
+{
+  set<int> eqs;
+
+  for (vector<pair<int, pair<string, string> > >::const_iterator it = equation_tags.begin();
+       it != equation_tags.end(); ++it)
+    if (it->second.first == "dynamic")
+      eqs.insert(it->first);
+
+  return eqs.size();
+}
+

preprocessor/DynamicModel.hh

@@ -31,6 +31,10 @@ using namespace std;
 class DynamicModel : public ModelTree
 {
 private:
+  //! Stores equations declared as [static]
+  /*! They will be used in toStatic() to replace equations marked as [dynamic] */
+  vector<BinaryOpNode *> static_only_equations;
+
   typedef map<pair<int, int>, int> deriv_id_table_t;
   //! Maps a pair (symbol_id, lag) to a deriv ID
   deriv_id_table_t deriv_id_table;
@@ -224,6 +228,15 @@ public:
   //! Replaces the model equations in dynamic_model with those in this model
   void replaceMyEquations(DynamicModel &dynamic_model) const;
 
+  //! Adds an equation marked as [static]
+  void addStaticOnlyEquation(expr_t eq);
+
+  //! Returns number of static only equations
+  size_t staticOnlyEquationsNbr() const;
+  
+  //! Returns number of dynamic only equations
+  size_t dynamicOnlyEquationsNbr() const;
+
   //! Writes LaTeX file with the equations of the dynamic model
   void writeLatexFile(const string &basename) const;
 

preprocessor/DynareBison.yy

@@ -603,6 +603,8 @@ tags_list : tags_list COMMA tag_pair
 
 tag_pair : NAME EQUAL QUOTED_STRING
            { driver.add_equation_tags($1, $3); }
+         | NAME
+           { driver.add_equation_tags($1, new string()); }
          ;
 
 hand_side : '(' hand_side ')'

preprocessor/ExprNode.cc

@@ -469,6 +469,12 @@ NumConstNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
   return const_cast<NumConstNode *>(this);
 }
 
+bool
+NumConstNode::isInStaticForm() const
+{
+  return true;
+}
+
 VariableNode::VariableNode(DataTree &datatree_arg, int symb_id_arg, int lag_arg) :
   ExprNode(datatree_arg),
   symb_id(symb_id_arg),
@@ -1314,6 +1320,12 @@ VariableNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
     }
 }
 
+bool
+VariableNode::isInStaticForm() const
+{
+  return lag == 0;
+}
+
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, int param1_symb_id_arg, int param2_symb_id_arg) :
   ExprNode(datatree_arg),
   arg(arg_arg),
@@ -2322,6 +2334,17 @@ UnaryOpNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
+bool
+UnaryOpNode::isInStaticForm() const
+{
+  if (op_code == oSteadyState || op_code == oSteadyStateParamDeriv
+      || op_code == oSteadyStateParam2ndDeriv
+      || op_code == oExpectation)
+    return false;
+  else
+    return arg->isInStaticForm();
+}
+
 BinaryOpNode::BinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
                            BinaryOpcode op_code_arg, const expr_t arg2_arg) :
   ExprNode(datatree_arg),
@@ -3558,6 +3581,12 @@ BinaryOpNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+bool
+BinaryOpNode::isInStaticForm() const
+{
+  return arg1->isInStaticForm() && arg2->isInStaticForm();
+}
+
 TrinaryOpNode::TrinaryOpNode(DataTree &datatree_arg, const expr_t arg1_arg,
                              TrinaryOpcode op_code_arg, const expr_t arg2_arg, const expr_t arg3_arg) :
   ExprNode(datatree_arg),
@@ -4155,6 +4184,12 @@ TrinaryOpNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+bool
+TrinaryOpNode::isInStaticForm() const
+{
+  return arg1->isInStaticForm() && arg2->isInStaticForm() && arg3->isInStaticForm();
+}
+
 ExternalFunctionNode::ExternalFunctionNode(DataTree &datatree_arg,
                                            int symb_id_arg,
                                            const vector<expr_t> &arguments_arg) :
@@ -4732,6 +4767,16 @@ ExternalFunctionNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) con
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+bool
+ExternalFunctionNode::isInStaticForm() const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); ++it)
+    if (!(*it)->isInStaticForm())
+      return false;
+  
+  return true;
+}
+
 FirstDerivExternalFunctionNode::FirstDerivExternalFunctionNode(DataTree &datatree_arg,
                                                                int top_level_symb_id_arg,
                                                                const vector<expr_t> &arguments_arg,

preprocessor/ExprNode.hh

 /*
- * Copyright (C) 2007-2012 Dynare Team
+ * Copyright (C) 2007-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -393,6 +393,9 @@ public:
 
   //! Move a trend variable with lag/lead to time t by dividing/multiplying by its growth factor
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const = 0;
+
+  //! Returns true if the expression is in static form (no lead, no lag, no expectation, no STEADY_STATE)
+  virtual bool isInStaticForm() const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -448,6 +451,7 @@ public:
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
+  virtual bool isInStaticForm() const;
 };
 
 //! Symbol or variable node
@@ -508,6 +512,7 @@ public:
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
+  virtual bool isInStaticForm() const;
 };
 
 //! Unary operator node
@@ -583,6 +588,7 @@ public:
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
+  virtual bool isInStaticForm() const;
 };
 
 //! Binary operator node
@@ -677,6 +683,7 @@ public:
   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;
+  virtual bool isInStaticForm() const;
 };
 
 //! Trinary operator node
@@ -739,6 +746,7 @@ public:
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
+  virtual bool isInStaticForm() const;
 };
 
 //! External function node
@@ -812,6 +820,7 @@ public:
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
+  virtual bool isInStaticForm() const;
 };
 
 class FirstDerivExternalFunctionNode : public ExternalFunctionNode

preprocessor/ModFile.cc

@@ -236,6 +236,19 @@ ModFile::checkPass()
              << "(deprecated) or the dsge_var option must be passed to the estimation statement (preferred)." << endl;
         exit(EXIT_FAILURE);
       }
+
+  if (dynamic_model.staticOnlyEquationsNbr() != dynamic_model.dynamicOnlyEquationsNbr())
+    {
+      cerr << "ERROR: the number of equations marked [static] must be equal to the number of equations marked [dynamic]" << endl;
+      exit(EXIT_FAILURE);
+    }
+  
+  if (dynamic_model.staticOnlyEquationsNbr() > 0 &&
+      (mod_file_struct.ramsey_policy_present || mod_file_struct.discretionary_policy_present))
+    {
+      cerr << "ERROR: marking equations as [static] or [dynamic] is not possible with ramsey_policy or discretionary_policy" << endl;
+      exit(EXIT_FAILURE);
+    }
 }
 
 void

preprocessor/ModelTree.cc

@@ -1385,9 +1385,12 @@ ModelTree::addEquation(expr_t eq)
 }
 
 void
-ModelTree::addEquationTags(int i, const string &key, const string &value)
+ModelTree::addEquation(expr_t eq, vector<pair<string, string> > &eq_tags)
 {
-  equation_tags.push_back(make_pair(i, make_pair(key, value)));
+  int n = equation_number();
+  for (size_t i = 0; i < eq_tags.size(); i++)
+    equation_tags.push_back(make_pair(n, eq_tags[i]));
+  addEquation(eq);
 }
 
 void

preprocessor/ModelTree.hh

@@ -297,8 +297,8 @@ public:
   int mfs;
   //! Declare a node as an equation of the model
   void addEquation(expr_t eq);
-  //! Adds tags to equation number i
-  void addEquationTags(int i, const string &key, const string &value);
+  //! Declare a node as an equation of the model, also giving its tags
+  void addEquation(expr_t eq, vector<pair<string, string> > &eq_tags);
   //! Declare a node as an auxiliary equation of the model, adding it at the end of the list of auxiliary equations
   void addAuxEquation(expr_t eq);
   //! Returns the number of equations in the model

preprocessor/ParsingDriver.cc

@@ -249,8 +249,7 @@ ParsingDriver::add_predetermined_variable(string *name)
 void
 ParsingDriver::add_equation_tags(string *key, string *value)
 {
-  int n = model_tree->equation_number();
-  model_tree->addEquationTags(n, *key, *value);
+  eq_tags.push_back(make_pair(*key, *value));
   delete key;
   delete value;
 }
@@ -1933,7 +1932,28 @@ expr_t
 ParsingDriver::add_model_equal(expr_t arg1, expr_t arg2)
 {
   expr_t id = model_tree->AddEqual(arg1, arg2);
-  model_tree->addEquation(id);
+
+  // Detect if the equation is tagged [static]
+  bool is_static_only = false;  
+  for (vector<pair<string, string> >::const_iterator it = eq_tags.begin();
+       it != eq_tags.end(); ++it)
+    if (it->first == "static")
+      {
+        is_static_only = true;
+        break;
+      }
+
+  if (is_static_only)
+    {
+      if (!id->isInStaticForm())
+        error("An equation tagged [static] cannot contain leads, lags, expectations or STEADY_STATE operators");
+      
+      dynamic_model->addStaticOnlyEquation(id);
+    }
+  else
+    model_tree->addEquation(id, eq_tags);
+
+  eq_tags.clear();
   return id;
 }
 

preprocessor/ParsingDriver.hh

@@ -201,6 +201,8 @@ private:
   expr_t add_model_variable(int symb_id, int lag);
   //! For parsing the graph_format option
   SymbolList graph_formats;
+  //! Temporary storage for equation tags
+  vector<pair<string, string> > eq_tags;
 
   //! The mod file representation constructed by this ParsingDriver
   ModFile *mod_file;

tests/Makefile.am

@@ -6,6 +6,7 @@ MODFILES = \
 	gsa/ls2003.mod \
 	ramst.mod \
 	ramst_a.mod \
+	ramst_static_tag.mod \
 	example1.mod \
 	example2.mod \
 	example1_use_dll.mod \

tests/ramst_static_tag.mod

+var c k;
+varexo x;
+
+parameters alph gam delt bet aa;
+alph=0.5;
+gam=0.5;
+delt=0.02;
+bet=0.05;
+aa=0.5;
+
+
+model;
+c + k - aa*x*k(-1)^alph - (1-delt)*k(-1);
+[dynamic] c^(-gam) - (1+bet)^(-1)*(aa*alph*x(+1)*k^(alph-1) + 1 - delt)*c(+1)^(-gam);
+[static] k = ((delt+bet)/(x*aa*alph))^(1/(alph-1));
+end;
+
+initval;
+x = 1;
+k = ((delt+bet)/(1.0*aa*alph))^(1/(alph-1));
+c = aa*k^alph-delt*k;
+end;
+
+steady;
+
+check;
+
+shocks;
+var x;
+periods 1;
+values 1.2;
+end;
+
+simul(periods=200);
+
+rplot c;
+rplot k;