diff --git a/ComputingTasks.cc b/ComputingTasks.cc
index a9745cf97ad0e7497293d4b5ecfea067cf0dfce2..9dcbe9e08e6a8bc882d0f2dff318e4827184cd0c 100644
--- a/ComputingTasks.cc
+++ b/ComputingTasks.cc
@@ -183,9 +183,9 @@ VarModelStatement::VarModelStatement(const SymbolList &symbol_list_arg,
 }
 
 void
-VarModelStatement::getVarModelNameAndVarList(map<string, SymbolList > &var_model_info)
+VarModelStatement::getVarModelNameAndVarList(map<string, pair<SymbolList, int> > &var_model_info)
 {
-  var_model_info[name] = symbol_list;
+  var_model_info[name] = make_pair(symbol_list, atoi(options_list.num_options.find("var.order")->second.c_str()));
 }
 
 void
diff --git a/ComputingTasks.hh b/ComputingTasks.hh
index 651b28ea71519f0fef3a4825b45f8f3688731ec4..1fd82e8ddb469a200fc620e54ad718487ddc7dc6 100644
--- a/ComputingTasks.hh
+++ b/ComputingTasks.hh
@@ -122,7 +122,7 @@ public:
                     const OptionsList &options_list_arg,
                     const string &name_arg,
                     const SymbolTable &symbol_table_arg);
-  void getVarModelNameAndVarList(map<string, SymbolList > &var_model_info);
+  void getVarModelNameAndVarList(map<string, pair<SymbolList, int> > &var_model_info);
   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 a01a2970c4b440ccb3e6260b913c530bfa8b64da..348a782c3a92ddf7e0ac7b1966f09c2f632785d0 100644
--- a/DynamicModel.cc
+++ b/DynamicModel.cc
@@ -3198,12 +3198,67 @@ DynamicModel::runTrendTest(const eval_context_t &eval_context)
 }
 
 void
-DynamicModel::setVarExpectationIndices(map<string, SymbolList > var_model_info)
+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)
+{
+  // 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;
+  for (map<string, pair<SymbolList, int> >::const_iterator it = var_model_info.begin();
+      it != var_model_info.end(); it++)
+    for (size_t i = 0; i < equations.size(); i++)
+      if (equations[i]->isVarModelReferenced(it->first))
+        {
+          vector<string> symbol_list = it->second.first.get_symbols();
+          int order = it->second.second;
+          for (vector<string>::const_iterator it1 = symbol_list.begin();
+               it1 != symbol_list.end(); it1++)
+            if (order > 2)
+              if (var_endos_and_lags.find(*it1) != var_endos_and_lags.end())
+                var_endos_and_lags[*it1] = min(var_endos_and_lags[*it1], -1*order);
+              else
+                var_endos_and_lags[*it1] = -1*order;
+          break;
+        }
+
+  if (var_endos_and_lags.empty())
+    return;
+
+  // Ensure that the minimum lag value exists in the model equations. If not, add an equation for it
+  for (size_t i = 0; i < equations.size(); i++)
+    equations[i]->getEndosAndMaxLags(model_endos_and_lags);
+
+  int count = 0;
+  for (map<string, int>::const_iterator it = var_endos_and_lags.begin();
+       it != var_endos_and_lags.end(); it++)
+    {
+      map<string, int>::const_iterator it1 = model_endos_and_lags.find(it->first);
+      if (it1 == model_endos_and_lags.end())
+        {
+          cerr << "ERROR: Variable used in var that is not used in the model: " << it->first << endl;
+          exit(EXIT_FAILURE);
+        }
+      else
+        if (it->second < it1->second)
+          {
+            int symb_id = symbol_table.getID(it->first);
+            expr_t newvar = AddVariable(symb_id, it->second);
+            expr_t auxvar = AddVariable(symbol_table.addVarModelEndoLagAuxiliaryVar(symb_id, it->second, newvar), 0);
+            addEquation(AddEqual(newvar, auxvar), -1);
+            addAuxEquation(AddEqual(newvar, auxvar));
+            count++;
+          }
+    }
+
+  if (count > 0)
+    cout << "Accounting for var_model lags not in model block: added " << count << " auxiliary variables and equations." << 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,
@@ -4594,8 +4649,8 @@ DynamicModel::substituteLeadLagInternal(aux_var_t type, bool deterministic_model
         case avDiffForward:
           cout << "forward vars";
           break;
-        case avMultiplier:
-          cerr << "avMultiplier encountered: impossible case" << endl;
+        default:
+          cerr << "DynamicModel::substituteLeadLagInternal: impossible case" << endl;
           exit(EXIT_FAILURE);
         }
       cout << ": added " << neweqs.size() << " auxiliary variables and equations." << endl;
diff --git a/DynamicModel.hh b/DynamicModel.hh
index ff9ddfd1a676a4d72a6e75b37943b8b1fc3dc909..ce81c3f933f3563aa7fd4c61006f5b19b92205db 100644
--- a/DynamicModel.hh
+++ b/DynamicModel.hh
@@ -248,7 +248,9 @@ public:
   void getNonZeroHessianEquations(map<int, string> &eqs) const;
 
   //! Set indices for var expectation in dynamic model file
-  void setVarExpectationIndices(map<string, SymbolList> 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);
 
   //! 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/ExprNode.cc b/ExprNode.cc
index c727124a09132a74463ac342c23b184f757ad555..cda13bcf147641e727601e3149f7f42f5504d188 100644
--- a/ExprNode.cc
+++ b/ExprNode.cc
@@ -278,6 +278,11 @@ ExprNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int lag_ar
   return false;
 }
 
+void
+ExprNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+}
+
 NumConstNode::NumConstNode(DataTree &datatree_arg, int id_arg) :
   ExprNode(datatree_arg),
   id(id_arg)
@@ -472,6 +477,11 @@ NumConstNode::isVariableNodeEqualTo(SymbolType type_arg, int variable_id, int la
   return false;
 }
 
+void
+NumConstNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+}
+
 bool
 NumConstNode::containsEndogenous(void) const
 {
@@ -509,10 +519,16 @@ NumConstNode::isInStaticForm() const
 }
 
 void
-NumConstNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+NumConstNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
 }
 
+bool
+NumConstNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return false;
+}
+
 expr_t
 NumConstNode::substituteStaticAuxiliaryVariable() const
 {
@@ -1508,8 +1524,25 @@ VariableNode::isInStaticForm() const
 }
 
 void
-VariableNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+VariableNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
+{
+}
+
+bool
+VariableNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return false;
+}
+
+void
+VariableNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
 {
+  string varname = datatree.symbol_table.getName(symb_id);
+  if (type == eEndogenous)
+    if (model_endos_and_lags.find(varname) == model_endos_and_lags.end())
+      model_endos_and_lags[varname] = min(model_endos_and_lags[varname], lag);
+    else
+      model_endos_and_lags[varname] = lag;
 }
 
 UnaryOpNode::UnaryOpNode(DataTree &datatree_arg, UnaryOpcode op_code_arg, const expr_t arg_arg, int expectation_information_set_arg, int param1_symb_id_arg, int param2_symb_id_arg) :
@@ -2585,11 +2618,23 @@ UnaryOpNode::isInStaticForm() const
 }
 
 void
-UnaryOpNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+UnaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
   arg->setVarExpectationIndex(var_model_info);
 }
 
+bool
+UnaryOpNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return arg->isVarModelReferenced(model_info_name);
+}
+
+void
+UnaryOpNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+  arg->getEndosAndMaxLags(model_endos_and_lags);
+}
+
 expr_t
 UnaryOpNode::substituteStaticAuxiliaryVariable() const
 {
@@ -3894,12 +3939,26 @@ BinaryOpNode::isInStaticForm() const
 }
 
 void
-BinaryOpNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+BinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
   arg1->setVarExpectationIndex(var_model_info);
   arg2->setVarExpectationIndex(var_model_info);
 }
 
+bool
+BinaryOpNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return arg1->isVarModelReferenced(model_info_name)
+    || arg2->isVarModelReferenced(model_info_name);
+}
+
+void
+BinaryOpNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+  arg1->getEndosAndMaxLags(model_endos_and_lags);
+  arg2->getEndosAndMaxLags(model_endos_and_lags);
+}
+
 expr_t
 BinaryOpNode::substituteStaticAuxiliaryVariable() const
 {
@@ -4573,13 +4632,29 @@ TrinaryOpNode::isInStaticForm() const
 }
 
 void
-TrinaryOpNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+TrinaryOpNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
   arg1->setVarExpectationIndex(var_model_info);
   arg2->setVarExpectationIndex(var_model_info);
   arg3->setVarExpectationIndex(var_model_info);
 }
 
+bool
+TrinaryOpNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return arg1->isVarModelReferenced(model_info_name)
+    || arg2->isVarModelReferenced(model_info_name)
+    || arg3->isVarModelReferenced(model_info_name);
+}
+
+void
+TrinaryOpNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+  arg1->getEndosAndMaxLags(model_endos_and_lags);
+  arg2->getEndosAndMaxLags(model_endos_and_lags);
+  arg3->getEndosAndMaxLags(model_endos_and_lags);
+}
+
 expr_t
 TrinaryOpNode::substituteStaticAuxiliaryVariable() const
 {
@@ -4888,12 +4963,28 @@ AbstractExternalFunctionNode::isInStaticForm() const
 }
 
 void
-AbstractExternalFunctionNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+AbstractExternalFunctionNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
   for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); it++)
     (*it)->setVarExpectationIndex(var_model_info);
 }
 
+bool
+AbstractExternalFunctionNode::isVarModelReferenced(const string &model_info_name) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); ++it)
+    if (!(*it)->isVarModelReferenced(model_info_name))
+      return true;
+  return false;
+}
+
+void
+AbstractExternalFunctionNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+  for (vector<expr_t>::const_iterator it = arguments.begin(); it != arguments.end(); ++it)
+    (*it)->getEndosAndMaxLags(model_endos_and_lags);
+}
+
 pair<int, expr_t>
 AbstractExternalFunctionNode::normalizeEquation(int var_endo, vector<pair<int, pair<expr_t, expr_t> > >  &List_of_Op_RHS) const
 {
@@ -6085,10 +6176,21 @@ VarExpectationNode::isInStaticForm() const
   return false;
 }
 
+bool
+VarExpectationNode::isVarModelReferenced(const string &model_info_name) const
+{
+  return model_name == model_info_name;
+}
+
+void
+VarExpectationNode::getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const
+{
+}
+
 void
-VarExpectationNode::setVarExpectationIndex(map<string, SymbolList> var_model_info)
+VarExpectationNode::setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info)
 {
-  vector<string> vs = var_model_info[model_name].get_symbols();;
+  vector<string> vs = var_model_info[model_name].first.get_symbols();
   yidx = find(vs.begin(), vs.end(), datatree.symbol_table.getName(symb_id)) - vs.begin();
 }
 
diff --git a/ExprNode.hh b/ExprNode.hh
index 6e57a3c7e400366937164533f4dd220cd3523259..cc566b7841b6852860a1020d6f729b9f16b73df5 100644
--- a/ExprNode.hh
+++ b/ExprNode.hh
@@ -452,8 +452,14 @@ public:
   //! Substitute auxiliary variables by their expression in static model
   virtual expr_t substituteStaticAuxiliaryVariable() const = 0;
 
-  // Add index information for var_model variables
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info) = 0;
+  //! Add index information for var_model variables
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, int> > &var_model_info) = 0;
+
+  //! Returns true if model_info_name is referenced by a VarExpectationNode
+  virtual bool isVarModelReferenced(const string &model_info_name) const = 0;
+
+  //! Fills map
+  virtual void getEndosAndMaxLags(map<string, int> &model_endos_and_lags) const = 0;
 };
 
 //! Object used to compare two nodes (using their indexes)
@@ -514,7 +520,9 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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;
 };
 
@@ -581,7 +589,9 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -668,7 +678,9 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -774,7 +786,9 @@ public:
   //! Returns the non-zero hand-side of an equation (that must have a hand side equal to zero)
   expr_t getNonZeroPartofEquation() const;
   virtual bool isInStaticForm() const;
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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
   virtual expr_t substituteStaticAuxiliaryVariable() const;
   //! Substitute auxiliary variables by their expression in static model auxiliary variable definition
@@ -850,7 +864,9 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -933,7 +949,9 @@ public:
   virtual expr_t cloneDynamic(DataTree &dynamic_datatree) const = 0;
   virtual expr_t removeTrendLeadLag(map<int, expr_t> trend_symbols_map) const;
   virtual bool isInStaticForm() const;
-  virtual void setVarExpectationIndex(map<string, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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
   virtual expr_t substituteStaticAuxiliaryVariable() const;
 };
@@ -1100,7 +1118,9 @@ public:
   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, SymbolList> var_model_info);
+  virtual void setVarExpectationIndex(map<string, pair<SymbolList, 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 310e6ce5285d5d3dd9ce40bd0c7e3bfc02c68c81..0b52a1b561c51bb9e3454610c048656a18ebcd5d 100644
--- a/ModFile.cc
+++ b/ModFile.cc
@@ -335,6 +335,23 @@ ModFile::transformPass(bool nostrict, bool compute_xrefs)
         }
     }
 
+  // Var Model
+  map<string, pair<SymbolList, int> > var_model_info;
+  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);
+    }
+
+  if (!var_model_info.empty())
+    {
+      dynamic_model.setVarExpectationIndices(var_model_info);
+      dynamic_model.addEquationsForVar(var_model_info);
+    }
+  dynamic_model.fillVarExpectationFunctionsToWrite();
+
   if (symbol_table.predeterminedNbr() > 0)
     dynamic_model.transformPredeterminedVariables();
 
@@ -411,19 +428,6 @@ ModFile::transformPass(bool nostrict, bool compute_xrefs)
         exit(EXIT_FAILURE);
       }
 
-  // Var Model
-  map<string, SymbolList > var_model_info;
-  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);
-    }
-  if (!var_model_info.empty())
-    dynamic_model.setVarExpectationIndices(var_model_info);
-  dynamic_model.fillVarExpectationFunctionsToWrite();
-
   // Freeze the symbol table
   symbol_table.freeze();
 
diff --git a/SymbolTable.cc b/SymbolTable.cc
index ce7fed3f8db817af1a41c23c71b3428f428a35aa..183c45c84872f70ec65d873aaf4db962a322919b 100644
--- a/SymbolTable.cc
+++ b/SymbolTable.cc
@@ -352,6 +352,7 @@ SymbolTable::writeOutput(ostream &output) const throw (NotYetFrozenException)
             break;
           case avEndoLag:
           case avExoLag:
+          case avVarModel:
             output << "M_.aux_vars(" << i+1 << ").orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id())+1 << ";" << endl
                    << "M_.aux_vars(" << i+1 << ").orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
             break;
@@ -472,6 +473,7 @@ SymbolTable::writeCOutput(ostream &output) const throw (NotYetFrozenException)
 	      break;
 	    case avEndoLag:
 	    case avExoLag:
+            case avVarModel:
 	      output << "av[" << i << "].orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
 		     << "av[" << i << "].orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
 	      break;
@@ -564,6 +566,7 @@ SymbolTable::writeCCOutput(ostream &output) const throw (NotYetFrozenException)
           break;
         case avEndoLag:
         case avExoLag:
+        case avVarModel:
           output << "av" << i << ".orig_index = " << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) << ";" << endl
                  << "av" << i << ".orig_lead_lag = " << aux_vars[i].get_orig_lead_lag() << ";" << endl;
           break;
@@ -683,6 +686,28 @@ SymbolTable::addExpectationAuxiliaryVar(int information_set, int index, expr_t e
   return symb_id;
 }
 
+int
+SymbolTable::addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException)
+{
+  int symb_id;
+  ostringstream varname;
+  varname << "AUX_VARMODEL_" << orig_symb_id << "_" << -orig_lead_lag;
+
+  try
+    {
+      symb_id = addSymbol(varname.str(), eEndogenous);
+    }
+  catch (AlreadyDeclaredException &e)
+    {
+      cerr << "ERROR: you should rename your variable called " << varname.str() << ", this name is internally used by Dynare" << endl;
+      exit(EXIT_FAILURE);
+    }
+
+  aux_vars.push_back(AuxVarInfo(symb_id, avVarModel, orig_symb_id, orig_lead_lag, 0, 0, expr_arg));
+
+  return symb_id;
+}
+
 int
 SymbolTable::addMultiplierAuxiliaryVar(int index) throw (FrozenException)
 {
@@ -980,6 +1005,7 @@ SymbolTable::writeJuliaOutput(ostream &output) const throw (NotYetFrozenExceptio
               break;
             case avEndoLag:
             case avExoLag:
+            case avVarModel:
               output << getTypeSpecificID(aux_vars[i].get_orig_symb_id()) + 1 << ", "
                      << aux_vars[i].get_orig_lead_lag() << ", NaN, NaN";
               break;
diff --git a/SymbolTable.hh b/SymbolTable.hh
index 37f95c181b07db5bc81bfe959fd90897e7f7a06b..6d55e2f45c590f47bfece55dae6dbb6237663b58 100644
--- a/SymbolTable.hh
+++ b/SymbolTable.hh
@@ -42,7 +42,8 @@ enum aux_var_t
     avExoLag = 3,         //!< Substitute for exo lags >= 2
     avExpectation = 4,    //!< Substitute for Expectation Operator
     avDiffForward = 5,    //!< Substitute for the differentiate of a forward variable
-    avMultiplier = 6      //!< Multipliers for FOC of Ramsey Problem
+    avMultiplier = 6,     //!< Multipliers for FOC of Ramsey Problem
+    avVarModel = 7        //!< Variable for var_model with order > abs(min_lag()) present in model
   };
 
 //! Information on some auxiliary variables
@@ -251,6 +252,9 @@ public:
     Throws an exception if match not found.
   */
   int searchAuxiliaryVars(int orig_symb_id, int orig_lead_lag) const throw (SearchFailedException);
+  //! Adds an auxiliary variable when var_model is used with an order that is greater in absolute value
+  //! than the largest lag present in the model.
+  int addVarModelEndoLagAuxiliaryVar(int orig_symb_id, int orig_lead_lag, expr_t expr_arg) throw (AlreadyDeclaredException, FrozenException);
   //! Returns the number of auxiliary variables
   int AuxVarsSize() const { return aux_vars.size(); };
   //! Retruns expr_node for an auxiliary variable