diff --git a/src/ComputingTasks.cc b/src/ComputingTasks.cc
index 01c2cf84c390253a664b91e39a68a26154c4be82..61f91926f05f4dab1bc730aee554d53ff6ed8964 100644
--- a/src/ComputingTasks.cc
+++ b/src/ComputingTasks.cc
@@ -262,12 +262,16 @@ PacModelStatement::PacModelStatement(string name_arg,
                                      string discount_arg,
                                      int growth_symb_id_arg,
                                      int growth_lag_arg,
+                                     double steady_state_growth_rate_number_arg,
+                                     int steady_state_growth_rate_symb_id_arg,
                                      const SymbolTable &symbol_table_arg) :
   name{move(name_arg)},
   aux_model_name{move(aux_model_name_arg)},
   discount{move(discount_arg)},
   growth_symb_id{growth_symb_id_arg},
   growth_lag{growth_lag_arg},
+  steady_state_growth_rate_number{steady_state_growth_rate_number_arg},
+  steady_state_growth_rate_symb_id{steady_state_growth_rate_symb_id_arg},
   symbol_table{symbol_table_arg}
 {
 }
@@ -306,7 +310,12 @@ void
 PacModelStatement::writeOutput(ostream &output, const string &basename, bool minimal_workspace) const
 {
   output << "M_.pac." << name << ".auxiliary_model_name = '" << aux_model_name << "';" << endl
-         << "M_.pac." << name << ".discount_index = " << symbol_table.getTypeSpecificID(discount) + 1 << ";" << endl;
+         << "M_.pac." << name << ".discount_index = " << symbol_table.getTypeSpecificID(discount) + 1 << ";" << endl
+         << "M_.pac." << name << ".steady_state_growth_rate = ";
+  if (steady_state_growth_rate_symb_id < 0)
+    output << steady_state_growth_rate_number << ";" << endl;
+  else
+    output << symbol_table.getTypeSpecificID(steady_state_growth_rate_symb_id) + 1 << ";" << endl;
 
   if (growth_symb_id >= 0)
     {
@@ -399,12 +408,6 @@ PacModelStatement::writeJsonOutput(ostream &output) const
   output << "}";
 }
 
-tuple<string, string, int, int>
-PacModelStatement::getPacModelInfoForPacExpectation() const
-{
-  return { name, aux_model_name, growth_symb_id, growth_lag };
-}
-
 VarEstimationStatement::VarEstimationStatement(OptionsList options_list_arg) :
   options_list{move(options_list_arg)}
 {
diff --git a/src/ComputingTasks.hh b/src/ComputingTasks.hh
index f61b52f307249af8637b48cac7912092970d97b3..5a37b8b7f8b07843a99b665572e71b9f87b8fc34 100644
--- a/src/ComputingTasks.hh
+++ b/src/ComputingTasks.hh
@@ -136,11 +136,12 @@ public:
 
 class PacModelStatement : public Statement
 {
-private:
-  const string name;
-  const string aux_model_name;
-  const string discount;
+public:
+  const string name, aux_model_name, discount;
   const int growth_symb_id, growth_lag;
+private:
+  const double steady_state_growth_rate_number;
+  const int steady_state_growth_rate_symb_id;
   const SymbolTable &symbol_table;
   vector<int> lhs;
 public:
@@ -149,12 +150,13 @@ public:
                     string discount_arg,
                     int growth_symb_id_arg,
                     int growth_lag_arg,
+                    double steady_state_growth_rate_number_arg,
+                    int steady_state_growth_rate_symb_id_arg,
                     const SymbolTable &symbol_table_arg);
   void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings) override;
   void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const override;
   void writeJsonOutput(ostream &output) const override;
   void fillUndiffedLHS(vector<int> &lhs);
-  tuple<string, string, int, int> getPacModelInfoForPacExpectation() const;
 };
 
 class VarRestrictionsStatement : public Statement
diff --git a/src/DynamicModel.cc b/src/DynamicModel.cc
index 3cb7a05876cd37044a263e9e4327f84137117016..0b2c82500186b6856f968822e6ec0c370c38cb97 100644
--- a/src/DynamicModel.cc
+++ b/src/DynamicModel.cc
@@ -3615,6 +3615,15 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
     output << "-1";
   output << "];" << endl;
 
+  // Write Pac Model Consistent Expectation parameter info
+  for (auto & it : pac_mce_alpha_symb_ids)
+    {
+      output << modstruct << "pac." << it.first << ".mce_alpha_idxs = [";
+      for (auto it : it.second)
+        output << symbol_table.getTypeSpecificID(it) + 1 << " ";
+      output << "];" << endl;
+    }
+
   // Write PacExpectationInfo
   for (auto it : pac_expectation_info)
     it->ExprNode::writeOutput(output, ExprNodeOutputType::matlabDynamicModel);
@@ -4312,8 +4321,128 @@ DynamicModel::getPacMaxLag(const string &pac_model_name) const
   return 0;
 }
 
+int
+DynamicModel::getPacTargetSymbId(const string &pac_model_name) const
+{
+  for (auto & equation : equations)
+    if (equation->containsPacExpectation(pac_model_name))
+      {
+        pair<int, int> lhs (-1, -1);
+        equation->getPacLHS(lhs);
+        int lhs_symb_id = lhs.first;
+        int lhs_orig_symb_id = lhs_symb_id;
+        if (symbol_table.isAuxiliaryVariable(lhs_symb_id))
+          try
+            {
+              lhs_orig_symb_id = symbol_table.getOrigSymbIdForAuxVar(lhs_symb_id);
+            }
+          catch (...)
+            {
+            }
+        return equation->arg2->getPacTargetSymbId(lhs_symb_id, lhs_orig_symb_id);
+      }
+  return -1;
+}
+
+int
+DynamicModel::addPacModelConsistentExpectationEquation(const string & name, int pac_target_symb_id,
+                                                       int discount_symb_id, int pac_max_lag_m,
+                                                       ExprNode::subst_table_t &diff_subst_table)
+{
+  int mce_symb_id;
+  try
+    {
+      mce_symb_id = symbol_table.addSymbol("mce_Z_" + name, SymbolType::endogenous);
+    }
+  catch (SymbolTable::AlreadyDeclaredException &e)
+    {
+      cerr << "Variable name needed by PAC (mce_Z_" << name << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  expr_t A = One;
+  expr_t fp = Zero;
+  expr_t beta = AddVariable(discount_symb_id);
+  for (int i = 1; i <= pac_max_lag_m; i++)
+    try
+      {
+        int alpha_i_symb_id = symbol_table.addSymbol("mce_alpha_" + name + "_" + to_string(i),
+                                                     SymbolType::parameter);
+        pac_mce_alpha_symb_ids[name].push_back(alpha_i_symb_id);
+        A = AddPlus(A, AddVariable(alpha_i_symb_id));
+        fp = AddPlus(fp,
+                     AddTimes(AddTimes(AddVariable(alpha_i_symb_id),
+                                       AddPower(beta, AddPossiblyNegativeConstant(i))),
+                              AddVariable(mce_symb_id, i)));
+
+      }
+    catch (SymbolTable::AlreadyDeclaredException &e)
+      {
+        cerr << "Variable name needed by PAC (mce_alpha_" << name << "_" << i << ")" << endl;
+        exit(EXIT_FAILURE);
+      }
+
+  // Add diff nodes and eqs for pac_target_symb_id
+  int neqs = 0;
+  const VariableNode *target_base_diff_node;
+  expr_t diff_node_to_search = AddDiff(AddVariable(pac_target_symb_id));
+  ExprNode::subst_table_t::const_iterator sit = diff_subst_table.find(diff_node_to_search);
+  if (sit != diff_subst_table.end())
+    target_base_diff_node = sit->second;
+  else
+    {
+      int symb_id = symbol_table.addDiffAuxiliaryVar(diff_node_to_search->idx, diff_node_to_search);
+      target_base_diff_node = AddVariable(symb_id);
+      addEquation(dynamic_cast<BinaryOpNode *>(AddEqual((expr_t) target_base_diff_node,
+                                                        AddMinus(AddVariable(pac_target_symb_id),
+                                                                 AddVariable(pac_target_symb_id, -1)))), -1);
+      neqs++;
+    }
+
+  map<int, VariableNode *> target_aux_var_to_add;
+  const VariableNode *last_aux_var = target_base_diff_node;
+  for (int i = 1; i <= pac_max_lag_m - 1; i++, neqs++)
+    {
+      expr_t this_diff_node = AddDiff(AddVariable(pac_target_symb_id, i));
+      int symb_id = symbol_table.addDiffLeadAuxiliaryVar(this_diff_node->idx, this_diff_node,
+                                                         last_aux_var->symb_id, last_aux_var->lag);
+      VariableNode *current_aux_var = AddVariable(symb_id);
+      addEquation(dynamic_cast<BinaryOpNode *>(AddEqual(current_aux_var,
+                                                        AddVariable(last_aux_var->symb_id, 1))), -1);
+      last_aux_var = current_aux_var;
+      target_aux_var_to_add[i] = current_aux_var;
+    }
+
+  expr_t fs = Zero;
+  for (int k = 1; k <= pac_max_lag_m - 1; k++)
+    {
+      expr_t ssum = Zero;
+      for (int j = k+1; j <= pac_max_lag_m; j++)
+        {
+          int alpha_j_symb_id = -1;
+          string varname = "mce_alpha_" + name + "_" + to_string(j);
+          try
+            {
+              alpha_j_symb_id = symbol_table.getID(varname);
+            }
+          catch (SymbolTable::UnknownSymbolNameException &e)
+            {
+              alpha_j_symb_id = symbol_table.addSymbol(varname, SymbolType::parameter);
+            }
+          ssum = AddPlus(ssum,
+                         AddTimes(AddVariable(alpha_j_symb_id), AddPower(beta, AddPossiblyNegativeConstant(j))));
+        }
+      fs = AddPlus(fs, AddTimes(ssum, target_aux_var_to_add[k]));
+    }
+  addEquation(AddEqual(AddVariable(mce_symb_id),
+                       AddMinus(AddTimes(A, AddMinus((expr_t) target_base_diff_node, fs)), fp)), -1);
+  neqs++;
+  cout << "Pac Model Consistent Expectation: added " << neqs << " auxiliary variables and equations." << endl;
+  return mce_symb_id;
+}
+
 void
-DynamicModel::fillPacExpectationVarInfo(string &pac_model_name,
+DynamicModel::fillPacExpectationVarInfo(const string &pac_model_name,
                                         vector<int> &lhs,
                                         int max_lag,
                                         int pac_max_lag,
@@ -4326,15 +4455,34 @@ DynamicModel::fillPacExpectationVarInfo(string &pac_model_name,
 }
 
 void
-DynamicModel::substitutePacExpectation()
+DynamicModel::substitutePacExpectation(const string & name, int model_consistent_expectation_symb_id)
+{
+  map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
+  for (auto & it : local_variables_table)
+    it.second = it.second->substitutePacExpectation(name, model_consistent_expectation_symb_id, subst_table);
+
+  for (auto & equation : equations)
+    {
+      auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(name, model_consistent_expectation_symb_id, subst_table));
+      assert(substeq != nullptr);
+      equation = substeq;
+    }
+
+  for (map<const PacExpectationNode *, const BinaryOpNode *>::const_iterator it = subst_table.begin();
+       it != subst_table.end(); it++)
+    pac_expectation_info.insert(const_cast<PacExpectationNode *>(it->first));
+}
+
+void
+DynamicModel::substitutePacExpectation(const string & name)
 {
   map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
   for (auto & it : local_variables_table)
-    it.second = it.second->substitutePacExpectation(subst_table);
+    it.second = it.second->substitutePacExpectation(name, subst_table);
 
   for (auto & equation : equations)
     {
-      auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(subst_table));
+      auto *substeq = dynamic_cast<BinaryOpNode *>(equation->substitutePacExpectation(name, subst_table));
       assert(substeq != nullptr);
       equation = substeq;
     }
diff --git a/src/DynamicModel.hh b/src/DynamicModel.hh
index 0f77e8f825236632c24c7b6bca5389601068922b..a4b665c13c6179ace5664faaf36e1dbefa57a105 100644
--- a/src/DynamicModel.hh
+++ b/src/DynamicModel.hh
@@ -346,15 +346,18 @@ public:
   //! Get Pac equation parameter info
   void walkPacParameters();
   //! Add var_model info to pac_expectation nodes
-  void fillPacExpectationVarInfo(string &pac_model_name,
+  void fillPacExpectationVarInfo(const string &pac_model_name,
                                  vector<int> &lhs,
                                  int max_lag,
                                  int pac_max_lag,
                                  vector<bool> &nonstationary,
                                  int growth_symb_id, int growth_lag);
 
-  //! Substitutes pac_expectation operator
-  void substitutePacExpectation();
+  //! Substitutes pac_expectation operator with expectation based on auxiliary model
+  void substitutePacExpectation(const string & name);
+
+  //! Substitute pac_expectation operator with model consistent expectation
+  void substitutePacExpectation(const string & name, int model_consistent_expectation_symb_id);
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &basename,
@@ -454,6 +457,15 @@ public:
   //! Return max lag of pac equation
   int getPacMaxLag(const string &pac_model_name) const;
 
+  //! Return target of the pac equation
+  int getPacTargetSymbId(const string &pac_model_name) const;
+
+  //! Add model consistent expectation equation for pac model
+  int addPacModelConsistentExpectationEquation(const string & name, int pac_target_symb_id, int discount, int pac_max_lag_m, ExprNode::subst_table_t &diff_subst_table);
+
+  //! store symb_ids for alphas created in addPacModelConsistentExpectationEquation
+  map<string, vector<int>> pac_mce_alpha_symb_ids;
+
   //! Table to undiff LHS variables for pac vector z
   vector<int> getUndiffLHSForPac(const string &aux_model_name,
                                  ExprNode::subst_table_t &diff_subst_table) const;
diff --git a/src/DynareBison.yy b/src/DynareBison.yy
index 9d3cfecb1212aff59b6a5567c44bd3998674d15a..75f7d105157f974d447d1e63b4247b74cca88f7e 100644
--- a/src/DynareBison.yy
+++ b/src/DynareBison.yy
@@ -163,7 +163,7 @@ class ParsingDriver;
 %token NUMBER_OF_POSTERIOR_DRAWS_AFTER_PERTURBATION MAX_NUMBER_OF_STAGES
 %token RANDOM_FUNCTION_CONVERGENCE_CRITERION RANDOM_PARAMETER_CONVERGENCE_CRITERION
 %token CENTERED_MOMENTS AUTOLAG RECURSIVE_ORDER_ESTIMATION BARTLETT_KERNEL_LAG WEIGHTING_MATRIX PENALIZED_ESTIMATOR VERBOSE
-%token SIMULATION_MULTIPLE SEED BOUNDED_SHOCK_SUPPORT EQTAGS
+%token SIMULATION_MULTIPLE SEED BOUNDED_SHOCK_SUPPORT EQTAGS STEADY_STATE_GROWTH
 %token ANALYTICAL_GIRF IRF_IN_PERCENT EMAS_GIRF EMAS_DROP EMAS_TOLF EMAS_MAX_ITER
 
 %token <vector<string>> SYMBOL_VEC
@@ -398,6 +398,7 @@ pac_model_options : o_pac_name
                   | o_pac_aux_model_name
                   | o_pac_discount
                   | o_pac_growth
+                  | o_pac_steady_state_growth
                   ;
 
 var_expectation_model : VAR_EXPECTATION_MODEL '(' var_expectation_model_options_list ')' ';'
@@ -3133,6 +3134,9 @@ o_pac_discount : DISCOUNT EQUAL symbol { driver.option_str("pac.discount", $3);
 o_pac_growth : GROWTH EQUAL symbol { driver.set_pac_growth($3, 0); }
              | GROWTH EQUAL symbol '(' MINUS INT_NUMBER ')' { driver.set_pac_growth($3, stoi($6)); }
              ;
+o_pac_steady_state_growth : STEADY_STATE_GROWTH EQUAL signed_number { driver.set_pac_steady_state_growth($3); }
+                          | STEADY_STATE_GROWTH EQUAL symbol { driver.set_pac_steady_state_growth($3); }
+                          ;
 o_var_name : MODEL_NAME EQUAL symbol { driver.option_str("var.model_name", $3); };
 o_var_order : ORDER EQUAL INT_NUMBER { driver.option_num("var.order", $3); };
 o_series : SERIES EQUAL symbol { driver.option_str("series", $3); };
diff --git a/src/DynareFlex.ll b/src/DynareFlex.ll
index c703738c3beec7def0e79691de18fd8581527c63..22b3f8c00b7eb069871e52212bfc84b5086fdad5 100644
--- a/src/DynareFlex.ll
+++ b/src/DynareFlex.ll
@@ -1,6 +1,6 @@
 /* -*- C++ -*- */
 /*
- * Copyright (C) 2003-2018 Dynare Team
+ * Copyright (C) 2003-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -830,6 +830,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
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
 <DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
 <DYNARE_STATEMENT>discount {return token::DISCOUNT;}
+<DYNARE_STATEMENT>steady_state_growth {return token::STEADY_STATE_GROWTH;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>nan {return token::NAN_CONSTANT;}
diff --git a/src/ExprNode.cc b/src/ExprNode.cc
index 04bd29f5fa854dad91b586e2c99aafa1db96d4ef..c19698abf2cbd0e62ea50c92c4213b07fe884007 100644
--- a/src/ExprNode.cc
+++ b/src/ExprNode.cc
@@ -506,6 +506,12 @@ NumConstNode::PacMaxLag(int lhs_symb_id) const
   return 0;
 }
 
+int
+NumConstNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  return -1;
+}
+
 expr_t
 NumConstNode::decreaseLeadsLags(int n) const
 {
@@ -589,7 +595,14 @@ NumConstNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &no
 }
 
 expr_t
-NumConstNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+NumConstNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  return const_cast<NumConstNode *>(this);
+}
+
+expr_t
+NumConstNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -687,7 +700,7 @@ NumConstNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pa
 }
 
 void
-NumConstNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+NumConstNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
 }
 
@@ -1575,6 +1588,12 @@ VariableNode::PacMaxLag(int lhs_symb_id) const
   return 0;
 }
 
+int
+VariableNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  return -1;
+}
+
 expr_t
 VariableNode::substituteAdl() const
 {
@@ -1617,7 +1636,14 @@ VariableNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &no
 }
 
 expr_t
-VariableNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+VariableNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  return const_cast<VariableNode *>(this);
+}
+
+expr_t
+VariableNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<VariableNode *>(this);
 }
@@ -1992,7 +2018,7 @@ VariableNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pa
 }
 
 void
-VariableNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+VariableNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
 }
 
@@ -3294,6 +3320,12 @@ UnaryOpNode::PacMaxLag(int lhs_symb_id) const
   return arg->PacMaxLag(lhs_symb_id);
 }
 
+int
+UnaryOpNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  return arg->getPacTargetSymbId(lhs_symb_id, undiff_lhs_symb_id);
+}
+
 expr_t
 UnaryOpNode::substituteAdl() const
 {
@@ -3610,9 +3642,17 @@ UnaryOpNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nod
 }
 
 expr_t
-UnaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+UnaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t argsubst = arg->substitutePacExpectation(subst_table);
+  expr_t argsubst = arg->substitutePacExpectation(name, subst_table);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
+expr_t
+UnaryOpNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                      map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  expr_t argsubst = arg->substitutePacExpectation(name, mce_symb_id, subst_table);
   return buildSimilarUnaryOpNode(argsubst, datatree);
 }
 
@@ -3817,7 +3857,7 @@ UnaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pai
 }
 
 void
-UnaryOpNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+UnaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   arg->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
 }
@@ -5211,6 +5251,39 @@ BinaryOpNode::PacMaxLag(int lhs_symb_id) const
   return max(arg1->PacMaxLag(lhs_symb_id), arg2->PacMaxLag(lhs_symb_id));
 }
 
+int
+BinaryOpNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  set<pair<int, int>> endogs;
+  arg1->collectDynamicVariables(SymbolType::endogenous, endogs);
+  arg2->collectDynamicVariables(SymbolType::endogenous, endogs);
+  int ret_symb_id = -1;
+  for (auto & it : endogs)
+    {
+      int id = it.first;
+      while (datatree.symbol_table.isAuxiliaryVariable(id))
+        try
+          {
+            id = datatree.symbol_table.getOrigSymbIdForAuxVar(id);
+          }
+        catch (...)
+          {
+            break;
+          }
+      if (id != lhs_symb_id && id != undiff_lhs_symb_id)
+        if (ret_symb_id < 0)
+          ret_symb_id = it.first;
+        else
+          {
+            cerr << "Error: Pac model is of wrong format: endogenous vars found: "
+                 << datatree.symbol_table.getName(ret_symb_id) << ", "
+                 << datatree.symbol_table.getName(it.first) << endl;
+            exit(EXIT_FAILURE);
+          }
+    }
+  return ret_symb_id;
+}
+
 expr_t
 BinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -5394,10 +5467,19 @@ BinaryOpNode::countDiffs() const
 }
 
 expr_t
-BinaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+BinaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  expr_t arg1subst = arg1->substitutePacExpectation(name, subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, subst_table);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
+expr_t
+BinaryOpNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
-  expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
+  expr_t arg1subst = arg1->substitutePacExpectation(name, mce_symb_id, subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, mce_symb_id, subst_table);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -5903,7 +5985,7 @@ BinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pa
 }
 
 void
-BinaryOpNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+BinaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
   arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
@@ -6562,6 +6644,20 @@ TrinaryOpNode::PacMaxLag(int lhs_symb_id) const
   return max(arg1->PacMaxLag(lhs_symb_id), max(arg2->PacMaxLag(lhs_symb_id), arg3->PacMaxLag(lhs_symb_id)));
 }
 
+int
+TrinaryOpNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  int symb_id = arg1->getPacTargetSymbId(lhs_symb_id, undiff_lhs_symb_id);
+  if (symb_id >= 0)
+    return symb_id;
+
+  symb_id = arg2->getPacTargetSymbId(lhs_symb_id, undiff_lhs_symb_id);
+  if (symb_id >= 0)
+    return symb_id;
+
+  return arg3->getPacTargetSymbId(lhs_symb_id, undiff_lhs_symb_id);
+}
+
 expr_t
 TrinaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -6711,11 +6807,21 @@ TrinaryOpNode::countDiffs() const
 }
 
 expr_t
-TrinaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+TrinaryOpNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
-  expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
-  expr_t arg3subst = arg3->substitutePacExpectation(subst_table);
+  expr_t arg1subst = arg1->substitutePacExpectation(name, subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, subst_table);
+  expr_t arg3subst = arg3->substitutePacExpectation(name, subst_table);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
+expr_t
+TrinaryOpNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                        map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  expr_t arg1subst = arg1->substitutePacExpectation(name, mce_symb_id, subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(name, mce_symb_id, subst_table);
+  expr_t arg3subst = arg2->substitutePacExpectation(name, mce_symb_id, subst_table);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -6827,7 +6933,7 @@ TrinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, p
 }
 
 void
-TrinaryOpNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+TrinaryOpNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   arg1->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
   arg2->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
@@ -7086,6 +7192,19 @@ AbstractExternalFunctionNode::PacMaxLag(int lhs_symb_id) const
   return val;
 }
 
+int
+AbstractExternalFunctionNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  int symb_id = -1;
+  for (auto argument : arguments)
+    {
+      symb_id = argument->getPacTargetSymbId(lhs_symb_id, undiff_lhs_symb_id);
+      if (symb_id >= 0)
+        return symb_id;
+    }
+  return -1;
+}
+
 expr_t
 AbstractExternalFunctionNode::decreaseLeadsLags(int n) const
 {
@@ -7222,11 +7341,21 @@ AbstractExternalFunctionNode::countDiffs() const
 }
 
 expr_t
-AbstractExternalFunctionNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+AbstractExternalFunctionNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   vector<expr_t> arguments_subst;
   for (auto argument : arguments)
-    arguments_subst.push_back(argument->substitutePacExpectation(subst_table));
+    arguments_subst.push_back(argument->substitutePacExpectation(name, subst_table));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
+expr_t
+AbstractExternalFunctionNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                                       map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+   vector<expr_t> arguments_subst;
+  for (auto argument : arguments)
+    arguments_subst.push_back(argument->substitutePacExpectation(name, mce_symb_id, subst_table));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -7395,7 +7524,7 @@ AbstractExternalFunctionNode::addParamInfoToPac(pair<int, int> &lhs_arg, int opt
 }
 
 void
-AbstractExternalFunctionNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+AbstractExternalFunctionNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   for (auto argument : arguments)
     argument->fillPacExpectationVarInfo(model_name_arg, lhs_arg, max_lag_arg, pac_max_lag_arg, nonstationary_arg, growth_symb_id_arg, growth_lag_arg, equation_number_arg);
@@ -8746,6 +8875,12 @@ VarExpectationNode::PacMaxLag(int lhs_symb_id) const
   exit(EXIT_FAILURE);
 }
 
+int
+VarExpectationNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  return -1;
+}
+
 expr_t
 VarExpectationNode::decreaseLeadsLags(int n) const
 {
@@ -8916,7 +9051,14 @@ VarExpectationNode::substituteUnaryOpNodes(DataTree &static_datatree, diff_table
 }
 
 expr_t
-VarExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+VarExpectationNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
+expr_t
+VarExpectationNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                             map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -9033,7 +9175,7 @@ VarExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_a
 }
 
 void
-VarExpectationNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+VarExpectationNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
 }
 
@@ -9319,6 +9461,12 @@ PacExpectationNode::PacMaxLag(int lhs_symb_id) const
   return 0;
 }
 
+int
+PacExpectationNode::getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const
+{
+  return -1;
+}
+
 expr_t
 PacExpectationNode::decreaseLeadsLags(int n) const
 {
@@ -9650,7 +9798,7 @@ PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_a
 
 
 void
-PacExpectationNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
+PacExpectationNode::fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag_arg, int equation_number_arg)
 {
   if (model_name != model_name_arg)
     return;
@@ -9675,8 +9823,11 @@ PacExpectationNode::fillPacExpectationVarInfo(string &model_name_arg, vector<int
 }
 
 expr_t
-PacExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+PacExpectationNode::substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
+  if (model_name != name)
+    return const_cast<PacExpectationNode *>(this);
+
   map<const PacExpectationNode *, const BinaryOpNode *>::const_iterator myit =
     subst_table.find(const_cast<PacExpectationNode *>(this));
   if (myit != subst_table.end())
@@ -9728,6 +9879,18 @@ PacExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, con
   return subExpr;
 }
 
+expr_t
+PacExpectationNode::substitutePacExpectation(const string & name, int mce_symb_id,
+                                             map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  if (model_name != name)
+    return const_cast<PacExpectationNode *>(this);
+
+  expr_t subExpr = datatree.AddVariable(mce_symb_id);
+  subst_table[const_cast<PacExpectationNode *>(this)] = dynamic_cast<BinaryOpNode *>(subExpr);
+  return subExpr;
+}
+
 void
 ExprNode::decomposeAdditiveTerms(vector<pair<expr_t, int>> &terms, int current_sign) const
 {
diff --git a/src/ExprNode.hh b/src/ExprNode.hh
index ff84420ce946795c2d876bc28aff1eb14f786612..6c7f4aff1e93f0396c1d7c56981e88112e356e60 100644
--- a/src/ExprNode.hh
+++ b/src/ExprNode.hh
@@ -425,6 +425,9 @@ class ExprNode
       //! Takes account of undiffed LHS variables in calculating the max lag
       virtual int PacMaxLag(int lhs_symb_id) const = 0;
 
+      //! Get the target variable of the PAC model
+      virtual int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const = 0;
+
       virtual expr_t undiff() const = 0;
 
       //! Returns a new expression where all the leads/lags have been shifted backwards by the same amount
@@ -563,7 +566,11 @@ class ExprNode
       virtual expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const = 0;
 
       //! Substitute pac_expectation operator
-      virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
+      virtual expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
+
+      //! Substitute pac_expectation operator with model consistent expectation
+      virtual expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                              map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) = 0;
 
       //! Add ExprNodes to the provided datatree
       virtual expr_t clone(DataTree &datatree) const = 0;
@@ -604,7 +611,7 @@ class ExprNode
       virtual void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) = 0;
 
       //! Fills var_model info for pac_expectation node
-      virtual void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) = 0;
+      virtual void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) = 0;
 
       //! Fills the AR matrix structure
       virtual void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const = 0;
@@ -702,6 +709,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -716,7 +724,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -730,7 +740,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -793,6 +803,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -807,7 +818,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -821,7 +834,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -909,6 +922,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -926,7 +940,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -940,7 +956,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1034,6 +1050,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1052,7 +1069,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -1072,7 +1091,7 @@ public:
   expr_t getNonZeroPartofEquation() const;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &ar_params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRowHelper(expr_t arg1, expr_t arg2,
                                    int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
@@ -1163,6 +1182,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1179,7 +1199,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   bool isNumConstNodeEqualTo(double value) const override;
@@ -1193,7 +1215,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1291,6 +1313,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const override;
@@ -1305,7 +1328,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   expr_t decreaseLeadsLagsPredeterminedVariables() const override;
   expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
@@ -1321,7 +1346,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1496,6 +1521,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   void prepareForDerivation() override;
@@ -1516,7 +1542,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<tuple<int, expr_t, expr_t>>  &List_of_Op_RHS) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1537,7 +1565,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
@@ -1598,6 +1626,7 @@ public:
   int VarMinLag() const override;
   int VarMaxLag(DataTree &static_datatree, set<expr_t> &static_lhs) const override;
   int PacMaxLag(int lhs_symb_id) const override;
+  int getPacTargetSymbId(int lhs_symb_id, int undiff_lhs_symb_id) const override;
   expr_t undiff() const override;
   expr_t decreaseLeadsLags(int n) const override;
   void prepareForDerivation() override;
@@ -1618,7 +1647,9 @@ public:
   int findTargetVariable(int lhs_symb_id) const override;
   expr_t substituteDiff(DataTree &static_datatree, diff_table_t &diff_table, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
   expr_t substituteUnaryOpNodes(DataTree &static_datatree, diff_table_t &nodes, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const override;
-  expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
+  expr_t substitutePacExpectation(const string & name, int mce_symb_id,
+                                  map<const PacExpectationNode *, const BinaryOpNode *> &subst_table) override;
   pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<tuple<int, expr_t, expr_t>>  &List_of_Op_RHS) const override;
   void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1639,7 +1670,7 @@ public:
   expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const override;
   bool isInStaticForm() const override;
   void addParamInfoToPac(pair<int, int> &lhs_arg, int optim_share_arg, pair<int, pair<vector<int>, vector<bool>>> &ec_params_and_vars_arg, set<pair<int, pair<int, int>>> &params_and_vars_arg, const vector<tuple<int, int, int, double>> &non_optim_vars_params_and_constants) override;
-  void fillPacExpectationVarInfo(string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
+  void fillPacExpectationVarInfo(const string &model_name_arg, vector<int> &lhs_arg, int max_lag_arg, int pac_max_lag_arg, vector<bool> &nonstationary_arg, int growth_symb_id_arg, int growth_lag, int equation_number_arg) override;
   void fillAutoregressiveRow(int eqn, const vector<int> &lhs, map<tuple<int, int, int>, expr_t> &AR) const override;
   void fillErrorCorrectionRow(int eqn, const vector<int> &nontrend_lhs, const vector<int> &trend_lhs, map<tuple<int, int, int>, expr_t> &EC) const override;
   void findConstantEquations(map<VariableNode *, NumConstNode *> &table) const override;
diff --git a/src/ModFile.cc b/src/ModFile.cc
index 7d83522349ce0cd7badeab236e5fa2db5862ec25..b23fb1d2156efdab2c0b6e58275dc5ad98c88eb1 100644
--- a/src/ModFile.cc
+++ b/src/ModFile.cc
@@ -418,38 +418,48 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
       auto pms = dynamic_cast<PacModelStatement *>(statement.get());
       if (pms != nullptr)
          {
+           int max_lag;
            vector<int> lhs;
            vector<bool> nonstationary;
-           int growth_symb_id, growth_lag, max_lag;
-           string aux_model_name, pac_model_name;
-           tie ( pac_model_name, aux_model_name, growth_symb_id, growth_lag ) =
-             pms->getPacModelInfoForPacExpectation();
-           if (trend_component_model_table.isExistingTrendComponentModelName(aux_model_name))
+           if (trend_component_model_table.isExistingTrendComponentModelName(pms->aux_model_name))
              {
-               max_lag = trend_component_model_table.getMaxLag(aux_model_name) + 1;
-               lhs = dynamic_model.getUndiffLHSForPac(aux_model_name, diff_subst_table);
+               max_lag = trend_component_model_table.getMaxLag(pms->aux_model_name) + 1;
+               lhs = dynamic_model.getUndiffLHSForPac(pms->aux_model_name, diff_subst_table);
                // All lhs variables in a trend component model are nonstationary
-               nonstationary.insert(nonstationary.end(), trend_component_model_table.getDiff(aux_model_name).size(), true);
+               nonstationary.insert(nonstationary.end(), trend_component_model_table.getDiff(pms->aux_model_name).size(), true);
              }
-           else if (var_model_table.isExistingVarModelName(aux_model_name))
+           else if (var_model_table.isExistingVarModelName(pms->aux_model_name))
              {
-               max_lag = var_model_table.getMaxLag(aux_model_name);
-               lhs = var_model_table.getLhs(aux_model_name);
+               max_lag = var_model_table.getMaxLag(pms->aux_model_name);
+               lhs = var_model_table.getLhs(pms->aux_model_name);
                // nonstationary variables in a VAR are those that are in diff
-               nonstationary = var_model_table.getDiff(aux_model_name);
+               nonstationary = var_model_table.getDiff(pms->aux_model_name);
              }
+           else if (pms->aux_model_name == "")
+             max_lag = 0;
            else
              {
-               cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model"
-                    << endl;
+               cerr << "Error: aux_model_name not recognized as VAR model or Trend Component model" << endl;
                exit(EXIT_FAILURE);
              }
            pms->fillUndiffedLHS(lhs);
            dynamic_model.walkPacParameters();
-           int pac_max_lag = original_model.getPacMaxLag(pac_model_name);
-           dynamic_model.fillPacExpectationVarInfo(pac_model_name, lhs, max_lag,
-                                                   pac_max_lag, nonstationary, growth_symb_id, growth_lag);
-           dynamic_model.substitutePacExpectation();
+           int pac_max_lag_m = original_model.getPacMaxLag(pms->name);
+           if (pms->aux_model_name == "")
+             {
+               int pac_target_symb_id = dynamic_model.getPacTargetSymbId(pms->name);
+               int model_consistent_expectation_symb_id =
+                 dynamic_model.addPacModelConsistentExpectationEquation(pms->name, pac_target_symb_id,
+                                                                        symbol_table.getID(pms->discount), pac_max_lag_m,
+                                                                        diff_subst_table);
+               dynamic_model.substitutePacExpectation(pms->name, model_consistent_expectation_symb_id);
+             }
+           else
+             {
+               dynamic_model.fillPacExpectationVarInfo(pms->name, lhs, max_lag,
+                                                       pac_max_lag_m, nonstationary, pms->growth_symb_id, pms->growth_lag);
+               dynamic_model.substitutePacExpectation(pms->name);
+             }
          }
      }
 
diff --git a/src/ParsingDriver.cc b/src/ParsingDriver.cc
index a53f1ffdc1f49b1869d33c566f004ca5100ced09..acc646a7ef574937a9a1976e13f71de07c438020 100644
--- a/src/ParsingDriver.cc
+++ b/src/ParsingDriver.cc
@@ -2615,10 +2615,27 @@ ParsingDriver::pac_model()
     error("You must pass the model_name option to the pac_model statement.");
   auto name = it->second;
 
+  if (pac_growth_symb_id < 0)
+    error("You must pass the growth option to the pac_model statement");
+
+  string aux_model_name = "";
   it = options_list.string_options.find("pac.aux_model_name");
-  if (it == options_list.string_options.end())
-    error("You must pass the auxiliary_model_name option to the pac_model statement.");
-  auto aux_model_name = it->second;
+  if (it != options_list.string_options.end())
+    if (pac_steady_state_growth_rate_number >= 0 || pac_steady_state_growth_rate_symb_id >= 0)
+      {
+        pac_steady_state_growth_rate_number = -1;
+        pac_steady_state_growth_rate_symb_id = -1;
+        warning("when aux_model_name is used in the pac_model statement, steady_state_growth is ignored");
+      }
+    else
+      aux_model_name = it->second;
+  else
+    if (pac_steady_state_growth_rate_number < 0 && pac_steady_state_growth_rate_symb_id < 0)
+      error("when aux_model_name is not passed to the pac_model statement, you must pass steady_state_growth_rate");
+    else if (mod_file->symbol_table.getType(pac_growth_symb_id) == SymbolType::parameter
+             && (pac_steady_state_growth_rate_number >= 0
+                 || pac_steady_state_growth_rate_symb_id != pac_growth_symb_id))
+      error("when aux_model_name is not passed to the pac_model statement, steady_state_growth must be a parameter equal to growth");
 
   it = options_list.string_options.find("pac.discount");
   if (it == options_list.string_options.end())
@@ -2627,9 +2644,13 @@ ParsingDriver::pac_model()
 
   mod_file->addStatement(make_unique<PacModelStatement>(name, aux_model_name, discount,
                                                         pac_growth_symb_id, pac_growth_lag,
+                                                        pac_steady_state_growth_rate_number,
+                                                        pac_steady_state_growth_rate_symb_id,
                                                         mod_file->symbol_table));
   pac_growth_symb_id = -1;
   pac_growth_lag = 0;
+  pac_steady_state_growth_rate_number = -1;
+  pac_steady_state_growth_rate_symb_id = -1;
   options_list.clear();
 }
 
@@ -2642,6 +2663,21 @@ ParsingDriver::set_pac_growth(const string &name, int lag)
   pac_growth_lag = -lag;
 }
 
+void
+ParsingDriver::set_pac_steady_state_growth(const string &name_or_number)
+{
+  try
+    {
+      pac_steady_state_growth_rate_number = stod(name_or_number);
+    }
+  catch (...)
+    {
+      if (!mod_file->symbol_table.exists(name_or_number))
+        error("Unknown symbol used in pac_steady_state_growth option: " + name_or_number + "\n");
+      pac_steady_state_growth_rate_symb_id = mod_file->symbol_table.getID(name_or_number);
+    }
+}
+
 expr_t
 ParsingDriver::add_exp(expr_t arg1)
 {
diff --git a/src/ParsingDriver.hh b/src/ParsingDriver.hh
index a806eb53ee8bedf626a6f5acf6fa728499d38188..e0fa48ae8df6fc2b146da7e5c13cea17ff30cd87 100644
--- a/src/ParsingDriver.hh
+++ b/src/ParsingDriver.hh
@@ -257,6 +257,8 @@ private:
   //! Temporary storage for growth declared in pac_model
   int pac_growth_symb_id = -1;
   int pac_growth_lag = 0;
+  double pac_steady_state_growth_rate_number = -1;
+  int pac_steady_state_growth_rate_symb_id = -1;
 
   bool nostrict;
 
@@ -733,6 +735,8 @@ public:
   void pac_model();
   //! Adds growth for pac
   void set_pac_growth(const string &name, int lag = 0);
+  //! Adds steady state growth for pac
+  void set_pac_steady_state_growth(const string &name_or_number);
   //! Writes token "diff(arg1)" to model tree
   expr_t add_diff(expr_t arg1);
   //! Writes token "adl(arg1, lag)" to model tree
diff --git a/src/SymbolTable.cc b/src/SymbolTable.cc
index e0338ba8d60d6a73daa19733fdd9e6cad58d5648..aa9f5beef22362ae2626a8d1896e4bce930401ea 100644
--- a/src/SymbolTable.cc
+++ b/src/SymbolTable.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2018 Dynare Team
+ * Copyright (C) 2003-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -379,6 +379,7 @@ SymbolTable::writeOutput(ostream &output) const noexcept(false)
             break;
           case AuxVarType::diff:
           case AuxVarType::diffLag:
+          case AuxVarType::diffLead:
             if (aux_vars[i].get_orig_symb_id() >= 0)
               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;
@@ -496,6 +497,7 @@ SymbolTable::writeCOutput(ostream &output) const noexcept(false)
               break;
             case AuxVarType::diff:
             case AuxVarType::diffLag:
+            case AuxVarType::diffLead:
               if (aux_vars[i].get_orig_symb_id() >= 0)
                 output << "av[" << i << "].orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                        << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -601,6 +603,7 @@ SymbolTable::writeCCOutput(ostream &output) const noexcept(false)
           break;
         case AuxVarType::diff:
         case AuxVarType::diffLag:
+        case AuxVarType::diffLead:
           if (aux_vars[i].get_orig_symb_id() >= 0)
             output << "av" << i << ".orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                    << "av" << i << ".orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
@@ -741,6 +744,29 @@ SymbolTable::addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id
   return symb_id;
 }
 
+int
+SymbolTable::addDiffLeadAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lead) noexcept(false)
+{
+  ostringstream varname;
+  int symb_id;
+
+  varname << "AUX_DIFF_LEAD_" << index;
+
+  try
+    {
+      symb_id = addSymbol(varname.str(), SymbolType::endogenous);
+    }
+  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.emplace_back(symb_id, AuxVarType::diffLead, orig_symb_id, orig_lead, 0, 0, expr_arg, "");
+
+  return symb_id;
+}
+
 int
 SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) noexcept(false)
 {
@@ -873,7 +899,8 @@ SymbolTable::getOrigSymbIdForAuxVar(int aux_var_symb_id) const noexcept(false)
     if ((aux_var.get_type() == AuxVarType::endoLag
          || aux_var.get_type() == AuxVarType::exoLag
          || aux_var.get_type() == AuxVarType::diff
-         || aux_var.get_type() == AuxVarType::diffLag)
+         || aux_var.get_type() == AuxVarType::diffLag
+         || aux_var.get_type() == AuxVarType::diffLead)
         && aux_var.get_symb_id() == aux_var_symb_id)
       return aux_var.get_orig_symb_id();
   throw UnknownSymbolIDException(aux_var_symb_id);
@@ -884,7 +911,7 @@ SymbolTable::getOrigLeadLagForDiffAuxVar(int diff_aux_var_symb_id) const noexcep
 {
   int lag = 0;
   for (const auto & aux_var : aux_vars)
-    if (aux_var.get_type() == AuxVarType::diffLag
+    if ((aux_var.get_type() == AuxVarType::diffLag || aux_var.get_type() == AuxVarType::diffLead)
         && aux_var.get_symb_id() == diff_aux_var_symb_id)
       lag += 1 + getOrigLeadLagForDiffAuxVar(aux_var.get_orig_symb_id());
   return lag;
@@ -898,7 +925,7 @@ SymbolTable::getOrigSymbIdForDiffAuxVar(int diff_aux_var_symb_id) const noexcept
     if (aux_var.get_symb_id() == diff_aux_var_symb_id)
       if (aux_var.get_type() == AuxVarType::diff)
         orig_symb_id = diff_aux_var_symb_id;
-      else if (aux_var.get_type() == AuxVarType::diffLag)
+      else if (aux_var.get_type() == AuxVarType::diffLag || aux_var.get_type() == AuxVarType::diffLead)
         orig_symb_id = getOrigSymbIdForDiffAuxVar(aux_var.get_orig_symb_id());
   return orig_symb_id;
 }
@@ -1154,6 +1181,7 @@ SymbolTable::writeJuliaOutput(ostream &output) const noexcept(false)
               break;
             case AuxVarType::diff:
             case AuxVarType::diffLag:
+            case AuxVarType::diffLead:
               if (aux_var.get_orig_symb_id() >= 0)
                 output << getTypeSpecificID(aux_var.get_orig_symb_id()) + 1 << ", "
                        << aux_var.get_orig_lead_lag() << ", typemin(Int), string(), string()";
diff --git a/src/SymbolTable.hh b/src/SymbolTable.hh
index 49027f48a1124df16e0bd879a1188f53c7eab9db..bde5e91685ffdfb5dd935eae1c90b4d11f28a87e 100644
--- a/src/SymbolTable.hh
+++ b/src/SymbolTable.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2018 Dynare Team
+ * Copyright (C) 2003-2019 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -47,7 +47,8 @@ enum class AuxVarType
     varModel = 7,    //!< Variable for var_model with order > abs(min_lag()) present in model
     diff = 8,        //!< Variable for Diff operator
     diffLag = 9,     //!< Variable for timing between Diff operators
-    unaryOp = 10     //!< Variable for allowing the undiff operator to work when diff was taken of unary op, eg diff(log(x))
+    diffLead = 10,   //!< Variable for timing between Diff operators
+    unaryOp = 11     //!< Variable for allowing the undiff operator to work when diff was taken of unary op, eg diff(log(x))
   };
 
 //! Information on some auxiliary variables
@@ -309,6 +310,8 @@ public:
   int addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) noexcept(false);
   //! Takes care of timing between diff statements
   int addDiffLagAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) noexcept(false);
+  //! Takes care of timing between diff statements
+  int addDiffLeadAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lead) noexcept(false);
   //! An Auxiliary variable for a unary op
   int addUnaryOpAuxiliaryVar(int index, expr_t expr_arg, string unary_op, int orig_symb_id = -1, int orig_lag = 0) noexcept(false);
   //! Returns the number of auxiliary variables