diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index f45a4e3177b0ee94eb2138d9bf92d50b8d1b6710..bf5e74bd9bbf27dd62aa3fec42d2c328372643c5 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -268,30 +268,28 @@ VarModelStatement::VarModelStatement(const SymbolList &symbol_list_arg,
 }
 
 void
-VarModelStatement::getVarModelNameAndVarList(map<string, pair<SymbolList, int> > &var_model_info)
+VarModelStatement::getVarModelNameAndVarList(map<string, pair<map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >, pair<SymbolList, int> > > &var_model_info) const
 {
-  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()));
+  set<pair<int, int> > empty_int_set;
+  map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > eqtagmap;
+  if (!symbol_list.empty())
+    {
+      OptionsList::num_options_t::const_iterator it = options_list.num_options.find("var.order");
+      vector<string> empty_str_vec;
+      var_model_info[name] = make_pair(eqtagmap, make_pair(symbol_list, atoi(it->second.c_str())));
+    }
+  else
+    {
+      OptionsList::vec_str_options_t::const_iterator it1 = options_list.vector_str_options.find("var.eqtags");
+      for (vector<string>::const_iterator it = it1->second.begin(); it != it1->second.end(); it++)
+        eqtagmap[make_pair(*it, -1)] = make_pair(make_pair(0, empty_int_set), empty_int_set);
+      var_model_info[name] = make_pair(eqtagmap, make_pair(symbol_list, 0));
+    }
 }
 
 void
 VarModelStatement::checkPass(ModFileStructure &mod_file_struct, WarningConsolidation &warnings)
 {
-  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 (name.empty())
-    {
-      cerr << "ERROR: You must provide the model_name option to the var_model statement." << endl;
-      exit(EXIT_FAILURE);
-    }
 }
 
 void
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index 8432be9fa28f6dc70168b503633b4405b56f71ba..1799f37da142bd10c4d8b2269c20e4b2115e0ed7 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.
  *
@@ -129,7 +129,7 @@ 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);
+  void getVarModelNameAndVarList(map<string, pair<map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >, pair<SymbolList, int> > > &var_model_info) 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..f2179a6ddae2809ec95e0ee8adfe4aaf9af89f76 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,16 @@ 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 expr_t discount, const expr_t growth)
+{
+  pac_expectation_node_map_t::iterator it = pac_expectation_node_map.find(make_pair(model_name, make_pair(discount, growth)));
+  if (it != pac_expectation_node_map.end())
+    return it->second;
+
+  return new PacExpectationNode(*this, model_name, discount, growth);
+}
+
 expr_t
 DataTree::AddEqual(expr_t iArg1, expr_t iArg2)
 {
diff --git a/DataTree.hh b/DataTree.hh
index 3cda0bf2a7524583e5d8c3492fc07e838ee73d04..61be6c0b8e594be6cf8bb626a3fc1b80b112c0d2 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<expr_t, expr_t> >, 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, expr_t discount, expr_t growth);
   //! 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 67b28c81542501be1e901e96b5f4cd66dae736ee..fb23c71d0182026764b8194684d6ee7d7f371f84 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3228,14 +3228,62 @@ DynamicModel::runTrendTest(const eval_context_t &eval_context)
 }
 
 void
-DynamicModel::setVarExpectationIndices(map<string, pair<SymbolList, int> > var_model_info)
+DynamicModel::getVarModelVariablesFromEqTags(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  for (map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > >::iterator it = var_model_info.begin(); it != var_model_info.end(); it++)
+    {
+      map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > final_map;
+      for (map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >::iterator itvareqs = it->second.begin();
+           itvareqs != it->second.end(); itvareqs++)
+        {
+          int eqnumber = -1;
+          set<pair<int, int> > lhs, rhs;
+          string eqtag (itvareqs->first.first);
+          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 = iteqtag->first;
+                break;
+              }
+
+          if (eqnumber == -1)
+            {
+              cerr << "ERROR: equation tag '" << eqtag << "' not found in var_model " << it->first << endl;
+              exit(EXIT_FAILURE);
+            }
+
+          int nonstationary = 0;
+          for (vector<pair<int, pair<string, string> > >::const_iterator iteqtag = equation_tags.begin();
+               iteqtag != equation_tags.end(); iteqtag++)
+            if (iteqtag->first == eqnumber)
+              if (iteqtag->second.first.compare("data_type") == 0
+                  && iteqtag->second.second.compare("nonstationary") == 0)
+                {
+                  nonstationary = 1;
+                  break;
+                }
+
+          equations[eqnumber]->get_arg1()->collectDynamicVariables(eEndogenous, lhs);
+          equations[eqnumber]->get_arg2()->collectDynamicVariables(eEndogenous, rhs);
+
+          final_map[make_pair(eqtag, eqnumber)] = make_pair(make_pair(nonstationary, lhs), rhs);
+        }
+
+      var_model_info[it->first] = final_map;
+    }
+}
+
+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 +3334,64 @@ 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::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  for (size_t i = 0; i < equations.size(); i++)
+    equations[i]->fillPacExpectationVarInfo(var_model_info);
+}
+
+void
+DynamicModel::substitutePacExpectation()
+{
+  map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > 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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > >::const_iterator it = subst_table.begin(); it != subst_table.end(); it++)
+    pac_expectation_info[it->second.second.first] = it->second.second.second;
+}
+
+void
+DynamicModel::writePacExpectationInfo(ostream &output) const
+{
+  int i = 1;
+  for (map<string, pair<int, pair<vector<int>, vector<int> > > >::const_iterator it = pac_expectation_info.begin();
+       it != pac_expectation_info.end(); it++, i++)
+    {
+      output << "M_.pac_expectation(" << i << ").var_model_name = '"
+             << it->first << "';" << endl
+             << "M_.pac_expectation(" << i << ").growth_param_index = "
+             << symbol_table.getTypeSpecificID(it->second.first) + 1 << ";" << endl
+             << "M_.pac_expectation(" << i << ").h0_param_indices = [";
+      for (vector<int>::const_iterator it1 = it->second.second.first.begin();
+           it1 != it->second.second.first.end(); it1++)
+        {
+          if (it1 != it->second.second.first.begin())
+            output << " ";
+          output << symbol_table.getTypeSpecificID(*it1) + 1;
+        }
+      output << "];" << endl
+           << "M_.pac_expectation(" << i << ").h1_param_indices = [";
+      for (vector<int>::const_iterator it1 = it->second.second.second.begin();
+           it1 != it->second.second.second.end(); it1++)
+        {
+          if (it1 != it->second.second.second.begin())
+            output << " ";
+          output << symbol_table.getTypeSpecificID(*it1) + 1;
+        }
+      output << "];" << endl;
+    }
+}
+
 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,
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 0d5c94b7e822f1be9bf9c92ef334b529e9ddbdeb..c2a3fa6466c4da7650032aea27ec42cd99a126ce 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -226,6 +226,10 @@ private:
   //! Used for var_expectation and var_model
   map<string, set<int> > var_expectation_functions_to_write;
 
+  //! Used for pac_expectation operator
+  // maps model_name to (growth_idx, (h0_indices, h1_indices))
+  map<string, pair<int, pair<vector<int>, vector<int> > > > pac_expectation_info;
+
   //!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 +284,19 @@ public:
   //! Set the equations that have non-zero second derivatives
   void setNonZeroHessianEquations(map<int, string> &eqs);
 
+  //! Fill var_model_info with variables associated with equation tags
+  void getVarModelVariablesFromEqTags(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   //! 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);
+  //! Add var_model info to pac_expectation nodes
+  void fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
+  //! Substitutes pac_expectation operator
+  void substitutePacExpectation();
+
+  //! Write Pac Expectation info
+  void writePacExpectationInfo(ostream &output) const;
 
   //! Adds informations for simulation in a binary file
   void Write_Inf_To_Bin_File_Block(const string &dynamic_basename, const string &bin_basename,
diff --git a/DynareBison.yy b/DynareBison.yy
index 8ec4ceeb8881c2e149e31869c776c8dc4430b3d7..2b59c8613787decafef852846996b6ad46eeb2c2 100644
--- a/DynareBison.yy
+++ b/DynareBison.yy
@@ -117,7 +117,7 @@ class ParsingDriver;
 %token USE_PENALIZED_OBJECTIVE_FOR_HESSIAN INIT_STATE RESCALE_PREDICTION_ERROR_COVARIANCE GENERATE_IRFS
 %token NAN_CONSTANT NO_STATIC NOBS NOCONSTANT NODISPLAY NOCORR NODIAGNOSTIC NOFUNCTIONS NO_HOMOTOPY
 %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 '(' MODEL_NAME EQUAL NAME COMMA DISCOUNT EQUAL hand_side COMMA GROWTH EQUAL hand_side')'
+            { $$ = driver.add_pac_expectation($5, $9, $13); }
           | MINUS hand_side %prec UMINUS
             { $$ = driver.add_uminus($2); }
           | PLUS hand_side
diff --git a/DynareFlex.ll b/DynareFlex.ll
index 1ebe6a5c74884bd1e4b0700589fdbd2d0d7a78bb..6248ed6aa5fe03ec45afab1ecd043a448b4bf928 100644
--- a/DynareFlex.ll
+++ b/DynareFlex.ll
@@ -676,6 +676,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;}
@@ -821,6 +822,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 2b149036d125cd0a7a7480efaa05f1d1655585dd..694c7f49369ded7e53526eacd4059dac0df2d9ce 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -493,6 +493,12 @@ NumConstNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
   return const_cast<NumConstNode *>(this);
 }
 
+expr_t
+NumConstNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  return const_cast<NumConstNode *>(this);
+}
+
 expr_t
 NumConstNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -560,6 +566,11 @@ NumConstNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
 {
 }
 
+void
+NumConstNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+}
+
 bool
 NumConstNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -1261,6 +1272,12 @@ VariableNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
   return const_cast<VariableNode *>(this);
 }
 
+expr_t
+VariableNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  return const_cast<VariableNode *>(this);
+}
+
 expr_t
 VariableNode::decreaseLeadsLags(int n) const
 {
@@ -1599,6 +1616,11 @@ VariableNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
 {
 }
 
+void
+VariableNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+}
+
 bool
 VariableNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -2798,6 +2820,13 @@ UnaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *> &
   return newAuxVar;
 }
 
+expr_t
+UnaryOpNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  expr_t argsubst = arg->substitutePacExpectation(subst_table);
+  return buildSimilarUnaryOpNode(argsubst, datatree);
+}
+
 expr_t
 UnaryOpNode::decreaseLeadsLags(int n) const
 {
@@ -2970,6 +2999,12 @@ UnaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mod
   arg->setVarExpectationIndex(var_model_info);
 }
 
+void
+UnaryOpNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  arg->fillPacExpectationVarInfo(var_model_info);
+}
+
 bool
 UnaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -4409,6 +4444,14 @@ BinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
   return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
 }
 
+expr_t
+BinaryOpNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  expr_t arg1subst = arg1->substitutePacExpectation(subst_table);
+  expr_t arg2subst = arg2->substitutePacExpectation(subst_table);
+  return buildSimilarBinaryOpNode(arg1subst, arg2subst, datatree);
+}
+
 expr_t
 BinaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -4486,6 +4529,13 @@ BinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
   arg2->setVarExpectationIndex(var_model_info);
 }
 
+void
+BinaryOpNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  arg1->fillPacExpectationVarInfo(var_model_info);
+  arg2->fillPacExpectationVarInfo(var_model_info);
+}
+
 bool
 BinaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -5166,6 +5216,15 @@ TrinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
   return buildSimilarTrinaryOpNode(arg1subst, arg2subst, arg3subst, datatree);
 }
 
+expr_t
+TrinaryOpNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  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);
+}
+
 expr_t
 TrinaryOpNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -5240,6 +5299,14 @@ TrinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_m
   arg3->setVarExpectationIndex(var_model_info);
 }
 
+void
+TrinaryOpNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  arg1->fillPacExpectationVarInfo(var_model_info);
+  arg2->fillPacExpectationVarInfo(var_model_info);
+  arg3->fillPacExpectationVarInfo(var_model_info);
+}
+
 bool
 TrinaryOpNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -5487,6 +5554,15 @@ AbstractExternalFunctionNode::substituteDiff(subst_table_t &subst_table, vector<
   return buildSimilarExternalFunctionNode(arguments_subst, datatree);
 }
 
+expr_t
+AbstractExternalFunctionNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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)->substitutePacExpectation(subst_table));
+  return buildSimilarExternalFunctionNode(arguments_subst, datatree);
+}
+
 expr_t
 AbstractExternalFunctionNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -5588,6 +5664,13 @@ AbstractExternalFunctionNode::setVarExpectationIndex(map<string, pair<SymbolList
     (*it)->setVarExpectationIndex(var_model_info);
 }
 
+void
+AbstractExternalFunctionNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
+    (*it)->fillPacExpectationVarInfo(var_model_info);
+}
+
 bool
 AbstractExternalFunctionNode::isVarModelReferenced(const string &model_info_name) const
 {
@@ -6855,7 +6938,7 @@ void
 VarExpectationNode::prepareForDerivation()
 {
   preparedForDerivation = true;
-  // All derivatives are null, so non_null_derivatives is left empty
+  // Come back
 }
 
 expr_t
@@ -6962,6 +7045,12 @@ VarExpectationNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNo
   return const_cast<VarExpectationNode *>(this);
 }
 
+expr_t
+VarExpectationNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  return const_cast<VarExpectationNode *>(this);
+}
+
 expr_t
 VarExpectationNode::differentiateForwardVars(const vector<string> &subset, subst_table_t &subst_table, vector<BinaryOpNode *> &neweqs) const
 {
@@ -7040,6 +7129,11 @@ VarExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &
   yidx = find(vs.begin(), vs.end(), datatree.symbol_table.getName(symb_id)) - vs.begin();
 }
 
+void
+VarExpectationNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+}
+
 expr_t
 VarExpectationNode::substituteStaticAuxiliaryVariable() const
 {
@@ -7059,3 +7153,450 @@ VarExpectationNode::writeJsonOutput(ostream &output,
          << ", yindex = " << yidx
          << ")";
 }
+
+PacExpectationNode::PacExpectationNode(DataTree &datatree_arg,
+                                       const string &model_name_arg,
+                                       const expr_t discount_arg,
+                                       const expr_t growth_arg) :
+  ExprNode(datatree_arg),
+  model_name(model_name_arg),
+  discount(discount_arg),
+  growth(growth_arg)
+{
+  datatree.pac_expectation_node_map[make_pair(model_name, make_pair(discount, growth))] = 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(model_name, discount, growth);
+}
+
+expr_t
+PacExpectationNode::cloneDynamic(DataTree &dynamic_datatree) const
+{
+  return dynamic_datatree.AddPacExpectation(model_name, discount, growth);
+}
+
+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 << ", ";
+      discount->writeOutput(output, output_type, temporary_terms, tef_terms);
+      output << ", ";
+      growth->writeOutput(output, output_type, temporary_terms, tef_terms);
+      output << RIGHT_PAR(output_type);
+      return;
+    }
+
+  // If current node is a temporary term
+  temporary_terms_t::const_iterator it = temporary_terms.find(const_cast<PacExpectationNode *>(this));
+  if (it != temporary_terms.end())
+    {
+      if (output_type == oMatlabDynamicModelSparse)
+        output << "T" << idx << "(it_)";
+      else
+        output << "T" << idx;
+      return;
+    }
+
+  output << "[h0, h1, d] = hVectors(params, H, ids, idns);";
+}
+
+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)
+{
+  cerr << "PacExpectationNode::eval: shouldn't arrive here." << endl;
+  exit(EXIT_FAILURE);
+}
+
+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);
+}
+
+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->writeJsonOutput(output, temporary_terms, tef_terms, isdynamic);
+  output << ", ";
+  growth->writeJsonOutput(output, temporary_terms, tef_terms, isdynamic);
+  output << ")";
+}
+
+void
+PacExpectationNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+{
+  map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > >::const_iterator it = var_model_info.find(model_name);
+  if (it == var_model_info.end())
+    {
+      cerr << "ERROR: could not find the declaration of " << model_name
+           << " referenced by pac_expectation operator" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  for (map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+    {
+      if (it1->second.first.second.size() != 1)
+        {
+          cerr << "ERROR " << model_name
+               << ": you may only have one variable on the LHS of a VAR" << endl;
+          exit(EXIT_FAILURE);
+        }
+      set<pair<int, int> >::const_iterator setit = it1->second.first.second.begin();
+      if (setit->second != 0)
+        {
+          cerr << "ERROR " << model_name
+               << ": the LHS variable of a VAR cannot appear with a lead or a lag" << endl;
+          exit(EXIT_FAILURE);
+        }
+      lhs.push_back(setit->first);
+      rhs.push_back(it1->second.second);
+      if (it1->second.first.first)
+        nonstationary_vars.push_back(lhs.back());
+    }
+}
+
+expr_t
+PacExpectationNode::substitutePacExpectation(map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table)
+{
+  map<const expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > >::iterator myit =
+    subst_table.find(const_cast<PacExpectationNode *>(this));
+  if (myit != subst_table.end())
+    return const_cast<BinaryOpNode *>(myit->second.first);
+
+  bool stationary_vars_present = true;
+  if (nonstationary_vars.size() == lhs.size())
+    stationary_vars_present = false;
+
+  bool nonstationary_vars_present = true;
+  if (nonstationary_vars.empty())
+    nonstationary_vars_present = false;
+
+  map<int, set<int > > all_rhs_vars; // lag -> set< symb_id > (all vars that appear at a given lag)
+  // Order RHS vars by time (already ordered by equation tag)
+  for (vector<set<pair<int, int> > >::const_iterator it = rhs.begin();
+       it != rhs.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())
+        {
+          cerr << "ERROR " << model_name << ": " << datatree.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 = all_rhs_vars.find(abs(it1->second));
+          if (mit == all_rhs_vars.end())
+            {
+              if (it1->second > 0)
+                {
+                  cerr << "ERROR " << model_name <<
+                    ": you cannot have a variable with a lead in a VAR" << endl;
+                  exit(EXIT_FAILURE);
+                }
+              set<int> si;
+              si.insert(it1->first);
+              all_rhs_vars[abs(it1->second)] = si;
+            }
+          else
+            mit->second.insert(it1->first);
+        }
+
+  vector<int> h0_indices, h1_indices;
+  expr_t subExpr = datatree.AddNonNegativeConstant("0");
+  if (stationary_vars_present)
+    for (map<int, set<int> >::const_iterator it = all_rhs_vars.begin();
+         it != all_rhs_vars.end(); it++)
+      for (set<int>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+        {
+          string param_name_h0("h0_" + model_name
+                               + "_var_" + datatree.symbol_table.getName(*it1)
+                               + "_lag_" + to_string(it->first));
+          int new_param_symb_id = datatree.symbol_table.addSymbol(param_name_h0, eParameter);
+          h0_indices.push_back(new_param_symb_id);
+          subExpr = datatree.AddPlus(subExpr,
+                                     datatree.AddTimes(datatree.AddVariable(new_param_symb_id),
+                                                       datatree.AddVariable(*it1, it->first)));
+        }
+
+  if (nonstationary_vars_present)
+    for (map<int, set<int > >::const_iterator it = all_rhs_vars.begin();
+         it != all_rhs_vars.end(); it++)
+      for (set<int>::const_iterator it1 = it->second.begin(); it1 != it->second.end(); it1++)
+        {
+          string param_name_h1("h1_" + model_name
+                               + "_var_" + datatree.symbol_table.getName(*it1)
+                               + "_lag_" + to_string(it->first));
+          int new_param_symb_id = datatree.symbol_table.addSymbol(param_name_h1, eParameter);
+          h1_indices.push_back(new_param_symb_id);
+          subExpr = datatree.AddPlus(subExpr,
+                                     datatree.AddTimes(datatree.AddVariable(new_param_symb_id),
+                                                       datatree.AddVariable(*it1, it->first)));
+        }
+
+  int pg_symb_id = datatree.symbol_table.addSymbol("pac_growth_parameter", eParameter);
+  subExpr = datatree.AddPlus(subExpr,
+                             datatree.AddTimes(datatree.AddVariable(pg_symb_id),
+                                               growth));
+
+  subst_table[const_cast<PacExpectationNode *>(this)] =
+    make_pair(dynamic_cast<BinaryOpNode *>(subExpr),
+              make_pair(model_name,
+                        make_pair(pg_symb_id,
+                                  make_pair(h0_indices, h1_indices))));
+
+  return subExpr;
+}
diff --git a/ExprNode.hh b/ExprNode.hh
index e2e20d6ae8b2d847a9664dbb1a9c4503055c2131..6ec139454c1663b1251057dedfc62b889a7c97c4 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -137,6 +137,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 */
@@ -463,6 +464,9 @@ enum ExprNodeOutputType
       //! 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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &subst_table) = 0;
+
       //! Add ExprNodes to the provided datatree
       virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const = 0;
 
@@ -481,6 +485,9 @@ enum ExprNodeOutputType
       //! Returns true if model_info_name is referenced by a VarExpectationNode
       virtual bool isVarModelReferenced(const string &model_info_name) const = 0;
 
+      //! Fills var_model info for pac_expectation node
+      virtual void fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info) = 0;
+
       //! Fills map
       virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const = 0;
     };
@@ -535,6 +542,7 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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;
@@ -547,6 +555,7 @@ 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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;
@@ -611,6 +620,7 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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;
@@ -623,6 +633,7 @@ 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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
@@ -710,6 +721,7 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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;
@@ -722,6 +734,7 @@ 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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
@@ -821,6 +834,7 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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;
@@ -839,6 +853,7 @@ public:
   expr_t getNonZeroPartofEquation() const;
   virtual bool isInStaticForm() const;
   virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info);
+  virtual void fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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
@@ -912,6 +927,7 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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;
@@ -924,6 +940,7 @@ 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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
@@ -1003,6 +1020,7 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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;
@@ -1017,6 +1035,7 @@ 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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
@@ -1185,6 +1204,73 @@ public:
   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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
+  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;
+  const expr_t discount, growth;
+  vector<int> lhs, nonstationary_vars;
+  vector<set<pair<int, int> > > rhs;
+public:
+  PacExpectationNode(DataTree &datatree_arg, const string &model_name, const expr_t discount_arg, const expr_t 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 expr_t, pair<const BinaryOpNode *, pair<string, pair<int, pair<vector<int>, vector<int> > > > > > &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,
@@ -1203,6 +1289,7 @@ 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 fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info);
   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 f81cb0963af1896aa687046adfbd6ecf1065a9a3..b6b909e769a985a3e62646c7c0110f0119f50731 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -362,7 +362,7 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
     }
 
   // Var Model
-  map<string, pair<SymbolList, int> > var_model_info;
+  map<string, pair<map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >, pair<SymbolList, int> > > var_model_info;
   for (vector<Statement *>::const_iterator it = statements.begin();
        it != statements.end(); it++)
     {
@@ -373,17 +373,28 @@ ModFile::transformPass(bool nostrict, bool stochastic, bool compute_xrefs, const
 
   if (!var_model_info.empty())
     {
-      dynamic_model.setVarExpectationIndices(var_model_info);
-      dynamic_model.addEquationsForVar(var_model_info);
+      map<string, pair<SymbolList, int> > var_model_info_var_expectation;
+      map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > var_model_info_pac;
+      for (map<string, pair<map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >, pair<SymbolList, int> > >::const_iterator it = var_model_info.begin(); it != var_model_info.end(); it++)
+        {
+          var_model_info_pac[it->first] = it->second.first;
+          var_model_info_var_expectation[it->first] = it->second.second;
+        }
+      dynamic_model.setVarExpectationIndices(var_model_info_var_expectation);
+      dynamic_model.addEquationsForVar(var_model_info_var_expectation);
+
+      dynamic_model.getVarModelVariablesFromEqTags(var_model_info_pac);
+      dynamic_model.fillPacExpectationVarInfo(var_model_info_pac);
+      dynamic_model.substitutePacExpectation();
     }
   dynamic_model.fillVarExpectationFunctionsToWrite();
 
-  if (symbol_table.predeterminedNbr() > 0)
-    dynamic_model.transformPredeterminedVariables();
-
   // Create auxiliary variable and equations for Diff operator
   dynamic_model.substituteDiff();
 
+  if (symbol_table.predeterminedNbr() > 0)
+    dynamic_model.transformPredeterminedVariables();
+
   // Create auxiliary vars for Expectation operator
   dynamic_model.substituteExpectation(mod_file_struct.partial_information);
 
@@ -750,6 +761,8 @@ ModFile::writeOutputFiles(const string &basename, bool clear_all, bool clear_glo
   // Initialize M_.det_shocks
   mOutputFile << "M_.det_shocks = [];" << endl;
 
+  dynamic_model.writePacExpectationInfo(mOutputFile);
+
   if (linear == 1)
     mOutputFile << "options_.linear = 1;" << endl;
 
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index cf8760113268d06665ffd2fd40efc8bb965b63af..0a8924a2aa5300dba938df1d26b4903b74ef3ca9 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -1506,6 +1506,17 @@ 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);
+
+  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));
   var_map[it->second] = symbol_list.getSymbols();
   symbol_list.clear();
@@ -2661,6 +2672,12 @@ ParsingDriver::add_var_expectation(string *arg1, string *arg2, string *arg3)
   return varExpectationNode;
 }
 
+expr_t
+ParsingDriver::add_pac_expectation(string *model_name, expr_t discount, expr_t growth)
+{
+  return data_tree->AddPacExpectation(*model_name, discount, growth);
+}
+
 expr_t
 ParsingDriver::add_exp(expr_t arg1)
 {
diff --git a/ParsingDriver.hh b/ParsingDriver.hh
index 872c391728f85ac41962731727f7eb353af02772..bf5c12b16c5ab61d1aa1e83498b5db053f1e10aa 100644
--- a/ParsingDriver.hh
+++ b/ParsingDriver.hh
@@ -686,6 +686,8 @@ 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(string *model_name, expr_t discount, expr_t growth);
   //! Writes token "diff(arg1)" to model tree
   expr_t add_diff(expr_t arg1);
   //! Writes token "adl(arg1, lag)" to model tree