diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index bf5e74bd9bbf27dd62aa3fec42d2c328372643c5..a50236a5c07d90efebd63e545900b1a02a9a36c4 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -260,31 +260,90 @@ 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<map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >, pair<SymbolList, int> > > &var_model_info) const
+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");
+  var_model_info[name] = make_pair(symbol_list, atoi(it->second.c_str()));
+}
+
+void
+VarModelStatement::getVarModelEqTags(vector<string> &var_model_eqtags) const
 {
-  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));
-    }
+    return;
+
+  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;
+  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_arg.begin();
+       it != rhs_arg.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
@@ -298,8 +357,79 @@ 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.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;
   output << "M_.var." << name << " = options_.var;" << endl
          << "clear options_.var;" << endl;
+
+  cout << "DONE writing varmodel" << endl;
 }
 
 void
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index 1799f37da142bd10c4d8b2269c20e4b2115e0ed7..b0090cf7703b9ba561410ae8bf43b91c6cc3e89b 100644
--- a/ComputingTasks.hh
+++ b/ComputingTasks.hh
@@ -125,11 +125,23 @@ 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<bool> nonstationary, diff;
 public:
   VarModelStatement(const SymbolList &symbol_list_arg,
                     const OptionsList &options_list_arg,
-                    const string &name_arg);
-  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;
+                    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/DynamicModel.cc b/DynamicModel.cc
index da197ef821168a481d887162c268d0210ce71fde..72494a77f92068073bb28030d9be7fa25d919d46 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3235,50 +3235,91 @@ DynamicModel::runTrendTest(const eval_context_t &eval_context)
 }
 
 void
-DynamicModel::getVarModelVariablesFromEqTags(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, 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 (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++)
+  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)
         {
-          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;
-              }
+          cerr << "ERROR: equation tag '" << eqtag << "' not found" << endl;
+          exit(EXIT_FAILURE);
+        }
 
-          if (eqnumber == -1)
+      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)
             {
-              cerr << "ERROR: equation tag '" << eqtag << "' not found in var_model " << it->first << endl;
-              exit(EXIT_FAILURE);
+              nonstationary_bool = true;
+              break;
             }
 
-          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_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);
 
-          equations[eqnumber]->get_arg1()->collectDynamicVariables(eEndogenous, lhs);
-          equations[eqnumber]->get_arg2()->collectDynamicVariables(eEndogenous, rhs);
+      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);
+        }
 
-          final_map[make_pair(eqtag, eqnumber)] = make_pair(make_pair(nonstationary, lhs), rhs);
+      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);
         }
 
-      var_model_info[it->first] = final_map;
+      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]->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);
     }
 }
 
@@ -3342,10 +3383,12 @@ DynamicModel::addEquationsForVar(map<string, pair<SymbolList, int> > &var_model_
 }
 
 void
-DynamicModel::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+DynamicModel::fillPacExpectationVarInfo(string &var_model_name,
+                                        map<int, set<int > > &rhs,
+                                        vector<bool> &nonstationary)
 {
   for (size_t i = 0; i < equations.size(); i++)
-    equations[i]->fillPacExpectationVarInfo(var_model_info);
+    equations[i]->fillPacExpectationVarInfo(var_model_name, rhs, nonstationary);
 }
 
 void
diff --git a/DynamicModel.hh b/DynamicModel.hh
index 676718a679a9a63fbcb138bd5191df9d01dff638..d87b412b99d3f2febc387dad51107057200f99f7 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -283,14 +283,24 @@ 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);
+  //! 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);
   //! 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);
   //! 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);
+  void fillPacExpectationVarInfo(string &var_model_name,
+                                 map<int, set<int > > &rhs,
+                                 vector<bool> &nonstationary);
   //! Substitutes pac_expectation operator
   void substitutePacExpectation();
 
diff --git a/ExprNode.cc b/ExprNode.cc
index 9b9407885feabf22d26c10aa841101d4b4679c60..655c07234f1fbde1c94e02565241d93dfb3348b0 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -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()
 {
@@ -567,7 +579,7 @@ 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)
+NumConstNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
 }
 
@@ -1560,6 +1572,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
 {
@@ -1617,7 +1635,7 @@ 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)
+VariableNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
 }
 
@@ -2792,6 +2810,14 @@ UnaryOpNode::substituteAdl() 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
 {
@@ -3000,9 +3026,9 @@ UnaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mod
 }
 
 void
-UnaryOpNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+UnaryOpNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
-  arg->fillPacExpectationVarInfo(var_model_info);
+  arg->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
 }
 
 bool
@@ -4444,6 +4470,12 @@ BinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
   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)
 {
@@ -4530,10 +4562,10 @@ BinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_mo
 }
 
 void
-BinaryOpNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+BinaryOpNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
-  arg1->fillPacExpectationVarInfo(var_model_info);
-  arg2->fillPacExpectationVarInfo(var_model_info);
+  arg1->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
+  arg2->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
 }
 
 bool
@@ -5216,6 +5248,12 @@ TrinaryOpNode::substituteDiff(subst_table_t &subst_table, vector<BinaryOpNode *>
   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)
 {
@@ -5300,11 +5338,11 @@ TrinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_m
 }
 
 void
-TrinaryOpNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+TrinaryOpNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
-  arg1->fillPacExpectationVarInfo(var_model_info);
-  arg2->fillPacExpectationVarInfo(var_model_info);
-  arg3->fillPacExpectationVarInfo(var_model_info);
+  arg1->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
+  arg2->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
+  arg3->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
 }
 
 bool
@@ -5554,6 +5592,15 @@ AbstractExternalFunctionNode::substituteDiff(subst_table_t &subst_table, vector<
   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)
 {
@@ -5665,10 +5712,10 @@ AbstractExternalFunctionNode::setVarExpectationIndex(map<string, pair<SymbolList
 }
 
 void
-AbstractExternalFunctionNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+AbstractExternalFunctionNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
-    (*it)->fillPacExpectationVarInfo(var_model_info);
+    (*it)->fillPacExpectationVarInfo(var_model_name, rhs_arg, nonstationary_arg);
 }
 
 bool
@@ -6969,6 +7016,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
 {
@@ -7130,7 +7183,7 @@ VarExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &
 }
 
 void
-VarExpectationNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+VarExpectationNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
 }
 
@@ -7162,8 +7215,8 @@ PacExpectationNode::PacExpectationNode(DataTree &datatree_arg,
   model_name(model_name_arg),
   discount_symb_id(discount_symb_id_arg),
   growth_symb_id(growth_symb_id_arg),
-  stationary_vars_present(true),
-  nonstationary_vars_present(true)
+  stationary_vars_present(false),
+  nonstationary_vars_present(false)
 {
   datatree.pac_expectation_node_map[make_pair(model_name, make_pair(discount_symb_id, growth_symb_id))] = this;
 }
@@ -7341,6 +7394,12 @@ PacExpectationNode::compile(ostream &CompileCode, unsigned int &instruction_numb
   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
 {
@@ -7486,74 +7545,23 @@ PacExpectationNode::writeJsonOutput(ostream &output,
 }
 
 void
-PacExpectationNode::fillPacExpectationVarInfo(map<string, map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > > > &var_model_info)
+PacExpectationNode::fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg)
 {
-  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);
-    }
-
-  vector<set<pair<int, int> > > rhs;
-  vector<int> lhs, nonstationary_vars;
-  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());
-    }
-
-  if (nonstationary_vars.size() == lhs.size())
-    stationary_vars_present = false;
+  if (model_name != var_model_name)
+    return;
 
-  if (nonstationary_vars.empty())
-    nonstationary_vars_present = false;
+  z_vec = rhs_arg;
 
-  // 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);
-        }
+  for (vector<bool>::const_iterator it = nonstationary_arg.begin();
+       it != nonstationary_arg.end(); it++)
+    {
+      if (*it)
+        nonstationary_vars_present = true;
       else
-        {
-          map<int, set<int> >::iterator mit = z_vec.find(abs(it1->second));
-          if (mit == z_vec.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);
-              z_vec[abs(it1->second)] = si;
-            }
-          else
-            mit->second.insert(it1->first);
-        }
+        stationary_vars_present = true;
+      if (nonstationary_vars_present && stationary_vars_present)
+        break;
+    }
 }
 
 expr_t
diff --git a/ExprNode.hh b/ExprNode.hh
index 80052b1bbc1f2814e97a8a5e105c5cdf205629a1..156df34af11a4ccef5d0991e34146b96ef448afb 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -440,6 +440,9 @@ class ExprNode
       //! 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)
@@ -487,7 +490,7 @@ class ExprNode
       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;
+      virtual void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_arg) = 0;
 
       //! Fills map
       virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const = 0;
@@ -549,6 +552,7 @@ public:
   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;
@@ -556,7 +560,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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;
@@ -627,6 +631,7 @@ public:
   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;
@@ -634,7 +639,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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
@@ -728,6 +733,7 @@ public:
   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;
@@ -735,7 +741,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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
@@ -841,6 +847,7 @@ public:
   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;
@@ -854,7 +861,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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
@@ -934,6 +941,7 @@ public:
   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;
@@ -941,7 +949,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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
@@ -1028,6 +1036,7 @@ public:
   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;
@@ -1036,7 +1045,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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
@@ -1215,6 +1224,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;
@@ -1224,7 +1234,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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;
@@ -1284,6 +1294,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;
@@ -1293,7 +1304,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 void fillPacExpectationVarInfo(string &var_model_name, map<int, set<int > > &rhs_arg, vector<bool> &nonstationary_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 ee1cc504dca72b5cd95266b357c2bd9fd6c4d761..7d1f757af12f2b5034a198a7af56f66aa2fc5de1 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -361,37 +361,47 @@ 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<map<pair<string, int>, pair<pair<int, set<pair<int, int> > >, set<pair<int, int> > > >, 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.fillPacExpectationVarInfo(var_model_name, rhs_pac, nonstationary);
+              dynamic_model.substitutePacExpectation();
+            }
+        }
     }
 
-  if (!var_model_info.empty())
+  if (!var_model_info_var_expectation.empty())
     {
-      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();
 
-  // Create auxiliary variable and equations for Diff operator
-  dynamic_model.substituteDiff();
-
   if (symbol_table.predeterminedNbr() > 0)
     dynamic_model.transformPredeterminedVariables();
 
diff --git a/ParsingDriver.cc b/ParsingDriver.cc
index bd771b853c24d92cef19a7ee17abd6b86fd15172..47eeba3c7e6c10658704578543314dccf1cc1642 100644
--- a/ParsingDriver.cc
+++ b/ParsingDriver.cc
@@ -1517,7 +1517,7 @@ ParsingDriver::var_model()
     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->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();