diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index e2f60fc99a667b80b4c8a4c462555f6a8702e944..1a70427b899881087cd1247a0317abed5247b7d7 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -260,38 +260,97 @@ PriorPosteriorFunctionStatement::writeJsonOutput(ostream &output) const
 
 VarModelStatement::VarModelStatement(const SymbolList &symbol_list_arg,
                                      const OptionsList &options_list_arg,
-                                     const string &name_arg) :
+                                     const string &name_arg,
+                                     const SymbolTable &symbol_table_arg) :
   symbol_list(symbol_list_arg),
   options_list(options_list_arg),
-  name(name_arg)
+  name(name_arg),
+  symbol_table(symbol_table_arg)
 {
 }
 
 void
-VarModelStatement::getVarModelNameAndVarList(map<string, pair<SymbolList, int> > &var_model_info)
+VarModelStatement::getVarModelInfoForVarExpectation(map<string, pair<SymbolList, int> > &var_model_info) const
 {
+  if (symbol_list.empty())
+    return;
+
   OptionsList::num_options_t::const_iterator it = options_list.num_options.find("var.order");
-  if (it != options_list.num_options.end())
-    var_model_info[name] = make_pair(symbol_list, atoi(it->second.c_str()));
+  var_model_info[name] = make_pair(symbol_list, atoi(it->second.c_str()));
 }
 
 void
-VarModelStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
+VarModelStatement::getVarModelEqTags(vector<string> &var_model_eqtags) const
 {
-  OptionsList::vec_str_options_t::const_iterator itvs = options_list.vector_str_options.find("var.eqtags");
-  OptionsList::num_options_t::const_iterator it = options_list.num_options.find("var.order");
-  if ((it == options_list.num_options.end() && itvs == options_list.vector_str_options.end())
-      || (it != options_list.num_options.end() && itvs != options_list.vector_str_options.end()))
-    {
-      cerr << "ERROR: You must provide either the order or eqtags option to the var_model statement, but not both." << endl;
-      exit(EXIT_FAILURE);
-    }
+  if (!symbol_list.empty())
+    return;
 
-  if (name.empty())
-    {
-      cerr << "ERROR: You must provide the model_name option to the var_model statement." << endl;
-      exit(EXIT_FAILURE);
-    }
+  OptionsList::vec_str_options_t::const_iterator it1 =
+    options_list.vector_str_options.find("var.eqtags");
+  var_model_eqtags = it1->second;
+}
+
+void
+VarModelStatement::fillVarModelInfoFromEquations(vector<int> &eqnumber_arg, vector<int> &lhs_arg,
+                                                 vector<set<pair<int, int> > > &rhs_arg,
+                                                 vector<bool> &nonstationary_arg,
+                                                 vector<bool> &diff_arg, vector<int> &orig_diff_var_arg)
+{
+  eqnumber = eqnumber_arg;
+  lhs = lhs_arg;
+  rhs_by_eq = rhs_arg;
+  nonstationary = nonstationary_arg;
+  diff = diff_arg;
+  orig_diff_var = orig_diff_var_arg;
+
+  // Order RHS vars by time (already ordered by equation tag)
+  for (vector<set<pair<int, int> > >::const_iterator it = rhs_by_eq.begin();
+       it != rhs_by_eq.end(); it++)
+    for (set<pair<int, int> >::const_iterator it1 = it->begin();
+         it1 != it->end(); it1++)
+      if (find(lhs.begin(), lhs.end(), it1->first) == lhs.end()
+          && find(orig_diff_var.begin(), orig_diff_var.end(), it1->first) == orig_diff_var.end())
+        {/*
+          cerr << "ERROR " << name << ": " << symbol_table.getName(it1->first)
+               << " cannot appear in the VAR because it does not appear on the LHS" << endl;
+          exit(EXIT_FAILURE);
+         */
+        }
+      else
+        {
+          map<int, set<int> >::iterator mit = rhs.find(abs(it1->second));
+          if (mit == rhs.end())
+            {
+              if (it1->second > 0)
+                {
+                  cerr << "ERROR " << name << ": you cannot have a variable with a lead in a VAR" << endl;
+                  exit(EXIT_FAILURE);
+                }
+              set<int> si;
+              si.insert(it1->first);
+              rhs[abs(it1->second)] = si;
+            }
+          else
+            mit->second.insert(it1->first);
+        }
+}
+
+void
+VarModelStatement::getVarModelName(string &var_model_name) const
+{
+  var_model_name = name;
+}
+
+
+void
+VarModelStatement::getVarModelRHS(map<int, set<int > > &rhs_arg) const
+{
+  rhs_arg = rhs;
+}
+
+void
+VarModelStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
+{
 }
 
 void
@@ -300,6 +359,108 @@ VarModelStatement::writeOutput(ostream &output, const string &basename, bool min
   options_list.writeOutput(output);
   if (!symbol_list.empty())
     symbol_list.writeOutput("options_.var.var_list_", output);
+
+  output << "options_.var.eqn = [";
+  for (vector<int>::const_iterator it = eqnumber.begin();
+       it != eqnumber.end(); it++)
+    {
+      if (it != eqnumber.begin())
+        output << " ";
+      output << *it + 1;
+    }
+  output << "];" << endl
+         << "options_.var.lhs = [";
+  for (vector<int>::const_iterator it = lhs.begin();
+       it != lhs.end(); it++)
+    {
+      if (it != lhs.begin())
+        output << " ";
+      output << symbol_table.getTypeSpecificID(*it) + 1;
+    }
+  output << "];" << endl
+         << "options_.var.rhs.lag = [";
+  for (map<int, set<int> >::const_iterator it = rhs.begin();
+       it != rhs.end(); it++)
+    {
+      if (it != rhs.begin())
+        output << " ";
+      output << it->first;
+    }
+  output << "];" << endl;
+  int i = 1;
+  for (map<int, set<int> >::const_iterator it = rhs.begin();
+       it != rhs.end(); it++, i++)
+    {
+      output << "options_.var.rhs.vars_at_lag{" << i << "} = [";
+      for (set<int>::const_iterator it1 = it->second.begin();
+           it1 != it->second.end(); it1++)
+        {
+          if (it1 != it->second.begin())
+            output << " ";
+          output << symbol_table.getTypeSpecificID(*it1) + 1;
+        }
+      output << "];" << endl;
+    }
+  output << "options_.var.nonstationary = logical([";
+  for (vector<bool>::const_iterator it = nonstationary.begin();
+       it != nonstationary.end(); it++)
+    {
+      if (it != nonstationary.begin())
+        output << " ";
+      if (*it)
+        output << "1";
+      else
+        output << "0";
+    }
+  output << "]);" << endl
+         << "options_.var.diff = logical([";
+  for (vector<bool>::const_iterator it = diff.begin();
+       it != diff.end(); it++)
+    {
+      if (it != diff.begin())
+        output << " ";
+      if (*it)
+        output << "1";
+      else
+        output << "0";
+    }
+  output << "]);" << endl
+         << "options_.var.orig_diff_var = [";
+  for (vector<int>::const_iterator it = orig_diff_var.begin();
+       it != orig_diff_var.end(); it++)
+    {
+      if (it != orig_diff_var.begin())
+        output << " ";
+      if (*it == -1)
+        output << -1;
+      else
+        output << symbol_table.getTypeSpecificID(*it) + 1;
+    }
+  output << "];" << endl;
+  i = 1;
+  for (vector<set<pair<int, int > > >::const_iterator it = rhs_by_eq.begin();
+       it != rhs_by_eq.end(); it++, i++)
+    {
+      output << "options_.var.rhs.vars_at_eq{" << i << "}.var = [";
+      for (set<pair<int, int> >::const_iterator it1 = it->begin();
+           it1 != it->end(); it1++)
+        {
+          if (it1 != it->begin())
+            output << " ";
+          output << symbol_table.getTypeSpecificID(it1->first) + 1;
+        }
+      output << "];" << endl
+             << "options_.var.rhs.vars_at_eq{" << i << "}.lag = [";
+      for (set<pair<int, int> >::const_iterator it1 = it->begin();
+           it1 != it->end(); it1++)
+        {
+          if (it1 != it->begin())
+            output << " ";
+          output << it1->second;
+        }
+      output << "];" << endl;
+
+    }
   output << "M_.var." << name << " = options_.var;" << endl
          << "clear options_.var;" << endl;
 }
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index 8432be9fa28f6dc70168b503633b4405b56f71ba..cb7052d2710656e28ae0eb85e111d698e8c7de2b 100644
--- a/ComputingTasks.hh
+++ b/ComputingTasks.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -125,11 +125,24 @@ private:
   const SymbolList symbol_list;
   const OptionsList options_list;
   const string &name;
+  const SymbolTable &symbol_table;
+  vector<int> eqnumber, lhs, orig_diff_var;
+  map<int, set<int > > rhs; // lag -> set< symb_id > (all vars that appear at a given lag)
+  vector<set<pair<int, int> > > rhs_by_eq; // rhs by equation
+  vector<bool> nonstationary, diff;
 public:
   VarModelStatement(const SymbolList &symbol_list_arg,
                     const OptionsList &options_list_arg,
-                    const string &name_arg);
-  void getVarModelNameAndVarList(map<string, pair<SymbolList, int> > &var_model_info);
+                    const string &name_arg,
+                    const SymbolTable &symbol_table_arg);
+  void getVarModelInfoForVarExpectation(map<string, pair<SymbolList, int> > &var_model_info) const;
+  void getVarModelEqTags(vector<string> &var_model_eqtags) const;
+  void fillVarModelInfoFromEquations(vector<int> &eqnumber_arg, vector<int> &lhs_arg,
+                                     vector<set<pair<int, int> > > &rhs_arg,
+                                     vector<bool> &nonstationary_arg,
+                                     vector<bool> &diff_arg, vector<int> &orig_diff_var_arg);
+  void getVarModelName(string &var_model_name) const;
+  void getVarModelRHS(map<int, set<int > > &rhs_arg) const;
   virtual void checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings);
   virtual void writeOutput(ostream &output, const string &basename, bool minimal_workspace) const;
   void createVarModelMFunction(ostream &output, const map<string, set<int> > &var_expectation_functions_to_write) const;
diff --git a/DataTree.cc b/DataTree.cc
index 50f668ae6babbea79cea4a8998c27339bc695333..493137224671443cf3ff7048a72c2e65dfd422f7 100644
--- a/DataTree.cc
+++ b/DataTree.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -520,6 +520,17 @@ DataTree::AddVarExpectation(const int symb_id, const int forecast_horizon, const
   return new VarExpectationNode(*this, symb_id, forecast_horizon, model_name);
 }
 
+expr_t
+DataTree::AddPacExpectation(const string &model_name, const string &var_model_name, const int discount_id, const int growth_id)
+{
+  pac_expectation_node_map_t::iterator it =
+    pac_expectation_node_map.find(make_pair(model_name, make_pair(var_model_name, make_pair(discount_id, growth_id))));
+  if (it != pac_expectation_node_map.end())
+    return it->second;
+
+  return new PacExpectationNode(*this, model_name, var_model_name, discount_id, growth_id);
+}
+
 expr_t
 DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
 {
diff --git a/DataTree.hh b/DataTree.hh
index 3cda0bf2a7524583e5d8c3492fc07e838ee73d04..e55de76c5b9ddacb6a8c12c034ff9cf8ebd1e7e2 100644
--- a/DataTree.hh
+++ b/DataTree.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -49,6 +49,7 @@ class DataTree
   friend class FirstDerivExternalFunctionNode;
   friend class SecondDerivExternalFunctionNode;
   friend class VarExpectationNode;
+  friend class PacExpectationNode;
 protected:
   //! A reference to the symbol table
   SymbolTable &symbol_table;
@@ -80,6 +81,10 @@ protected:
   typedef map<pair<string, pair<int, int> >, VarExpectationNode *> var_expectation_node_map_t;
   var_expectation_node_map_t var_expectation_node_map;
 
+  // (model_name, (discount, growth)) -> PacExpectationNode
+  typedef map<pair<string, pair<string, pair<int, int> > >, PacExpectationNode *> pac_expectation_node_map_t;
+  pac_expectation_node_map_t pac_expectation_node_map;
+
   // ((arguments, deriv_idx), symb_id) -> FirstDerivExternalFunctionNode
   typedef map<pair<pair<vector<expr_t>, int>, int>, FirstDerivExternalFunctionNode *> first_deriv_external_function_node_map_t;
   first_deriv_external_function_node_map_t first_deriv_external_function_node_map;
@@ -226,6 +231,8 @@ public:
   expr_t AddEqual(expr_t iArg1, expr_t iArg2);
   //! Adds "var_expectation(arg1, arg2, model_name=arg3)" to model tree
   expr_t AddVarExpectation(const int symb_id, const int forecast_horizon, const string &model_name);
+  //! Adds pac_expectation command to model tree
+  expr_t AddPacExpectation(const string &model_name, const string &var_model_name, const int discount_id, const int growth_id);
   //! Adds a model local variable with its value
   void AddLocalVariable(int symb_id, expr_t value) throw (LocalVariableException);
   //! Adds an external function node
diff --git a/DynamicModel.cc b/DynamicModel.cc
index 054482adc0a247ee91002af316c8c028a9c4a23b..10024fc69f0932b67c03a101edbd884fd151625d 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3200,6 +3200,13 @@ DynamicModel::writeOutput(ostream &output, const string &basename, bool block_de
   else
     output << "-1";
   output << "];" << endl;
+
+  // Write PacExpectationInfo
+  deriv_node_temp_terms_t tef_terms;
+  temporary_terms_t temp_terms_empty;
+  for (set<const PacExpectationNode *>::const_iterator it = pac_expectation_info.begin();
+       it != pac_expectation_info.end(); it++)
+    (*it)->writeOutput(output, oMatlabDynamicModel, temp_terms_empty, tef_terms);
 }
 
 map<pair<int, pair<int, int > >, expr_t>
@@ -3228,14 +3235,103 @@ DynamicModel::runTrendTest(const eval_context_t &eval_context)
 }
 
 void
-DynamicModel::setVarExpectationIndices(map<string, pair<SymbolList, int> > var_model_info)
+DynamicModel::getVarModelVariablesFromEqTags(vector<string> &var_model_eqtags,
+                                             vector<int> &eqnumber,
+                                             vector<int> &lhs,
+                                             vector<set<pair<int, int> > > &rhs,
+                                             vector<bool> &nonstationary) const
+{
+  for (vector<string>::const_iterator itvareqs = var_model_eqtags.begin();
+       itvareqs != var_model_eqtags.end(); itvareqs++)
+    {
+      int eqnumber_int = -1;
+      set<pair<int, int> > lhs_set, lhs_tmp_set, rhs_set;
+      string eqtag (*itvareqs);
+      for (vector<pair<int, pair<string, string> > >::const_iterator iteqtag =
+             equation_tags.begin(); iteqtag != equation_tags.end(); iteqtag++)
+        if (iteqtag->second.first.compare("name") == 0
+            && iteqtag->second.second.compare(eqtag) == 0)
+          {
+            eqnumber_int = iteqtag->first;
+            break;
+          }
+
+      if (eqnumber_int == -1)
+        {
+          cerr << "ERROR: equation tag '" << eqtag << "' not found" << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      bool nonstationary_bool = false;
+      for (vector<pair<int, pair<string, string> > >::const_iterator iteqtag =
+             equation_tags.begin(); iteqtag != equation_tags.end(); iteqtag++)
+        if (iteqtag->first == eqnumber_int)
+          if (iteqtag->second.first.compare("data_type") == 0
+              && iteqtag->second.second.compare("nonstationary") == 0)
+            {
+              nonstationary_bool = true;
+              break;
+            }
+
+      equations[eqnumber_int]->get_arg1()->collectDynamicVariables(eEndogenous, lhs_set);
+      equations[eqnumber_int]->get_arg1()->collectDynamicVariables(eExogenous, lhs_tmp_set);
+      equations[eqnumber_int]->get_arg1()->collectDynamicVariables(eParameter, lhs_tmp_set);
+
+      if (lhs_set.size() != 1 || !lhs_tmp_set.empty())
+        {
+          cerr << "ERROR: A VAR may only have one endogenous variable on the LHS" << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      set<pair<int, int> >::const_iterator it = lhs_set.begin();
+      if (it->second != 0)
+        {
+          cerr << "ERROR: The variable on the LHS of a VAR may not appear with a lead or a lag" << endl;
+          exit(EXIT_FAILURE);
+        }
+
+      eqnumber.push_back(eqnumber_int);
+      lhs.push_back(it->first);
+      nonstationary.push_back(nonstationary_bool);
+
+      equations[eqnumber_int]->get_arg2()->collectDynamicVariables(eEndogenous, rhs_set);
+      rhs.push_back(rhs_set);
+    }
+}
+
+void
+DynamicModel::getDiffInfo(vector<int> &eqnumber, vector<bool> &diff, vector<int> &orig_diff_var) const
+{
+  for (vector<int>::const_iterator it = eqnumber.begin();
+       it != eqnumber.end(); it++)
+    {
+      diff.push_back(equations[*it]->get_arg1()->isDiffPresent());
+      if (diff.back())
+        {
+          set<pair<int, int> > diff_set;
+          equations[*it]->get_arg1()->collectDynamicVariables(eEndogenous, diff_set);
+          if (diff_set.empty() || diff_set.size() != 1)
+            {
+              cerr << "ERROR: problem getting variable for diff operator in equation " << *it << endl;
+              exit(EXIT_FAILURE);
+            }
+          set<pair<int, int> >::const_iterator it1 = diff_set.begin();
+          orig_diff_var.push_back(it1->first);
+        }
+      else
+        orig_diff_var.push_back(-1);
+    }
+}
+
+void
+DynamicModel::setVarExpectationIndices(map<string, pair<SymbolList, int> > &var_model_info)
 {
   for (size_t i = 0; i < equations.size(); i++)
     equations[i]->setVarExpectationIndex(var_model_info);
 }
 
 void
-DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > var_model_info)
+DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > &var_model_info)
 {
   // List of endogenous variables and the minimum lag value that must exist in the model equations
   map<string, int> var_endos_and_lags, model_endos_and_lags;
@@ -3286,6 +3382,50 @@ DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > var_model_i
     cout << "Accounting for var_model lags not in model block: added " << count << " auxiliary variables and equations." << endl;
 }
 
+void
+DynamicModel::walkPacParameters()
+{
+  for (size_t i = 0; i < equations.size(); i++)
+    {
+      bool pac_encountered = false;
+      pair<int, int> lhs (-1, -1);
+      set<pair<int, pair<int, int> > > params_and_vals;
+      equations[i]->walkPacParameters(pac_encountered, lhs, params_and_vals);
+      if (pac_encountered)
+        equations[i]->addParamInfoToPac(lhs, params_and_vals);
+    }
+}
+
+void
+DynamicModel::fillPacExpectationVarInfo(string &var_model_name,
+                                        vector<int> &lhs,
+                                        map<int, set<int > > &rhs,
+                                        vector<bool> &nonstationary)
+{
+  for (size_t i = 0; i < equations.size(); i++)
+    equations[i]->fillPacExpectationVarInfo(var_model_name, lhs, rhs, nonstationary, i);
+}
+
+void
+DynamicModel::substitutePacExpectation()
+{
+  map<const PacExpectationNode *, const BinaryOpNode *> subst_table;
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second = it->second->substitutePacExpectation(subst_table);
+
+  for (size_t i = 0; i < equations.size(); i++)
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substitutePacExpectation(subst_table));
+      assert(substeq != NULL);
+      equations[i] = 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::computingPass(bool jacobianExo, bool hessian, bool thirdDerivatives, int paramsDerivsOrder,
                             const eval_context_t &eval_context, bool no_tmp_terms, bool block, bool use_dll,
@@ -4750,10 +4890,37 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
 }
 
 void
-DynamicModel::substituteAdlAndDiff()
+DynamicModel::substituteAdl()
+{
+  for (int i = 0; i < (int) equations.size(); i++)
+    equations[i] = dynamic_cast<BinaryOpNode *>(equations[i]->substituteAdl());
+}
+
+void
+DynamicModel::substituteDiff()
 {
+  ExprNode::subst_table_t subst_table;
+  vector<BinaryOpNode *> neweqs;
+
+  // Substitute in model local variables
+  for (map<int, expr_t>::iterator it = local_variables_table.begin();
+       it != local_variables_table.end(); it++)
+    it->second = it->second->substituteDiff(subst_table, neweqs);
+
+  // Substitute in equations
   for (int i = 0; i < (int) equations.size(); i++)
-    equations[i] = dynamic_cast<BinaryOpNode *>(equations[i]->substituteAdlAndDiff());
+    {
+      BinaryOpNode *substeq = dynamic_cast<BinaryOpNode *>(equations[i]->substituteDiff(subst_table, neweqs));
+      assert(substeq != NULL);
+      equations[i] = substeq;
+    }
+
+  // Add new equations
+  for (int i = 0; i < (int) neweqs.size(); i++)
+    addEquation(neweqs[i], -1);
+
+  if (subst_table.size() > 0)
+    cout << "Substitution of Diff operator: added " << neweqs.size() << " auxiliary variables and equations." << endl;
 }
 
 void
diff --git a/DynamicModel.hh b/DynamicModel.hh
index edd506d25fe71f539de65cf2f1bf2849eb86971a..afcd72b9bb34ecc2aadb03360f2445ac85dba70e 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -226,6 +226,9 @@ private:
   //! Used for var_expectation and var_model
   map<string, set<int> > var_expectation_functions_to_write;
 
+  //! Used for pac_expectation operator
+  set<const PacExpectationNode *> pac_expectation_info; // PacExpectationNode pointers
+
   //!Maximum lead and lag for each block on endogenous of the block, endogenous of the previous blocks, exogenous and deterministic exogenous
   vector<pair<int, int> > endo_max_leadlag_block, other_endo_max_leadlag_block, exo_max_leadlag_block, exo_det_max_leadlag_block, max_leadlag_block;
 
@@ -280,10 +283,29 @@ public:
   //! Set the equations that have non-zero second derivatives
   void setNonZeroHessianEquations(map<int, string> &eqs);
 
+  //! Get equation info associated with equation tags from var_model
+  void getVarModelVariablesFromEqTags(vector<string> &var_model_eqtags,
+                                      vector<int> &eqnumber,
+                                      vector<int> &lhs,
+                                      vector<set<pair<int, int> > > &rhs,
+                                      vector<bool> &nonstationary) const;
+
+  // Get equtaino information on diff operator
+  void getDiffInfo(vector<int> &eqnumber, vector<bool> &diff, vector<int> &orig_diff_var) const;
+
   //! Set indices for var expectation in dynamic model file
-  void setVarExpectationIndices(map<string, pair<SymbolList, int> > var_model_info);
+  void setVarExpectationIndices(map<string, pair<SymbolList, int> > &var_model_info);
   //! Add aux equations (and aux variables) for variables declared in var_model at max order if they don't already exist
-  void addEquationsForVar(map<string, pair<SymbolList, int> > var_model_info);
+  void addEquationsForVar(map<string, pair<SymbolList, int> > &var_model_info);
+  //! Get Pac equation parameter info
+  void walkPacParameters();
+  //! Add var_model info to pac_expectation nodes
+  void fillPacExpectationVarInfo(string &var_model_name,
+                                 vector<int> &lhs,
+                                 map<int, set<int > > &rhs,
+                                 vector<bool> &nonstationary);
+  //! Substitutes pac_expectation operator
+  void substitutePacExpectation();
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename,
@@ -366,8 +388,11 @@ public:
   //! Transforms the model by removing trends specified by the user
   void detrendEquations();
 
-  //! Substitutes adl and diff operators
-  void substituteAdlAndDiff();
+  //! Substitutes adl operator
+  void substituteAdl();
+
+  //! Substitutes diff operator
+  void substituteDiff();
 
   //! Fill var_expectation_functions_to_write
   void fillVarExpectationFunctionsToWrite();
diff --git a/DynareBison.yy b/DynareBison.yy
index 53a9ea4267090ad0eaec7286882a3ee67fc7e472..97e59db1b9009e4a8b5f10e95109d6347701acc7 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -115,9 +115,9 @@ class ParsingDriver;
 %token NONLINEAR_FILTER_INITIALIZATION FILTER_ALGORITHM PROPOSAL_APPROXIMATION CUBATURE UNSCENTED MONTECARLO DISTRIBUTION_APPROXIMATION
 %token <string_val> NAME
 %token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN INIT_STATE FAST_REALTIME RESCALE_PREDICTION_ERROR_COVARIANCE GENERATE_IRFS
-%token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
+%token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY VAR_MODEL_NAME
 %token NOGRAPH POSTERIOR_NOGRAPH POSTERIOR_GRAPH NOMOMENTS NOPRINT NORMAL_PDF SAVE_DRAWS MODEL_NAME STDERR_MULTIPLES DIAGONAL_ONLY
-%token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE
+%token OBSERVATION_TRENDS OPTIM OPTIM_WEIGHTS ORDER OSR OSR_PARAMS MAX_DIM_COVA_GROUP ADVANCED OUTFILE OUTVARS OVERWRITE DISCOUNT
 %token PARALLEL_LOCAL_FILES PARAMETERS PARAMETER_SET PARTIAL_INFORMATION PERIODS PERIOD PLANNER_OBJECTIVE PLOT_CONDITIONAL_FORECAST PLOT_PRIORS PREFILTER PRESAMPLE
 %token PERFECT_FORESIGHT_SETUP PERFECT_FORESIGHT_SOLVER NO_POSTERIOR_KERNEL_DENSITY FUNCTION
 %token PRINT PRIOR_MC PRIOR_TRUNC PRIOR_MODE PRIOR_MEAN POSTERIOR_MODE POSTERIOR_MEAN POSTERIOR_MEDIAN MLE_MODE PRUNING
@@ -132,7 +132,7 @@ class ParsingDriver;
 %token UNIFORM_PDF UNIT_ROOT_VARS USE_DLL USEAUTOCORR GSA_SAMPLE_FILE USE_UNIVARIATE_FILTERS_IF_SINGULARITY_IS_DETECTED
 %token VALUES VAR VAREXO VAREXO_DET VAROBS VAREXOBS PREDETERMINED_VARIABLES VAR_EXPECTATION PLOT_SHOCK_DECOMPOSITION MODEL_LOCAL_VARIABLE
 %token WRITE_LATEX_DYNAMIC_MODEL WRITE_LATEX_STATIC_MODEL WRITE_LATEX_ORIGINAL_MODEL CROSSEQUATIONS COVARIANCE WRITE_LATEX_STEADY_STATE_MODEL
-%token XLS_SHEET XLS_RANGE LMMCP OCCBIN BANDPASS_FILTER COLORMAP VAR_MODEL QOQ YOY AOA
+%token XLS_SHEET XLS_RANGE LMMCP OCCBIN BANDPASS_FILTER COLORMAP VAR_MODEL QOQ YOY AOA PAC_EXPECTATION
 %left COMMA
 %left EQUAL_EQUAL EXCLAMATION_EQUAL
 %left LESS GREATER LESS_EQUAL GREATER_EQUAL
@@ -158,7 +158,7 @@ class ParsingDriver;
 %token INDXESTIMA INDXGDLS EQ_MS FILTER_COVARIANCE FILTER_DECOMPOSITION SMOOTHED_STATE_UNCERTAINTY
 %token EQ_CMS TLINDX TLNUMBER BANACT RESTRICTIONS POSTERIOR_SAMPLER_OPTIONS
 %token OUTPUT_FILE_TAG DRAWS_NBR_BURN_IN_1 DRAWS_NBR_BURN_IN_2 HORIZON
-%token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION
+%token SBVAR TREND_VAR DEFLATOR GROWTH_FACTOR MS_IRF MS_VARIANCE_DECOMPOSITION GROWTH
 %token MS_ESTIMATION MS_SIMULATION MS_COMPUTE_MDD MS_COMPUTE_PROBABILITIES MS_FORECAST
 %token SVAR_IDENTIFICATION EQUATION EXCLUSION LAG UPPER_CHOLESKY LOWER_CHOLESKY MONTHLY QUARTERLY
 %token MARKOV_SWITCHING CHAIN DURATION NUMBER_OF_REGIMES NUMBER_OF_LAGS
@@ -908,6 +908,8 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_var_expectation($3, new string("1"), $7); }
           | VAR_EXPECTATION '(' symbol COMMA INT_NUMBER COMMA MODEL_NAME EQUAL NAME ')'
             { $$ = driver.add_var_expectation($3, $5, $9); }
+          | PAC_EXPECTATION '(' pac_expectation_options_list ')'
+            { $$ = driver.add_pac_expectation(); }
           | MINUS hand_side %prec UMINUS
             { $$ = driver.add_uminus($2); }
           | PLUS hand_side
@@ -966,6 +968,21 @@ hand_side : '(' hand_side ')'
             { $$ = driver.add_steady_state($3); }
           ;
 
+
+pac_expectation_options_list : pac_expectation_options_list COMMA pac_expectation_options
+                             | pac_expectation_options
+                             ;
+
+pac_expectation_options : VAR_MODEL_NAME EQUAL NAME
+                          { driver.add_pac_expectation_var_model_name($3); }
+                        | MODEL_NAME EQUAL NAME
+                          { driver.add_pac_expectation_model_name($3); }
+                        | DISCOUNT EQUAL symbol
+                          { driver.add_pac_expectation_discount($3); }
+                        | GROWTH EQUAL symbol
+                          { driver.add_pac_expectation_growth($3); }
+                        ;
+
 comma_hand_side : hand_side
                   { driver.add_external_function_arg($1); }
                 | comma_hand_side COMMA hand_side
diff --git a/DynareFlex.ll b/DynareFlex.ll
index 86f0075685e319d01121a6a3d125cf598f458067..51a2b48a7fcd881845aa09323ddd19cbbc7036a9 100644
--- a/DynareFlex.ll
+++ b/DynareFlex.ll
@@ -342,6 +342,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_STATEMENT>optim		{return token::OPTIM;}
 <DYNARE_STATEMENT>periods	{return token::PERIODS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>model_name	{return token::MODEL_NAME;}
+<DYNARE_BLOCK>var_model_name    {return token::VAR_MODEL_NAME;}
 <DYNARE_STATEMENT>endogenous_terminal_period 	{return token::ENDOGENOUS_TERMINAL_PERIOD;}
 <DYNARE_STATEMENT>sub_draws	{return token::SUB_DRAWS;}
 <DYNARE_STATEMENT>minimal_solving_periods {return token::MINIMAL_SOLVING_PERIODS;}
@@ -677,6 +678,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
 }
 
  /* Inside a Dynare block */
+<DYNARE_BLOCK>growth {return token::GROWTH;}
 <DYNARE_BLOCK>var {return token::VAR;}
 <DYNARE_BLOCK>stderr {return token::STDERR;}
 <DYNARE_BLOCK>values {return token::VALUES;}
@@ -822,6 +824,8 @@ 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_STATEMENT,DYNARE_BLOCK>steady_state {return token::STEADY_STATE;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>expectation {return token::EXPECTATION;}
 <DYNARE_BLOCK>var_expectation {return token::VAR_EXPECTATION;}
+<DYNARE_BLOCK>pac_expectation {return token::PAC_EXPECTATION;}
+<DYNARE_BLOCK>discount {return token::DISCOUNT;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varobs {return token::VAROBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>varexobs {return token::VAREXOBS;}
 <DYNARE_STATEMENT,DYNARE_BLOCK>full {return token::FULL;}
diff --git a/ExprNode.cc b/ExprNode.cc
index fbdf400da0b4305f5bb443fd7ceeba96d9f6637e..687b6bde75aa68a2c39243b806672dae8d52ae45 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2017 Dynare Team
+ * Copyright (C) 2007-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -294,6 +294,12 @@ ExprNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_ar
   return false;
 }
 
+bool
+ExprNode::isDiffPresent() const
+{
+  return false;
+}
+
 void
 ExprNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
 {
@@ -307,6 +313,12 @@ NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
   datatree.num_const_node_map[id] = this;
 }
 
+bool
+NumConstNode::isDiffPresent() const
+{
+  return false;
+}
+
 void
 NumConstNode::prepareForDerivation()
 {
@@ -482,7 +494,19 @@ NumConstNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 }
 
 expr_t
-NumConstNode::substituteAdlAndDiff() const
+NumConstNode::substituteAdl() const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
+expr_t
+NumConstNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<NumConstNode *>(this);
+}
+
+expr_t
+NumConstNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<NumConstNode *>(this);
 }
@@ -554,6 +578,21 @@ NumConstNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
 {
 }
 
+void
+NumConstNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+}
+
+void
+NumConstNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+}
+
+void
+NumConstNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+}
+
 bool
 NumConstNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -1244,7 +1283,19 @@ VariableNode::maxLead() const
 }
 
 expr_t
-VariableNode::substituteAdlAndDiff() const
+VariableNode::substituteAdl() const
+{
+  return const_cast<VariableNode *>(this);
+}
+
+expr_t
+VariableNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<VariableNode *>(this);
+}
+
+expr_t
+VariableNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<VariableNode *>(this);
 }
@@ -1531,6 +1582,12 @@ VariableNode::detrend(int symb_id, bool log_trend, expr_t trend) const
     }
 }
 
+bool
+VariableNode::isDiffPresent() const
+{
+  return false;
+}
+
 expr_t
 VariableNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
 {
@@ -1587,6 +1644,21 @@ VariableNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
 {
 }
 
+void
+VariableNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+}
+
+void
+VariableNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+}
+
+void
+VariableNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+}
+
 bool
 VariableNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -2709,23 +2781,16 @@ UnaryOpNode::maxLead() const
 }
 
 expr_t
-UnaryOpNode::substituteAdlAndDiff() const
+UnaryOpNode::substituteAdl() const
 {
-  if (op_code != oDiff && op_code != oAdl)
+  if (op_code != oAdl)
     {
-      expr_t argsubst = arg->substituteAdlAndDiff();
+      expr_t argsubst = arg->substituteAdl();
       return buildSimilarUnaryOpNode(argsubst, datatree);
     }
 
-  if (op_code == oDiff)
-    {
-      expr_t argsubst = arg->substituteAdlAndDiff();
-      return datatree.AddMinus(argsubst,
-                               argsubst->decreaseLeadsLags(1));
-    }
-
-  expr_t arg1subst = arg->substituteAdlAndDiff();
-  expr_t retval;
+  expr_t arg1subst = arg->substituteAdl();
+  expr_t retval = NULL;
   ostringstream inttostr;
   if (adl_param_lag >= 0)
     {
@@ -2765,6 +2830,55 @@ UnaryOpNode::substituteAdlAndDiff() const
   return retval;
 }
 
+bool
+UnaryOpNode::isDiffPresent() const
+{
+  if (op_code == oDiff)
+    return true;
+  return arg->isDiffPresent();
+}
+
+expr_t
+UnaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  if (op_code != oDiff)
+    {
+      expr_t argsubst = arg->substituteDiff(subst_table, neweqs);
+      return buildSimilarUnaryOpNode(argsubst, datatree);
+    }
+
+  subst_table_t::iterator it = subst_table.find(const_cast<UnaryOpNode *>(this));
+  if (it != subst_table.end())
+    return const_cast<VariableNode *>(it->second);
+
+  expr_t argsubst = arg->substituteDiff(subst_table, neweqs);
+  assert(argsubst != NULL);
+
+  int symb_id = -1;
+  VariableNode *vn = dynamic_cast<VariableNode *>(argsubst);
+  if (vn != NULL)
+    symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst, vn->get_symb_id(), vn->get_lag());
+  else
+    symb_id = datatree.symbol_table.addDiffAuxiliaryVar(argsubst->idx, argsubst);
+
+  expr_t newAuxVar = datatree.AddVariable(symb_id, 0);
+
+  //AUX_DIFF_IDX = argsubst - argsubst(-1)
+  neweqs.push_back(dynamic_cast<BinaryOpNode *>(datatree.AddEqual(newAuxVar,
+                                                                  datatree.AddMinus(argsubst, argsubst->decreaseLeadsLags(1)))));
+
+  assert(dynamic_cast<VariableNode *>(newAuxVar) != NULL);
+  subst_table[this] = dynamic_cast<VariableNode *>(newAuxVar);
+  return newAuxVar;
+}
+
+expr_t
+UnaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  expr_t argsubst = arg->substitutePacExpectation(subst_table);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 expr_t
 UnaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -2937,6 +3051,24 @@ UnaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mod
   arg->setVarExpectationIndex(var_model_info);
 }
 
+void
+UnaryOpNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+  arg->walkPacParameters(pac_encountered, lhs, params_and_vals);
+}
+
+void
+UnaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+  arg->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+}
+
+void
+UnaryOpNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+  arg->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+}
+
 bool
 UnaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -4361,10 +4493,32 @@ BinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpN
 }
 
 expr_t
-BinaryOpNode::substituteAdlAndDiff() const
+BinaryOpNode::substituteAdl() const
+{
+  expr_t arg1subst = arg1->substituteAdl();
+  expr_t arg2subst = arg2->substituteAdl();
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
+expr_t
+BinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteDiff(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(subst_table, neweqs);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
+bool
+BinaryOpNode::isDiffPresent() const
+{
+  return arg1->isDiffPresent() || arg2->isDiffPresent();
+}
+
+expr_t
+BinaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t arg1subst = arg1->substituteAdlAndDiff();
-  expr_t arg2subst = arg2->substituteAdlAndDiff();
+  expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
@@ -4445,6 +4599,91 @@ BinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
   arg2->setVarExpectationIndex(var_model_info);
 }
 
+void
+BinaryOpNode::walkPacParametersHelper(const expr_t arg1, const expr_t arg2,
+                                      pair<int, int> &lhs,
+                                      set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+  set<int> params;
+  set<pair<int, int> > endogs;
+  arg1->collectVariables(eParameter, params);
+  arg2->collectDynamicVariables(eEndogenous, endogs);
+  if (params.size() == 1)
+    if (endogs.size() == 1)
+      params_and_vals.insert(make_pair(*(params.begin()), *(endogs.begin())));
+    else
+      if (endogs.size() == 2)
+        {
+          BinaryOpNode *testarg2 = dynamic_cast<BinaryOpNode *>(arg2);
+          VariableNode *test_arg1 = dynamic_cast<VariableNode *>(testarg2->get_arg1());
+          VariableNode *test_arg2 = dynamic_cast<VariableNode *>(testarg2->get_arg2());
+          if (testarg2 != NULL && testarg2->get_op_code() == oMinus
+              && test_arg1 != NULL &&test_arg2 != NULL
+              && lhs.first != -1)
+            {
+              int find_symb_id = -1;
+              try
+                {
+                  // lhs is an aux var (diff)
+                  find_symb_id = datatree.symbol_table.getOrigSymbIdForAuxVar(lhs.first);
+                }
+              catch (...)
+                {
+                  //lhs is not an aux var
+                  find_symb_id = lhs.first;
+                }
+              endogs.clear();
+
+              if (test_arg1->get_symb_id() == find_symb_id)
+                {
+                  test_arg1->collectDynamicVariables(eEndogenous, endogs);
+                  params_and_vals.insert(make_pair(*(params.begin()), *(endogs.begin())));
+                }
+              else if (test_arg2->get_symb_id() == find_symb_id)
+                {
+                  test_arg2->collectDynamicVariables(eEndogenous, endogs);
+                  params_and_vals.insert(make_pair(*(params.begin()), *(endogs.begin())));
+                }
+            }
+        }
+}
+
+void
+BinaryOpNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+  if (op_code == oTimes)
+    {
+      int orig_params_and_vals_size = params_and_vals.size();
+      walkPacParametersHelper(arg1, arg2, lhs, params_and_vals);
+      if (params_and_vals.size() == orig_params_and_vals_size)
+        walkPacParametersHelper(arg2, arg1, lhs, params_and_vals);
+    }
+  else if (op_code == oEqual)
+    {
+      set<pair<int, int> > general_lhs;
+      arg1->collectDynamicVariables(eEndogenous, general_lhs);
+      if (general_lhs.size() == 1)
+        lhs = *(general_lhs.begin());
+    }
+
+  arg1->walkPacParameters(pac_encountered, lhs, params_and_vals);
+  arg2->walkPacParameters(pac_encountered, lhs, params_and_vals);
+}
+
+void
+BinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+  arg1->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+  arg2->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+}
+
+void
+BinaryOpNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+  arg1->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+  arg2->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+}
+
 bool
 BinaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -5108,11 +5347,35 @@ TrinaryOpNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOp
 }
 
 expr_t
-TrinaryOpNode::substituteAdlAndDiff() const
+TrinaryOpNode::substituteAdl() const
+{
+  expr_t arg1subst = arg1->substituteAdl();
+  expr_t arg2subst = arg2->substituteAdl();
+  expr_t arg3subst = arg3->substituteAdl();
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
+expr_t
+TrinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  expr_t arg1subst = arg1->substituteDiff(subst_table, neweqs);
+  expr_t arg2subst = arg2->substituteDiff(subst_table, neweqs);
+  expr_t arg3subst = arg3->substituteDiff(subst_table, neweqs);
+  return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
+}
+
+bool
+TrinaryOpNode::isDiffPresent() const
+{
+  return arg1->isDiffPresent() || arg2->isDiffPresent() || arg3->isDiffPresent();
+}
+
+expr_t
+TrinaryOpNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
-  expr_t arg1subst = arg1->substituteAdlAndDiff();
-  expr_t arg2subst = arg2->substituteAdlAndDiff();
-  expr_t arg3subst = arg3->substituteAdlAndDiff();
+  expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
+  expr_t arg3subst = arg3->substitutePacExpectation(subst_table);
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
@@ -5190,6 +5453,30 @@ TrinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_m
   arg3->setVarExpectationIndex(var_model_info);
 }
 
+void
+TrinaryOpNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+  arg1->walkPacParameters(pac_encountered, lhs, params_and_vals);
+  arg2->walkPacParameters(pac_encountered, lhs, params_and_vals);
+  arg3->walkPacParameters(pac_encountered, lhs, params_and_vals);
+}
+
+void
+TrinaryOpNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+  arg1->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+  arg2->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+  arg3->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+}
+
+void
+TrinaryOpNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+  arg1->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+  arg2->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+  arg3->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+}
+
 bool
 TrinaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -5420,11 +5707,38 @@ AbstractExternalFunctionNode::substituteExpectation(subst_table_t &subst_table,
 }
 
 expr_t
-AbstractExternalFunctionNode::substituteAdlAndDiff() const
+AbstractExternalFunctionNode::substituteAdl() const
+{
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    arguments_subst.push_back((*it)->substituteAdl());
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
+expr_t
+AbstractExternalFunctionNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  vector<expr_t> arguments_subst;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    arguments_subst.push_back((*it)->substituteDiff(subst_table, neweqs));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
+bool
+AbstractExternalFunctionNode::isDiffPresent() const
+{
+  bool result = false;
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    result = result || (*it)->isDiffPresent();
+  return result;
+}
+
+expr_t
+AbstractExternalFunctionNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   vector<expr_t> arguments_subst;
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    arguments_subst.push_back((*it)->substituteAdlAndDiff());
+    arguments_subst.push_back((*it)->substitutePacExpectation(subst_table));
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
@@ -5529,6 +5843,27 @@ AbstractExternalFunctionNode::setVarExpectationIndex(map<string, pair<SymbolList
     (*it)->setVarExpectationIndex(var_model_info);
 }
 
+void
+AbstractExternalFunctionNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->walkPacParameters(pac_encountered, lhs, params_and_vals);
+}
+
+void
+AbstractExternalFunctionNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->addParamInfoToPac(lhs_arg, params_and_vals_arg);
+}
+
+void
+AbstractExternalFunctionNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->fillPacExpectationVarInfo(var_model_name_arg, lhs_arg, rhs_arg, nonstationary_arg, equation_number_arg);
+}
+
 bool
 AbstractExternalFunctionNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -6796,7 +7131,7 @@ void
 VarExpectationNode::prepareForDerivation()
 {
   preparedForDerivation = true;
-  // All derivatives are null, so non_null_derivatives is left empty
+  // Come back
 }
 
 expr_t
@@ -6827,6 +7162,12 @@ VarExpectationNode::eval(const eval_context_t &eval_context) const throw (EvalEx
   return it->second;
 }
 
+bool
+VarExpectationNode::isDiffPresent() const
+{
+  return false;
+}
+
 void
 VarExpectationNode::computeXrefs(EquationInfo &ei) const
 {
@@ -6892,7 +7233,19 @@ VarExpectationNode::substituteExpectation(subst_table_t &subst_table, vector<Bin
 }
 
 expr_t
-VarExpectationNode::substituteAdlAndDiff() const
+VarExpectationNode::substituteAdl() const
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
+expr_t
+VarExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
+expr_t
+VarExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
 {
   return const_cast<VarExpectationNode *>(this);
 }
@@ -6975,6 +7328,21 @@ VarExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &
   yidx = find(vs.begin(), vs.end(), datatree.symbol_table.getName(symb_id)) - vs.begin();
 }
 
+void
+VarExpectationNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+}
+
+void
+VarExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+}
+
+void
+VarExpectationNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+}
+
 expr_t
 VarExpectationNode::substituteStaticAuxiliaryVariable() const
 {
@@ -6994,3 +7362,500 @@ VarExpectationNode::writeJsonOutput(ostream &output,
          << ", yindex = " << yidx
          << ")";
 }
+
+PacExpectationNode::PacExpectationNode(DataTree &datatree_arg,
+                                       const string &model_name_arg,
+                                       const string &var_model_name_arg,
+                                       const int discount_symb_id_arg,
+                                       const int growth_symb_id_arg) :
+  ExprNode(datatree_arg),
+  model_name(model_name_arg),
+  var_model_name(var_model_name_arg),
+  discount_symb_id(discount_symb_id_arg),
+  growth_symb_id(growth_symb_id_arg),
+  stationary_vars_present(false),
+  nonstationary_vars_present(false)
+{
+  datatree.pac_expectation_node_map[make_pair(model_name, make_pair(var_model_name, make_pair(discount_symb_id, growth_symb_id)))] = this;
+}
+
+void
+PacExpectationNode::computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                          map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
+                                          bool is_matlab, NodeTreeReference tr) const
+{
+  temp_terms_map[tr].insert(const_cast<PacExpectationNode *>(this));
+}
+
+void
+PacExpectationNode::computeTemporaryTerms(map<expr_t, int> &reference_count,
+                                          temporary_terms_t &temporary_terms,
+                                          map<expr_t, pair<int, int> > &first_occurence,
+                                          int Curr_block,
+                                          vector< vector<temporary_terms_t> > &v_temporary_terms,
+                                          int equation) const
+{
+  expr_t this2 = const_cast<PacExpectationNode *>(this);
+  temporary_terms.insert(this2);
+  first_occurence[this2] = make_pair(Curr_block, equation);
+  v_temporary_terms[Curr_block][equation].insert(this2);
+}
+
+expr_t
+PacExpectationNode::toStatic(DataTree &static_datatree) const
+{
+  return static_datatree.AddPacExpectation(string(model_name), string(var_model_name), discount_symb_id, growth_symb_id);
+}
+
+expr_t
+PacExpectationNode::cloneDynamic(DataTree &dynamic_datatree) const
+{
+  return dynamic_datatree.AddPacExpectation(string(model_name), string(var_model_name), discount_symb_id, growth_symb_id);
+}
+
+void
+PacExpectationNode::writeOutput(ostream &output, ExprNodeOutputType output_type,
+                                const temporary_terms_t &temporary_terms,
+                                deriv_node_temp_terms_t &tef_terms) const
+{
+  assert(output_type != oMatlabOutsideModel);
+
+  if (IS_LATEX(output_type))
+    {
+      output << "PAC_EXPECTATION" << LEFT_PAR(output_type) << model_name << ", "
+             << var_model_name << ", " << discount_symb_id;
+      if (growth_symb_id >= 0)
+        output << ", " << growth_symb_id;
+      output << RIGHT_PAR(output_type);
+      return;
+    }
+
+  output <<"M_.pac_expectation." << model_name << ".var_model_name = '" << var_model_name << "';" << endl
+         << "M_.pac_expectation." << model_name << ".discount_index = "
+         << datatree.symbol_table.getTypeSpecificID(discount_symb_id) + 1 << ";" << endl
+         << "M_.pac_expectation." << model_name << ".equation_number = " << equation_number + 1 << ";" << endl
+         << "M_.pac_expectation." << model_name << ".lhs_var = "
+         << datatree.symbol_table.getTypeSpecificID(lhs_pac_var.first) + 1 << ";" << endl
+         << "M_.pac_expectation." << model_name << ".lhs_lag = " << lhs_pac_var.second << ";" << endl;
+
+  if (growth_symb_id >= 0)
+    {
+      output << "M_.pac_expectation." << model_name << ".growth_neutrality_param_index = "
+             << datatree.symbol_table.getTypeSpecificID(growth_param_index) + 1 << ";" << endl
+             << "M_.pac_expectation." << model_name << ".growth_index = "
+             << datatree.symbol_table.getTypeSpecificID(growth_symb_id) + 1 << ";" << endl
+             << "M_.pac_expectation." << model_name << ".growth_type = ";
+      switch(datatree.symbol_table.getType(growth_symb_id))
+        {
+        case eEndogenous:
+          output << "'endogenous';" << endl;
+          break;
+        case eExogenous:
+          output << "'exogenous';" << endl;
+          break;
+        case eParameter:
+          output << "'exogenous';" << endl;
+          break;
+        default:
+          cerr << "pac_expectation: error encountered in growth type" << endl;
+          exit(EXIT_FAILURE);
+        }
+    }
+
+  output << "M_.pac_expectation." << model_name << ".equation_params = [";
+  for (set<pair<int, pair<int, int> > >::const_iterator it = params_and_vals.begin();
+       it != params_and_vals.end(); it++)
+    {
+      if (it != params_and_vals.begin())
+        output << " ";
+      output << datatree.symbol_table.getTypeSpecificID(it->first) + 1;
+    }
+  output << "];" << endl
+         << "M_.pac_expectation." << model_name << ".equation_vars = [";
+  for (set<pair<int, pair<int, int> > >::const_iterator it = params_and_vals.begin();
+       it != params_and_vals.end(); it++)
+    {
+      if (it != params_and_vals.begin())
+        output << " ";
+      output << datatree.symbol_table.getTypeSpecificID(it->second.first) + 1;
+    }
+  output << "];" << endl
+         << "M_.pac_expectation." << model_name << ".equation_var_lags = [";
+  for (set<pair<int, pair<int, int> > >::const_iterator it = params_and_vals.begin();
+       it != params_and_vals.end(); it++)
+    {
+      if (it != params_and_vals.begin())
+        output << " ";
+      output << it->second.second;
+    }
+  output << "];" << endl
+         << "M_.pac_expectation." << model_name << ".h0_param_indices = [";
+  for (vector<int>::const_iterator it = h0_indices.begin();
+       it != h0_indices.end(); it++)
+    {
+      if (it != h0_indices.begin())
+        output << " ";
+      output << datatree.symbol_table.getTypeSpecificID(*it) + 1;
+    }
+  output << "];" << endl
+         << "M_.pac_expectation." << model_name << ".h1_param_indices = [";
+  for (vector<int>::const_iterator it = h1_indices.begin();
+       it != h1_indices.end(); it++)
+    {
+      if (it != h1_indices.begin())
+        output << " ";
+      output << datatree.symbol_table.getTypeSpecificID(*it) + 1;
+    }
+  output << "];" << endl;
+}
+
+int
+PacExpectationNode::maxEndoLead() const
+{
+  return 0;
+}
+
+int
+PacExpectationNode::maxExoLead() const
+{
+  return 0;
+}
+
+int
+PacExpectationNode::maxEndoLag() const
+{
+  return 0;
+}
+
+int
+PacExpectationNode::maxExoLag() const
+{
+  return 0;
+}
+
+int
+PacExpectationNode::maxLead() const
+{
+  return 0;
+}
+
+expr_t
+PacExpectationNode::decreaseLeadsLags(int n) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+void
+PacExpectationNode::prepareForDerivation()
+{
+  cerr << "PacExpectationNode::prepareForDerivation: shouldn't arrive here." << endl;
+  exit(EXIT_FAILURE);
+}
+
+expr_t
+PacExpectationNode::computeDerivative(int deriv_id)
+{
+  cerr << "PacExpectationNode::computeDerivative: shouldn't arrive here." << endl;
+  exit(EXIT_FAILURE);
+}
+
+expr_t
+PacExpectationNode::getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables)
+{
+  cerr << "PacExpectationNode::getChainRuleDerivative: shouldn't arrive here." << endl;
+  exit(EXIT_FAILURE);
+}
+
+bool
+PacExpectationNode::containsExternalFunction() const
+{
+  return false;
+}
+
+double
+PacExpectationNode::eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException)
+{
+  throw EvalException();
+}
+
+void
+PacExpectationNode::computeXrefs(EquationInfo &ei) const
+{
+}
+
+void
+PacExpectationNode::collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const
+{
+}
+
+void
+PacExpectationNode::collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const
+{
+  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<PacExpectationNode *>(this));
+  if (it != temporary_terms.end())
+    temporary_terms_inuse.insert(idx);
+}
+
+void
+PacExpectationNode::compile(ostream &CompileCode, unsigned int &instruction_number,
+                            bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                            const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                            deriv_node_temp_terms_t &tef_terms) const
+{
+  cerr << "PacExpectationNode::compile not implemented." << endl;
+  exit(EXIT_FAILURE);
+}
+
+bool
+PacExpectationNode::isDiffPresent() const
+{
+  return false;
+}
+
+pair<int, expr_t >
+PacExpectationNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > > &List_of_Op_RHS) const
+{
+  //COME BACK
+  return make_pair(0, const_cast<PacExpectationNode *>(this));
+}
+
+expr_t
+PacExpectationNode::substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::substituteAdl() const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+bool
+PacExpectationNode::containsEndogenous(void) const
+{
+  return true;
+}
+
+bool
+PacExpectationNode::containsExogenous() const
+{
+  return false;
+}
+
+bool
+PacExpectationNode::isNumConstNodeEqualTo(double value) const
+{
+  return false;
+}
+
+expr_t
+PacExpectationNode::decreaseLeadsLagsPredeterminedVariables() const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+bool
+PacExpectationNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const
+{
+  return false;
+}
+
+expr_t
+PacExpectationNode::replaceTrendVar() const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::detrend(int symb_id, bool log_trend, expr_t trend) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+expr_t
+PacExpectationNode::removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+bool
+PacExpectationNode::isInStaticForm() const
+{
+  return false;
+}
+
+bool
+PacExpectationNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return model_name == model_info_name;
+}
+
+void
+PacExpectationNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+}
+
+void
+PacExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
+{
+}
+
+expr_t
+PacExpectationNode::substituteStaticAuxiliaryVariable() const
+{
+  return const_cast<PacExpectationNode *>(this);
+}
+
+void
+PacExpectationNode::writeJsonOutput(ostream &output,
+                                    const temporary_terms_t &temporary_terms,
+                                    deriv_node_temp_terms_t &tef_terms,
+                                    const bool isdynamic) const
+{
+  output << "pac_expectation("
+         << ", model_name = " << model_name
+         << ", " << discount_symb_id;
+  if (growth_symb_id >= 0)
+    output << ", " << growth_symb_id;
+  output << ")";
+}
+
+void
+PacExpectationNode::walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const
+{
+  pac_encountered = true;
+}
+
+void
+PacExpectationNode::addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg)
+{
+  if (lhs_arg.first == -1)
+    {
+      cerr << "Pac Expectation: error in obtaining LHS varibale." << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  if (params_and_vals_arg.size() != 3)
+    {
+      cerr << "Pac Expectation: error in obtaining RHS parameters." << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  lhs_pac_var = lhs_arg;
+  params_and_vals = params_and_vals_arg;
+}
+
+
+void
+PacExpectationNode::fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg)
+{
+  if (var_model_name != var_model_name_arg)
+    return;
+
+  lhs = lhs_arg;
+  z_vec = rhs_arg;
+  equation_number = equation_number_arg;
+
+  for (vector<bool>::const_iterator it = nonstationary_arg.begin();
+       it != nonstationary_arg.end(); it++)
+    {
+      if (*it)
+        nonstationary_vars_present = true;
+      else
+        stationary_vars_present = true;
+      if (nonstationary_vars_present && stationary_vars_present)
+        break;
+    }
+}
+
+expr_t
+PacExpectationNode::substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table)
+{
+  map<const PacExpectationNode *, const BinaryOpNode *>::const_iterator myit =
+    subst_table.find(const_cast<PacExpectationNode *>(this));
+  if (myit != subst_table.end())
+    return const_cast<BinaryOpNode *>(myit->second);
+
+  int maxlag = z_vec.size();
+  expr_t subExpr = datatree.AddNonNegativeConstant("0");
+  if (stationary_vars_present)
+    for (int i = 1; i < maxlag + 1; i++)
+      for (vector<int>::const_iterator it = lhs.begin(); it != lhs.end(); it++)
+        {
+          stringstream param_name_h0;
+          param_name_h0 << "h0_" << model_name
+                        << "_var_" << datatree.symbol_table.getName(*it)
+                        << "_lag_" << i;
+          int new_param_symb_id = datatree.symbol_table.addSymbol(param_name_h0.str(), eParameter);
+          h0_indices.push_back(new_param_symb_id);
+          subExpr = datatree.AddPlus(subExpr,
+                                     datatree.AddTimes(datatree.AddVariable(new_param_symb_id),
+                                                       datatree.AddVariable(*it, -i)));
+        }
+
+  if (nonstationary_vars_present)
+    for (int i = 1; i < maxlag + 1; i++)
+      for (vector<int>::const_iterator it = lhs.begin(); it != lhs.end(); it++)
+        {
+          stringstream param_name_h1;
+          param_name_h1 << "h1_" << model_name
+                        << "_var_" << datatree.symbol_table.getName(*it)
+                        << "_lag_" << i;
+          int new_param_symb_id = datatree.symbol_table.addSymbol(param_name_h1.str(), eParameter);
+          h1_indices.push_back(new_param_symb_id);
+          subExpr = datatree.AddPlus(subExpr,
+                                     datatree.AddTimes(datatree.AddVariable(new_param_symb_id),
+                                                       datatree.AddVariable(*it, -i)));
+        }
+
+  if (growth_symb_id >= 0)
+    {
+      growth_param_index = datatree.symbol_table.addSymbol(model_name +
+                                                           "_pac_growth_neutrality_correction",
+                                                           eParameter);
+      subExpr = datatree.AddPlus(subExpr,
+                                 datatree.AddTimes(datatree.AddVariable(growth_param_index),
+                                                   datatree.AddVariable(growth_symb_id)));
+    }
+
+  subst_table[const_cast<PacExpectationNode *>(this)] = dynamic_cast<BinaryOpNode *>(subExpr);
+
+  return subExpr;
+}
diff --git a/ExprNode.hh b/ExprNode.hh
index 9fa28e05c513dd4f363f1cf7fba59b6966a9984e..cc8409adc536417a1f3947332c5af25dddcb4313 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2017 Dynare Team
+ * Copyright (C) 2007-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -35,6 +35,7 @@ using namespace std;
 class DataTree;
 class VariableNode;
 class BinaryOpNode;
+class PacExpectationNode;
 
 typedef class ExprNode *expr_t;
 
@@ -123,7 +124,7 @@ enum ExprNodeOutputType
 #define MIN_COST(is_matlab) ((is_matlab) ? MIN_COST_MATLAB : MIN_COST_C)
 
 //! Base class for expression nodes
-    class ExprNode
+class ExprNode
     {
       friend class DataTree;
       friend class DynamicModel;
@@ -137,6 +138,7 @@ enum ExprNodeOutputType
       friend class TrinaryOpNode;
       friend class AbstractExternalFunctionNode;
       friend class VarExpectationNode;
+      friend class PacExpectationNode;
     private:
       //! Computes derivative w.r. to a derivation ID (but doesn't store it in derivatives map)
       /*! You shoud use getDerivative() to get the benefit of symbolic a priori and of caching */
@@ -438,6 +440,9 @@ enum ExprNodeOutputType
       //! Returns true if the expression contains one or several exogenous variable
       virtual bool containsExogenous() const = 0;
 
+      //! Returns true if the expression contains a diff operator
+      virtual bool isDiffPresent(void) const = 0;
+
       //! Return true if the nodeID is a variable withe a type equal to type_arg, a specific variable id aqual to varfiable_id and a lag equal to lag_arg and false otherwise
       /*!
         \param[in] the type (type_arg), specifique variable id (variable_id and the lag (lag_arg)
@@ -457,8 +462,14 @@ enum ExprNodeOutputType
       */
       virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const = 0;
 
-      //! Substitute adl and diff operators
-      virtual expr_t substituteAdlAndDiff() const = 0;
+      //! Substitute adl operator
+      virtual expr_t substituteAdl() const = 0;
+
+      //! Substitute diff operator
+      virtual expr_t substituteDiff(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;
 
       //! Add ExprNodes to the provided datatree
       virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const = 0;
@@ -478,6 +489,15 @@ enum ExprNodeOutputType
       //! Returns true if model_info_name is referenced by a VarExpectationNode
       virtual bool isVarModelReferenced(const string &model_info_name) const = 0;
 
+      //! Fills parameter information related to PAC equation
+      virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const = 0;
+
+      //! Adds PAC equation param info to pac_expectation
+      virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg) = 0;
+
+      //! Fills var_model info for pac_expectation node
+      virtual void fillPacExpectationVarInfo(string &var_model_name, vector<int> &lhs, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg) = 0;
+
       //! Fills map
       virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const = 0;
     };
@@ -530,12 +550,15 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
   virtual expr_t replaceTrendVar() const;
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
@@ -543,6 +566,9 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   virtual expr_t substituteStaticAuxiliaryVariable() const;
@@ -605,12 +631,15 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
   virtual expr_t replaceTrendVar() const;
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
@@ -618,6 +647,9 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   //! Substitute auxiliary variables by their expression in static model
@@ -703,12 +735,15 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
   virtual expr_t replaceTrendVar() const;
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
@@ -716,6 +751,9 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   //! Substitute auxiliary variables by their expression in static model
@@ -796,6 +834,9 @@ public:
   {
     return powerDerivOrder;
   }
+  void walkPacParametersHelper(const expr_t arg1, const expr_t arg2,
+                               pair<int, int> &lhs,
+                               set<pair<int, pair<int, int> > > &params_and_vals) const;
   virtual expr_t toStatic(DataTree &static_datatree) const;
   virtual void computeXrefs(EquationInfo &ei) const;
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
@@ -813,12 +854,15 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
   virtual expr_t replaceTrendVar() const;
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
@@ -832,6 +876,9 @@ public:
   expr_t getNonZeroPartofEquation() const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   //! Substitute auxiliary variables by their expression in static model
@@ -903,12 +950,15 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
   virtual expr_t replaceTrendVar() const;
   virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
@@ -916,6 +966,9 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   //! Substitute auxiliary variables by their expression in static model
@@ -993,13 +1046,16 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual expr_t buildSimilarExternalFunctionNode(vector<expr_t> &alt_args, DataTree &alt_datatree) const = 0;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
   virtual void writePrhs(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const string &ending) const;
   virtual expr_t replaceTrendVar() const;
@@ -1008,6 +1064,9 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   //! Substitute auxiliary variables by their expression in static model
@@ -1174,7 +1233,84 @@ public:
   virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
   virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
-  virtual expr_t substituteAdlAndDiff() const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
+  virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
+  virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
+                       bool lhs_rhs, const temporary_terms_t &temporary_terms,
+                       const map_idx_t &map_idx, bool dynamic, bool steady_dynamic,
+                       deriv_node_temp_terms_t &tef_terms) const;
+  virtual void collectTemporary_terms(const temporary_terms_t &temporary_terms, temporary_terms_inuse_t &temporary_terms_inuse, int Curr_Block) const;
+  virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
+  virtual bool containsEndogenous(void) const;
+  virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
+  virtual bool isNumConstNodeEqualTo(double value) const;
+  virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
+  virtual bool isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_arg) const;
+  virtual expr_t replaceTrendVar() const;
+  virtual expr_t detrend(int symb_id, bool log_trend, expr_t trend) const;
+  virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
+  virtual bool isInStaticForm() const;
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
+  virtual bool isVarModelReferenced(const string &model_info_name) const;
+  virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
+  virtual expr_t substituteStaticAuxiliaryVariable() const;
+  virtual void writeJsonOutput(ostream &output, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms, const bool isdynamic) const;
+};
+
+class PacExpectationNode : public ExprNode
+{
+private:
+  const string model_name, var_model_name;
+  const int discount_symb_id, growth_symb_id;
+  bool stationary_vars_present, nonstationary_vars_present;
+  vector<int> lhs;
+  pair<int, int> lhs_pac_var;
+  map<int, set<int > > z_vec; // lag -> set< symb_id > (all vars that appear at a given lag)
+  vector<int> h0_indices, h1_indices;
+  int growth_param_index, equation_number;
+  set<pair<int, pair<int, int> > > params_and_vals;
+public:
+  PacExpectationNode(DataTree &datatree_arg, const string &model_name, const string &var_model_name,
+                     const int discount_arg, const int growth_arg);
+  virtual void computeTemporaryTerms(map<expr_t, pair<int, NodeTreeReference> > &reference_count,
+                                     map<NodeTreeReference, temporary_terms_t> &temp_terms_map,
+                                     bool is_matlab, NodeTreeReference tr) const;
+  virtual void writeOutput(ostream &output, ExprNodeOutputType output_type, const temporary_terms_t &temporary_terms, deriv_node_temp_terms_t &tef_terms) const;
+  virtual void computeTemporaryTerms(map<expr_t, int> &reference_count,
+                                     temporary_terms_t &temporary_terms,
+                                     map<expr_t, pair<int, int> > &first_occurence,
+                                     int Curr_block,
+                                     vector< vector<temporary_terms_t> > &v_temporary_terms,
+                                     int equation) const;
+  virtual expr_t toStatic(DataTree &static_datatree) const;
+  virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
+  virtual int maxEndoLead() const;
+  virtual int maxExoLead() const;
+  virtual int maxEndoLag() const;
+  virtual int maxExoLag() const;
+  virtual int maxLead() const;
+  virtual expr_t decreaseLeadsLags(int n) const;
+  virtual void prepareForDerivation();
+  virtual expr_t computeDerivative(int deriv_id);
+  virtual expr_t getChainRuleDerivative(int deriv_id, const map<int, expr_t> &recursive_variables);
+  virtual bool containsExternalFunction() const;
+  virtual double eval(const eval_context_t &eval_context) const throw (EvalException, EvalExternalFunctionException);
+  virtual void computeXrefs(EquationInfo &ei) const;
+  virtual expr_t substituteEndoLeadGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
+  virtual expr_t substituteEndoLagGreaterThanTwo(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExoLead(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool deterministic_model) const;
+  virtual expr_t substituteExoLag(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substituteExpectation(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs, bool partial_information_model) const;
+  virtual expr_t substituteAdl() const;
+  virtual expr_t substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
+  virtual expr_t substitutePacExpectation(map<const PacExpectationNode *, const BinaryOpNode *> &subst_table);
   virtual pair<int, expr_t> normalizeEquation(int symb_id_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const;
   virtual void compile(ostream &CompileCode, unsigned int &instruction_number,
                        bool lhs_rhs, const temporary_terms_t &temporary_terms,
@@ -1184,6 +1320,7 @@ public:
   virtual void collectDynamicVariables(SymbolType type_arg, set<pair<int, int> > &result) const;
   virtual bool containsEndogenous(void) const;
   virtual bool containsExogenous() const;
+  virtual bool isDiffPresent(void) const;
   virtual bool isNumConstNodeEqualTo(double value) const;
   virtual expr_t differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const;
   virtual expr_t decreaseLeadsLagsPredeterminedVariables() const;
@@ -1193,6 +1330,9 @@ public:
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void walkPacParameters(bool &pac_encountered, pair<int, int> &lhs, set<pair<int, pair<int, int> > > &params_and_vals) const;
+  virtual void addParamInfoToPac(pair<int, int> &lhs_arg, set<pair<int, pair<int, int> > > &params_and_vals_arg);
+  virtual void fillPacExpectationVarInfo(string &var_model_name_arg, vector<int> &lhs_arg, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg, int equation_number_arg);
   virtual bool isVarModelReferenced(const string &model_info_name) const;
   virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const;
   virtual expr_t substituteStaticAuxiliaryVariable() const;
diff --git a/ModFile.cc b/ModFile.cc
index 1de813fcde3c47fc60359e16bfad9dd21930651c..eb4a374b2b73e53277f46c15769ee4030dbff03e 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -346,7 +346,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
 {
   // Save the original model (must be done before any model transformations by preprocessor)
   // - except adl and diff which we always want expanded
-  dynamic_model.substituteAdlAndDiff();
+  dynamic_model.substituteAdl();
   dynamic_model.setLeadsLagsOrig();
   dynamic_model.cloneDynamic(original_model);
 
@@ -361,20 +361,45 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
         }
     }
 
+  // Create auxiliary variable and equations for Diff operator
+  dynamic_model.substituteDiff();
+
   // Var Model
-  map<string, pair<SymbolList, int> > var_model_info;
+  map<string, pair<SymbolList, int> > var_model_info_var_expectation;
   for (vector<Statement *>::const_iterator it = statements.begin();
        it != statements.end(); it++)
     {
       VarModelStatement *vms = dynamic_cast<VarModelStatement *>(*it);
       if (vms != NULL)
-        vms->getVarModelNameAndVarList(var_model_info);
+        {
+          vms->getVarModelInfoForVarExpectation(var_model_info_var_expectation);
+
+          vector<string> var_model_eqtags;
+          vms->getVarModelEqTags(var_model_eqtags);
+          if (!var_model_eqtags.empty())
+            {
+              vector<int> eqnumber, lhs, orig_diff_var;
+              vector<set<pair<int, int> > > rhs;
+              vector<bool> nonstationary, diff;
+              dynamic_model.getVarModelVariablesFromEqTags(var_model_eqtags,
+                                                           eqnumber, lhs, rhs, nonstationary);
+              original_model.getDiffInfo(eqnumber, diff, orig_diff_var);
+              vms->fillVarModelInfoFromEquations(eqnumber, lhs, rhs, nonstationary, diff, orig_diff_var);
+              string var_model_name;
+              map<int, set<int > > rhs_pac;
+              vms->getVarModelName(var_model_name);
+              vms->getVarModelRHS(rhs_pac);
+              dynamic_model.walkPacParameters();
+              dynamic_model.fillPacExpectationVarInfo(var_model_name, lhs, rhs_pac, nonstationary);
+              dynamic_model.substitutePacExpectation();
+            }
+        }
     }
 
-  if (!var_model_info.empty())
+  if (!var_model_info_var_expectation.empty())
     {
-      dynamic_model.setVarExpectationIndices(var_model_info);
-      dynamic_model.addEquationsForVar(var_model_info);
+      dynamic_model.setVarExpectationIndices(var_model_info_var_expectation);
+      dynamic_model.addEquationsForVar(var_model_info_var_expectation);
     }
   dynamic_model.fillVarExpectationFunctionsToWrite();
 
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index cf8760113268d06665ffd2fd40efc8bb965b63af..4e2ede104be0e99a735ab9c036406ccd1bc2c6f9 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -1506,7 +1506,18 @@ ParsingDriver::var_model()
   if (it == options_list.string_options.end())
     error("You must pass the model_name option to the var_model statement.");
   const string *name = new string(it->second);
-  mod_file->addStatement(new VarModelStatement(symbol_list, options_list, *name));
+
+  if (options_list.vector_str_options.find("var.eqtags") != options_list.vector_str_options.end())
+    if (!symbol_list.empty())
+      error("You cannot pass a symbol list when passing equation tags to the var_model statement");
+    else if (options_list.num_options.find("var.order") != options_list.num_options.end())
+      error("You cannot pass the order option when passing equation tags to the var_model statement");
+
+  if (!symbol_list.empty())
+    if (options_list.num_options.find("var.order") == options_list.num_options.end())
+      error("You must pass the order option when passing a symbol list to the var_model statement");
+
+  mod_file->addStatement(new VarModelStatement(symbol_list, options_list, *name, mod_file->symbol_table));
   var_map[it->second] = symbol_list.getSymbols();
   symbol_list.clear();
   options_list.clear();
@@ -2661,6 +2672,76 @@ ParsingDriver::add_var_expectation(string *arg1, string *arg2, string *arg3)
   return varExpectationNode;
 }
 
+expr_t
+ParsingDriver::add_pac_expectation()
+{
+  if (pac_expectation_model_name.empty())
+    error("pac_expectation: you must pass the model_name option");
+
+  if (pac_expectation_var_model_name.empty())
+    error("pac_expectation: you must pass the var_model_name option");
+
+  if (pac_expectation_discount.empty())
+    error("pac_expectation: you must pass the discount option");
+
+  int pac_expectation_discount_id =
+    mod_file->symbol_table.getID(pac_expectation_discount);
+
+  int pac_expectation_growth_id = -1;
+  if (!pac_expectation_growth.empty())
+    pac_expectation_growth_id = mod_file->symbol_table.getID(pac_expectation_growth);
+
+  expr_t pac_exp_node = data_tree->AddPacExpectation(pac_expectation_model_name,
+                                                     pac_expectation_var_model_name,
+                                                     pac_expectation_discount_id,
+                                                     pac_expectation_growth_id);
+
+  pac_expectation_model_name = pac_expectation_discount = pac_expectation_growth = "";
+
+  return pac_exp_node;
+}
+
+void
+ParsingDriver::add_pac_expectation_model_name(string *arg)
+{
+  if (!pac_expectation_model_name.empty())
+    error("pac_expectation: you can only pass the model_name option once");
+  pac_expectation_model_name = *arg;
+  delete arg;
+}
+
+void
+ParsingDriver::add_pac_expectation_var_model_name(string *arg)
+{
+  if (!pac_expectation_var_model_name.empty())
+    error("pac_expectation: you can only pass the var_model_name option once");
+  pac_expectation_var_model_name = *arg;
+  delete arg;
+}
+
+void
+ParsingDriver::add_pac_expectation_discount(string *arg)
+{
+  if (!pac_expectation_discount.empty())
+    error("pac_expectation: you can only pass the discount option once");
+  check_symbol_is_parameter(arg);
+  pac_expectation_discount = *arg;
+  delete arg;
+}
+
+void
+ParsingDriver::add_pac_expectation_growth(string *arg)
+{
+  if (!pac_expectation_growth.empty())
+    error("pac_expectation: you can only pass the growth option once");
+  check_symbol_existence(*arg);
+  SymbolType type = mod_file->symbol_table.getType(mod_file->symbol_table.getID(*arg));
+  if (type != eParameter && type != eEndogenous && type != eExogenous)
+    error("pac_expectation growth argument must either be a parameter or an endogenous or exogenous variable.");
+  pac_expectation_growth = *arg;
+  delete arg;
+}
+
 expr_t
 ParsingDriver::add_exp(expr_t arg1)
 {
diff --git a/ParsingDriver.hh b/ParsingDriver.hh
index 872c391728f85ac41962731727f7eb353af02772..3b4d28ab9ca963eb244882fe9a4c43839f320e17 100644
--- a/ParsingDriver.hh
+++ b/ParsingDriver.hh
@@ -251,6 +251,10 @@ private:
 
   //! Used by VAR restrictions
   void clear_VAR_storage();
+
+  //! Used by pac_expectation
+  string pac_expectation_model_name, pac_expectation_var_model_name, pac_expectation_discount, pac_expectation_growth;
+
 public:
   ParsingDriver(WarningConsolidation &warnings_arg, bool nostrict_arg) : warnings(warnings_arg), nostrict(nostrict_arg), model_error_encountered(false)
   {
@@ -686,6 +690,13 @@ public:
   expr_t add_expectation(string *arg1,  expr_t arg2);
   //! Writes token "VAR_EXPECTATION(arg1, arg2, arg3)" to model tree
   expr_t add_var_expectation(string *arg1,  string *arg2, string *arg3);
+  //! Writes token "PAC_EXPECTATION(model_name, discount, growth)" to model tree
+  expr_t add_pac_expectation();
+  //! Adds arguments for pac_expectation
+  void add_pac_expectation_model_name(string *arg);
+  void add_pac_expectation_var_model_name(string *arg);
+  void add_pac_expectation_discount(string *arg);
+  void add_pac_expectation_growth(string *arg);
   //! 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/SymbolTable.cc b/SymbolTable.cc
index dde3d141bc047ed7d1202920b61538721ab0cd92..ad61e3749c2b089a88d80af2b86d7bfd2cd7365d 100644
--- a/SymbolTable.cc
+++ b/SymbolTable.cc
@@ -374,6 +374,11 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
             aux_vars[i].get_expr_node()->writeOutput(output, oLatexDynamicModel);
             output << ")';" << endl;
             break;
+          case avDiff:
+            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;
+            break;
           }
       }
 
@@ -482,6 +487,11 @@ SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
               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;
               break;
+            case avDiff:
+              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;
+              break;
             }
         }
     }
@@ -575,6 +585,11 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
           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;
           break;
+        case avDiff:
+          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;
+          break;
         }
       output << "aux_vars.push_back(" << "av" << i << ");" << endl;
     }
@@ -691,6 +706,35 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t e
   return symb_id;
 }
 
+int
+SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException)
+{
+  ostringstream varname;
+  int symb_id;
+
+  varname << "AUX_DIFF_" << index;
+
+  try
+    {
+      symb_id = addSymbol(varname.str(), eEndogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_vars.push_back(AuxVarInfo(symb_id, avDiff, orig_symb_id, orig_lag, 0, 0, expr_arg));
+
+  return symb_id;
+}
+
+int
+SymbolTable::addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException)
+{
+  return addDiffAuxiliaryVar(index, expr_arg, 0, 0);
+}
+
 int
 SymbolTable::addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException)
 {
@@ -765,6 +809,16 @@ SymbolTable::searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const thro
   throw SearchFailedException(orig_symb_id, orig_lead_lag);
 }
 
+int
+SymbolTable::getOrigSymbIdForAuxVar(int aux_var_symb_id) const throw (UnknownSymbolIDException)
+{
+  for (size_t i = 0; i < aux_vars.size(); i++)
+    if ((aux_vars[i].get_type() == avEndoLag || aux_vars[i].get_type() == avExoLag || aux_vars[i].get_type() == avDiff)
+        && aux_vars[i].get_symb_id() == aux_var_symb_id)
+      return aux_vars[i].get_orig_symb_id();
+  throw UnknownSymbolIDException(aux_var_symb_id);
+}
+
 expr_t
 SymbolTable::getAuxiliaryVarsExprNode(int symb_id) const throw (SearchFailedException)
 // throw exception if it is a Lagrange multiplier
@@ -1008,6 +1062,11 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
               output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
                      << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
               break;
+            case avDiff:
+              if (aux_vars[i].get_orig_symb_id() >= 0)
+                output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
+                       << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
+              break;
             case avMultiplier:
               output << "NaN, NaN, " << aux_vars[i].get_equation_number_for_multiplier() + 1
                      << ", NaN";
diff --git a/SymbolTable.hh b/SymbolTable.hh
index 76aa3896141782b1533d21bad56da4adaa58953b..62804a59b644dc148c5a31749b44e1a7897f47e4 100644
--- a/SymbolTable.hh
+++ b/SymbolTable.hh
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2003-2017 Dynare Team
+ * Copyright (C) 2003-2018 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -43,7 +43,8 @@ enum aux_var_t
     avExpectation = 4,    //!< Substitute for Expectation Operator
     avDiffForward = 5,    //!< Substitute for the differentiate of a forward variable
     avMultiplier = 6,     //!< Multipliers for FOC of Ramsey Problem
-    avVarModel = 7        //!< Variable for var_model with order > abs(min_lag()) present in model
+    avVarModel = 7,       //!< Variable for var_model with order > abs(min_lag()) present in model
+    avDiff = 8            //!< Variable for Diff operator
   };
 
 //! Information on some auxiliary variables
@@ -280,9 +281,14 @@ public:
     Throws an exception if match not found.
   */
   int searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const throw (SearchFailedException);
+  //! Serches aux_vars for the aux var represented by aux_var_symb_id and returns its associated orig_symb_id
+  int getOrigSymbIdForAuxVar(int aux_var_symb_id) const throw (UnknownSymbolIDException);
   //! Adds an auxiliary variable when var_model is used with an order that is greater in absolute value
   //! than the largest lag present in the model.
   int addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException);
+  //! Adds an auxiliary variable when the diff operator is encountered
+  int addDiffAuxiliaryVar(int index, expr_t expr_arg) throw (FrozenException);
+  int addDiffAuxiliaryVar(int index, expr_t expr_arg, int orig_symb_id, int orig_lag) throw (FrozenException);
   //! Returns the number of auxiliary variables
   int
   AuxVarsSize() const